text:
# Phenomenological Schwinger-Dyson approach to 𝜂-𝜂' mass matrix Talk delivered by D. Klabučar.

11institutetext: aRudjer Bošković Institute, P.O.B. 180, 10002 Zagreb, Croatia
bPhysics Department, Faculty of Science, Zagreb University

## Abstract

It is reviewed pedagogically how a very successful description of the $`\eta `$$`\eta ^{}`$ mass matrix can be achieved in the consistently coupled Schwinger-Dyson and Bethe-Salpeter approach in spite of the limitations of the ladder approximation. This description is in agreement with both phenomenology and lattice results.

The Schwinger-Dyson and Bethe-Salpeter (SD-BS) approach is the bound-state approach that is chirally well-behaved. (For example, see Ref. Roberts:2000hi ; Alkofer:2000wg ; Roberts:2000aa  for recent reviews, and references therein for various phenomenological and other applications of SD-BS approach.) Therefore, among the bound-state approaches it is probably the most suitable one to treat the light pseudoscalar mesons. One solves the Schwinger-Dyson (SD) equation for dressed propagators of the light $`u,d`$ and $`s`$ quarks using a phenomenologically successful interaction. These light-quark propagators are then employed in consistent solving of Bethe-Salpeter (BS) equations for various quark-antiquark ($`q\overline{q}`$) relativistic bound states. Namely, in both SD and BS equations we employ the same, ladder approximation, and the same interaction due to an effective, dressed gluon exchange. In the chiral limit (and close to it), light pseudoscalar ($`P`$) meson $`q\overline{q}`$ bound states ($`P=\pi ^{0,\pm },K^{0,\pm },\eta `$) then simultaneously manifest themselves also as (quasi-)Goldstone bosons of dynamically broken chiral symmetry. This resolves the dichotomy “$`q\overline{q}`$ bound state vs. Goldstone boson", enabling one to work with the mesons as explicit $`q\overline{q}`$ bound states (for example, in Refs. Kekez:1996az ; Klabucar:1997zi ; Kekez:1998xr ; Kekez:1998rw ; Kekez:2000aw ; Klabucar:2000me ; Kekez:2001ph ; Bistrovic:1999dy ) while reproducing (even analytically, in the chiral limit) the famous results of the axial anomaly for the light pseudoscalar mesons, namely the amplitudes for $`P\gamma \gamma `$ and $`\gamma ^{}P^0P^+P^{}`$ Alkofer:1995jx . This is unique among the bound state approaches – for example, see Refs. Roberts:2000hi ; Kekez:1998xr ; Alkofer:1995jx  and references therein. Nevertheless, one keeps the advantage of bound state approaches that from the $`q\overline{q}`$ substructure one can calculate many important quantities (such as the pion decay constant $`f_\pi `$) which are just parameters in most of other chiral approaches to the light-quark sector. The description Klabucar:1997zi ; Kekez:2000aw ; Klabucar:2000me ; Kekez:2001ph  of $`\eta `$$`\eta ^{}`$ complex is especially noteworthy, as it is very successful in spite of the limitations of the SD-BS approach in the ladder approximation.

For the description of $`\eta `$ and $`\eta ^{}`$, the crucial issues are the meson mixing and construction of physical meson states. They are formulated in Refs. Klabucar:1997zi ; Kekez:2000aw ; Klabucar:2000me  for the SD-BS approach, where solving of appropriate BS equations yields the eigenvalues of the squared masses, $`M_{u\overline{u}}^2,M_{d\overline{d}}^2,M_{s\overline{s}}^2`$ and $`M_{u\overline{s}}^2`$, of the respective quark-antiquark bound states $`|u\overline{u},|d\overline{d},|s\overline{s}`$ and $`|u\overline{s}`$. The last one is simply the kaon, and $`M_{u\overline{s}}`$ is its mass $`M_K`$. Nevertheless, the first three do not correspond to any physical pseudoscalar mesons. Thus, $`M_{u\overline{u}}^2,M_{d\overline{d}}^2,M_{s\overline{s}}^2`$ do not automatically represent any physical masses, although the mass matrix (actually, to be precise, the non-anomalous part of the mass matrix) is simply $`\widehat{M}_{NA}^2=\text{diag}(M_{u\overline{u}}^2,M_{d\overline{d}}^2,M_{s\overline{s}}^2)`$ in the basis $`|q\overline{q},(q=u,d,s)`$. However, the flavor SU(3) quark model leads one to recouple these states into the familiar octet-singlet basis of the zero-charge subspace of the light unflavored pseudoscalar mesons:

