text:
# Rigidity of the elastic domain structure near the boundary of its existence in thin epitaxial films

\[

## Abstract

We consider an interesting and practically important case of elastic domain structure, which is the analogue of c/a domain pattern with 90 walls in perovskites, and is solvable analytically for arbitrary misfit strain. There is no critical thickness, below which the domain structure cannot exist, when the “extrinsic” misfit is zero and the domains are of equal width. At the boundary of polydomain-monodomain transition the period of the pattern diverges, as does the dynamic stiffness of the domain structure. It is unlikely, therefore, that one can achieve a softness of the dielectric response of the c/a elastic domains in ferroelectric-ferroelastic thin films.

\]

Equilibrium domain structures in epitaxial ferroelectric films may appear even in the case of complete compensation of the depolarizing electric fields by the electrodes or finite conductivity of the film. This takes place if the film also behaves as a ferroelastic. Formation of a ferroelastic domain structure was first considered a while ago as a mechanism to relax the misfit imposed by a substrate . In widely studied perovskite ferroelectrics, which are improper ferroelastics, several elastic domain (and heterophase) structures were predicted . One of the typical structures, so-called $`c/a/c/a\mathrm{}`$ domain pattern, is of a special fundamental and practical interest. Firstly, the ratio of the widths of the $`c`$\- and $`a`$-domains changes with external conditions (film thickness and misfit strain). At some parameters the domain pattern does not exist, i.e. there is a phase transition between multi- and monodomain states. Although the corresponding phase diagrams have been discussed by several authors (see  and references therein), the behavior close to the phase transition was not studied. The model calculations have predicted the response of the $`c/a`$ domain structure to be large, even giant for some specific values of the (“extrinsic”, see below) misfit $`w`$ and the film thickness $`l`$, and some data was even interpreted in this way. The enhancement of the response was expected near the border of existence of the $`c/a`$ structure in the $`wl`$ plane, i.e. near the phase transition. Previous analyses, mainly numerical, have not actually specified the system parameters needed to realize the giant susceptibility. We attempt, therefore, an exact solution of a simplest relevant model, analogous to the actual $`c/a`$ structures. The present results put strict constraints on the system parameters where the softness of the dielectric response might be expected. It becomes clear that it would be very difficult to realize the domain structure with such a property. As a corollary, the existing data  cannot be interpreted as being due to a large contribution of the elastic domain walls. We reveal in the present paper some unusual features of the polydomain-monodomain phase transition and show that the response of the system becomes more ”rigid” in the vicinity of this transition.

We present here the first analytical calculations of the energy of a domain pattern analogous to the standard $`c/a`$ structure. The structure can be viewed, without loss of generality, as a result of a ferroelastic transition in an epitaxial film, which breaks the symmetry of the parent phase and the substrate. All the strains are considered in the reference frame of the high symmetry phase, and the $`z`$-axis is selected to be perpendicular to the plane of the film. We consider the system far from the transition, where the pattern consists of the domains having the spontaneous strain, $`u_{xx}^0u_{zz}^0,`$ of opposite signs, separated by the domain walls inclined at 45 with respect to the $`z`$-axis, as dictated by the elastic compatibility conditions and the lattice symmetry, Fig. 1. In ferroelectric perovskites this would correspond to the 90 walls separating the domains with the polarization parallel and perpendicular to the plane. The film consists of $`a`$ and $`c`$type domains with the widths $`a_1`$ and $`a_2,`$ respectively. From the exact expressions for the free energy $`F`$ \[see Eqs.(6),(10)\], we calculate the period $`2a=a_1+a_2`$ of the structure, the parameter of asymmetry $`\delta =`$ $`\left(a_1a_2\right)/2a,`$ and the corresponding stiffness $`\kappa `$. The stiffness $`\kappa `$ is closely related to the “mechanical force constant” discussed in Ref. . We shall consider two stiffnesses: the usually measured dynamic stiffness $`\kappa _{\mathrm{}}`$, calculated at fixed period of the structure, and the static stiffness $`\kappa _0,`$ corresponding to the situation where the period of domain structure is allowed to relax via creation or annihilation of the domain walls. The latter corresponds to very slow processes, making its observation very challenging. One certainly deals with the dynamic stiffness $`\kappa _{\mathrm{}}`$ in applications. We shall show that the two stiffnesses exhibit quite different behavior near the transition: the static one diminishes, while the dynamic one diverges. We also study the change of the stiffness with the film thickness and discuss the possibilities of obtaining the films with small dynamic stiffness.

We consider first the case where the sole origin of the misfit is the spontaneous strain, $`u_{xx}^0u_{zz}^0,`$ appearing at the phase transition. We shall call this misfit “intrinsic”, to distinguish it from “extrinsic” misfit of any other origin. The latter includes e.g. the misfit of parent phase with the substrate, and other, “noncritical”, strain components, which might appear during the phase transition. The effect of extrinsic strain is a focus of the present paper. Similar to previous authors  we neglect the uncompensated electric fields, which would only increase the stiffness.

The method of calculating the energy of the domain structure is the same as in our previous paper  but with different equations of state for the stress tensor components $`\sigma _{ik}`$:

$`\sigma _{zz}`$ $`=`$ $`\left(\lambda +2\mu \right)e_{zz}+\lambda \left(e_{xx}+e_{yy}\right),`$ (1)
$`\sigma _{xz}`$ $`=`$ $`2\mu e_{xz},`$ (2)
with $`\sigma _{xx}`$ and $`\sigma _{yy}`$ components obtained from (1) by cyclic permutation of $`x,y,z`$. Here $`\lambda ,\mu `$ are the Lamé coefficients, and $`e_{ik}=u_{ik}u_{ik}^0`$ is the elastic strain, with $`u_{ik}^0`$ the components of the spontaneous strain. In contact with substrate, the film with symmetry-breaking misfit with the substrate must split into domains in such a way that the average in-plane strain will be zero, since the homogeneous strain in the substrate would cost an infinite energy. Without the “extrinsic” strain the spontaneous strain would alternate from domain to domain (which all would have the same width) as $`u_{xx}^0=u_{zz}^0\pm u_0,`$ $`u_{yy}^0=0.`$ On the other hand, with the “extrinsic” misfit strain $`w`$ the distribution of the spontaneous strain in the domains would be

$$u_{xx}^0(u_{zz}^0)\pm ()u_0+w,u_{yy}^0=w.$$
(3)
We assume for the substrate the same equations of state, as (1),(2), but with $`u_{xx}^0=u_{yy}^0=u_{zz}^0=0.`$

To find the energy of the domain structure, we have to determine the energy of the homogeneous and inhomogeneous (stray) stresses. The energy of the homogeneous stresses $`\left(F_h\right),`$ which appears for domains with non-equal widths (i.e. for $`\delta 0),`$ can be easily found for an epitaxial film with any domain structure, provided that the elastic moduli in all domains are the same, and the non-linear effects can be neglected. We can readily find $`F_h`$ from the general expression for the elastic energy $`F_{el}=\frac{1}{2}𝑑V\sigma _{ik}u_{ik}^0`$  as

$`{\displaystyle \frac{F_h}{𝒜}}`$ $`=`$ $`{\displaystyle \frac{l}{2}}(\overline{\sigma }_{xx}\overline{u}_{xx}^0+\overline{\sigma }_{yy}\overline{u}_{yy}^0+\overline{\sigma }_{zz}\overline{u}_{zz}^0)`$ (4)
$`=`$ $`Ml\left[\left(\overline{u}_{xx}^0\right)^2+\left(\overline{u}_{yy}^0\right)^2+2\nu \overline{u}_{xx}^0\overline{u}_{yy}^0\right]`$ (5)
$`=`$ $`Mu_0^2l\left[(\delta \delta _0)^2+w^2(1\nu ^2)/u_0^2\right]`$ (6)
where $`\delta _0=w(1+\nu )/u_0`$ is the relative extrinsic misfit, $`\nu =\lambda /2(\lambda +\mu )`$ the Poisson’s ratio, and $`M=2\mu \left(\lambda +\mu \right)/\left(\lambda +2\mu \right)E/2(1\nu ^2)`$, where $`E`$ is the Young’s modulus. The overbar here and below marks the averaging over the film (domain pattern). The condition $`\overline{\sigma }_{zz}=0`$ was used together with Eq.(1) to obtain for the planar stresses $`\overline{\sigma }_{xx}=2M(\overline{e}_{xx}+\nu \overline{e}_{yy}),`$ $`\overline{\sigma }_{yy}=\overline{\sigma }_{xx}(xy),`$ where $`\overline{e}_{xx}\overline{u}_{xx}\overline{u}_{xx}^0=\overline{u}_{xx}^0=u_0\delta w,`$ $`\overline{e}_{yy}=\overline{u}_{yy}^0=w,`$ Eq.(3). We have also used the absence of the in-plane strains, $`\overline{u}_{xx}=\overline{u}_{yy}=0`$, imposed by the substrate, and an obvious relation $`\overline{u}_{xx}^0=\frac{1}{2}(1+\delta )(u_0+w)+\frac{1}{2}(1\delta )(u_0+w)=u_0\delta +w.`$ Note that Eq.(5) is rather general and can be applied to other kinds of domain patterns in the epitaxial films (cf. ).

The stray energy of inhomogeneous stresses in a stripe-like domain structure periodic in $`x`$direction, Fig. 1, is to be found from exact solutions for the strain field produced by the domains. The pattern is defined by the distribution of spontaneous strains $`u_{xx}^0\left(x\right)`$ and $`u_{zz}^0\left(x\right).`$ The condition of the local equilibrium, $`\sigma _{ik}/x_k=0,`$ gives two sets of two equations for the displacement components $`u_x,`$ $`u_z`$ in the film and the substrate with the use of a standard relation $`2u_{ik}=u_i/x_k+u_k/x_i`$. The boundary conditions are given by the absence of stresses, $`\sigma _{zz}=`$ $`\sigma _{xz}=0`$ at the free surface of the film ($`z=l`$) and at $`z\mathrm{}`$. Both the strain tensor and the displacement vector are continuous at the film-substrate interface. The original system of partial differential equations is reduced to a system of the ordinary differential equations with the use of the Fourier transform and then solved, as described in detail in Refs. , and we obtain the following simple expression for the stray energy per area $`𝒜:`$

$$\frac{F_{stray}}{𝒜}=\frac{2Mu_0^2a}{\pi ^3}\underset{n=1}{\overset{\mathrm{}}{}}\left[1\left(1\right)^n\mathrm{cos}\pi n\delta \right]\frac{1\left(1+2k^2\right)e^{2k}}{n^3},$$
(7)
where $`k=\pi nl/a`$. The series is calculated with the result

$`{\displaystyle \frac{F_{stray}}{𝒜}}`$ $`=`$ $`{\displaystyle \frac{2}{\pi ^3}}Mu_0^2a[\zeta \left(3\right)ReLi_3(e^{i\pi \delta })`$ (10)
$`Li_3\left(e^b\right)+ReLi_3\left(e^{b+i\pi \delta }\right)`$
$`+{\displaystyle \frac{b^2}{2}}Re\mathrm{ln}{\displaystyle \frac{1e^b}{1+e^{b+i\pi \delta }}}],`$
where $`Li_n\left(z\right)_{s=1}^{\mathrm{}}z^s/s^n=z\mathrm{\Phi }(z,n,1)`$, $`b=2\pi l/a`$. For $`\delta =0`$ this formula gives the total elastic energy of a symmetric domain structure (a pattern of opposite domains of equal widths). We have found earlier a somewhat similar formula for another symmetric domain structure of $`a_1|a_2`$ type. Let us consider the following general cases.

Zero extrinsic misfit $`\left(w=0\right).`$ In this case the energy of the homogeneous elastic field is simply $`F_h/𝒜=Mu_0^2l\delta ^2.`$ The stray energy has a simple form for the following two limiting cases. In the standard case of narrow domains $`\left(al\right)`$ and $`\delta 1`$

$$\frac{F_{stray}}{𝒜}=\frac{2}{\pi ^3}Mu_0^2a\left[\frac{7}{4}\zeta \left(3\right)\frac{\pi ^2\mathrm{ln}2}{2}\delta ^2+\frac{\pi ^4}{96}\delta ^4\right].$$
(11)
One has to add the surface energy of the domain walls to find the equilibrium domain width

$$F_{dw}/𝒜=\sqrt{2}\gamma l/a,$$
(12)
where $`\gamma Mu_0^2\mathrm{\Delta }/\sqrt{2}`$ is the energy of the unit surface of the domain walls, and $`\mathrm{\Delta }`$ is a characteristic microscopic length . The total free energy, $`F=F_h+F_{stray}+F_{dw},`$ is minimal for symmetric domain structure ($`\delta =0`$ and $`a_1=a_2=a)`$, with the standard (Kittel) domain width (cf.)

$$a=a_K\left(\frac{2\pi ^3}{7\zeta (3)}l\mathrm{\Delta }\right)^{1/2}\sqrt{l\mathrm{\Delta }}l.$$
(13)
Using this result, we see immediately that the high-frequency response of the domain pattern (at a fixed domain width $`a)`$ is

$$\kappa _{\mathrm{}}=\frac{^2}{\delta ^2}\left(\frac{F}{𝒜}\right)=2Mu_0^2l\left(1\frac{\mathrm{ln}2}{\pi }\frac{a}{l}+\frac{\pi a}{8l}\delta ^2\right)>0,$$
(14)
so the pattern is stable, although it softens because of the negative contribution to the stray energy from the terms $`\delta ^2`$ in (11). Note also that $`\kappa _{\mathrm{}}`$ increases with $`\delta ,`$ and this, as we shall see, is a general result.

If the elastic modulus $`\mu `$ were soft, the domains can become wide, $`al`$ $`\left(b1\right)`$. There

$$\frac{F_{stray}}{𝒜}=Mu_0^2l\left[1\delta ^2\left(1\frac{\pi l}{a}\right)\frac{8l}{\pi a}\mathrm{ln}\frac{e^{3/4}a}{4\pi l}\right],$$
(15)
and one finds a symmetric structure with the large domain width

$$a=4\pi e^{1/4}l\mathrm{exp}\left(\pi \mathrm{\Delta }/8l\right),$$
(16)
which would be $`l`$ if we had a substrate with $`\mathrm{\Delta }d_{at}`$ (cf. the answer for a ferroelectric capacitor ). This is contrary to a previously considered Kittel case with narrow domains $`a_Kl.`$ In the case of the wide domains the response softens considerably, but remains positive,

$$\kappa _{\mathrm{}}=2\pi Mu_0^2l\frac{l}{a}>0.$$
(17)
It is unlikely, however, that the softness of the present domain structure is of any practical importance. The problem is that even a small extrinsic misfit in very thin films will push the domain structure to the boundary of its existence , where it becomes rigid for experimentally accessible frequencies of external field, as discussed below.

Non-zero extrinsic misfit $`\left(w0\right).`$ The energy of homogeneous stresses is given by Eq. (6). To find the domain structure close to the boundary of its stability we need all characteristics in the limit $`\delta 1,`$ and $`al`$ ($`b1),`$ i.e. close to a monodomain state. The stray energy in this limit is

$$\frac{F_{stray}}{𝒜}=\frac{1}{4}Mu_0^2\left(1\delta ^2\right)^2\left[\mathrm{ln}\frac{4el}{a\left(1\delta ^2\right)}\frac{\pi l}{a}\right].$$
(18)
Note that Roytburd has approximated the numerically computed stray energy with the functional dependence $`\left(1\delta ^2\right)^2,`$ i.e. without the important log term. Minimizing the total energy, $`F_{tot}=F_h+F_{stray}+F_{dw},`$ with respect to the half-period of the structure $`a`$, we obtain the equation

$$\frac{\left(1\delta ^2\right)^2}{4\pi }\mathrm{ln}\frac{4l}{a\left(1\delta ^2\right)}=\frac{l\mathrm{\Delta }}{a^2},$$
(19)
which has a solution

$$\frac{a}{l}=\frac{1}{\left(1\delta ^2\right)\mathrm{ln}^{1/2}\left(4l/\pi \mathrm{\Delta }\right)}\sqrt{\frac{8\pi \mathrm{\Delta }}{l}},$$
(20)
where we have omitted the terms $`\mathrm{ln}\mathrm{ln}\left(4l/\pi \mathrm{\Delta }\right)\mathrm{ln}\left(4l/\pi \mathrm{\Delta }\right).`$ It shows that the domain period diverges as $`1/(1\delta ^2)`$ close to the phase boundary with the monodomain state. The total energy then takes the form

$$\frac{F_{tot}}{Mu_0^2l𝒜}=(\delta \delta _0)^2+\varphi \left(1\delta ^2\right)\frac{1}{4}\left(1\delta ^2\right)^2,$$
(21)
where from the equilibrium value of $`\delta `$ is to be found (subject to the constraint $`\delta ^21)`$, with the parameter $`\varphi \left(\mathrm{\Delta }/2\pi l\right)^{1/2}1,`$ and $`=\mathrm{ln}\left(4l/\pi \mathrm{\Delta }\right)`$ the logarithm of a large number $`l/\mathrm{\Delta }1.`$

The transition into a monodomain state $`\left(\delta ^2=1\right)`$ occurs as a function of the extrinsic relative misfit at $`\left|\delta _0\right|=1\varphi .`$ This condition suggests a critical thickness of the film $`l_c,`$ where the polydomain-monodomain transition takes place at a given misfit $`\delta _0`$ (i.e. the phase boundary in $`\delta _0l`$ plane),

$$\frac{l_c}{\mathrm{\Delta }}=\frac{1}{\pi \left(1\delta _0\right)^2}\mathrm{ln}\frac{1}{1\delta _0}.$$
(22)
Thus, the domain structure may exist only at $`l>l_c\left(\delta _0\right).`$ This formula is obtained for $`l_c/\mathrm{\Delta }1`$, but it also correctly gives an absence of the critical thickness of the film when the parameter $`\delta _0`$ is zero ($`l_c=0`$ when $`\delta _0=0).`$ The absence of the critical film thickness has been suggested earlier for another domain pattern, $`a_1|a_2,`$ from numerical computations. This is contrary to the speculations by Roytburd, who obtained $`l_c0`$ for $`\delta _0=0,`$ apparently as an artifact of the employed approximations. The critical point is approached linearly with $`\delta _0`$

$$\delta =1\left(1\varphi \delta _0\right)/\varphi .$$
(23)
We see that the slope is $`d\delta /d\delta _0=1/\varphi 1,`$ so that the approach to the critical point is very steep, it looks almost discontinuous.

We readily obtain the high-frequency (measurable) dynamical stiffness close to the phase boundary,

$`{\displaystyle \frac{\kappa _{\mathrm{}}}{Mu_0^2l}}`$ $`=`$ $`{\displaystyle \frac{1}{1\delta ^2}}\sqrt{{\displaystyle \frac{8\mathrm{\Delta }}{\pi l}}}\mathrm{ln}{\displaystyle \frac{2e^2l}{\pi \mathrm{\Delta }}}`$ (24)
$``$ $`{\displaystyle \frac{1}{1\delta ^2}}\sqrt{{\displaystyle \frac{8\mathrm{\Delta }}{\pi l}}}\mathrm{}(!).`$ (25)
which diverges when $`\delta ^21,`$ i.e. close to the phase boundary. This is in striking disagreement with the results of Pertsev et al. who claimed a softness of the $`c/a`$ domain structure close to the phase boundary, which actually does not materialize. The static stiffness, which has no practical significance since it requires very long time to relax the pattern in order to optimize the number of domain walls, sharply vanishes with increasing relative misfit $`\delta _0`$,

$$\frac{\kappa _0}{Mu_0^2l}1\delta ^22\left(1\varphi \delta _0\right)/\varphi 0,$$
(26)
see the exact result in Fig. 2. Interestingly, close to the transition into monodomain state the system splits into two groups of wide ($`a_1)`$ and narrow ($`a_2)`$ domains with

$`a_1`$ $`=`$ $`{\displaystyle \frac{1}{1\delta }}\sqrt{{\displaystyle \frac{8\pi l\mathrm{\Delta }}{}}}{\displaystyle \frac{a_K}{(1\delta )\sqrt{}}}\mathrm{},`$ (27)
$`a_2`$ $`=`$ $`\sqrt{{\displaystyle \frac{2\pi l\mathrm{\Delta }}{}}}{\displaystyle \frac{a_K}{\sqrt{}}}<a_{K.}`$ (28)
We see that the width of the wide domains diverges, the density of the domains walls decreases, but the width of the narrow domains, $`a_2,`$ remains small. The narrow domains are somewhat compressed compared to the standard Kittel width $`a_K`$ (13).

It becomes obvious from the present analysis that it is unlikely that one can succeed in making a “soft” domain pattern with the small $`\kappa _{\mathrm{}}`$ in epitaxial thin ferroelectric films with elastic domains, suitable for the applications. In fact, the only way to do this would be to avoid any extrinsic misfit, i.e. to keep $`\delta _0=0,`$ and to reduce the film thickness. However, in this case the interval of the extrinsic misfit strain, which allows the very existence of the domains, is very narrow. Then, the proximity to the phase boundary in systems with two kinds of inequivalent domains means that the stiffness rapidly increases, oppositely to what one actually desired. We have tacitly assumed above that the stiffness of the domain structure determines its dielectric response. Indeed, one can assume that there is a spontaneous polarization parallel to the film plane in one domain and perpendicular to it in the other in the present $`c/a`$ domain structure (as in $`c`$-domains in perovskites). One can also assume, as in Ref. , that $`\delta \delta _{eq}=lP_sE_{ext}/\kappa _{\mathrm{}}`$ in external field, where $`P_s`$ is the spontaneous polarization, $`E_{ext}`$ the external electric field, and $`\delta _{eq}`$ is the equilibrium value of $`\delta `$ in zero field. We have shown earlier  that this relation may not hold, and in some cases it strongly overestimates the dielectric response. Additionally, the neglected contribution of uncompensated electric fields would only make the response stiffer. However, even this overestimate leads one to a conclusion that the dielectric response of the pattern with inequivalent elastic domains in the epitaxial thin films is actually suppressed.


format: text

text:
# Remarks on the Bergmann–Thomson Expression on Angular Momentum in General Relativity

## I Bergmann–Thomson expression on angular momentum in GR

One can most easily obtain the canonical Bergmann–Thomson expression on total angular momentum density, matter and gravitation, in the following way. At first, let us transform the Einstein equations written in mixed form and multiplied by $`\sqrt{|g|}`$

$$\sqrt{|g|}G_i^k=\beta \sqrt{|g|}T_i^k$$
(1)
to the so–called superpotential form

$${}_{E}{}^{}K_{i}^{k}={}_{F}{}^{}U_{i}^{[kl]}{}_{,l}{}^{},$$
(2)
where

$${}_{E}{}^{}K_{i}^{k}:=\sqrt{|g|}\left(T_i^k+_Et_i^k\right)$$
(3)
is the canonical, Einstein’s energy–momentum complex, matter and gravitation, and $`{}_{F}{}^{}U_{i}^{[kl]}`$ mean von Freud superpotentials.

From (2) we get, after series operations, the Bergmann–Thomson energy–momentum complex $`{}_{BT}{}^{}K_{}^{jk}`$, matter and gravitation. At first, we form

$$g^{ij}{}_{E}{}^{}K_{i}^{k}=g^{ij}{}_{F}{}^{}U_{i}^{[kl]}{}_{,l}{}^{}=\left(g^{ij}{}_{F}{}^{}U_{i}^{[kl]}\right)_{,l}{}_{F}{}^{}U_{i}^{[kl]}g_{,l}^{ij}$$
(4)
or

$$g^{ij}{}_{E}{}^{}K_{i}^{k}+{}_{F}{}^{}U_{i}^{[kl]}g_{,l}^{ij}={}_{F}{}^{}U_{}^{j[kl]}{}_{,l}{}^{}.$$
(5)
Then, we write (5) in the form

$${}_{BT}{}^{}K_{}^{jk}={}_{F}{}^{}U_{}^{j[kl]}{}_{,l}{}^{},$$
(6)
where

$${}_{BT}{}^{}K_{}^{jk}:=_EK^{jk}+{}_{F}{}^{}U_{i}^{[kl]}g_{,l}^{ij}=:\sqrt{|g|}\left(T^{jk}+_{BT}t^{jk}\right)$$
(7)
is the Bergmann–Thomson energy–momentum complex  of matter and gravitation which satisfies local conservation laws

$${}_{BT}{}^{}K_{}^{jk}{}_{,k}{}^{}=0.$$
(8)
Here $`{}_{BT}{}^{}t_{}^{jk}_{BT}t^{kj}`$ mean the components of the so–called Bergmann–Thomson energy–momentum pseudotensor of the gravitational field . Finally, from (6) we get

$`x^i{}_{BT}{}^{}K_{}^{jk}x^j{}_{BT}{}^{}K_{}^{ik}`$ $`=`$ $`x^i{}_{F}{}^{}U_{}^{j[kl]}{}_{,l}{}^{}x^j{}_{F}{}^{}U_{}^{i[kl]}{}_{,l}{}^{}`$ (9)
$`=`$ $`\left(x^i{}_{F}{}^{}U_{}^{j[kl]}x^j{}_{F}{}^{}U_{}^{j[kl]}\right)_{,l}_FU^{j[ki]}+_FU^{i[kj]},`$ (10)
or

$$x^i{}_{BT}{}^{}K_{}^{jk}x^j{}_{BT}{}^{}K_{}^{ik}+_{BT}S^{ijk}=M_{,l}^{[ij][kl]},$$
(11)
where

$`{}_{BT}{}^{}S_{}^{ijk}`$ $`:=`$ $`{}_{F}{}^{}U_{}^{i[jk]}_FU^{j[ik]}={\displaystyle \frac{\alpha }{\sqrt{|g|}}}\left[(g)\left(g^{kj}g^{il}g^{ki}g^{jl}\right)\right]_{,l}`$ (12)
$`=:`$ $`{\displaystyle \frac{\alpha }{\sqrt{|g|}}}g_{,l}^{[ij][kl]}`$ (13)
and

$$M^{[ij][kl]}:=x^i{}_{F}{}^{}U_{}^{j[kl]}x^j{}_{F}{}^{}U_{}^{i[kl]}.$$
(14)
The expression

$$x^i{}_{BT}{}^{}K_{}^{jk}x^j{}_{BT}{}^{}K_{}^{ik}+_{BT}S^{ijk}=:_{BT}M^{ijk}$$
(15)
is the Bergmann–Thomson angular momentum complex matter and gravitation and the quantities $`M^{[ij][kl]}`$ are superpotentials .

The complex $`{}_{BT}{}^{}M_{}^{ijk}`$ satisfies local conservation laws

$${}_{BT}{}^{}M_{,k}^{ijk}=0.$$
(16)
One can interpret physically the angular momentum complex (13) as a sum of the orbital part

$$O^{ijk}:=x^i{}_{BT}{}^{}K_{}^{jk}x^j{}_{BT}{}^{}K_{}^{ik}=\sqrt{|g|}\left(x^iT^{jk}x^jT^{ik}\right)+\sqrt{|g|}\left(x^i{}_{BT}{}^{}t_{}^{jk}x^j{}_{BT}{}^{}T_{}^{ik}\right)$$
(17)
of the angular momentum density of matter and gravitation (The matter part includes also spin density ) and a spinorial part

$${}_{BT}{}^{}S_{}^{ijk}=_FU^{i[jk]}_FU^{j[ik]}=\frac{\alpha }{\sqrt{|g|}}g_{,l}^{[ij][kl]}$$
(18)
of the gravitational angular momentum density.

We have from (6)

$$2_{BT}K^{[ij]}=\sqrt{|g|}(_{BT}t^{ij}_{BT}t^{ji})={}_{BT}{}^{}S_{}^{ijk}{}_{,k}{}^{}=\left(\frac{\alpha }{\sqrt{|g|}}\right)_{,k}g_{,l}^{[ij][kl]},$$
(19)
because the dynamical energy–momentum tensor of matter $`T^{ik}`$ is symmetric: $`T^{ik}=T^{ki}`$.

The equality (17) justifies the above proposal of the physical interpretation of the pseudotensor $`{}_{BT}{}^{}S_{}^{ijk}=()_{BT}S^{jik}`$ as a quantity describing canonical spin density of the gravitational field ***In Special Relativity the antisymmetric part of an energy–momentum tensor is proportional to the ordinary divergence of a quantity which describes canonical spin densities . Here we follow this line..

It is very interesting that formal application of the special-relativistic symmetrization procedure given by Belinfante (see, e.g., ) to the complex $`{}_{BT}{}^{}K_{}^{ij}`$ with $`S^{ijk}=_{BT}S^{ijk}`$ leads us immediately to the new, symmetric complex $`{}_{BT}{}^{}K_{}^{Nij}=\sqrt{|g|}\left(T^{ij}\frac{c^4}{8\pi G}G^{ij}\right)`$ which trivially vanishes, i.e., it leads us to the Lorentz and Levi-Civita solution of the energy-momentum problem for gravitational field: $`{}_{g}{}^{}T_{}^{ik}=()\frac{c^4}{8\pi G}G^{ik}`$.

From (10) we get the following expression on the components $`M^{ik}=()M^{ki}`$ of the global angular momentum of matter and gravitation for an isolated system endowed with an asymtotically Lorentzian coordinates The Bergmann–Thomson expression (10), likely as the canonical Einstein energy–momentum complex, can be reasonably use only in the case of an isolated system endowed with an asymptotically flat coordinates (see, e.g., ).

$${}_{BT}{}^{}M_{}^{ik}=\frac{1}{c}\underset{overspherehavingR\mathrm{}}{}\left(x^i{}_{F}{}^{}U_{}^{j[0\alpha ]}x^j{}_{F}{}^{}U_{}^{i[0\alpha ]}\right)n_\alpha r^2𝑑\mathrm{\Omega }.$$
(20)
$`r^2=x^2+y^2+z^2`$; $`n_\alpha `$ are the components of the unit (exterior) normal to the sphere and $`d\mathrm{\Omega }=\mathrm{sin}\theta d\theta d\phi `$.

For the Schwarzschild spacetime equipped with an asymptotically Lorentzian coordinates $`(x^0=ct,x^1=x,x^2=y,x^3=z)`$ we get from (18) the expected result

$${}_{BT}{}^{}M_{}^{ik}=0.$$
(21)
On the other hand, in the case of the stationary and axially symmetric Kerr’s spacetime (see, e.g., ) we obtain, in an asymptotically Lorentzian coordinates also, that the only one component

$${}_{BT}{}^{}M_{}^{12}=()_{BT}M^{21}=mac\frac{1}{3}mac=\frac{2}{3}mac$$
(22)
of the $`{}_{BT}{}^{}M_{}^{ik}`$ is different from zero.

One should interpret the result (20) as giving the global angular momentum, spinorial and orbital of matter and gravitation, for Kerr spacetime.

## II Landau–Lifschitz expression on angular momentum density in General Relativity and its comparison with Bergmann–Thomson expression

One can easily obtain the Landau–Lifschitz expression on angular momentum density for matter and gravitation in the following way . At first, one should transform the Einstein equations multiplied by $`(g)`$

$$(g)G^{ik}=\beta (g)T^{ik}$$
(23)
to the superpotential form

$$(g)\left(T^{ik}+_{LL}t^{ik}\right)=h_{,l}^{ikl},$$
(24)
where

$$h^{ikl}=()h^{ilk}=\lambda _{,m}^{iklm}$$
(25)
and

$$\lambda ^{iklm}:=\alpha (g)\left(g^{ik}g^{lm}g^{il}g^{km}\right).$$
(26)
$`\alpha =\frac{1}{2\beta }=\frac{c^4}{16\pi G}`$.

$`{}_{LL}{}^{}t_{}^{ik}=_{LL}t^{ki}`$ are components of the so–called Landau–Lifschitz gravitational energy–momentum pseudotensor and the sum

$$(g)\left(T^{ik}+_{LL}t^{ik}\right)=:_{LL}K^{ik}$$
(27)
forms the so–called Landau–Lifschitz energy–momentum complex The complex (25) has worse transformational properties than the Einstein canonical energy–momentum complex . of matter and gravitation. $`h^{ikl}=()h^{ilk}`$ form the Landau–Lifschitz superpotentials.

Then, from the (22)–(23), one can obtain easily

$`x^i{}_{LL}{}^{}K_{}^{kl}x^k{}_{LL}{}^{}K_{}^{il}`$ $`=`$ $`x^ih_{,m}^{klm}x^kh_{,m}^{ilm}`$ (28)
$`=`$ $`\left(x^ih^{klm}x^kh^{ilm}\right)_{,m}h^{kli}+h^{ilk},`$ (29)
or

$$x^i{}_{LL}{}^{}K_{}^{kl}x^k{}_{LL}{}^{}K_{}^{il}+_{LL}S^{ikl}=L_{,m}^{[ik][lm]},$$
(30)
where

$${}_{LL}{}^{}S_{}^{ikl}:=h^{ikl}h^{kil}0$$
(31)
and

$$L^{[ik][lm]}:=x^ih^{klm}x^kh^{ilm}.$$
(32)
The expression (27) exactly corresponds to the Bergmann–Thomson expression (10). Namely, $`x^i{}_{LL}{}^{}K_{}^{kl}x^k{}_{LL}{}^{}K_{}^{il}`$ gives orbital part of the angular momentum density of matter and gravitation,§§§Material part $`(g)\left(x^iT^{kl}x^kT^{il}\right)`$ includes material spin also. $`{}_{LL}{}^{}S_{}^{ikl}=h^{ikl}h^{kil}`$ gives the non–vanishing spinorial part of the gravitational angular momentum density and $`L^{[ik][lm]}=x^ih^{klm}x^kh^{ilm}`$ are superpotentials for the total angular momentum density of matter and gravitation.

The total angular momentum density

$${}_{LL}{}^{}M_{}^{ikl}:=x_{LL}^iK^{kl}x_{LL}^kK^{il}+_{LL}S^{ikl}$$
(33)
satisfies local conservation laws

$${}_{LL}{}^{}M_{,l}^{ikl}=0.$$
(34)
From (27) one can obtain the following integral expression on the global angular momentum of an isolated system endowed with an asymptotically Lorentzian coordinates $`(x^0=ct,x^1=x,x^2=y,x^3=z)`$

$${}_{LL}{}^{}M_{}^{ik}=\frac{1}{c}\underset{oversphereS^2havingradiusgoingto\mathrm{}}{}L^{[ik][0\alpha ]}n_\alpha r^2𝑑\mathrm{\Omega }.$$
(35)
This expression gives the same values of the components $`M^{ik}=()M^{ki}`$ as the Bergmann–Thomson expression (18) gives. Especially, for the Kerr’s spacetime it gives

$${}_{LL}{}^{}M_{}^{12}=()_{LL}M^{21}=mac\frac{1}{3}mac=\frac{2}{3}mac;$$
(36)
other components vanish.

However, in the case we have

$${}_{LL}{}^{}S_{}^{ikl}=\lambda _{,m}^{klmi}$$
(37)
from which it follows

$${}_{LL}{}^{}S_{}^{ikl}{}_{,l}{}^{}=0.$$
(38)
The last equality guarantees symmetry of the Landau–Lifschitz energy–momentum complex:

$${}_{LL}{}^{}K_{}^{ik}_{LL}K^{ki}=S_{,l}^{ikl}=0.$$
(39)
Landau–Lifschitz, in their book , modify the expression (32) on global angular momentum of an isolated system equipped with an asymptotically Lorentzian coordinates by using (34). Namely, they subtract the gravitational spin angular density $`{}_{LL}{}^{}S_{}^{ikl}=()_{LL}S^{kil}`$ from the both sides of the (27) and obtain the new expression on angular momentum density of matter and gravitation

$$x^i{}_{LL}{}^{}K_{}^{kl}x^k{}_{LL}{}^{}K_{}^{il}=\left(L^{[ik][lm]}\lambda ^{klmi}\right)_{,m}.$$
(40)
This new expression includes the total matter angular momentum density, spinorial and orbital, but its gravitational part consists of the orbital gravitational angular momentum density only.

The new expression (37) leads to the following integrals on global quantities of an isolated system equipped with an asymptotically Lorentzian coordinates $`(x^0=ct,x^1=x,x^2=y,x^3=z)`$

$$M^{ik}=\frac{1}{c}\underset{oversphereS^2havingradiusR\mathrm{}}{}\left(L^{[ik][0\alpha ]}+\lambda ^{i0\alpha k}\right)n_\alpha r^2𝑑\mathrm{\Omega }.$$
(41)
In the case of the Schwarzschild spacetime equipped with an asymptotically Lorentzian coordinates $`(x^0=ct,x,y,z)`$ the expression (38) gives $`M^{ik}=()M^{ki}=0`$ and in the case of the Kerr’s spacetime it gives

$$M^{12}=()M^{21}=mac\frac{1}{3}mac+\frac{1}{3}mac=mac;$$
(42)
other components of the $`M^{ik}`$ vanish.

An analogical modification one can do with the Bergmann–Thomson expression (18). Namely, by using (11), one can rewrite (10) in the following form

$`x^i{}_{BT}{}^{}K_{}^{jk}x^j{}_{BT}{}^{}K_{}^{ik}`$ $``$ $`\left({\displaystyle \frac{\alpha }{\sqrt{|g|}}}\right)_{,l}g^{[ij][kl]}+\left({\displaystyle \frac{\alpha g^{[ij][kl]}}{\sqrt{|g|}}}\right)_{,l}`$ (43)
$`=`$ $`M_{,l}^{[ij][kl]}.`$ (44)
By moving the term $`\left(\frac{\alpha g^{[ij][kl]}}{\sqrt{|g|}}\right)_{,l}`$ onto the right hand side of the (40) one gets a new form of the Bergmann–Thomson expression (10) with the new angular momentum complex $`{}_{BT}{}^{}\overline{M}_{}^{ikl}=()_{BT}\overline{M}^{kil}`$ and with the new superpotentials $`\overline{M}^{[ij][kl]}`$.

Namely, we get

$`x^i{}_{BT}{}^{}K_{}^{jk}x^j{}_{BT}{}^{}K_{}^{ik}`$ $``$ $`\left({\displaystyle \frac{\alpha }{\sqrt{|g|}}}\right)_{,l}g^{[ij][kl]}`$ (45)
$`=`$ $`(M^{[ij][kl]}{\displaystyle \frac{\alpha }{\sqrt{|g|}}}g^{[ij][kl]})_{,l}=:\overline{M}_{,l}^{[ij][kl]}.`$ (46)
The last expression leads us to the following integrals on global quantities $`\overline{M}^{ik}=()\overline{M}^{ki}`$ in the case of an isolated system equipped with an asymptotically Lorentzian coordinates $`(x^0=ct,x,y,z)`$

$`\overline{M}^{ik}`$ $`=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle \underset{oversphereS^2havingradiusR\mathrm{}}{}}\left(M^{[ij][0\alpha ]}\alpha {\displaystyle \frac{g^{[ij][0\alpha ]}}{\sqrt{|g|}}}\right)n_\alpha r^2𝑑\mathrm{\Omega }`$ (47)
$`:=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle \underset{oversphereS^2havingradiusR\mathrm{}}{}}\overline{M}^{[ij][0\alpha ]}n_\alpha r^2𝑑\mathrm{\Omega }.`$ (48)
The integrals (42) correspond to the Landau–Lifschitz integrals (38) and, especially, for the Kerr’s spacetime they give the same values as the Landau–Lifschitz integrals (38) give, i.e., they give

$$\overline{M}^{12}=()\overline{M}^{21}=mac\frac{1}{3}mac+\frac{1}{3}mac=mac;$$
(49)
other components vanish.

So, the Bergmann–Thomson and Landau–Lifschitz integral expressions on global angular momentum are fully equivalent in the case of an isolated system equipped with an asymptotically Lorentzian coordinates. But we favorize the Bergmann–Thomson expression on angular momentum density because it has better local transformational properties in comparison with the Landau–Lifschitz expression and it is more closely related to the canonical energy–momentum complex $`{}_{E}{}^{}K_{i}^{k}=\sqrt{|g|}\left(T_i^k+_Et_i^k\right)`$.

In the case of the Kerr’s spacetime one should interpret the value $`M^{12}=()M^{21}=\overline{M}^{12}=()\overline{M}^{21}=mac`$ of the integrals (38) or (42) as referring only to the material part of the angular momentum complex. It is justified by the following fact : the integral including the energy–momentum tensor of matter $`T^{ik}=T^{ki}`$ only

$`M_2^1`$ $`=`$ $`()M^{12}={\displaystyle \frac{1}{c}}{\displaystyle \underset{x^0=ct=const}{}}T_i^0\xi ^i\sqrt{|g|}𝑑x𝑑y𝑑z`$ (50)
$`=`$ $`{\displaystyle \frac{c^3}{16\pi G}}{\displaystyle \underset{oversphereS^2havingradiusR\mathrm{}}{}}\left(^0\xi ^\alpha ^\alpha \xi ^0\right)\sqrt{|g|}n_\alpha r^2𝑑\mathrm{\Omega },`$ (51)
where $`\xi ^i=x\delta _2^iy\delta _1^i`$ are the components of the spatial Killing vector field $`\stackrel{}{\xi }=\xi ^i_i`$ which exists in the case, gives

$$M_2^1=()M^{12}=()mac.$$
(52)
So, the matter integral (44) has the same value in the case as the integrals (38) and (42). In consequence, we conclude that only the material parts $`()g\left(x^iT^{kl}x^kT^{il}\right),\sqrt{|g|}\left(x^iT^{kl}x^kT^{il}\right)`$ of the suitable angular momentum complex (Landau-Lifschitz or Bergmann-Thomson) give contribution to the integrals (38) and (42). Gravitational parts give no contribution in the used coordinates.


text:
# Resonance Effects in Topological Discrete sine-Gordon System

## 1. Introduction

Nonlinear systems which have soliton solutions have played an increasingly important role in physics. One such area in physics is condensed matter theory where the soliton has been used to model crystal discolations, charge density waves, magnetic domain walls, etc . In many cases it has been necessary to model interactions of the soliton with other objects such as impurities, defects, phonons, antisolitons, or other solitons. In (1+1)-dimensional field theories, these interactions have sometimes been quite spectacular, particularly if the 1-dimensional soliton, or the kink, supports extra modes of vibration. Such a mode can be excited during an interaction and this leads to a “resonance effect” . This effect has been seen in both continuous and discrete systems \[1-8\]. These resonance effects give rise to “windows” in the space of impact velocity, i.e. there are certain intervals of velocities in this space for which the kinks form a bound state (a bion) and others for which the kinks scatter off each other.

This paper deals with kink-antikink collisions in the topological discrete sine-Gordon system (TDSG) introduced by Speight and Ward in 1994. This system is particularly interesting since it is a discrete system but with many of the features of the continuum sine-Gordon. In particular, it maintains the Bogomoln’yi bound. The resulting Bogomoln’yi equation is a first-order difference equation which can be easily solved. The solution is analogous to the well-known continuum sine-Gordon static kink solution. Moreover, as we show in section 4, the TDSG kink supports an internal shape mode.

The paper is divided into 5 sections. Section 2 contains a brief summary of kink dynamics in the TDSG system , . In section 3 we consider the kink-antikink collisions in the TDSG. In section 4 we show that TDSG kink supports an internal mode of vibration and section 5 contains some concluding remarks.

## 2. Kink dynamics in TDSG

The TDSG system is defined by the Lagrangian density

$$=\frac{h}{4}\underset{nZ}{}(\dot{\psi _n}^2(D_n^2+F_n^2)),$$
(1)
where the kinetic energy is $`E_K=\frac{h}{4}\dot{\psi }_n^2`$ and the potential energy is $`E_P=\frac{h}{4}(D_n^2+F_n^2)`$.

The quantities $`D_n`$ and $`F_n`$ are given by

$`D_n={\displaystyle \frac{2}{h}}\mathrm{sin}{\displaystyle \frac{1}{2}}(\psi _{n+1}\psi _n)\text{and}F_n=\mathrm{sin}{\displaystyle \frac{1}{2}}(\psi _{n+1}+\psi _n),`$ (2)
giving

$$E_P=\frac{h}{4}\frac{4}{h^2}\mathrm{sin}^2\frac{1}{2}(\psi _{n+1}\psi _n)+\mathrm{sin}^2\frac{1}{2}(\psi _{n+1}+\psi _n).$$
(3)
$`\psi _n=\psi (x)`$ at the $`n^{th}`$ lattice site and $`\psi _{n+1}=\psi (x+h)`$, where $`h`$ is the lattice spacing. The first term in equation (3) represents an attractive force analogous to the Hooke force. The second term is the substrate potential, found by taking an average of the potential of the two nearest neighbours.

In the continuum limit, as $`h0`$, $`D_n\psi _x`$, and $`F_n\mathrm{sin}\psi `$, and one recovers the standard expression of the continuum potential

$$E_P^{cont}(\psi )=\frac{1}{4}_{\mathrm{}}^{\mathrm{}}(\psi _x^2+\mathrm{sin}^2\psi )𝑑x.$$
(4)
The key feature of this model is the choice of the derivative on the lattice. This choice is motivated by the Bogomol’nyi argument of the continuum system. Rather than choose the standard forward difference, Speight and Ward choose,

$$\frac{\mathrm{cos}\psi _{n+1}\mathrm{cos}\psi _n}{h}=D_nF_n,$$
(5)
with the kink boundary conditions

$$\underset{n\mathrm{}}{lim}\psi _n=0\text{and}\underset{n\mathrm{}}{lim}\psi _n=\pi .$$

This leads to the following discrete Bogomol’nyi argument

$`0`$ $`{\displaystyle \frac{h}{4}}{\displaystyle \underset{nZ}{}}(D_nF_n)^2`$ (6)
$`=E_P+{\displaystyle \frac{1}{2}}{\displaystyle \underset{nZ}{}}(\mathrm{cos}\psi _{n+1}\mathrm{cos}\psi _n)`$ (7)
$`=E_P1,`$ (8)
so $`E_P1`$, with equality if $`D_n=F_n`$, i.e. if

$`{\displaystyle \frac{2}{h}}\mathrm{sin}{\displaystyle \frac{1}{2}}(\psi _{n+1}\psi _n)`$ $`=\mathrm{sin}{\displaystyle \frac{1}{2}}(\psi _{n+1}+\psi _n)`$ (9)
$`\text{or,}\psi _n`$ $`=2\mathrm{tan}^1\mathrm{exp}a(nhb).`$ (10)
The constant $`a=h^1\mathrm{ln}[\frac{2+h}{2h}]`$ is the kink’s slope and $`bR`$ is the position. Equation (10) is the static kink solution. So the moduli space of static solutions for the TDSG is isomorphic to $``$ rather than $``$, which is what one usually expects for a discrete system. Moreover, $`E_P`$ is independent of $`b`$ suggesting that the model has no Peierls-Nabarro barrier. This means that the kink can move along the lattice arbitrarily slowly without getting pinned.

The Euler-Lagrange equation gives the following equation of motion,

$$\ddot{\psi _n}=\frac{4h^2}{4h^2}\mathrm{cos}\psi _n(\mathrm{sin}\psi _{n+1}+\mathrm{sin}\psi _{n1})\frac{4+h^2}{4h^2}\mathrm{sin}\psi _n(\mathrm{cos}\psi _{n+1}+\mathrm{cos}\psi _{n1}).$$
(11)
Equation (11) can then be used to study the behaviour of a moving kink. Speight and Ward have solved (11) using a 4th-order Runge-Kutta scheme. Their initial condition was a static kink Galilean boosted to a velocity $`v`$. Their results show that for small initial velocities the kink wobbles with a period $`h/v`$ as it moves along the lattice. This is a purely dynamical effect since there is no Peierls-Nabarro barrier.

For fast moving kinks it is found that the kink radiates as it moves along the lattice. This is observed through a gradual decrease in the kink velocity. But it is found that there is certain velocity threshold below which the radiation from the kink is significantly reduced.

The existence of such a velocity threshold can be understood by considering the linearised equation of motion .

$$\ddot{\psi _n}=\frac{4h^2}{4h^2}(\psi _{n+1}+\psi _{n1})\frac{4+h^2}{2h^2}\psi _n.$$
(12)
Equation (12) can be used to derive a dispersion relation for small amplitude travelling waves, given by

$$\omega ^2=\frac{4+h^2}{2h^2}\frac{4h^2}{2h^2}\mathrm{cos}kh,$$
(13)
where $`k`$ is the wavenumber. Since $`0<h<2`$, we have $`1<\omega <2/h`$. The case $`\omega =1`$ corresponds to the threshold velocity. The frequency at which the kink hits the lattice sites is $`v/h`$ per unit time. So provided $`vh/2\pi `$, the kink will radiate (at $`\omega =2\pi v/h`$). But if $`v<h/2\pi `$ then $`\omega <1`$ and the kink does not radiate. Hence there is a ”preferred” velocity for the kink. Since this velocity depends only on $`h`$, it can be regarded as a feature of the lattice.

The choice of the functions $`D_n`$ and $`F_n`$ is not unique. For instance, $`D_n`$ could be multiplied by a function $`f(x,h)`$ and $`F_n`$ by $`f(x,h)^1`$ and the product is independent of $`f`$ so one would still have a topologically stable kink solution. The only condition on $`f(x,h)`$ is that $`lim_{h0}f=1`$. Similarly there is freedom in the choice of $`\ddot{\psi }`$. Zakrzewski  has investigated the effects on the kink motion when this choice is exploited.

The $`\ddot{\psi }`$ in (11) is replaced by

$$\ddot{\psi }\frac{\ddot{\psi }}{1+\alpha \dot{\psi }^2}$$
(14)
where $`\alpha =\alpha (h)`$ and $`\alpha (0)=1`$. The equation of motion (11) is now changed to

$$\frac{\ddot{\psi }_n}{(1+\alpha \dot{\psi _n}^2)}=\frac{f^2}{h^2}(\mathrm{sin}(\psi _{n+1}\psi _n)\mathrm{sin}(\psi _n\psi _{n1}))\frac{1}{4f^2}(\mathrm{sin}(\psi _{n+1}+\psi _n)+\mathrm{sin}(\psi _n+\psi _{n1})).$$
(15)
The choices of $`f`$ and $`\alpha `$ are motivated by the fact that equation (15) still needs to admit the Galilean boosted sine-Gordon kink solution. One choice is the following

$$f^2=\frac{\frac{v_1^2}{1v_1^2}+\sqrt{\left(\frac{v_1^2}{1v_1^2}\right)^2+\frac{4\mathrm{sinh}^22\beta }{h^2}}}{\frac{8\mathrm{sinh}^2\beta }{h^2}}$$
(16)
where $`\beta =h/2\sqrt{1v_1^2}`$, and

$$\alpha (h)=\frac{1v_1^2}{v_1^2}\mathrm{sinh}^2(2\beta ).$$
(17)
Numerical simulations of (15) using a 4th-order Runge-Kutta method has shown that the kink propagates without emitting any radiation if the initial velocity of the kink is chosen to be the velocity $`v_1`$ used in the expression of $`f^2`$. For initial velocities $`vv_1`$, there is virtually no emission of radiation so the kink velocity hardly decreases. For $`v>v_1`$, the kink initially propagates at velocity greater than $`v_1`$ but then slows down to a velocity just less than $`v_1`$ after which it’s virtually constant with little emission of radiation. The initial increase in velocity is due to the readjustment of the kink configuration to match the lattice field configuration.

## 3. Kink-antikink Collisions in the TDSG System

### 3.1 Preliminaries

In this section we present results of kink-antikink interactions in the Speight-Ward model. The results are obtained by solving (11) with the initial condition corresponding to Galilean boosted static kink and antikink solutions at each end of the lattice grid, i.e.

$`\psi (x,0)`$ $`=2\mathrm{tan}^1\mathrm{exp}a(x+b)2\mathrm{tan}^1\mathrm{exp}a(xb),`$ (18)
$`\dot{\psi }(x,0)`$ $`=av(\text{sech}a(x+b)+\text{sech}a(xb)).`$ (19)
The boundary condition is taken to be $`\psi (x_{min})=\psi (x_{max})=0`$, i.e. the ends of the kink and anti-kink are held fixed at all times.

The kink and antikink move towards each other with relative velocity $`2v`$. Equation (11) is solved using a 4th-order Runge-Kutta algorithm. The programme is run with various lattice spacings $`h`$. The time step is chosen to be 0.05. This conserves the energy to within 0.05% in a simulation which runs for 1000 units. We use the same time step for all values of $`h`$.

Our choice of the boundary condition means that any radiation emitted by the kink will reflect off the edges and eventually interact with the kinks. This is easily avoided by making the grid sufficiently large.

As mentioned already, the amount of radiation emitted by the kink depends on the kink velocity and also on the lattice spacing. For large velocities and coarser lattices the amount of radiation is significantly greater. In these cases we have damped the first few lattice sites at each end of the grid. That is, the value of $`\dot{\psi }`$ at these sites is decreased by a constant amount (10%) for the entire duration of the simulation. The total energy of the system is of course no longer conserved. The energy loss depends on the initial velocity of the kink and on the lattice spacing, but it is usually within 10% of the initial energy, for a simulation running for 1000 units.

### 3.2 Simulation Results

The simulations showed that the outcome of an interaction depends on the impact velocity of the kinks. The velocity of the kink (or the antikink) is defined to be the rate of change of $`X(t)`$, where $`X(t)`$ is the “average” position of the kink (or the antikink). $`X(t)`$ for the kink is given by

$$X(t)=\frac{{\displaystyle \underset{i=N}{\overset{0}{}}}x_ie_i}{{\displaystyle \underset{i=N}{\overset{0}{}}}e_i}.$$
(20)
Here $`e_i`$ is the total energy of the field and the radiation at the $`i^{th}`$ lattice site, where $`x_i=x_{min}+ih`$. The quantity $`_N^0e_i`$ is the total energy of the field and the radiation for the single kink. We found that for cases when the phonon radiation from the kinks is not significant, $`X(t)`$ makes a good approximation to the kink position. We define the velocity of the kink to be $`\dot{X}(t)`$, which is computed using the forward difference, $`\dot{X}(t)=\frac{1}{\mathrm{\Delta }t}(X(t+\mathrm{\Delta }t)X(t))`$.

Since the kink wobbles as it moves along the lattice, $`\dot{X}(t)`$ is an oscillating function. The period of oscillation is the wobble period, $`h/v`$, where $`v`$ is the velocity used in the initial condition (19). This wobble can be understood by considering the geodesic motion on the kink moduli space. This moduli space has a periodic metric .

Figure 1 shows various plots of kink velocity against time for kinks with different initial positions. The kinks are on a lattice with $`h=1.4`$ and the initial velocity of the kink is $`v=0.197`$. In all cases the average velocity of the kink is significantly decreased. This is due to phonon radiation from the kinks. The plots also show that the kink velocity depends on the initial position of the kink, modulo $`nh`$, where $`n`$.

The outcome of an interaction has sensitive dependence on the velocity of the kink. In particular, there are three different outcomes:

*(a)* If the kink impact velocity is above a certain critical velocity, $`v_c`$, then the kinks simply pass through each other, i.e. there is a smallest number $`v_c`$ such that the outcome of the interaction for all $`vv_c`$ is always a passing through behaviour. The value of $`v_c`$ depends on the lattice spacing and also on $`b`$. Figure 2 shows the dependence of $`v_c`$ against $`h`$. The dotted curve represents the quadratic $`0.073h^2+0.036h`$. This curve is an empirical fit to the data points and is not derived from the theory.

The dependence of $`v_c`$ on the position of the kinks is due to dynamical dressing and the kink wobble. For small initial velocities (figure 1), the kink velocity oscillates with a period $`(h/v)`$. The amplitude of oscillation is of the order of $`10^2`$. For large velocities, the phonon radiation make the oscillations erratic and also significantly decrease the overall velocity of the kink. The initial velocity also changes due to the readjustment of the kink field to suit the lattice distribution (dynamical dressing). The initial propagating velocity of the kink is not always $`v`$, the velocity used in expression (19). In some cases it is larger than $`v`$ and in others it is less than $`v`$. This change depends on the magnitude of $`v`$ and on the lattice spacing. The effect can be seen in figure 1.

For initial kink velocities which give impact velocities below $`v_c`$, we found that the outcome is either (b) a bounce or *(c)* a trapped breather-like state (a bion). The bounce is essentially a kink-antikink scattering, where the kinks after the first collision pass through each other but fail to escape to infinity. Instead, they travel a small distance and then turn around for another collision. After the second collision the kinks either escape to infinity or turn around again for more collisions. For a 2-bounce event (meaning the kinks pass each other twice) the ultimate outcome of this scattering would be a reflection, and for a 3-bounce a transmission.

The velocities for the bounces occur in small intervals. For velocities outside of these intervals the outcome is always a bound state. This is similar to the results for the $`\varphi ^4`$ model , except in this case the outcome does not depend fractally on the impact velocity. We have only seen 2, 3 and 4-bounce events. But the velocities at which the bounces occur are not easy to find and the intervals are small ($`10^7`$ for a 4-bounce), so it may be that higher bounce events occur but we have not seen them.

Figure 3 shows an example of the 2-bounce event. In this figure we have plotted the field $`\psi (0,t)`$ against $`t`$. $`\psi (0,t)`$ is the field at the centre of mass of the kink-antikink system.The kinks are on a lattice of unit spacing with an initial separation of 10 units. $`v_c`$ for this combination is found to be 0.100. The kinks bounce for all initial velocities in the range $`0.0969v0.0972`$ (2-bounce window).

Figures 4 and 5 show examples of three and four bounce events. In figure 4 the kinks transmit whereas in 5 they reflect off each other. As can be seen from the figure for the 2-bounce case, the kink and antikink collide, pass through each other and travel a small distance before stopping, turning around and travelling in the other direction. The kinks then collide again, pass through each other and this time escape to infinity.

Figure 6 is a graph of $`\psi (0,t)`$ against $`t`$ representing the oscillations at the centre of mass of a bion. This is essentially an $`n`$-bounce event for $`n`$ very large. The erratic oscillations suggest that the system might be chaotic. The period of oscillations become constant after about 10 cycles and the bion then oscillates for a long time (i.e for hundreds of cycles). We have calculated the maximum Lyapunov exponent of the time series $`\psi (0,t)`$, using the algorithm of Wolf *et al* . Its value is approximately 0.18, suggesting that the oscillations are mildly chaotic. In general, a positive maximum Lyapunov exponent is taken as a formal definition of chaos.

The results of the simulation can be understood using the resonance energy exchange mechanism. This requires the TDSG kink to support shape modes. We show next that the TDSG kink does support such a mode. The frequency of this mode is calculated in two ways, by using a collective coordinate approximation and also numerically, by Fourier analysing the deformations of a perturbed static kink. Note that we could also have used a matrix-based approach similar to the one used by Speight in , where it was used to calculate the Casimir energy of the kink.

## 4. Kink internal shape mode

#### 4.1 Collective coordinate approximation

We can obtain a collective coordinate approximation of the TDSG system by treating the scale factor $`a`$ in the kink as a dynamical variable. This reduces the number of degrees of freedom from infinity to one. The Lagrangian of the reduced system is given by

$$=f(a,b)\dot{a}^2g(a,b),$$
(21)
where $`f(a,b)`$ and $`g(a,b)`$ are functions given by

$`f(a,b)`$ $`={\displaystyle \frac{h}{4}}{\displaystyle \underset{x}{}}(xb)^2\text{sech}^2a(xb),`$ (22)
$`g(a,b)`$ $`={\displaystyle \frac{1}{h}}\mathrm{tanh}({\displaystyle \frac{1}{2}}ah)+{\displaystyle \frac{h}{4}}\mathrm{coth}({\displaystyle \frac{1}{2}}ah).`$ (23)
The corresponding equation of motion is

$$2f\ddot{a}+f^{}\dot{a}^2g^{}=0.$$
(24)
where the prime denotes differentiation with respect to $`a`$. $`g(a,b)`$ has a minimum (stable equilibrium) at $`a=s:=\frac{2}{h}\mathrm{tanh}^1\frac{h}{2}`$. Small amplitude oscillations $`a(t)=s+ϵ(t)`$ satisfy the linearised equation,

$$2f(s,b)\ddot{ϵ}=g^{\prime \prime }(s,b)ϵ.$$
(25)
Hence the waves oscillate with frequency

$$\nu =\frac{1}{2\pi }\sqrt{\frac{g^{\prime \prime }(s)}{2f(s,b)}}.$$
(26)
$`g^{\prime \prime }(s)=(1\frac{1}{4}h^2)^2`$, so the frequency of the shape mode depends on the lattice spacing. For $`h=1`$, the shape mode frequency is 0.152. Table 1 shows the dependence of frequency on the lattice spacing $`h`$.

| h | frequency |
| --- | --- |
| 0.1 | 0.175 |
| 0.2 | 0.175 |
| 0.5 | 0.170 |
| 0.8 | 0.161 |
| 0.9 | 0.157 |
| 1.0 | 0.152 |
| 1.2 | 0.143 |

As a check on the accuracy of the collective coordinate approach, we have used it to calculate the frequency of the shape mode of the $`\varphi ^4`$ kink. The value found is $`\sqrt{(}1.55)`$, which is fairly close to the actual value $`\sqrt{(}1.5)`$.

#### 4.2 Numerical method

Numerically, the shape mode frequency is obtained by Fourier analysing the deformations of the perturbed static kink. We set the initial condition for equation (11) to be

$`\psi (x,0)`$ $`=\psi _K(x,0)+ϵ`$ (27)
$`\dot{\psi }(x,0)`$ $`=0,`$
where $`\psi _K(x)`$ is the static kink solution and $`ϵ`$ is the pertubation. The value of $`ϵ`$ is chosen to be -0.03. The field $`\psi (x,t)`$ is sampled at various values of $`x`$. From this we obtain the variations in the field, $`\delta \psi (x,t)=\psi (x,t)\psi (x,0)`$, for each of the sampled values of $`x`$. $`\delta \psi (x,t)`$ were then Fourier analysed using the MATLAB FFT routine.

In figure 7 we have plotted the power spectrum for the data sampled at $`x=0`$ for a kink on an $`h=1`$ lattice. The power spectrum of the data for other values of $`x`$ give identical power spectrums.

Figure 7 was constructed using $`N=2^{14}=16384`$ data points. Aliasing of the power does not occur since the amplitude at the Nyquist frequency is essentially zero. There is a peak at $`\nu =0`$ corresponding to the Goldstone mode, and a peak at $`\nu =0.161`$ corresponding to the shape mode.

This procedure is repeated for other values of $`h`$. In contrast to the results suggested by the collective coordinate approach, we find that the frequency of the shape mode is more or less independent of $`h`$. For $`0.5h1.2`$, we found $`0.158<\nu <0.161`$. The value of the frequency (for $`h=1`$) however is within the range of frequencies found using the collective coordinate method.

In the continuum limit the TDSG shape mode frequency is precisely the quasimode frequency of the continuum sG system found by Rice, $`\omega _R=0.175`$ , although it is unlikely that such a quasimode exists

It is interesting to compare the situation with other discrete sine-Gordon systems. The well-known Frenkel-Kontorrova model does not support a kink with an internal shape mode for all values of the discreteness parameters. The sine-Lattice model (s-L) however does support one, as long as the system is only weakly discrete . The shape mode for the s-L system lies just above the lower phonon band suggesting that it is a genuine kink shape mode rather than just a phonon resonance. Moreover, in the continuum limit this shape mode converges to the sine-Gordon quasimode frequency. The existence of a shape mode in a discrete sine-Gordon system depends on the anharmonicity of the potential term.

## 5. Conclusion

We have found that the kink-antikink interactions in the TDSG system exhibit resonance phenomena similar to those in the continuum $`\varphi ^4`$ system. These are found to be due to the excitation of an intermal mode of the TDSG kink. A collective coordinate analysis shows that the frequency of this mode depends on the lattice spacing, and in the continuum limit $`(h0)`$ the frequency seems to approach the quasimode frequency predicted by Rice. The recent analysis done in  however says that there cannot be a quasimode for the continuum sine-Gordon kink, but they do predict the existence of a quasimode for a lattice sine-Gordon kink. This is not the quasimode we have found, so this is something genuinely new. Further numerical simulations are required to determine what is happening in the limit $`h0`$.

Also, as we have mentioned already, there is a modified TDSG system  which admits an exact travelling-kink solution (with a fixed velocity). It would be interesting to perform kink-antikink interactions for this system. Since the expressions we have used to perform the interactions are only an approximation to the equation of motion, it would be interesting to see the difference.

We have also done kink-antikink collision simulations for the TD$`\varphi ^4`$ system . The results are found to be analoguos to the conventional discrete $`\varphi ^4`$ system , though in the “windows” structure in the space of impact velocity is not a fractal.

*Acknowledgements* I would like to thank Richard Ward for helpful discussions and EPSRC for a research studentship.