$`|\pi ^0`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|u\overline{u}|d\overline{d}),`$ (1)
$`|\eta _8`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{6}}}(|u\overline{u}+|d\overline{d}2|s\overline{s}),`$ (2)
$`|\eta _0`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}(|u\overline{u}+|d\overline{d}+|s\overline{s}).`$ (3)
With $`|u\overline{u}`$ and $`|d\overline{d}`$ being practically chiral states as opposed to a significantly heavier $`|s\overline{s}`$, Eqs. (13) do not define the octet and singlet states of the exact SU(3) flavor symmetry, but the effective octet and singlet states. However, as pointed out by Gilman and Kauffman Gilman:1987ax  (following Chanowitz, their Ref. ), in spite of this flavor symmetry breaking by the $`s`$ quark, these equations still implicitly assume nonet symmetry in the sense that the same states $`|q\overline{q}`$ ($`q=u,d,s`$) appear in both the octet member $`\eta _8`$ (2) and the singlet $`\eta _0`$ (3). Nevertheless, in order to avoid the UA(1) problem, this symmetry must ultimately be broken by gluon anomaly at least at the level of the masses of pseudoscalar mesons.

In the basis (13), the non-anomalous part of the (squared-)mass matrix of $`\pi ^0`$ and etas is

$$\widehat{M}_{NA}^2=\left(\begin{array}{ccc}M_\pi ^2& 0& 0\hfill \\ 0& M_{88}^2& M_{80}^2\hfill \\ 0& M_{08}^2& M_{00}^2\hfill \end{array}\right).$$
(4)
The $`\eta _8`$ “mass" $`M_{88}M_{\eta _8}`$ can be related to the kaon mass through the Gell-Mann–Okubo (GMO) relation, although the kaon does not appear in this scheme as it obviously cannot mix with $`\pi ^0`$ and etas, since it is strangely flavored. Equation (4) shows that also the isospin $`I=1`$ state $`\pi ^0`$ decouples from any mixing with the $`I=0`$ states $`\eta _8`$ and $`\eta _0`$, thanks to our working in the isospin limit throughout. Therefore, we are concerned only with the diagonalization of the $`2\times 2`$ submatrix in the subspace of etas in order to find their physical masses and corresponding $`q\overline{q}`$ content. In the isospin limit, obviously $`M_{u\overline{u}}=M_{d\overline{d}}`$, which we then can strictly identify with our model $`\pi ^0`$ mass $`M_\pi `$. Since in this model we can also calculate $`M_{s\overline{s}}^2=s\overline{s}|\widehat{M}_{NA}^2|s\overline{s}`$, this gives us our calculated entries in the mass matrix:

$$M_{88}^2\eta _8|\widehat{M}_{NA}^2|\eta _8M_{\eta _8}^2=\frac{2}{3}(M_{s\overline{s}}^2+\frac{1}{2}M_\pi ^2),$$
(5)
$$M_{80}^2\eta _8|\widehat{M}_{NA}^2|\eta _0=M_{08}^2=\frac{\sqrt{2}}{3}(M_\pi ^2M_{s\overline{s}}^2),$$
(6)
$$M_{00}^2\eta _0|\widehat{M}_{NA}^2|\eta _0=\frac{2}{3}(\frac{1}{2}M_{s\overline{s}}^2+M_\pi ^2).$$
(7)
The last one, $`M_{00}`$, is the non-anomalous part of the $`\eta _0`$ “mass" $`M_{\eta _0}`$. Namely, all the model masses $`M_{q\overline{q}^{}}`$ ($`q,q^{}=u,d,s`$) and corresponding $`q\overline{q}^{}`$ bound state amplitudes are obtained in the ladder approximation, and thus (irrespective of what the concrete model could be) with an interaction kernel which cannot possibly capture the effects of gluon anomaly. Fortunately, the large $`N_c`$ expansion indicates that the leading approximation in that expansion describes the bulk of main features of QCD. The gluon anomaly is suppressed as $`1/N_c`$ and is viewed as a perturbation in the large $`N_c`$ expansion. It is coupled only to the singlet combination $`\eta _0`$ (3); only the $`\eta _0`$ mass receives, from the gluon anomaly, a contribution which, unlike quasi-Goldstone masses $`M_{q\overline{q}^{}}`$’s comprising $`\widehat{M}_{NA}^2`$, does not vanish in the chiral limit. As discussed in detail in Sec. V of Ref. Klabucar:1997zi , in the present bound-state context it is thus best to adopt the standard way (see, e.g., Refs. Miransky ; Donoghue:dd ) to parameterize the anomaly effect. We thus break the $`U_A(1)`$ symmetry, and avoid the $`U_A(1)`$ problem, by shifting the $`\eta _0`$ (squared) mass by an amount denoted by $`3\beta `$ (in the notation of Refs. Kekez:2000aw ; Klabucar:2000me ). The complete mass matrix is then $`\widehat{M}^2=\widehat{M}_{NA}^2+\widehat{M}_A^2`$, where

$$\widehat{M}_A^2=\left(\begin{array}{ccc}0& \mathrm{\hspace{0.17em}0}& 0\hfill \\ 0& \mathrm{\hspace{0.17em}0}& 0\hfill \\ 0& \mathrm{\hspace{0.17em}0}& 3\beta \hfill \end{array}\right),$$
(8)
and where the value of the anomalous $`\eta _0`$ mass shift $`3\beta `$ is related to the topological susceptibility of the vacuum, but in the present approach must be treated as a parameter to be determined outside of our bound-state model, i.e., fixed by phenomenology or taken from the lattice calculations Alles:1996nm .

We could now go straight to the nonstrange-strange (NS-S) basis, but before doing this, it may be instructive to rewrite for a moment the matrix (8) in the flavor, $`|q\overline{q}`$ basis, where

$$\widehat{M}_A^2=\beta \left(\begin{array}{ccc}1& 1& 1\hfill \\ 1& 1& 1\hfill \\ 1& 1& 1\hfill \end{array}\right),$$
(9)
since for some readers it may be the best place to introduce the effect of flavor symmetry breaking. Namely, Eq. (9) tells us that due to the gluon anomaly, there are transitions $`|q\overline{q}|q^{}\overline{q}^{}`$; $`q,q^{}=u,d,s`$. However, the amplitudes for the transition from, and into, light $`u\overline{u}`$ and $`d\overline{d}`$ pairs need not be the same as those for the significantly more massive $`s\overline{s}`$. The modification of the anomalous mass matrix (9) which allows for possible effects of the breaking of the SU(3) flavor symmetry is then

$$\widehat{M}_A^2=\beta \left(\begin{array}{ccc}1& 1& X\hfill \\ 1& 1& X\hfill \\ X& X& X^2\hfill \end{array}\right).$$
(10)
There are arguments Kekez:2000aw ; Klabucar:2000me , supported by phenomenology, that the transition suppression is estimated well by the nonstrange-to-strange ratio of respective constituent masses, $`X_u/_s`$, or, as commented below, of respective decay constants, $`Xf_\pi /f_{s\overline{s}}`$.

After adding the anomalous contribution (8) to Eq. (4), pion still remains decoupled and we obviously still can restrict ourselves to $`2\times 2`$ submatrix in the subspace of etas. However, when dealing with quark degrees of freedom when the symmetry between the nonstrange (NS) and strange (S) sectors is broken as described above, the most suitable basis for that subspace is the so-called NS-S basis:

$`|\eta _{\text{NS}}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|u\overline{u}+|d\overline{d})={\displaystyle \frac{1}{\sqrt{3}}}|\eta _8+\sqrt{{\displaystyle \frac{2}{3}}}|\eta _0,`$ (11)
$`|\eta _\text{S}`$ $`=`$ $`|s\overline{s}=\sqrt{{\displaystyle \frac{2}{3}}}|\eta _8+{\displaystyle \frac{1}{\sqrt{3}}}|\eta _0.`$ (12)
The $`\eta `$$`\eta ^{}`$ mass matrix in this basis is

$$\widehat{M}^2=\left(\begin{array}{cc}M_{\eta _{NS}}^2\hfill & M_{\eta _S\eta _{NS}}^2\hfill \\ M_{\eta _{NS}\eta _S}^2\hfill & M_{\eta _S}^2\hfill \end{array}\right)=\left(\begin{array}{cc}M_{u\overline{u}}^2+2\beta \hfill & \sqrt{2}\beta X\hfill \\ \sqrt{2}\beta X\hfill & M_{s\overline{s}}^2+\beta X^2\hfill \end{array}\right)\begin{array}{c}\\ \varphi \end{array}\left(\begin{array}{cc}M_\eta ^2\hfill & 0\hfill \\ 0\hfill & M_\eta ^{}^2\hfill \end{array}\right),$$
(13)
where the indicated diagonalization is given by the NS–S mixing relations111The effective-singlet-octet mixing angle $`\theta `$, defined by analogous relations where $`\eta _{\text{NS}}\eta _8,\eta _\text{S}\eta _0,\varphi \theta `$, is related to the NS–S mixing angle $`\varphi `$ as $`\theta =\varphi \mathrm{arctan}\sqrt{2}=\varphi 54.74^{}`$. The relation between our approach and the two-mixing-angle scheme is clarified in the Appendix of Ref. Kekez:2000aw .

$$|\eta =\mathrm{cos}\varphi |\eta _{\text{NS}}\mathrm{sin}\varphi |\eta _\text{S},|\eta ^{}=\mathrm{sin}\varphi |\eta _{\text{NS}}+\mathrm{cos}\varphi |\eta _\text{S},$$
(14)
rotating $`\eta _{\text{NS}},\eta _\text{S}`$ to the mass eigenstates $`\eta ,\eta ^{}`$. Now the NS–S mass matrix (13) tells us that due to the gluon anomaly, there are transitions $`|\eta _{\text{NS}}|\eta _\text{S}`$. However, the amplitude for the transition from, and into, $`\eta _{\text{NS}}`$, need not be the same as those for the more massive $`\eta _\text{S}`$. The role of the flavor-symmetry-breaking factor $`X`$ is to allow for that possibility. As remarked little earlier, there are arguments Kekez:2000aw ; Klabucar:2000me , supported by phenomenology, that the transition suppression is estimated well by the nonstrange-to-strange ratio of respective quark constituent masses, $`_u`$ and $`_s`$. Due to the Goldberger-Treiman relation, this ratio is essentially equal Klabucar:1997zi ; Kekez:2000aw ; Klabucar:2000me  to the ratio of $`\eta _{\text{NS}}`$ and $`\eta _\text{S}`$ pseudoscalar decay constants $`f_{\eta _{\text{NS}}}=f_\pi `$ and $`f_{\eta _\text{S}}=f_{s\overline{s}}`$ which are calculable in the SD-BS approach. Same as $`_u`$ and $`_s`$, they were found in our earlier papers, especially Klabucar:1997zi ; Kekez:2000aw ; Klabucar:2000me . In other words, we can estimate the flavor-symmetry-breaking suppression factor as $`X_u/_s`$, or equivalently, as $`Xf_\pi /f_{s\overline{s}}`$. Our model results $`_u/_s=0.622`$ and $`f_\pi /f_{s\overline{s}}=0.689`$ are in both cases reasonably close to $`X_{exp}0.78`$ extracted phenomenologically Kekez:2000aw ; Klabucar:2000me  from the empirical mass matrix $`\widehat{m}_{exp}^2`$ featuring experimental pion and kaon masses, or, after diagonalization, $`\eta `$ and $`\eta ^{}`$ masses – see Eq. (7) in Ref. Klabucar:2000me . (In our model calculations below, we use $`X=f_\pi /f_{s\overline{s}}=0.689`$ to show the robustness of our approach to slight variations. Namely, our earlier works Kekez:2000aw ; Klabucar:2000me ; Klabucar:2001gr  mostly presented results based on slightly different values of $`X`$ obtained with the help of ratios of $`\gamma \gamma `$ amplitudes, which is yet another, but again related way of obtaining $`X`$.)

In our present notation, capital $`M_a`$’s denote the calculated, model pseudoscalar masses, whereas lowercase $`m_a`$’s denote the corresponding empirical masses. The empirical mass matrix $`\widehat{m}_{exp}^2`$ can be obtained from the calculated, model one in Eq. (13) by i) obvious substitutions $`M_{u\overline{u}}M_\pi m_\pi `$ and $`M_{s\overline{s}}m_{s\overline{s}}`$, and ii) by noting that $`m_{s\overline{s}}`$, the “empirical" mass of the unphysical $`s\overline{s}`$ pseudoscalar bound state, is given in terms of masses of physical particles as $`m_{s\overline{s}}^2=2m_K^2m_\pi ^2`$ due to the Gell-Mann-Oakes-Renner (GMOR) relation. Since $`M_{u\overline{u}}`$, obtained by solving the BS equation, is identical to our model pion mass $`M_\pi `$, it was fitted to the empirical pion mass $`m_\pi `$, e.g., in Ref. Kekez:2000aw . Similarly, $`M_{u\overline{s}}M_K`$ is fitted to the empirical kaon mass $`m_K`$. Therefore we also have $`M_{s\overline{s}}^22m_K^2m_\pi ^2`$, since our model has good chiral behavior and also satisfies the GMOR relation, thanks to which we have $`M_{s\overline{s}}^2=2M_K^2M_\pi ^2`$ in a very good approximation. We thus see that in our model mass matrix, the parts stemming from its non-anomalous part $`\widehat{M}_{NA}^2`$ (4) are already close to the corresponding parts in $`\widehat{m}_{exp}^2`$. We can thus expect a good overall description of the masses in $`\eta `$ and $`\eta ^{}`$ complex. We now proceed to verify this expectation.

The anomalous entry $`\beta `$ is fixed phenomenologically to be $`\beta _{exp}0.28`$ GeV2, along with $`X_{exp}0.78`$, by requiring that trace and determinant of $`\widehat{m}_{exp}^2`$ have their experimental values. But, this can be done also with our calculated, model mass matrix $`\widehat{M}^2`$. Requiring that the empirical value of the trace $`m_\eta ^2+m_\eta ^{}^21.22`$ GeV2 be fixed by Eq. (13), yields

$$\beta =\frac{1}{2+X^2}[(m_\eta ^2+m_\eta ^{}^2)_{exp}(M_{u\overline{u}}^2+M_{s\overline{s}}^2)]$$
(15)
where $`X=f_\pi /f_{s\overline{s}}=0.689`$ and $`M_{u\overline{u}}=0.1373`$ GeV and $`M_{s\overline{s}}=0.7007`$ GeV are now our model-calculated Kekez:2000aw  quantities, giving us $`\beta =0.286`$ GeV2. Since $`(M_{u\overline{u}}^2+M_{s\overline{s}}^2)=2m_K^2`$ holds to a very good approximation, our approach satisfies well the first equality (from the matrix trace) in

$$2\beta +\beta X^2=m_\eta ^2+m_\eta ^{}^22m_K^2=\frac{2N_f}{f_\pi ^2}\chi ,$$
(16)
where the second equality is the Witten-Veneziano (WV) formula WV , with $`\chi `$ being the topological susceptibility of the pure Yang-Mills gauge theory. Our model values of $`X`$ and $`\beta `$ ($`f_\pi `$ is fitted to its experimental value) thus imply $`\chi =(178\mathrm{MeV})^4`$, in excellent agreement with the lattice result $`\chi =(175\pm 5\mathrm{MeV})^4`$ of Alles et al. Alles:1996nm .

The mixing angle is then determined to be $`\varphi =43.2^{}`$ (or equivalently, $`\theta =11.5^{}`$), for example through the relation

$$\mathrm{tan}2\varphi =\frac{2\sqrt{2}\beta X}{M_{\eta _S}^2M_{\eta _{NS}}^2},$$
(17)
where

$$M_{\eta _{NS}}^2=M_{u\overline{u}}^2+2\beta =M_\pi ^2+2\beta =0.592\mathrm{GeV}^2=(769\mathrm{MeV})^2$$
(18)
and

$$M_{\eta _S}^2=M_{s\overline{s}}^2+\beta X^2=0.627\mathrm{GeV}^2=(792\mathrm{MeV})^2$$
(19)
are our calculated $`\eta _{NS}`$ and $`\eta _S`$ masses. They have reasonable values, in a good agreement with, e.g., $`\eta _{NS}`$ and $`\eta _S`$ masses calculated in the dynamical SU(3) linear $`\sigma `$ model Klabucar:2001gr .

The diagonalization of the $`NS`$-$`S`$ mass matrix gives us the $`\eta `$ and $`\eta ^{}`$ masses:

$`M_\eta ^2`$ $`=`$ $`\mathrm{cos}^2\varphi M_{\eta _{NS}}^2\sqrt{2}\beta X\mathrm{sin}2\varphi +\mathrm{sin}^2\varphi M_{\eta _S}^2`$ (20)
$`M_\eta ^{}^2`$ $`=`$ $`\mathrm{sin}^2\varphi M_{\eta _{NS}}^2+\sqrt{2}\beta X\mathrm{sin}2\varphi +\mathrm{cos}^2\varphi M_{\eta _S}^2.`$ (21)
Plugging in the above predictions for $`\beta ,X,M_{\eta _{NS}}`$ and $`M_{\eta _S}`$, our model $`\eta `$ and $`\eta ^{}`$ masses then turn out to be $`M_\eta =575`$ MeV and $`M_\eta ^{}=943`$ MeV. This is in good agreement with the respective empirical values of 547 MeV and 958 MeV.

However, the above is not all that can be said about agreement with experiment and other approaches. The second thing we may point out is the reasonable agreement we find if we insert our values of $`\beta `$, $`X`$ and $`M_{q\overline{q}^{}}`$’s into our model mass matrix and compare it with the $`\eta `$-$`\eta ^{}`$ mass matrix obtained on lattice by UKQCD collaboration McNeile:2000hf .

Third, Ref. Kekez:2000aw  clearly shows that our approach and results are not in conflict, but in fact agree very well with results in the two-mixing-angle scheme (reviewed and discussed in, e.g, Ref. Feldmann99IJMPA ). Actually, our results can also be given Kekez:2000aw ; KlKeFB  in the two-mixing-angle scheme.

Fourth, what we found from the mass matrix is consistent with what we found in the same SD-BS approach through another route, i.e. from $`\eta ,\eta ^{}\gamma \gamma `$ processes Klabucar:1997zi ; Kekez:2000aw ; Klabucar:2000me ; Kekez:2001ph .

The above shows that the consistently coupled SD-BS approach provides a surprisingly satisfactory description of the $`\eta `$-$`\eta ^{}`$ complex, especially if one recalls that $`\beta `$, parameterizing the anomalous $`\eta _0`$ mass shift, was the only new parameter. Namely, all other model parameters were fixed already by Ref. jain93b  providing the model we used in Refs. Kekez:1996az ; Klabucar:1997zi ; Kekez:1998xr ; Kekez:1998rw ; Kekez:2000aw ; Klabucar:2000me ; Kekez:2001ph . Nevertheless, we would like to point out that even this one parameter, $`\beta `$, can be fixed beforehand. Instead of being a parameter, $`\beta `$ can be obtained through WV formula (16) from the lattice results on the topological susceptibility $`\chi `$. The central value of widely accepted $`\chi =(175\pm 5\mathrm{MeV})^4`$ Alles:1996nm  would lead to $`\beta `$ less than 7% below our model value, which is within upper error bar anyway. Thus, we can eliminate $`\beta `$ as a free parameter and still achieve almost as satisfactory description as the one given above.

Acknowledgment: D. Klabučar thanks the organizers, M. Rosina, B. Golli and S. Širca, for their hospitality and for the partial support which made possible his participation at Mini-Workshop Bled 2002.


format: text

text:
# A single saddle model for the 𝛽-relaxation in supercooled liquids

## I Introduction

In recent years much attention has been devoted in the field of glass-forming liquids to the so-called ’energy landscape’ paradigm. The main background of this approach is the existence of a strong relation between the topological properties of the potential energy landscape (i.e. the potential energy as a function of the system’s coordinates) and the dynamical behaviour of the system.

Historically this idea goes back to Goldstein gold , who, in 1969, argued how at low enough temperatures supercooled liquids spend a long time in local minima of the potential energy, with rare activated jumps among minima stillinger ; zwanzig ; glotzer . According to Goldstein this low-$`T`$ description breaks down at a temperature $`T_c`$ above which activation is no longer the main mechanism of diffusion. This temperature $`T_c`$ has been subsequently identified with the critical temperature where idealized Mode Coupling Theory (MCT) locates a divergence of the relaxation time bengtzelius ; leutheusser ; goetze , supporting the view that $`T_c`$ marks a crossover from a low-$`T`$ activated hopping dynamics, to a continuous flow dynamics at higher $`T`$, which is well described by MCT angell ; ullo ; cummins ; sokolov ; doliwa .

According to this interpretation, from a topological point of view, it is clear that at low temperature the crucial role is played by the minima of the potential energy function, since the system spends most of the time close to them. This observation can be translated in quantitative terms and various approaches have been elaborated to use information available on the minima and their structure, to compute (thermo)dynamical quantities stillinger ; rahman ; cotterill ; keyes-1 .

On the other hand, above $`T_c`$ the system is no longer confined around a minimum, but rather ’flows’ in configuration space. From a topological point of view one may imagine that the system explores regions of the configuration space which are rich in negative modes. A confirmation of this perspective comes from the Instantaneous Normal Modes (INM) approach rahman ; cotterill : the INM spectrum, i.e. the average density of states of the Hessian matrix of the potential energy, displays above $`T_c`$ a finite fraction of imaginary modes, indicating that the system at equilibrium explores regions of the landscape which have negative curvature seeley ; madan ; keyes-1 . There have been various attempts to use the INM analysis to identify diffusive modes cho ; bembenek-1 ; sastry ; sciortino ; lanave , but it still remains unclear under what conditions this can be done, and how far the procedure can be pushed.

From a theoretical point of view, a crucial question is whether also above $`T_c`$ the behaviour of the system can be characterized by well defined topological entities, as it happens at low temperature with minima. This point has been addressed in detail in mean-field disordered systems in the context of spin-glasses kurchan-1 ; cavagna-1 ; cavagna-2 ; franz ; cavagna-2.2 ; thalmann , as well as in finite dimensional ordered models pettini . Given the numerous analogies between the phenomenology of certain spin-glass models and that of structural glasses thirumalai , and exploiting the results achieved in that field, a new topological interpretation of the dynamics above $`T_c`$ has been developed by various authors. Within this scenario, above $`T_c`$ it is still possible to describe the behaviour of the system in terms of well defined topological entities. These are not minima, but rather saddle points of the potential energy function, that is stationary points with a non zero fraction of unstable directions. More precisely, we can summarize this ’saddles interpretation’ in a few main points:

* For $`T>T_c`$ at equilibrium the system is ’close’ to a saddle point. The typical topological properties of this saddle point, for example its energy $`E`$ and the number of unstable directions $`K`$, depend on the temperature $`T`$.
* Since $`K0`$ the system explores a region where negative (non-confining) modes are available. Diffusion can then in principle occur in two different ways: either via barrier hopping along the positive eigenvalues, or exploiting the presence of these free directions (see cavagna-3  for a discussion of this point). When $`K(T)`$ is large this second mechanism will be faster than the first one, and will reasonably prevail. However, upon lowering the temperature $`K(T)`$ decreases and it will become less and less efficient.
* As the temperature approaches $`T_c`$ the typical saddles become less unstable, and at $`T_c`$ they finally cease to have negative modes and turn into minima. Thus, as $`TT_c`$ the system finds less and less free directions to diffuse and the second mechanism of diffusion previously described will freeze. If barriers at $`T_c`$ are relevant, as numerical simulations indicate for fragile systems grigera , this implies a dramatic increase of the relaxation time and a slowing down of the dynamics as $`T_c`$ is approached. Around $`T_c`$ a crossover occurs to a different regime where, being negative modes extremely rare, the main mechanism of diffusion becomes barrier crossing and the Goldstein scenario is thus recovered cavagna-3 .

This scenario has been tested in different ways and up to now various indications in its favour have been gathered: Numerical experiments on glass forming fragile models angelani ; kurt ; shah-1 ; grigera ; keyes-3 ; shah-2 ; wales ; dawson  have shown that it is possible to identify relevant saddles at a given temperature $`T`$ and classify them through their index $`K(T)`$. Also, and more importantly, $`K(T)0`$ as $`TT_c`$, $`T_c`$ being the Mode Coupling temperature. Besides, analytical computations in simple models cavagna-4  show that the system at equilibrium is truly close to a saddle according to a rigorous well defined notion of distance in the configuration space. As in numerical experiments, the index of the closest saddle goes to zero as the critical temperature is approached.

Despite these encouraging results the saddle scenario is not immune of some criticisms. The results of numerical simulations kurt ; grigera ; angelani  highly rely on numerical algorithms to locate the saddles, and the sampling procedure has been deeply questioned (see for example wales ). Also the whole energy landscape approach, even if supported by great part of the community, has recently received many objections (see for example berthier  where a description of the supercooled behaviour in terms of heterogeneities is preferred). Besides, the saddle scenario still remains from a theoretical point of view rather vague and incomplete. In particular, there are several questions which should be addressed, in our opinion, in a more formal and quantitative way:

* From an intuitive point of view and in the light of the numerical results, it is reasonable that saddles do play a crucial role in the dynamical behaviour of the system above $`T_c`$. However, is this just a useful qualitative interpretation, or can we push this approach further? In other words, is it possible to cast the ’saddle dynamics’ into a formal description?
* What is the precise relation between the vanishing of negative modes and the increase of the relaxation time as $`TT_c`$? In particular, what is the link between $`K(T)`$ and the relaxation time $`\tau `$?
* Is $`K(T)`$ really the topological quantity most relevant to the dynamics?
* Is a saddle description compatible with Mode Coupling Theory?

In this paper we will address these questions in a formal way by considering a very simple model of harmonic relaxation in a saddle. In our analysis we only assume as starting hypotheses two facts which have been numerically or analytically tested up to now: 1) the system above $`T_c`$ is close to a saddle point with non zero index $`K(T)`$, and 2) the number of negative modes of this saddle goes to zero as $`TT_c`$. The first assumption will enable us to perform an harmonic expansion of the potential energy around a saddle point, while the second will be used as a constraint on the behaviour of the saddle spectrum close to $`T_c`$.

Our aim is to find an explicit and quantitative link between the topological properties of the energy landscape and the dynamical behaviour of the system. The input of our model is then the spectrum of the saddle, while our results concern the relation between this spectrum and the relaxation time of the system. We will leave the shape of this spectrum unspecified as long as we can (apart from the above-mentioned constraints) in order to make our statements as general as possible. Only after will we discuss specific cases and their physical relevance. We will then show how even a simple relaxational dynamics is able to provide precise answers to the previous questions, and in certain cases, even approximately reproduce the results of more complicated dynamical theories such as Mode Coupling.

The paper is organized in the following way. In Sec. II we introduce our model and the main quantity to be considered, namely the mean square distance (MSD) (in configuration space) of the system from the typical reference saddle point at temperature $`T=(1+ϵ)T_c`$. The MSD will be expressed in terms of the topological properties of the saddle, more precisely its spectrum. In Sec. III we use the Laplace method to find a relation between the asymptotic dynamical behaviour of the MSD and the spectrum of the saddle point. In Sec. IV we analyze the consequences of our result, and discuss how the shape of the spectrum may determine qualitatively and quantitatively the long time dynamics. We look at different classes of spectra and focus on the conditions under which a two-step relaxation with a diverging time-scale is obtained. In Sec. V we compare out results with the prediction of the Mode Coupling Theory. Sec. VI is devoted to our conclusions.

## II The harmonic saddle model

Let us consider a glass forming system with $`N`$ particles and potential energy $`V(𝐫)`$, where $`𝐫=\{r_i\}`$ is the (mass-weighted) vector of the particles positions and $`i=1\mathrm{}N`$ is a particle index nota1 . As explained in the Introduction, the stationary points of $`V`$ are, from a topological point of view, the quantities of interest to describe the behaviour of the system.

At equilibrium in the low temperature phase $`T<T_c`$, the system (i.e. its representative point $`\underset{¯}{r}`$) explores for most of the time a region close to a minimum of $`V`$, only occasionally jumping via barrier crossing into the basin of another minimum. This separation of time scales between intra and inter basin motion has been exploited by both the INM and the Inherent Structures stillinger  approaches to compute dynamic (velocity autocorrelation) and thermodynamic quantities (equilibrium energy, entropy etc.) via an harmonic expansion of the potential energy around a reference configuration. On the other hand, above $`T_c`$ the system is typically close to a stationary point which is not a minimum, but a saddle with a certain number of negative modes. One can wonder whether also in this case it is possible to use this information to simplify the computations with a reasonable approximation, as the harmonic one adopted in the low temperature phase. If the relevant saddles have a large instability index $`K(T)`$, the system very easily finds escape directions to flow from one region of the phase to another one: even if the closest stationary point will typically be a saddle with the same statistical properties (index, energy etc.), it will not be the same saddle even for relatively short times. However, we know from the numerical and analytical works quoted in the Introduction, that when $`TT_c`$, $`K(T)0`$, thus close enough to $`T_c`$ the relevant saddles have a very small instability index. Then, if at some time the system is close to a given saddle it will take a long time before finding one of the few available escape directions and flow away. In this case the system remains close to the same saddle (we could say in its ’basin’) for some time, and an harmonic expansion seems justified, at least on time-scales shorter than the time when the saddle is definitely left.

Given that, we shall now consider our system at a temperature $`T=(1+ϵ)T_c`$ where $`ϵ`$ is a small parameter, and assume that at time it is $`t=0`$ close to a saddle point whose coordinates are $`r_0`$ in configuration space. We wish to examine the dynamical behaviour of the system and to do that we resort to a simplified treatment where we use an harmonic expansion of the potential $`V`$ around $`r_0`$. Our starting point is then the following Langevin equation:

$$\dot{x}_i=\underset{k=1}{\overset{N}{}}M_{ik}x_k+\xi _i(t),$$
(1)
where $`x_i(t)=r_i(t)r_i^0`$ are the displacements from the saddle, the matrix $`M_{ik}`$ is the second derivative of the potential $`V`$ at the saddle point and $`\xi `$ is a $`\delta `$-correlated noise, $`\xi _i(t)\xi _k(t^{})=2T\delta _{ik}\delta (tt^{})`$ foot-3 .

This equation has formally a very simple form of the Ornstein-Ulhenbeck type in $`3N`$ dimensions. Of course it is not as trivial as it may seem at a first sight in that the matrix $`M_{ik}`$ is not known. This is not only a consequence of having left the potential $`V`$ unspecified, since even when an explicit energy function is assumed, still the exact position $`𝐫_\mathrm{𝟎}`$ of a typical saddle at temperature $`T`$ is not determinable. This is the glassy nature of the system we are considering: even if the Hamiltonian is deterministic, as the temperature is lowered particles tend to arrange in disordered configurations. Thus $`M_{ik}`$ has rather to be considered a disordered matrix whose distribution is determined in a complicated way by the equilibrium Boltzmann measure (see e.g. cavagna-4 ; euclidian ). Thus in the following the Hessian matrix $`𝐌`$ will be treated as a random variable whose statistical properties are in principle accessible. In particular, we will express the dynamics in terms of the density of eigenvalues (spectrum) of $`𝐌`$:

$$\rho (\lambda ;r_0)=1/3N\underset{\alpha }{}\delta (\lambda \lambda _\alpha ),$$
(2)
where $`\lambda _\alpha `$ are the eigenvalues of $`M_{ik}`$. Of course, also $`\rho (\lambda ;r_0)`$ depends on the precise saddle where it is evaluated, and it is thus a stochastic quantity, as the Hessian is. However, its statistical properties are in general more easily computed via numerical simulations angelani ; kurt  or analytical procedures noi-inm ; grigera-2 ; stefano  and we have more information on its typical behaviour. For example, we know from previous works that the instability degree of a typical saddle decreases as the critical temperature is approached, which means that the spectrum evaluated in a typical saddle has less and less negative modes. Consequently, we will assume that the average spectrum $`\rho (\lambda )`$ depends parametrically on $`ϵ(TT_c)/T_c`$, and that its negative support vanishes for $`TT_c`$, i.e.

$$\rho _ϵ(0)0\text{for }ϵ0.$$
(3)
From eq. (1) we can easily compute the mean square displacement (MSD),

$$d_2(t)\frac{1}{3N}\underset{i=1}{\overset{N}{}}\underset{a=1}{\overset{3}{}}x_i^a(t)^2=\frac{T}{3N}\underset{\alpha =1}{\overset{3N}{}}\frac{1e^{2\lambda _\alpha t}}{\lambda _\alpha }=T_{\mathrm{}}^+\mathrm{}𝑑\lambda \rho _ϵ(\lambda )\frac{1e^{2\lambda t}}{\lambda },$$
(4)
where we have only assumed that in the thermodynamic limit fluctuations of the spectrum around its average value are negligible (as indicated by numerical simulations), and we have therefore substituted $`\rho (\lambda ;r_0)`$ with $`\rho _ϵ(\lambda )`$. Also, we have considered as initial conditions $`x_i(0)=0`$, that is we have assumed the system to start exactly on the top of the saddle. As we shall see, the initial condition is not important as long as the system does not start too far from the saddle (in Sec. III we will say exactly how far). Indeed the additional term arising in equation (4) when $`x_i(0)0`$ can be shown in this case to be exponentially small at large times.

This equation is the starting point of our analysis. All the physics is clearly encoded in the behaviour of the spectrum, both in its shape as a function of $`\lambda `$ and in its dependence on the temperature. We immediately see that if $`\rho _ϵ(\lambda )`$ has some negative support the integral will diverge for $`t\mathrm{}`$, meaning that the system asymptotically leaves the saddle if there are some negative eigenvalues, as expected. As already said, since we want our harmonic approximation to be reasonable, we are working in the regime $`ϵ1`$, when the support of $`\rho _ϵ(\lambda )`$ is almost entirely positive. In this case the systems remains a long time close to the saddle and we can study what happens in the large time limit, but before the saddle is left.

## III The asymptotic dynamics

For mathematical convenience it is better to look at the time derivative of $`d_2(t,ϵ)`$,

$$\dot{d}_2(t,ϵ)=2T_{\mathrm{}}^+\mathrm{}𝑑\lambda \rho _ϵ(\lambda )e^{2\lambda t}.$$
(5)
Since the MSD is a well defined physical quantity at any given time $`t`$, this integral should converge for any value of $`t`$. Convergence for $`\lambda \mathrm{}`$ requires that

$$\mathrm{log}\rho _ϵ(\lambda )<2\lambda t\text{for }\lambda \mathrm{}.$$
(6)
If we exclude oscillating functions this implies a concavity condition on $`\rho (\lambda )`$ (this is obvious graphically): it must exist a $`\lambda _0(ϵ)`$ such that

$$\frac{d^2}{d\lambda ^2}\mathrm{log}\rho _ϵ(\lambda )<0\text{for }\lambda <\lambda _0(ϵ).$$
(7)
The support of the spectrum must be entirely positive for $`ϵ=0`$, and this implies that $`\lambda _0(ϵ)>0`$ for $`ϵ`$ small enough. We can split the integral in (5) separating the domains $`\lambda <\lambda _0`$ and $`\lambda >\lambda _0`$. We have

$$2T_{\lambda _0}^+\mathrm{}𝑑\lambda \rho _ϵ(\lambda )e^{2\lambda t}\rho _{max}\frac{e^{2\lambda _0t}}{t},$$
(8)
where $`\rho _{max}`$ is the maximum of $`\rho _ϵ(\lambda )`$. As we shall see the remaining part of the integral in $`\dot{d}_2(t,ϵ)`$ is much larger than $``$ for $`t1`$. We can thus disregard $``$ and write

$$\dot{d}_2(t,ϵ)=2T_{\mathrm{}}^{\lambda _0}𝑑\lambda \rho _ϵ(\lambda )e^{2\lambda t}=2T_{\mathrm{}}^{\lambda _0}𝑑\lambda e^{tS_ϵ(\lambda ,t)},$$

with $`S_ϵ(\lambda ,t)=2\lambda \mathrm{log}\rho _ϵ(\lambda )/t`$. Thus, only the left tail of the spectrum contributes to the behaviour of $`\dot{d}_2(t,ϵ)`$.

To evaluate this integral in the regime $`t1`$ we can use the Laplace (saddle point) method bender  to find:

$$\dot{d}_2(t,ϵ)=d_0T\frac{e^{tS_ϵ(\widehat{\lambda }_ϵ(t),t)}}{\sqrt{tS_ϵ^{\prime \prime }(\widehat{\lambda }_ϵ(t),t)}},$$
(9)
where $`d_0`$ is a constant independent of $`t`$ and $`ϵ`$, the prime indicates the derivative with respect to $`\lambda `$ and $`\widehat{\lambda }_ϵ(t)`$ is the solution of the saddle point equation $`S_ϵ^{}(\widehat{\lambda }_ϵ,t)=0`$, namely

$$\frac{\rho _ϵ^{}(\widehat{\lambda }_ϵ)}{\rho _ϵ(\widehat{\lambda }_ϵ)}=2t.$$
(10)
From the convergence condition (7) we see that $`\rho _ϵ^{}(\lambda )/\rho _ϵ(\lambda )`$ is a monotonically decreasing function of $`\lambda `$, and therefore eq. (10) has a unique solution $`\widehat{\lambda }_ϵ(t)`$. Besides, we see that $`\widehat{\lambda }_ϵ(t)`$ changes sign at a well defined time: if we define

$$t_ϵ=\frac{1}{2}\frac{\rho _ϵ^{}(0)}{\rho _ϵ(0)},$$
(11)
we have that

$`\widehat{\lambda }_ϵ(t)>0\text{for }t<t_ϵ`$ (12)
$`\widehat{\lambda }_ϵ(t)<0\text{for }t>t_ϵ`$ (13)
If we interpret the saddle value $`\widehat{\lambda }(t)`$ as the relevant mode at time $`t`$, this result translates the intuitive fact that while at short times the system does not realize the presence of escape directions and ’relaxes’ in the pseudo-basin around the saddle given by the positive modes, as time increases it finally finds the unstable modes and exploits them to flow away. This argument already tells us that $`t_ϵ`$ represents a crucial dynamical time scale in our problem. This can be appreciated in a more explicit and quantitative way if we analyze in detail the behaviour of $`\dot{d}_2(t,ϵ)`$.

As we show in Appendix A, it is possible to rewrite (9) as

$$\dot{d}_2(t,ϵ)=D_0T\sqrt{\dot{\widehat{\lambda }_ϵ}(t)}\mathrm{exp}\left(2_1^t𝑑t^{}\widehat{\lambda }_ϵ(t^{})\right),$$
(14)
where $`D_0`$ is a constant. This equation shows that the sign of $`\widehat{\lambda }_ϵ(t)`$ is the key factor determining the asymptotic behaviour of $`\dot{d}_2(t,ϵ)`$ and thus, depending on the value of $`t`$ compared to $`t_ϵ`$, we find two time regimes:

* Early time region ($`1t<t_ϵ`$): In this regime $`\widehat{\lambda }_ϵ(t)>0`$ and thus $`\dot{d}_2(t,ϵ)`$ is a decreasing function of $`t`$, reaching its minimum at $`t=t_ϵ`$. From (14) we see that $`\dot{d}_2(t_ϵ,ϵ)0`$ for $`t_ϵ\mathrm{}`$. This means that if the time scale $`t_ϵ`$ diverges at $`T_c`$, then $`d_2(t,ϵ)`$ develops a plateau, whose length is of order $`t_ϵ`$ MSD . The value of $`d_2`$ at the plateau is given by

$$d_p(ϵ)=T_{1/t_ϵ}^{\mathrm{}}\frac{d\lambda }{\lambda }\rho _ϵ(\lambda ).$$
(15)
The demonstration of this expression is sketched in Appendix B. Physically, it is rather intuitive from eq. (4) if one remembers that at $`t=t_ϵ`$ only positive modes contribute, and makes the additional approximation that those with $`\lambda >1/t_ϵ`$ have reached their asymptotic value, while the rest do not contribute at all. Since the length of the plateau is proportional to $`t_ϵ`$, it is clear that the plateau itself is better defined the closer we are to $`T_c`$. Note also that for large $`t_ϵ`$, the expression for the plateau is just the value of $`d_2(t,ϵ)`$ which would be obtained for a strictly positive spectrum by taking the limit $`t\mathrm{}`$ in (4). This means that for $`1t<t_ϵ`$ the system “thermalizes” in the saddle using only the positive eigenvalues, and $`d_2`$ reaches the value it would have reached in a purely harmonic well. This justifies our initial condition $`x_i(0)=0`$: any other choice with $`d_2(t,ϵ)d_p`$ would have been the same. The average energy density is given by $`E(t)=E_0+\frac{1}{2}T[1\dot{d}_2(t,ϵ)]`$, where $`E_0`$ is the bare energy of the saddle. Therefore, at the plateau, that is for $`tt_ϵ`$, the corrections to the harmonic thermal energy $`T/2`$ are small if $`t_ϵ1`$.
* Late time region ($`tt_ϵ`$): In this regime we have that $`\widehat{\lambda }_ϵ(t)<0`$. Now the system has finally found the unstable directions and its dynamics is ruled by the escape from the saddle basin. The integral in the exponential of (14) changes sign at a time $`\tau `$ defined by $`_1^\tau 𝑑t^{}\widehat{\lambda }_ϵ(t^{})=0`$. Besides, in Appendix C we show that the prefactor in the square root does not go to zero exponentially for $`t\mathrm{}`$. Thus, $`\dot{d}_2(t,ϵ)\mathrm{}`$ for $`t\tau `$, corresponding to $`d_2(t,ϵ)`$ leaving the plateau, i.e. to the system leaving the saddle. We could in principle identify $`\tau `$ as an ’escape’ time and consider it as a second relevant dynamical time scale. However, as we have stressed in the Introduction, our harmonic expansion is meaningful as long as the system remains close to the saddle, therefore it is not clear whether $`\tau `$ has a genuine physical meaning or whether other terms in the potential expansion should already be considered on these time scales.

We stress that these results are valid for a general $`\rho _ϵ(\lambda )`$, with minimal requirements on its behaviour dictated by physical consistency. To summarize, the most important result is the existence of a relevant time-scale $`t_ϵ`$, directly expressed in terms of topological properties of the landscape (the spectrum at the saddle point). This result answers one of the questions discussed in the Introduction, namely the relation between a topology and a relaxation time. Contrarily to the naive expectation, the relaxation time $`t_ϵ`$ is not naturally related to the fraction of negative eigenvalues $`k_ϵ=_{\mathrm{}}^0\rho _ϵ(\lambda )𝑑\lambda `$. Rather, if there is a degree of universality in saddle dynamics, the correct scaling variable seems to be $`\rho _ϵ^{}(0)/\rho _ϵ(0)`$ (see (eq. 11)). A second general result is that the MSD exhibits a two step relaxation close to the critical temperature, provided the relevant time scale $`t_ϵ`$ diverges at $`T_c`$. Given the expression of $`t_ϵ`$ and given that $`\rho _ϵ(0)0`$ for $`ϵ0`$, it is likely that this is actually the most general scenario for most of reasonable spectra shapes (even if there are simple counter-examples, as the case $`\rho _ϵ(\lambda )=ϵf(\lambda )`$).

On the other hand, the way the system approaches and leaves the plateau, which is encoded in expression (14), depends on the specific form of $`\rho _ϵ(\lambda )`$. Therefore, to extend further our analysis we need at this point to specify the spectrum in some detail. Thus in the next section we shall consider different spectral shapes.

## IV Specific spectra

Let us now consider some specific forms for the spectrum $`\rho (\lambda )`$ and compute explicitly the behaviour of $`d_2(t)`$. Since we are interested in the asymptotic behaviour we shall specify the spectrum only in the left tail, that is for $`\lambda <\lambda _0`$.

### IV.1 Exponential Case

Let us consider

$$\rho _ϵ(\lambda )\rho _ϵ(0)\mathrm{exp}(a\lambda )$$
(16)
where $`lim_{ϵ0}\rho _ϵ(0)=0`$ in order to ensure a positive definite spectrum at $`T_c`$. This spectrum is of the form $`\rho _ϵ(\lambda )=ϵf(\lambda )`$ and, as noted in the previous section, it has a non-diverging time scale $`t_ϵ=a/2`$. Besides, from expression (4) we see that $`d_2(t)`$ diverges at a finite time $`t=t_ϵ`$. Equation (5) can be easily integrated and we get

$$\dot{d}_2(t,ϵ)2T\rho _ϵ(0)\mathrm{erfc}(\lambda _0(a2t))+.$$
(17)
Therefore in this case we do not get any diverging time scale as a function of $`ϵ`$ nor a two step dynamics. Rather what happens is that $`d_2(t,ϵ)`$ diverges at the same time for any temperature above $`T_c`$, but in a steeper and steeper way as $`T_c`$ is approached.

### IV.2 The Gaussian case

The Gaussian tail is less trivial. In this case we consider

$$\rho _ϵ(\lambda )e^{\frac{b}{2}(\lambda \overline{\lambda }_ϵ)^2},$$
(18)
with $`lim_{ϵ0}\overline{\lambda }_ϵ^1=0`$. By applying the Laplace method we find $`\widehat{\lambda }_ϵ(t)=\overline{\lambda }_ϵ2t/b`$ and $`t_ϵ=b\overline{\lambda }_ϵ/2`$. Then $`\dot{d}_2`$ is given by

$$\dot{d}_2(t,ϵ)D_0T\mathrm{exp}\left[\frac{2t^2}{b}2\overline{\lambda }_ϵt\right],$$
(19)
and

$$d_2(t,ϵ)d_p(ϵ)+D_0Te^{\overline{\lambda }_ϵt_ϵ}_0^{tt_ϵ}e^{\frac{2\tau ^2}{b}}𝑑\tau .$$
(20)
(The same expressions can of course be obtained by exactly solving the Gaussian integral eq. (5)). Therefore for a Gaussian tail we get a two-step dynamics with a plateau and a diverging time-scale. The plateau is approached exponentially, as we see from eq. (20).

### IV.3 The power law case

This is actually the most interesting case, since it gives predictions compatible with MCT. Besides, the $`p`$-spin spherical model which has been widely used as a mean-field description of glassy physics and exactly obeys MCT equations, has a spectrum belonging to this class (see Sec. V-C).

Let us consider

$$\rho _ϵ(\lambda )=(ϵ^\mu +\lambda )^\eta ,$$
(21)
with $`\mu \eta >0`$, such that $`\rho _ϵ(0)=ϵ^{\mu \eta }0`$, for $`ϵ0`$. From eq. (10) we have $`\widehat{\lambda }_ϵ(t)=ϵ^\mu +\eta /2t`$, while eq. (14) gives

$$\dot{d}_2(t,ϵ)=D_0T\frac{e^{2ϵ^\mu t}}{t^{\eta +1}}.$$
(22)
By integrating from $`t_ϵ`$ (given by (11)) up to $`t`$, we find,

$$d_2(t,ϵ)=d_p(ϵ)+ϵ^{\mu \eta }h(t/t_ϵ),t_ϵ=\frac{\eta }{2}ϵ^\mu ,$$
(23)
where $`h(x)`$ is a scaling function obtained by the expansion of the exponential in (22). The approach to the plateau can be found by noting that for $`x1`$ we have $`h(x)1/x^\eta `$, and thus

$$d_2(t,ϵ)d_p(ϵ)t^\eta .$$
(24)
These three general shapes of the the spectrum give rise to different dynamical behaviours. In particular, only a power spectrum gives rise to a diverging time-scale and a power-law approach to the plateau. This fact is particularly important when comparing the predictions of this simple model with those of the MCT, as we shall do in the next section.

To conclude, we note that a useful quantity to characterize the approach to the plateau is the time-dependent effective exponent crisanti ,

$$\nu _{\mathrm{eff}}(t,ϵ)=1+\frac{d\mathrm{log}\dot{d}_2(t,ϵ)}{d\mathrm{log}t}=1+\frac{t}{2}\frac{\ddot{\widehat{\lambda }}_ϵ(t)}{\dot{\widehat{\lambda }}_ϵ(t)}\mathrm{\hspace{0.17em}2}t\widehat{\lambda }_ϵ(t).$$

If in some time regime the MSD has a power law dependence on time, $`d_2(t,ϵ)=d_p\pm t^\nu `$, this must show up as a constant contribution to the effective exponent in some extended region of time, that is $`\nu _{\mathrm{eff}}(t)=\nu `$. If we apply this definition to the previous cases, we see easily that the only one exhibiting a constant region for $`\nu _{\mathrm{eff}}(t)`$ is the third one (see Fig. 1), for which find $`\nu _{\mathrm{eff}}(t)=\eta (1t/t_ϵ)`$.

## V Comparison with MCT

Mode Coupling Theory (MCT) bengtzelius ; leutheusser ; goetze  is a dynamical approach which has been widely tested on a great range of glass forming systems. It is not an exact theory, in that it makes some assumptions on the memory functions entering the exact dynamical equations, but it deals with the complete dynamics of the system and takes into account in detail its structural properties. For structural fragile glasses it is now commonly accepted that MCT represents an excellent description of the supercooled liquid system for $`TT_c`$ (see for ex. kob ).

Our single saddle model is, on the other hand, a much more simplified approach and we cannot expect it to reproduce the whole set of MCT results. However, if we believe that close to $`T_c`$ the main mechanism driving the dynamical slowing down is the vanishing of the negative modes available to the system, our model should describe reasonably the behaviour of the system in this temperature range and should therefore reproduce at least the main MCT predictions.

### V.1 MCT predictions

The MCT predicts the existence of a critical temperature $`T_c`$ where the relaxation time of the system diverges. In real systems a real dynamical transition is not observed, however $`T_c`$ still represents a well defined relevant temperature where a slowing down of the dynamics takes place and a clear crossover to a different dynamical regime arrives.

A general prediction of MCT is the presence, as $`T`$ approaches $`T_c`$ from above, of a two-step relaxation of the relevant dynamical quantities, such as, for example correlation, response functions and the MSD. More than this, MCT provides with detailed predictions about the approach and escape from the plateau and the final relaxation to the equilibrium asymptotic value. It distinguishes two dynamical regimes, both of which become critical at the transition: the so-called $`\beta `$-relaxation regime, which concerns the dynamics around the plateau, and the $`\alpha `$-relaxation regime, which is related to structural relaxation to the ergodic values. Here, we will be concerned with the $`\beta `$ relaxation, since, as already underlined, our model is a reasonable approximation only until when the plateau is left.

According to MCT, in the $`\beta `$-relaxation regime the correlation function $`C(t)`$ has the following scaling form (as usual $`t1`$),

$$C(t)=C_p+ϵ^{1/2}g(t/t_ϵ),t_ϵ=ϵ^{\frac{1}{2a}}.$$
(25)
The correlation function has a plateau $`C_p`$ of length $`t_ϵ`$, and $`t_ϵ`$ diverges as a power law for $`ϵ0`$. Moreover, MCT predicts a power law approach to the plateau (early $`\beta `$ regime), that is,

$$C(t)C_pt^a,1tt_ϵ,$$
(26)
whereas for times larger than $`t_ϵ`$ the correlation function leaves the plateau as (late $`\beta `$ regime)

$$C(t)C_p(t/t_0)^b,t_ϵ<t\tau ,$$
(27)
where $`t_0=ϵ^{\frac{a+b}{2ab}}`$, and $`\tau `$ is the $`\alpha `$-relaxation time.

For the supercooled system studied in kob  the value of the exponent $`a`$ is $`a=0.28`$.

### V.2 Some general requirements

Among the spectral shapes we have considered, only the power case gives a scaling compatible with MCT. Looking at equations (23) and (24), we conclude that the dynamical scaling of MCT is reproduced by spectra of the form

$$\rho _ϵ(\lambda )=(ϵ^\mu +\lambda )^\eta ,\mu \eta =1/2,$$
(28)
the predicted dynamical exponent being $`\alpha =\eta `$.

Equation (28), also implies a simple condition on the behaviour of the spectrum in $`\lambda =0`$: only spectra that behave as

$$\rho _ϵ(0)ϵ^{1/2}$$
(29)
are compatible with MCT. This seems to be a quite strong topological requirement, very easy to check in specific cases.

Finally, we note that the important point is that the power law behaviour is obeyed in the main range of the negative support as $`T_c`$ is approached. The presence of small correcting tails around the lower band edge $`\lambda =ϵ^\mu `$ is not important if they arise in a vanishing (as $`ϵ0`$) region. Indeed, let us imagine that the power law behaviour (28) holds down to a certain value $`\overline{\lambda }>ϵ^\mu `$. It is simple to show that the dynamical behaviour predicted by a power law spectrum still holds up to time-scales $`t\overline{t}`$ with $`\overline{t}=\rho _ϵ^{}(\overline{\lambda })/2\rho _ϵ(\overline{\lambda })`$, while for greater times the tail effect starts dominating. Therefore if $`t_ϵ/\overline{t}=(ϵ^\mu +\overline{\lambda })/ϵ^\mu `$ vanishes as $`ϵ0`$, the dynamical scaling still holds in a diverging time window. In finite dimensional systems one usually gets small tails in the spectrum where localized modes are concentrated. The previous observation suggests that for a supercooled liquid, where we may expect small tails to smoothen the spectrum close to the lower band edge, only the extended modes representing the main negative support are truly relevant for diffusion.

### V.3 The $`p`$-spin spherical model

To perform a real comparison between MCT and our simplified harmonic model we need to consider specific cases. In particular, an excellent test would be a system where both approaches can be fully applied. That is, we need a model where i) the MCT dynamical phenomenology is reproduced and ii) information on the local topology is analytically available.

The $`p`$-spin spherical model (PSM) crisanti  satisfies both these requirements. Indeed, it can be shown that the exact dynamical equations for this model have a MCT structure bouchaud , that is for the PSM the MCT is *exact* (contrarily to structural glasses where it represents an approximation). Besides, for this model saddles have been shown analytically to play a relevant role cavagna-4  and the spectrum of a typical saddle at temperature $`T=(1+ϵ)T_c`$ is exactly known kurchan-2 ; kurchan-3 . Its left tail is

$$\rho _ϵ(\lambda )=(ϵ+\lambda )^{1/2}.$$
(30)
Since in the PSM $`d_2(t)=2[1C(t)]`$, the MCT scaling forms of eqs. (25), (26) and (27) directly apply also to the mean-square displacement. The exact resolution of the dynamical MCT equations for $`C(t)`$ can be found in crisanti  and gives $`a_{PSM}0.4`$.

On the other hand the behaviour predicted by our harmonic model can be simply inferred by the results of the previous sub-section with $`\mu =1`$ and $`\eta =1/2`$. We get

$$d_2(t,ϵ)=d_p(ϵ)+ϵ^{1/2}h(t/t_ϵ),t_ϵ=\frac{1}{4ϵ},$$
(31)
and a power law approach to the plateau:

$$d_2(t,ϵ)d_p(ϵ)t^{1/2}.$$
(32)
This behaviour is consistent with a MCT scaling form with an exponent $`a=1/2`$, to be compared with the exact value $`a_{\mathrm{PSM}}0.4`$. Our conclusion is then that for the PSM the single saddle model correctly reproduces the general dynamical scaling of the complete MCT dynamics and the power law approach to the plateau of the early $`\beta `$ regime. The power law approach to the plateau can also be seen by setting $`tt_ϵ`$ in the effective exponent. This gives, $`\nu _{\mathrm{eff}}(t,ϵ)=1/2(1t/t_ϵ)`$. From this expression we see that for $`t>t_ϵ`$ the model gives no power law departure from the plateau (see Fig. 1) psm . This must not be regarded as significant since the model itself looses validity at large times (for this reason we have only focused on the early $`\beta `$ regime).

Unfortunately, our model does not exactly reproduce the exponent $`a`$ (or, in other terms, it does not give the correct form of the scaling function $`g(x)`$ (25), which is responsible for the value of the exponent $`a`$). Of course, this leaves us partly unsatisfied. However, it is known that anharmonicities are very strong in the PSM, and it is reasonable that they are already relevant close to the saddle. If this is the case one should be able to show that taking into account anharmonicities, for example via a perturbative approach, modifies the exponent in the correct direction. We are now working in this direction.

## VI Conclusions

In this paper we have analyzed a simple model of relaxational dynamics around an harmonic saddle. The physical problem we wanted to address was the description of supercooled liquids close to the crossover temperature $`T_c`$, and our main purpose was to model in a formal way the connection between topological properties of the potential energy and dynamical behaviour which is at the basis of the energy landscape approach.

Despite the extreme simplicity of our model and the generality of our assumptions our results are not trivial. We have obtained a general expression which relates the spectrum of the saddle and a relevant time-scale, and we have outlined the conditions under which such a time-scale diverges giving rise to a two-step dynamics with a well defined plateau at $`T_c`$. This time-scale is related to unstable modes, however it is not proportional to the instability index, as naive expectation would have suggested. Rather, it is related to the behaviour of the spectrum close to zero: since negative modes tend to disappear as $`T_c`$ is approached it is the way they turn into soft modes that determines the long time dynamics.

The precise way the system relaxes to the plateau depends on the shape of the saddle spectrum. We have analyzed different spectral shapes and we have shown that power law spectra give rise to the same dynamical scaling as the one predicted by MCT in the $`\beta `$-regime. For this kind of spectra then our model seems to be consistent with MCT.

Beside these general qualitative predictions we have also looked at the specific case of the $`p`$-spherical model, where both the spectrum and the dynamics can be analytically computed. Here the saddle model reproduces the dynamical scaling predicted by the MCT, which is in this case exact. The $`p`$-spin spherical model is interesting because it represents a concrete model where the two approaches, MCT and single-saddle relaxation, can be independently performed and ultimately compared. Our analysis indicates that, despite its brutal disregarding of anharmonic contributions, the single saddle model already reproduces the main features of the complete exact dynamics, namely the dynamical scaling and power law approach to the plateau. Our interpretation of this result has been sketched in the Introduction: the presence of the plateau and the slowing down of the dynamics described by the MCT equations are mainly due to the vanishing of the negative modes close to the transition temperature $`T_c`$, and these effects are already taken into account at the harmonic level in the single saddle model. In the PSM we have explicitly demonstrated this, however we expect the same to hold whenever MCT is a good description of the problem, as for structural fragile glasses. To support this conclusion we should in principle proceed as we did for the PSM: compute the dynamical quantities using the single saddle model and compare them with the MCT predictions. The problem is that, contrarily to the $`p`$-spin case, the spectrum of a typical saddle is not in general known for models of glass forming systems. We may expect that numerical simulations or analytic computations will provide this topological information and our analysis be completed in the future. On the other hand, there is something important we have already done at this stage: we have outlined under what general conditions for $`\rho (\lambda )`$ our single saddle model is compatible with the MCT predictions. If the spectra of real glass systems will turn out to satisfy these conditions we will then have very strong support for our arguments.

###### Acknowledgements.

We acknowledge rewarding discussions with S. Ciliberti, G. Parisi, F. Ricci-Tersenghi and P. Verrocchio.

## Appendix A Computation of the MSD

Here we show how to derive expression (14) for the MSD. For convenience we omit in the following the sub-index $`ϵ`$. We remind that the dot indicates a derivative with respect to time, while the prime a derivative with respect to $`\lambda `$. Let us call $`N(t)=e^{tS(\widehat{\lambda })}`$. By deriving $`N(t)`$ with respect to time the numerator of eq. (9) and using the saddle point equation (10) we obtain

$$\dot{N}(t)=2\widehat{\lambda }(t)N(t),$$
(33)
and thus

$$N(t)=N(t_0)\mathrm{exp}\left(2_{t_0}^t𝑑t^{}\widehat{\lambda }(t^{})\right).$$
(34)
For what concerns the denominator of (9) we have to evaluate

$$tS^{\prime \prime }(\widehat{\lambda })=\left(\frac{\rho ^2}{\rho ^2}\frac{\rho ^{\prime \prime }}{\rho }\right)=t^2\frac{\rho ^{\prime \prime }(\widehat{\lambda })}{\rho (\widehat{\lambda })},$$
(35)
where we have used again the saddle point equation. We define $`f(t)=\rho (\widehat{\lambda }(t))`$. Derivatives of $`f(t)`$ with respect to time give

$`\dot{f}(t)`$ $`=`$ $`t\dot{\widehat{\lambda }}f(t),`$ (36)
$`{\displaystyle \frac{1}{f(t)}}\ddot{f}(t)`$ $`=`$ $`{\displaystyle \frac{\rho ^{\prime \prime }(\widehat{\lambda })}{\rho (\widehat{\lambda })}}\dot{\widehat{\lambda }}^2+t\ddot{\widehat{\lambda }},`$ (37)
$`{\displaystyle \frac{1}{f(t)}}\ddot{f}(t)`$ $`=`$ $`\dot{\widehat{\lambda }}+t\ddot{\widehat{\lambda }}+t^2\dot{\widehat{\lambda }}^2.`$ (38)
Putting together these relations we have

$$tS^{\prime \prime }(\widehat{\lambda })=\frac{1}{\dot{\widehat{\lambda }}}.$$
(39)
So finally eq.(9) becomes

$$\dot{d}_2(t,ϵ)=d_0\sqrt{\dot{\widehat{\lambda }}(t)}\mathrm{exp}\left(2_1^t𝑑t^{}\widehat{\lambda }(t^{})\right),t1,$$
(40)
where we have fixed $`t_0=1`$ for convenience, and $`d_0=N(1)`$.

We recall that this formula is valid for $`t1`$, and on the assumption that the integral for $`\lambda >\lambda _0`$ is sub-dominant. We can see from (4) that this is actually true. Indeed, since the saddle point solution $`\widehat{\lambda }`$ is a decreasing function of time, it will surely become smaller that $`\lambda _0`$, for $`t`$ large enough, therefore

$$_1^t𝑑t^{}\widehat{\lambda }(t^{})\lambda _0t$$
(41)
at large times and $`(t)`$ is exponentially sub-dominant with respect to (40).

## Appendix B The plateau of the MSD

We now want to determine the value of the plateau. First, we need to appropriately define what we mean with plateau. We have seen that $`\dot{d}_2(t)`$ has a minimum for $`tt_ϵ`$. This means that $`d_2(t_ϵ)`$ is as close as possible to the limiting value where $`d_2(t)`$ would go if there were no negative eigenvalues. Thus, we may define an effective plateau for $`ϵ0`$ (but very small) as

$$q_ϵ=d_2(t_ϵ).$$
(42)
In this way for $`ϵ0`$ $`q_ϵ`$ approaches a well defined value $`q(T_c)`$ related to relaxation in a harmonic well. We now show that

$$q_ϵ=_{1/t_ϵ}^{\mathrm{}}\frac{\rho (\lambda )}{\lambda }+\mathrm{O}\left[\rho (1/t_ϵ)\right]+\mathrm{O}[\rho (1/t_ϵ)],$$
(43)
and thus

$$q_ϵq(T_c)=_0^{\mathrm{}}\frac{\rho (\lambda )}{\lambda }\text{for }ϵ0,$$
(44)
as expected. To prove (43) we proceed as follows. We split $`d_2(t)`$ into different parts

$`d_2(t)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\lambda \rho (\lambda ){\displaystyle \frac{1e^{2\lambda t}}{\lambda }}A+B+C+D+E+F=`$ (45)
$`{\displaystyle _{1/t_ϵ}^{\mathrm{}}}{\displaystyle \frac{\rho (\lambda )}{\lambda }}𝑑\lambda {\displaystyle _{1/t_ϵ}^{\mathrm{}}}{\displaystyle \frac{\rho (\lambda )e^{2\lambda t}}{\lambda }}𝑑\lambda +{\displaystyle _0^{1/t_ϵ}}\rho (\lambda ){\displaystyle \frac{1e^{2\lambda t}}{\lambda }}𝑑\lambda +`$
$`{\displaystyle _{\mathrm{}}^{1/t_ϵ}}{\displaystyle \frac{\rho (\lambda )}{\lambda }}𝑑\lambda {\displaystyle _{\mathrm{}}^{1/t_ϵ}}{\displaystyle \frac{\rho (\lambda )e^{\lambda t}}{\lambda }}𝑑\lambda +{\displaystyle _{1/t_ϵ}^0}\rho (\lambda ){\displaystyle \frac{1e^{2\lambda t}}{\lambda }}𝑑\lambda ,`$
and we evaluate the different integrals. By using the saddle point method for $`t>t_ϵ`$ we find

$$B\frac{t}{t_ϵ}\rho (1/t_ϵ)e^{2t/t_ϵ}B\rho (1/t_ϵ)\text{for }tt_ϵ.$$
(46)
C is easy:

$$C\rho (1/t_ϵ)_0^{t/t_ϵ}𝑑x\frac{1e^{2x}}{x}Cc\rho (1/t_ϵ)\text{for }tt_ϵ,$$
(47)
where $`c`$ is a constant prefactor. To evaluate D we note that to ensure the convergence of $`d_2(t)`$ the function $`f(\lambda )=\rho (\lambda )e^{2\lambda t}`$ should go to zero for $`\lambda \mathrm{}`$. We can thus estimate D with the saddle point method:

$$D=_{\mathrm{}}^{1/t_ϵ}f(\lambda )\frac{e^{2\lambda t}}{\lambda }𝑑\lambda \frac{t_ϵ}{t}f(1/t_ϵ)e^{2t/t_ϵ}D\rho (1/t_ϵ)\text{for }tt_ϵ.$$
(48)
The integral E is trickier. Once again we use the saddle point method. The procedure is the same as the one used for $`\dot{d}_2(t)`$. Since the factor $`1/\lambda `$ is algebraic it can be discarded in the evaluation of the maximum so that the maximum of the integrand still occurs for a value $`\widehat{\lambda }(t)`$ as given by equation (10). But, by definition, $`\widehat{\lambda }(t_ϵ)=0`$. So for $`t`$ close enough to $`t_ϵ`$ the real maximum lies outside the integration range and the integral is dominated by the upper integration limit. Thus we have

$$E\frac{t_ϵ}{t+\frac{\rho ^{}(1/t_ϵ)}{\rho (1/t_ϵ)}}e^{2t/t_ϵ}\rho (1/t_ϵ)E\rho (1/t_ϵ)\text{for }tt_ϵ.$$
(49)
Finally, exactly as for C we have

$$F\rho (1/t_ϵ)_{t/t_ϵ}^0𝑑x\frac{1e^x}{x}F\rho (1/t_ϵ)\text{for }tt_ϵ.$$
(50)
## Appendix C The prefactor of eq. (14)

We now argue that the prefactor of the square root of eq. (14) does not go to zero exponentially for $`t\mathrm{}`$, at least for physical spectra shapes. Let us suppose that for $`t\mathrm{}`$ we have

$$\dot{\widehat{\lambda }}(t)Ae^{at}.$$
(51)
This implies that

$$\widehat{\lambda }(t)\lambda _{\mathrm{}}t\mathrm{},$$
(52)
that is, if we define $`g(\lambda )=\frac{\rho ^{}(\lambda )}{2\rho (\lambda )}`$,

$$g(\lambda )\mathrm{}\lambda \lambda _{\mathrm{}},$$
(53)
with $`\lambda _{\mathrm{}}<0`$. This is what happens when $`\rho `$ has a cut in the left tail (semicircle). Condition (51) implies (we integrate between $`t`$ and $`\mathrm{}`$)

$$\widehat{\lambda }(t)\lambda _{\mathrm{}}+A/ae^{at}.$$
(54)
Given that $`g(\lambda )`$ is a decreasing function of $`\lambda `$, using the saddle point equation we find

$$2t=g(\widehat{\lambda })g\left(\lambda _{\mathrm{}}+A/ae^{at}\right).$$
(55)
Calling $`x\lambda _{\mathrm{}}+A/ae^{at}`$, and for $`x\lambda _{\mathrm{}}`$, this relation gives

$$g(x)\frac{1}{a}\mathrm{log}(x\lambda _{\mathrm{}}).$$
(56)
If we now recall that $`g(x)=\rho ^{}(x)/\rho (x)`$, and we integrate between $`\lambda `$ and $`0`$ and take the limit $`\lambda \lambda _{\mathrm{}}`$, we find

$$\rho (\lambda _{\mathrm{}})>0.$$
(57)
Thus, if equation (51) holds, the spectrum must be nonzero at the left cut (the converse is not true). Therefore, if we discard this unphysical case, we are sure that for $`t\mathrm{}`$ the Gaussian fluctuations do not go to zero exponentially. It is possible that for spectra of the kind we exclude, the exponential factor increases faster than a simple exponential and therefore still kills the prefactor contribution, but we have not been able to prove under what conditions this happens.


text:
# Tip Oscillation of Dendritic Patterns in a Phase Field Model

## 1 Introduction

Crystal growth has been intensively studied as a problem of pattern formation far from equilibrium. Dendrites are typical growth patterns found in crystal growth exhibited by most metal alloys and some plastic crystals. The phase field model, which has the form of the Ginzburg-Landau equation coupled with a diffusion equation, is a useful model for such growth patterns. For melt growth, the Ginzburg-Landau equation represents the dynamics of a phase transition from a liquid phase to a solid phase, and it is coupled with a heat conduction equation for the latent heat generated at the growing interface.

The standard theory of crystal growth predicts that a needle crystal selected by the solvability condition is stable, and side branches are produced owing to thermal noise. However, there are some experiments suggesting that side branches may appear as a result of regular tip oscillations. Liu and Goldenfeld investigated the linear stability of a steady needle crystal using the boundary layer model and found an instability in a certain range of parameter values of the surface tension and kinetic coefficients. Ihle numerically investigated the competition between the kinetic and surface tension anisotropy in dendritic growth using a fully dynamical front-tracking method and found oscillatory growth. However, the numerical results do not seem to exhibit clear oscillation, probably owing to numerical noise resulting from the discretization. Such numerical noise effects seem to be small in the phase field model. We have numerically investigated the competition between the kinetic and surface tension anisotropy using a phase field model and found that the growth direction changes from $`10`$ to $`11`$ as the kinetic coefficient is changed. This is a transition from a surface tension dendrite to a kinetic dendrite. Such transitions have been observed in experiments of viscous fingering and NH4Cl crystal. In a previous work, we found that complicated patterns grow from small round seeds in an intermediate parameter region where the two anisotropies compete. However, we could not find clear limit cycle oscillation of the tip growth. In this work, we consider the situation in which the anisotropy parameter is changed slowly, in a stepwise manner, and find a transition from a steady needle pattern to a tip-oscillating dendrite.

## 2 Tip oscillation of dendritic patterns in the phase field model

The phase field model in two dimensions is written

$`\tau (\theta )_tp`$ $`=`$ $`\{p\lambda u(1p^2)\}(1p^2)`$
$`+_x\{W(\theta )^2_xpW(\theta )W^{}(\theta )_yp\}`$
$`+_y\{W(\theta )^2_yp+W(\theta )W^{}(\theta )_xp\},`$
$`_tu`$ $`=`$ $`D^2u+_tp/2,`$ (1)
where $`p(x,y,t)`$ is the order parameter, $`\lambda `$ is a coupling constant, and $`\tau (\theta )`$ is the time constant for the order parameter, with $`\theta =\mathrm{arctan}(_yp/_xp)`$ the angle between the normal direction of the contour of constant $`p`$ and the $`x`$ axis. The solid and liquid phase correspond, respectively, to $`p=1`$ and $`p=1`$. The normalized temperature is denoted by $`u(x,y,t)=(TT_M)/(L/C_p)`$ where $`T(x,y,t),T_M,L`$ and $`C_p`$ are, respectively, the temperature, the melting temperature, the latent heat and the specific heat. The diffusion constant for $`u`$ is denoted by $`D`$, and $`W^2(\theta )`$ is an anisotropic diffusion constant for the order parameter. The generalized Gibbs-Thomson condition (the Wilson-Frenkel formula),

$$u_i=d_0(\theta )\kappa \beta (\theta )v_n,$$
(2)
is assumed for the sharp interface model, where $`u_i`$ is the normalized temperature at the interface, $`d_0(\theta ),\kappa ,\beta (\theta )`$ and $`v_n`$ denote, respectively, the anisotropic capillary length, the interface curvature, the anisotropic kinetic coefficient and the normal interface velocity. Karma and Rappel derived the values of $`d_0(\theta )`$ and $`\beta (\theta )`$ in the sharp-interface limit to the phase field model as

$`d_0(\theta )`$ $`=`$ $`{\displaystyle \frac{I}{\lambda J}}\{W(\theta )+W\mathrm{"}(\theta )\},`$
$`\beta (\theta )`$ $`=`$ $`{\displaystyle \frac{I}{\lambda J}}{\displaystyle \frac{\tau (\theta )}{W(\theta )}}\left[1\lambda {\displaystyle \frac{W^2(\theta )}{2D\tau (\theta )}}{\displaystyle \frac{K+JF}{I}}\right],`$ (3)
where $`I=2\sqrt{2}/3,J=16/15,K=0.13604`$ and $`F=\sqrt{2}\mathrm{ln}2`$. If $`\tau (\theta )=W^2(\theta )`$ and the coupling constant $`\lambda `$ is assumed to be $`\lambda =(2ID)/(K+JF)`$, the anisotropic kinetic coefficient becomes zero, and the standard Gibbs-Thomson condition is realized. Karma and Rappel performed numerical simulations for this case, and obtained dendrites with tip velocities and tip shapes that agree well with the theory based on the solvability condition. The growth law $`v\rho ^2=`$const., where $`v`$ is the steady growth velocity of the tip and $`\rho `$ is the tip radius, is then satisfied.

We performed numerical simulations of the phase field model with an anisotropic kinetic coefficient. We found that if the kinetic effect is included, the tip radius does not become sharper and sharper, in contrast to the case in which there is no kinetic effect and the relation $`v\rho ^2`$=const. is not satisfied, when the supercooling $`\mathrm{\Delta }`$ is close to 1. In this paper, we study the competition of the surface tension anisotropy and the kinetic coefficient anisotropy. Four-fold rotational symmetry is assumed for the anisotropy: $`W(\theta )=1+e_4\mathrm{cos}(4\theta )`$ and $`\tau (\theta )=W(\theta )\{1+e_k\mathrm{cos}(4\theta )\}`$, where $`e_4`$ and $`e_k`$ denote the anisotropy parameters of the surface tension and the kinetic coefficient. The coupling constant and the diffusion constant are chosen as $`\lambda =1.9ID/(K+JF)`$ and $`D=2`$. This choice of parameter values corresponds to the anisotropy coefficients $`d_0(\theta )115e_4\mathrm{cos}(4\theta )`$ and $`\beta (\theta )0.05+(e_k0.95e_4)\mathrm{cos}(4\theta )`$ in the sharp interface limit, owing to Eq. (3).

When $`e_k`$ is larger than 0.95$`e_4`$, the $`11`$ growth becomes favorable, owing to the kinetic effect. Numerical simulations were performed using the finite difference method with $`\mathrm{\Delta }x=0.4`$ and $`\mathrm{\Delta }t=0.015`$. The system size was $`L_x\times L_y=640\times 160`$. The uniform undercooling $`u=\mathrm{\Delta }`$ was employed as the initial condition. Mirror-symmetric boundary conditions were used at $`x=0`$ and $`y=0,L_y`$ and the fixed boundary conditions $`p=1`$ and $`u=\mathrm{\Delta }`$ were used at $`x=L_x`$. The parameters $`\mathrm{\Delta }`$ and $`e_4`$ were fixed as $`\mathrm{\Delta }=0.7`$ and $`e_4=0.06`$, and the parameter $`e_k`$ for the kinetic anisotropy was changed step by step.

The process used in our numerical simulation is as follows. If the tip of the dendritic pattern reaches $`x=600`$, the growth for the parameter $`e_k`$ is stopped, and the numerical data for the order parameter and the temperature are saved. The order parameter and the temperature profiles in the tip region ($`\stackrel{~}{p}(x,y),\stackrel{~}{u}(x,y)`$ for $`x>520`$) are used as the initial conditions for the parameter $`e_k`$ at the next step, that is, $`p(x,y,0)=\stackrel{~}{p}(x520,y),u(x,y,0)=\stackrel{~}{u}(x520,y)`$.

Figure 1 displays a snapshot of a needle pattern for $`e_k=0.07`$. This needle pattern is growing in the $`10`$ direction, owing to the surface tension anisotropy. As the parameter $`e_k`$ is increased, the growth direction changes from the $`10`$ direction to the $`11`$ direction due to the kinetic effect. In region of intermediate values of the parameter $`e_k`$, a dendritic pattern with side branches appears. Figure 2(a) displays a dendritic pattern for $`e_k=0.08`$. Regular side branching is seen. The interface temperature $`u(y,t)`$ is evaluated for each $`y`$ as the temperature at the position $`(x,y)`$ where $`p(x,y)=0`$ is satisfied. For this parameter value, the average tip velocity is approximately $`v=0.265`$, and the diffusion length $`l=2D/v15`$. The system size is sufficiently larger than the diffusion length. Figure 2(a) displays a snapshot pattern at the instant that the tip reaches the position $`x=600`$, and the effect of the boundary at $`x=L_x=640`$ is not yet important. Figure 2(b) displays the interface temperature at $`y=20`$. The interface temperature and the curvature oscillate regularly at this value of the parameter. Figure 2(c) displays the amplitude of the oscillation $`|\mathrm{\Delta }u(y)|=\mathrm{max}u(y,t)\mathrm{min}u(y,t)`$ as a function of $`y`$, where $`\mathrm{max}u(y,t)`$ and $`\mathrm{min}u(y,t)`$ are, respectively, the maximum and minimum values in the oscillatory time evolution of $`u(y,t)`$. This amplitude increases monotonically, but, it does not seem to fit a curve such as $`\mathrm{exp}(cs^{1/4})`$ or $`\mathrm{exp}(cs)`$, where $`s`$ is the arc length along the interface and $`c`$ is a constant, as is expected from some theories. It seems, rather, to fit a linear curve, $`cs`$, although we do not understand the reason for this. Figure 3 displays the amplitude of the oscillation $`|\mathrm{\Delta }u(y)|`$ at $`y=20`$ as a function of $`e_k`$. As seen, there is a transition from a stationary state to oscillatory growth, and the amplitude increases continuously from zero; that is, there is no jump at the transition. The transition is therefore concluded to be a supercritical Hopf bifurcation. This limit cycle process can be interpreted as follows. The surface tension dendrite grows in the $`10`$ direction. The side branches grow in the $`11`$ direction, owing to the kinetic effect. Then, the dendritic tip becomes round, the radius of curvature of the tip region increases, and the growth velocity in the $`10`$ direction decreases. As the velocity is decreased, the kinetic effect, which is proportional to the growth velocity, becomes weak, and the surface tension effect become relatively large. Then, growth in the $`10`$ becomes dominant, and the velocity in the $`10`$ direction increases. As the growth velocity increases, the kinetic effect becomes relatively large, and the side branches grow in the $`11`$ direction again. In this way, a negative feedback acts between the tip velocity in the $`10`$ direction and the kinetic effect, and a stable limit cycle appears.

Next, we consider the transition for kinetic dendrites. If $`e_k<0`$, $`e_4<0`$ and $`e_k`$ is sufficiently smaller than $`e_4`$, a steady needle pattern grows in the $`10`$ direction, owing to the kinetic anisotropy. As $`e_k`$ remaining less than 0 is increased, the kinetic anisotropy becomes relatively weak, and the surface tension becomes dominant. Then a transition from $`10`$ growth to $`11`$ growth occurs. In the simulations, the parameters $`\mathrm{\Delta }`$ and $`e_4`$ were fixed as $`\mathrm{\Delta }=0.7`$ and $`e_4=0.06`$, and the parameter $`e_k`$ was changed stepwise. The system size was $`L_x\times L_y=640\times 120`$. Figure 4(a) displays a needle pattern for $`e_k=0.088`$. The needle pattern is marginally stable for $`e_k0.087`$, and it becomes unstable for $`e_k>0.086`$. Figure 4(b) displays the snapshot of a pattern for $`e_k=0.085`$. Side branches develop in the tip region, and the growth of the side branches does not saturate. Figure 4(c) displays the temporal evolution of the interface temperature $`u(y,t)`$ at $`y=20`$. The amplitude of the oscillation of the interface temperature grows but does not saturate. The side branches develop and they become main branches if $`L_y`$ is sufficiently large. The negative feedback does not act in this case. If the side branches grow, owing to the surface tension effect, the dendritic tip becomes more round and the growth velocity in the $`10`$ direction decreases. Then the kinetic effect becomes further weakened, the surface tension effect becomes more dominant, and the growth in the $`11`$ direction does not saturate.

Elucidation of the oscillatory growth caused by the competition between two types of anisotropies is the main result of this paper. However, we also demonstrate a transition to the tip oscillation of dendritic growth, where the surface tension anisotropy has eight-fold rotational symmetry. In a previous work, we showed the existence of tip oscillation in such a system, but have not demonstrated the transition from the needle pattern to the oscillation. The parameter dependence of the oscillatory growth is observed when $`\tau (\theta )=W^2(\theta )`$, where $`W(\theta )=1+e_8\mathrm{cos}8\theta `$. The coupling constant $`\lambda `$ is chosen to be $`\lambda =(2ID)/(K+JF)`$, so that the kinetic effect can be neglected. The numerical simulation was performed using the finite difference method with $`\mathrm{\Delta }x=0.3`$ and $`\mathrm{\Delta }t=0.01`$. The system size was $`L_x\times L_y=600\times 90`$ and the undercooling was fixed as $`\mathrm{\Delta }=0.65`$. Figure 5(a) displays a regularly oscillating dendrite for $`e_8=0.01`$. The interface temperature is seen to exhibit limit cycle oscillation. Figure 5(b) displays the amplitude of the oscillation of the interface temperature $`|\mathrm{\Delta }u(y)|`$ at $`y=18`$ as a function of $`e_8`$. A supercritical bifurcation occurs at $`e_80.045`$ for this value of the undercooling. There are few real crystals that possess eight-fold rotational symmetry. However, it is not so unusual for the surface tension to be expressed as $`d(\theta )=1+e_4\mathrm{cos}4\theta +e_8\mathrm{cos}8\theta `$, in which case, the rotational symmetry is four-fold, but the crystal has a component of eight-fold symmetry. As the parameter $`e_4(>0)`$ is increased for a fixed value of $`e_8>0`$, $`10`$ growth becomes dominant in comparison with $`11`$ growth. Tip oscillation appears for a small value of $`e_4`$, but, its amplitude decreases as $`e_4`$ is increased. Tip oscillation is not observed for $`e_4>0.016`$ when $`e_8`$ is fixed to 0.01. (The transition to non-oscillatory growth seems to be discontinuous, with a jump, although this is not entirely clear.) It appears that the competition between $`10`$ growth and $`11`$ growth may be important in the realization of tip oscillation.

## 3 Summary and discussion

We have found a deterministic transition from a needle pattern to a regularly oscillating dendrite in the phase field model. This was found in surface tension dendrites with four-fold rotational symmetry and competitive kinetic anisotropy. A similar supercritical transition was found for surface tension dendrites with eight-fold rotational symmetry without the kinetic effect. Regularly oscillating surface tension dendrites with four-fold or six-fold rotational symmetry without the kinetic effect have not yet been found in this model. Regularly oscillating dendrites with side branches possess mirror symmetry as a results of the method used in our numerical simulation. The boundary condition is not, however, essential for the oscillation. We have confirmed that a mirror-symmetric dendrite grows naturally in a simulation in which a small needle pattern including a small antisymmetric perturbation is placed in the middle of the left side of the $`L_x\times 2L_y`$ box as the initial conditions. These results are not consistent with the linear stability analysis of Liu and Goldenfeld and the simulation of Ihle. They predict oscillatory growth with antisymmetric side branches for both surface tension dendrites and kinetic dendrites. The detailed mechanism of the oscillatory growth is not yet well understood, and its elucidation is left as a future study.

## Acknowledgements

We would like to thank Professors H. Honjo and S. Ohta for valuable discussions.


text:
# 1 A fragmentation process, which is leading twist

## Acknowledgments

C.E.C. thanks the NSF for support under Grant PHY-9900657. Both of us thank the DOE for support under contract DE-AC05-84ER40150 under which the Southeastern Universities Research Association (SURA) operates the Thomas Jefferson National Accelerator Facility.