text:
# Quasi-saddles as relevant points of the potential energy surface in the dynamics of supercooled liquids

## I Introduction

The supercooled liquid regime is the interesting state of the matter which precedes the glass formation (once the crystal state is avoided). Many efforts have been devoted to better understand the relevant physical processes taking place in the supercooled state, the most important being the enormous increase of both relaxation times and inverse diffusivity by many order of magnitude upon decreasing the temperature meta\_liq . In recent years the numerical investigation of simple model liquids has provided a very useful microscopic description of the supercooled regime and a detailed picture of the structural rearrangement of atoms during the dynamic evolution.

One of the most powerful frameworks for the study of supercooled state has been the so called Potential Energy Surface (PES) description of the system deb\_nature ; wales\_sci . This approach focuses on the properties of the multi-dimensional surface of the total potential energy sampled during the time evolution of the representative point in the $`3N`$-dimensional configuration space ($`N`$ is the total number of particles). Different landscape features may control the motion of the particles angell ; sastry ; miller ; keyes\_cm ; kur ; thomas ; lanave ; heuer , and the aim is to find those characteristics of the PES that have a direct relationship with (if possible the main causes of) the emergent behavior of the relevant physical quantities. The challenge is to find the “good” landscape features that one supposes to be responsible for the interesting phenomena, or at least useful to give a clear interpretation of them. Among others, two different landscape approaches have been widely used in the last decades in this context: the Instantaneous Normal Modes (INM) approach keyes , based on the investigation of the PES very close to the instantaneous point in the configuration space during the molecular dynamics evolution of the system, and the Inherent Structures (IS) approach stillinger , based on the analysis of the minima of PES visited by the system during its dynamic evolution. In this paper we describe and revisit in a detailed way an additional approach, that has been very useful to give new insight in the analysis of the relevant processes taking place in the supercooled liquid regime. This approach focuses on the minima of the square gradient of the total potential energy, “closest” to the instantaneous points of the molecular dynamics trajectory. This approach allows one to obtain a microscopic interpretation of the relevant processes, the main result being the characterization of the dynamics above and below the Mode Coupling Temperature $`T_c`$ mct . Moreover, the analysis of these points allows one to obtain information about some relevant characteristics of the PES, of great importance to construct simplified models of the landscape.

Lets us briefly resume the two main PES approaches, the INM and th IS approach. The INM method is based on the investigation of the PES around the instantaneous configurations $`𝐫`$ ($`𝐫`$ represents the $`3N`$-dimesional vector of the representative point in the configuration space) during the molecular dynamics evolution of the system. The diffusive quantities are supposed to be related to the shape of the energy surface at $`𝐫`$, that is to say to the eigenvalues and eigenvectors of the Hessian matrix (the second derivative of the potential energy). The main hypothesis of the INM approach is that the relevant diffusive directions have to be searched among the downward curvatures (eigenvectors with negative eigenvalues). Many attempts have been devoted to extend this approach to different liquid systems and to develop a theory of the supercooled liquid state based on INM concepts. Moreover many efforts have been spent to recognize the true diffusive directions among all the downward ones, as it is known there are downward curvatures that do not correspond to diffusive directions, notably in the crystalline state. The INM negative curvature directions are classified as shoulder modes (related to anharmonicities of the PES) and double wells (with a double well shaped one-dimensional profile). Diffusive directions are finally identified as those leading to different minima gezelter  (this analysis involves the search of the minima of the PES). The temperature dependence of the diffusion coefficient has been shown to be related to the temperature dependence of the fraction of diffusive directions lanave\_1 . On the basis of this analogy, it has been conjectured a structural interpretation of the Mode Coupling Temperature $`T_c`$ as the temperature at which the number of diffusive directions goes to zero fs\_DW .

The second IS landscape approach is based on the analysis of the inherent structures, the minima of the PES. To each instantaneous configuration $`𝐫`$ one associates a inherent structure $`𝐫_{_{IS}}`$ (using, for example, a steepest descent path starting from $`𝐫`$), $`𝐫𝐫_{_{IS}}`$, partitioning the whole PES in the basins of attraction of minima. This kind of analysis has been very useful not only in the study of the diffusive directions in the supercooled regime (see the previous discussion about the INM), but also in the study of thermodynamic quantities (for example to evaluate the configurational entropy fs\_entropy ; colu ) and in the study of the out-of-equilibrium dynamics kob ; fs\_aging .

In the present work we adopt a new approach in the analysis of the PES, in some way intermediate between the previous two, in that it maintains information about the possible diffusive degrees of freedom (as the INM approach) and uses a mapping in the PES that associates instantaneous configurations $`𝐫`$ to new points in the PES (as the IS approach). In this approach noi\_selle ; cav\_selle  the minima $`𝐫__S`$ of the square gradient of the potential energy $`W=|V|^2`$ (saddles and some inflection points of the PES), reached starting from equilibrated configurations $`𝐫`$, are calculated. This procedure partitions the PES in basins of attraction of this new points: $`𝐫𝐫__S`$. The choice of $`W`$, suggested long time ago by Weber and Stillinger weber , was motivated by the fact that saddles points of $`V`$ are absolute minima of $`W`$. As pointed out recently by Doye and Wales doye\_wales  and discussed in more details in the following, the numerical minimization of $`W`$ locates mostly local minima of $`W`$, i.e. points which are not real saddles of $`V`$. For this reason in the following we refer to the local minima of $`W`$ as Quasi Saddle Points (QSP) nota . In Ref.s noi\_selle  and cav\_selle  the topological and metrical properties of such points have been used for a description of the dynamics of supercooled liquids. The investigated systems were monatomic Lennard-Jones and/or Binary Mixture Lennard-Jones. After equilibration of the system at a given temperature $`T`$, corresponding to a given instantaneous potential energy $`e`$, the minima of $`W`$ were searched and associated to the stationary points of $`V`$ (saddles). Then the energy $`e__S`$ and the number of negative curvatures $`n__S`$ at these saddle points were measured. The dependence of $`e__S`$ and $`n__S`$ on the temperature noi\_selle  or on $`e`$ cav\_selle  show interesting properties: (i) the quantity $`n__S`$ extrapolates to zero at $`T_c`$, demonstrating the validity of the conjecture that $`T_c`$ marks the transition between a dynamics among minima at low $`T`$ and a dynamics among saddles at high $`T`$; (ii) the aspects of the energy landscape seen by a given minimum is highly regular (as demonstrated by the linear dependence of $`n__S`$ on $`e__Se_{_{IS}}`$ noi\_selle  or on $`e__S`$ cav\_selle ). The previous observation leads to the conclusion that the knowledge of the saddle points properties can be used to predict the supercooled dynamics and points out the relevance of these peculiar PES points.

Following these papers, the concept of saddles \- as obtained by a minimization of $`|V|^2`$ \- has been further used in the investigation of the supercooled liquid dynamics shah ; cav\_nuovo  and aging dynamics noi\_aging .

As we anticipated above, there is a drawback in the approach used in Ref.s noi\_selle  and cav\_selle . In fact, while it is true that in a saddle point (stationary point) of the $`3N`$-dimensional surface $`V(𝐫)`$, the function $`W`$ has a minimum (actually, $`W=0`$ there and $`W`$ has a global minimum), the reverse is not guaranteed to be true. There are local minima of $`W`$, with $`W0`$, that are NOT stationary points of $`V`$. In these points, along one or more directions $`𝐫(\xi )`$, parametrized by $`\xi `$, the function $`V(\xi )`$ has an inflection point. In these inflection points (QSP), the Hessian of $`V`$ has a number of vanishing eigenvalues (excluding the three translation directions) equal to the number of directions $`𝐫(\xi )`$ where $`V`$ has an inflection point. This property has been used in noi\_selle  in order to identify the QSP: a threshold was arbitrary chosen, and all the minimizations of $`W`$ that lead to points where the Hessian of $`V`$ has at least one eigenvalue smaller (in modulus) than the threshold were discarded from the analysis. In noi\_selle  we claimed that only a small fraction of minimizations ended in a QSP points, and we estimated in $`2\%`$ the error introduced in $`n__S`$ by the choice of the threshold.

Doye and Wales doye\_wales  have shown - by the analysis of the same system investigated in noi\_selle  and cav\_selle , the BMLJ with periodic boundary conditions - that almost all (around $`95\%`$ in their work) the minimizations of $`W`$ ended in a QSP point, and not in a true saddle as claimed in noi\_selle  and cav\_selle . Motivated by such contradictory results, we reanalyzed the minima of $`W`$ found in noi\_selle . In this work we show the result of such an analysis. We conclude that the results of Doye and Wales are correct and the error in the analysis in noi\_selle  originates from the “imperfect” minimization of $`W`$: in a point close to (but not exactly at) QSP, the Hessian of $`V`$ no longer has “small” eigenvalues, and the criterion used to identify false saddles fails.

One puzzling question, however, still remains open: why the quantity $`n__S`$ \- i.e. the number of negative eigenvalues at the minima of $`W`$ \- is related to the critical temperature of the system ? It remains true that $`n__S(T)`$ extrapolates to zero at $`T_c`$. In the present paper we address this question, and we will further stress the importance of the quantity $`n__S`$; we will show that also the transport properties (namely the diffusion coefficient) can be determined by the knowledge of $`n__S`$.

The paper is organized as follow: In section II we describe the new approach, defining the saddles of the PES and studying very carefully their operative definition and the problems related to the presence of false saddles. Then in section III we expose the relevant results obtained with the new PES approach, giving a topological interpretation of the mode coupling dynamic transition. In section IV the relationship between diffusivity and saddle properties is evidenced, giving a simple interpretation of the diffusion processes in terms of saddle order. In section V we describe some landscape features as emerged from saddle analysis. In section VI a simplified model of the PES is described, together with its relevant properties compared to the features of model liquids investigated above. At the end (section VII) a brief summary of the main results and conclusions are reported.

## II Quasi Saddle Points

The analyzed system is a simple model liquid, a monatomic Modified Lennard-Jones 6-12 (MLJ). The model is able to support strong supercooling without the occurrence of crystallization, due to the presence of a small perturbation term in the Hamiltonian that inhibits ordering (this term is a function of the static structure factor, see Ref. noi\_fdt  for more details). The system is composed of $`N=256`$ particles enclosed in a cubic box with periodic boundary conditions. Truncated ($`R_c=2.6`$) and shifted LJ potential are used. Equilibrium configuration are prepared by standard microcanonical molecular dynamics simulations at constant density $`\rho =1`$ (standard LJ units are used hereafter) and at temperatures ranging from the normal liquid phase ($`T1.6`$) down to the mode coupling temperature ($`T_c=0.475`$, estimated from diffusivity).

The aim of our approach is to use a method of investigation of the PES that allows us to reject all the non-relevant degrees of freedom for the description of the slow supercooled dynamics. The idea is that there are particular points of the PES associated in some way to the instantaneous configurations, and that the properties of these points contain all the important information for the determination of the slow dynamics. Possible and natural candidates are the saddle points of the PES. In principle one can think to partition the whole configuration space in basins of attraction of saddles (the definition of the partition is obviously not unique, a possible choice could be that obtained using the Voronoi polyhedra in the $`3N`$-dimensional space). However a useful partition in basins of saddles (that is local, i.e. allows us to associate a saddle point to each instantaneous configuration using only local information of the PES) is still not available. A possible partition is that obtained using the pseudo-potential weber  $`W=|V|^2`$, associating to each instantaneous configuration during the molecular dynamics a minimum of $`W`$:

instantaneous configuration $``$ minimum of $`W`$ .

However, as already pointed out, the relation between all the minima of $`W`$ and saddles of $`V`$ is not one-to-one: absolute minima of $`W`$ (with $`W=0`$) are all saddles of $`V`$

absolute minima of $`W`$ $``$ saddles of $`V`$ ,

but local minima (with $`W>0`$) are inflection points of $`V`$ (more precisely they are inflection points with the first derivative of $`V`$ of the same sign of the first derivative of the curvature in the inflection direction):

local minima of $`W`$ $``$ inflection points of $`V`$ .

In Fig.s 1 and 2 examples of two possible one-dimensional profiles of $`V`$ and corresponding $`W`$ are shown as a function of a generic one-dimensional coordinate. In Fig. 1 the two absolute minima of $`W`$ correspond to two saddle points of $`V`$ (in one dimension a minimum and a maximum). The arrows define the basins of attraction of minima of $`W`$, and the boundaries of these basins are indicated with open symbols. In Fig. 2 the local minimum of $`W`$ corresponds to an inflection point of $`V`$.

We note that at inflection points the eigenvalues of the Hessian vanish, so in principle one can recognize them simply analyzing the eigenvalues of the Hessian of $`V`$. In a previous work noi\_selle  we used this eigenvalues criterion (choosing a threshold to discriminate inflection directions) and we found that local minima of $`W`$ were found very rarely. A closer inspection, however, has revealed that, due to “imperfect” minimizations of $`W`$, the eigenvalues associated to a point “close” (due to the not perfect minimization) to an inflection point were of the same order of magnitude of the eigenvalues associated to a saddle. Refining the analysis of the minima of $`W`$, we find that the number of points that are true saddles is very small (about $`1`$% of the total), so almost all the minima of $`W`$ found are inflection points. However there is another interesting characteristic that has to be considered, the number of inflection directions (number of zero eigenvalues of Hessian, excluding the tranlations) at the minima of $`W`$: if this number is small one can think of the minima of $`W`$ as a “quasi” saddle point QSP, as they are true saddle points in the subspace orthogonal to the small number of inflection directions. One can conjecture that the properties of a QSP are very similar to that of a true saddle point associate to an instantaneous configuration; however, a relationship between them has still to be proved. In the simulated MLJ system the number of inflection directions $`n_0`$ with respect to the number of negative curvatures $`n__S`$ are reported in Table 1 for the different investigated temperatures. The values of $`n_0`$ (from $`1`$ to $`4`$) indicate that they correspond to few directions in the configuration space; a small value that allows us to speak about QSP as good candidates to approximate the properties of the true saddle points. We note that the fraction $`f=n_0/n__S`$ is higher for the low temperature data (where the quantity $`n__S`$ is small), evidencing how these points are more influenced by inflection directions. However, also in this case it is possible to estimate the error in the calculation of the order (see Fig. 5 in the next section), evidencing the robustness of the obtained results.

A natural question now arises: what are the relevant features of the found minima of $`W`$ that are useful to describe the long time dynamics of the supercooled liquids analyzed ? In the spirit of our early simple conjecture that the relevant quantity is related to the number of negative curvatures at saddle points, we speculate that the relevant information is in the negative curvatures and not in the few inflection directions. So, no matter if they are true or quasi saddle points, we are going to check if their properties are able to describe and to give a microscopic interpretation of the dynamic processes in the supercooled regime.

## III Dynamic transition

We firstly analyze the temperature behavior of some quantities related to quasi saddle points, and compare them to the corresponding behavior of the same quantities calculated at instantaneous configurations and at inherent structures.

For each temperature (ranging from $`T1.6`$ to $`TT_c`$), we analyze $`20`$ independent equilibrium configurations. For each configuration we calculate the associated IS, implementing a steepest descent algorithm which moves in the direction of $`V=𝐅`$, and the associated minimum of $`W`$, moving in the direction of $`W=𝐅`$, where $``$ is the Hessian matrix. So, at each temperature, we have three different points of the configuration space: the instantaneous configuration $`𝐫`$, the inherent structure $`𝐫_{_{IS}}`$, and the quasi saddle configuration $`𝐫__S`$.

In Fig. 3 the temperature dependence of the energies of the different points are shown: $`e`$, $`e_{_{IS}}`$ and $`e__S`$.

The value of the mode coupling temperature $`T_c=0.475`$ (obtained by diffusivity data, see section IV) is also indicated. The energy of the inherent structures shows a constant behavior down to $`T0.8`$, below which it starts to decrease abruptly (note however the small variation range with respect to that of the $`e__S`$). Below this temperature the system starts to visit basins of minima of lower and lower energy. The quantity $`e__S`$ shows a different behavior, decreasing gradually with temperature and approaching $`e_{_{IS}}`$ close to $`T_c`$. The fact that $`e__S`$ lies well below $`e`$ indicates that the process of minimization of the pseudo-potential $`W`$ consists in a downward path in the PES along the majority of the degrees of freedom. The coincidence of $`e_{_{IS}}`$ and $`e__S`$ at $`T_c`$ suggests the hypothesis that the sampled saddles start to be mainly minima at this temperature. The answer to this question is in the behavior of the number of negative curvatures as a function of temperature. In Fig. 4 the fraction of negative curvatures (number of negative curvatures over the total degrees of freedom $`3N`$) is shown for the analyzed configurations $`𝐫`$ ($`n`$) and $`𝐫__S`$ ($`n__S`$) (obviously this number is zero for the inherent structures $`𝐫_{_{IS}}`$).

The first thing to observe is that the quantity $`n__S`$ lies well below the instantaneous order $`n`$, indicating that the process of minimization of $`W`$ leads to a point well below in the PES, with a number of negative curvatures less than that at the equilibrium starting point. We have then found a direct method to reject many negative eigenvalues (corresponding to not diffusive directions) of the Hessian at the instantaneous configurations. If the remaining negative curvature directions are those useful to describe the slow dynamics of the system, will be soon clear. A first good indication of the relevance of the quantity $`n__S`$ for the description of the dynamics is the fact that its well defined temperature behavior extrapolates to zero very close to the $`T_c`$ value. In Fig. 5 the temperature behavior of $`n__S(T)`$ is shown in a larger scale. It is shown the $`n__S`$ calculated in a direct way (as the fraction number of negative curvatures at the minima of $`W`$) and an upper and lower estimation of the “true” saddle order, considering all the inflection directions as if they contributed to the saddle order (upper points) and as if they didn’t contribute to it (lower points). Fitting the three different data with a power law $`(TT_0)^\gamma `$, we obtain very similar value of $`T_0T_c`$, indicating the robustness of the analysis.

The sampled minima of $`W`$ are then QSP with a number of negative curvatures that decreases lowering the temperature, and approaches zero at the temperature of the dynamic mode coupling transition. At $`T_c`$ and below it, $`T<T_c`$, the sampled saddles are in average minima of $`V`$, i.e. saddles with order zero. The emergent scenario for the dynamics is then the following:

* for $`T>T_c`$ the system lies with high probability close to borders of basins of attraction of inherent structures,
* for $`T<T_c`$ the system spends most of the time trapped in the basins of inherent structures.

These findings confirm the conjecture that $`T_c`$ marks a dynamics crossover, evidencing a possible structural PES interpretation of this change in the dynamics.

We conclude this section showing the distributions of the quasi saddle order for different temperatures. In Fig. 6 the distribution $`P_n`$ of saddle order $`n`$ (now we indicate with $`n`$ the absolute saddle order, i.e. the total number of negative curvatures, not the fraction) is shown for four different temperatures. A simple conjecture allow us to obtain a prediction for these distributions. Assuming the independence of the relevant slow degrees of freedom and a simple one dimensional profile of the independent coordinate (with only a minimum and a maximum) we have that the distribution of the saddle order $`n`$ of $`N`$ degrees of freedom is simply the binomial distribution $`P_n=\left(\genfrac{}{}{0pt}{}{N}{n}\right)p^n(1p)^{Nn}`$, where $`p`$ is the probability that the single degree of freedom is close to (inside the basin of attraction of) the maximum. In the large $`N`$ and small $`p`$ limit we obtain a Poisson distribution

$$P_n\frac{<n>^n}{\mathrm{\Gamma }(n+1)}e^{<n>},$$
(1)
written in term of the $`\mathrm{\Gamma }`$ function ($`\mathrm{\Gamma }(n+1)=n!`$) and of the mean value of the saddle order $`<n>=pN`$ (the temperature dependence is in the $`p`$ parameter - small $`p`$ means low temperature). Lines in Fig. 6 are the distributions obtained from Eq. 1, with the $`<n>`$ value fixed by the data. They well describe the measured distributions, suggesting the correctness of the simple viewpoint of independent degrees of freedom (a simple model with this ingredient will be study in detail in the last section).

## IV Diffusion

As demonstrated in the previous section, some features of the saddles probed during the dynamic evolution of the system are able to capture the dynamic crossover at $`T_c`$. In order to corroborate this finding we now go to analyze the diffusion properties in the supercooled regime and to evidence a possible relationship with the saddle order $`n__S`$.

For different temperature above $`T_c`$ we have calculated the mean square displacement:

$$R^2(t)=\frac{1}{N}\underset{i}{}<|\stackrel{}{r}_i(t)\stackrel{}{r}_i(0)|^2>,$$
(2)
and, from its long time behavior, the diffusion coefficient $`D`$:

$$D=\underset{t\mathrm{}}{lim}\frac{R^2(t)}{6t}.$$
(3)
In Fig. 7 the diffusion coefficient $`D`$ is shown as a function of temperature. The mode coupling theory predicts a power law divergence of the inverse diffusivity close to $`T_c`$, so the behavior of the diffusion coefficient is a method to estimate $`T_c`$. Fitting the data with power law (line in Fig. 7)

$$D(T)(TT_c)^\gamma ,$$
(4)
we obtain the value $`T_c=0.475`$.

The fact that also the saddle order $`n__S`$ extrapolates to zero at $`T_c`$, suggests the possibility of a relationship between $`n__S`$ and the diffusivity $`D`$. A direct prediction of a relation between the two quantities arises from a simple interpretation of the dynamic processes.

Supposing that the slow degrees of freedom responsible of the diffusion processes are in some way related to the number of negative curvatures of the probed saddle points at a given temperature, one can try to predict a relationship between the saddle order $`n__S`$ and the diffusion coefficient. Assuming that the number of diffusive directions at a given temperature are proportional to the saddle order at the same temperature, we can think of the diffusive process as a random walk process in a space of dimension equal (or proportional) the the saddle order $`n__S`$. In this case one expects that the diffusion coefficient is proportional to the temperature times the saddle order:

$$D(T)Tn__S(T).$$
(5)
We note that this hypothesis doesn’t mean that in general all the negative curvatures at a generic saddle point are related to a slow diffusive process, as there are saddles related to small local rearrangements of atoms. It means only that the process of minimization of $`W`$, starting from an equilibrated configuration, leads in average to a point (saddles or quasi-saddle doesn’t matter) with negative curvature directions in some way related to the true diffusive directions.

In Fig. 8 the diffusivity over the temperature calculated through molecular dynamics simulations and the saddle order $`n__S`$ are shown as a function of temperature: it seems that the simple random walk model on a reduced space works well. Eliminating the temperature and plotting the two quantities one against the other in a double log scale (inset of Fig. 8) we see that the Eq. 5 is well suited, except the last low temperature point, which is however the most influenced by a possible not good equilibration procedure. This kind of analysis deserves surely further investigations, for example analyzing in a more direct way the paths during the molecular dynamics evolution and their relation with the negative curvature directions at saddle points. However one result emerges in a strong way: the quantity $`n__S`$ seems to be very useful to reproduce dynamic processes and to give a possible PES interpretation of them.

## V Landscape features

The previous analysis of the PES allows us to infer some topological characteristics of the landscape itself. The first feature of the PES we analyze is the energy organization of the saddles, i.e. a possible relation between energy and order of saddle points. Plotting (see Fig. 9a) the saddle energy $`e__S`$ as a function of the saddle order $`n__S`$ (using all the saddles found at different temperatures), as done in cav\_selle , i.e.

$$e__S(n__S)\mathrm{\Delta }e^{}n__S,$$
(6)
the data can be reproduced with $`\mathrm{\Delta }e^{}=13.5`$. However, data at low energy (and order) deviate from the linear fit. Observing that these points are those sampled in the temperature range where the inherent structures start to decrease in their energy values, one is attempted to relate the deviation from the linear fit of the energy-order relation to the inherent structure behavior.

In Fig. 9b the quantity plotted against the saddle order $`n__S`$ is the elevation energy $`e__Se_{_{IS}}`$ of saddles with respect to the corresponding local minima (minima visited at the same temperatures), obtaining now a remarkable linear relationship in the full temperature range:

$$e__S(n__S)e_{_{IS}}=\mathrm{\Delta }en__S,$$
(7)
with $`\mathrm{\Delta }e=13.3`$. This suggests that the energy landscape above a given minimum is organized in families of equally spaced energy saddle points, with only a single energy barrier parameter (in average) $`\mathrm{\Delta }e`$, that represents the energy gap between a saddle of order $`n__S`$ and a saddle of order $`n__S+1`$. A possible explanation of this result is the following: saddles of a given order are obtained from a combination of independent saddles of order $`1`$, that lie a fixed amount over minima. One can think of the excitations as as a “gas” of non-interacting degrees of freedom (in the next section we develop this viewpoint in a deeper way, introducing a simple landscape model). Similar results are found in the out-of-equilibrium regime noi\_aging , suggesting that the PES properties obtained are not influenced by the kind of dynamics the system use to explore its PES.

Another interesting topological information about the PES is obtained calculating the distance relations between adjacent inherent structures. For each saddle point we have perturbed the system along a randomly chosen negative curvature direction and then started a minimization procedure of the potential $`V`$ in order to find the underlying minimum. The same minimization has been performed perturbing the system along the previous negative curvature direction but in the opposite versus. In this way we obtain couples of inherent structures associated to each negative direction (we call them adjacent inherent structures). We then evaluated the Euclidean distance between adjacent IS, defined as

$$d_{_{IS}}^2=\frac{1}{N}\underset{i}{}<|\stackrel{}{r}_i(IS_1)\stackrel{}{r}_i(IS_2)|^2>,$$
(8)
where $`IS_1`$ and $`IS_2`$ are two adjacent minima and the average is over different couples of adjacent IS.

In Fig. 10 the quantity $`d_{_{IS}}^2`$ is shown as a function of the order $`n__S`$ of the starting saddle. The observed linear relation between the two quantities

$$d_{_{IS}}^2(n__S)=d_0^2n__S,$$
(9)
(with $`d_0^2=9.6`$) suggests a simple topological interpretation: the descent path from a saddle of order $`n__S`$ towards the underlying minima can be view as a sequence of independent random steps, each of them decreasing the saddle order by $`1`$ and increasing the mean square distance between the underlying local minima by a constant amount $`d_0^2`$.

The above results seem to indicate that in some aspects the PES exhibits a very simple and organized structure. This suggests the possible use of simple solvable mathematical models in order to capture some of the relevant features of the PES explored in the supercooled regime.

## VI Trigonometric Model

A very simplified model that is able to capture some of the regular characteristics of the PES of liquids evidenced above, is the so called Trigonometric Model (TM), introduced by Madan and Keyes madan\_keyes . In this section we calculate for the TM the behavior of the quantities previously analyzed for the Lennard-Jones system.

The TM is a model for $`N`$ independent degrees of freedom with Hamiltonian:

$$H_{_{TM}}=\mathrm{\Delta }\underset{i}{}[1\mathrm{cos}(\phi _i)],$$
(10)
where $`\{\phi _i\}`$ are angular variables: $`\phi _i[0,2\pi )`$. The PES of TM reproduce the regularity of the average saddle properties of the Lennard-Jones PES: saddles of order $`n__S`$ are $`2\mathrm{\Delta }`$ over (in energy) and $`\pi `$ near (in distance) saddles of order $`n__S1`$. The thermodynamics is easy computable, as the partition function is factorized:

$$Z_{_{TM}}(\beta )=Z_0^N(\beta ),$$
(11)
where $`\beta =T^1`$ (we use the unit $`K_B=1`$) and

$$Z_0(\beta )=_0^{2\pi }𝑑\phi e^{\beta \mathrm{\Delta }[1cos(\phi )]}=2\pi e^{\beta \mathrm{\Delta }}I_0(\beta \mathrm{\Delta }),$$
(12)
where $`I_0`$ is the Bessel function of order zero. The energy density $`e_{_{TM}}=N^1_\beta \mathrm{log}(Z_{_{TM}})`$ is easy written in term of $`I_0`$ and $`I_1`$ (the Bessel function of order $`1`$, $`I_1(x)=I_0^{}(x)`$):

$$e_{_{TM}}(\beta )=\mathrm{\Delta }\left[1\frac{I_1(\beta \mathrm{\Delta })}{I_0(\beta \mathrm{\Delta })}\right].$$
(13)
In order to test the reliability of TM in the description of liquid behavior, we now calculate the diffusion coefficient and the saddle order as a function of temperature. The diffusion coefficient $`D_{_{TM}}`$ is exactly computable calculating the mobility in linear response risken , with the introduction of a Langevin dynamics:

$$\gamma \dot{\phi _i}(t)=F_iH_{_{TM}}+\eta (t),$$
(14)
where $`F`$ is an external force, $`\gamma `$ is the friction, $`\eta `$ is a random variable with mean zero and $`\delta `$-correlate in time: $`<\eta (t)\eta (t^{})>=(2\gamma /\beta )\delta (tt^{})`$ (we choose $`m=1`$ mass unit). Following risken  we obtain:

$$D_{_{TM}}=\frac{1}{\gamma \beta I_0^2(\mathrm{\Delta }\beta )}.$$
(15)
Using the expansion of $`I_0`$ for large argument ( $`I_0(x)\mathrm{exp}(x)/\sqrt{x}`$, for large $`x`$) we obtain the Arrhenius behavior for the diffusivity at low temperature:

$$D_{_{TM}}e^{2\beta \mathrm{\Delta }}\mathrm{for}\mathrm{low}T$$
(16)
(we note the value $`2\mathrm{\Delta }`$ is the energy barrier value for the single degree of freedom). The dimensionless quantity

$$\frac{D_{_{TM}}}{D_0}=I_0^2(\mathrm{\Delta }\beta ),$$
(17)
where $`D_0=(\gamma \beta )^1`$ is the free Brownian diffusion coefficient, is shown in Fig. 11 as a function of temperature (line).

An interesting observation arises from the comparison between this theoretical canonical result and a numerical iso-potential (dynamics at constant potential energy) calculation. The points in Fig. 11 are obtained from a simulation of $`N=1000`$ variables constrained to move at constant potential energy (the dynamics is a random walk dynamics in which the forces $`_iH`$ are projected to the constant potential energy surface). The coincidence between analytical and numerical data evidences how canonical and iso-potential calculations lead to the same equilibrium dynamic properties. In Fig. 12 the diffusion coefficient $`D_{_{TM}}`$ is shown as a function of temperature in a semi-logarithmic scale. We note that, though $`D_{_{TM}}`$ is a smooth function, one can define a change in the behavior extrapolating a critical temperature $`T_c`$, using a power law fit in the low temperature range (from $`0.5`$ to $`1.5`$):

$$D_{_{TM}}(T)(TT_c)^\gamma ,$$
(18)
obtaining the values: $`T_c/\mathrm{\Delta }=0.38`$ and $`\gamma =1.2`$. Below $`T_c`$ the diffusion coefficient can be approximated by Eq. 16, a low T Arrhenius behavior with energy barrier value $`2\mathrm{\Delta }`$ (see the inset in Fig. 12).

The saddle order $`n_{_{TM}}`$ is also computable in an exact way. Defining it as the negative curvatures at minima of $`W=|H|^2`$ (we note that, in the TM, all the minima of $`W`$ are true saddles of $`H`$ and, moreover, the instantaneous order coincides with the saddle order) we can calculate the equilibrium mean value of $`n_{_{TM}}`$ as the probability that the variable $`\phi `$ is in the range $`[\pi /2,3\pi /2]`$, where the $`H`$ profile as a negative curvature:

$`n_{_{TM}}`$ $`=`$ $`Z_0^1{\displaystyle _{\pi /2}^{3\pi /2}}𝑑\phi e^{\beta \mathrm{\Delta }[1cos(\phi )]}`$ (19)
$`=`$ $`[2\pi I_0(\beta \mathrm{\Delta })]^1{\displaystyle _{\pi /2}^{3\pi /2}}𝑑\phi e^{\beta \mathrm{\Delta }cos(\phi )}.`$
The temperature behavior of $`n_{_{TM}}`$ is shown in Fig. 13. Also in this case it is possible to fit the data with a power law fit close to $`T_c`$ (from $`0.5`$ to $`1.5`$), obtaining a value of $`T_c/\mathrm{\Delta }=0.41`$, close to that obtained from diffusivity ($`T_c/\mathrm{\Delta }=0.38`$).

Having calculated the temperature behavior of diffusivity $`D_{_{TM}}(\beta )`$ and saddle order $`n_{_{TM}}(\beta )`$, we can now try to study the relationship between them. Due to the fact that a direct analytical expression of $`D_{_{TM}}`$ as a function of $`n_{_{TM}}`$ is not available, we study it numerically eliminating the temperature parameter. The relationship is shown in the inset of Fig. 13, where the quantity $`D_{_{TM}}/D_0`$ is plotted against $`n_{_{TM}}`$ in a logarithmic scale (this plot corresponds to that of the Lennard-Jones system \- inset of Fig. 8 \- as the quantity $`D_{_{TM}}/D_0`$ is proportional to $`\beta D_{_{TM}}`$). The investigated $`T`$ range is close to the critical temperature $`T_c`$ (as in the liquid system) and also in this case we find an approximatively linear relationship (the dashed line in the inset of Fig. 13 corresponds to a power law with exponent $`1`$):

$$D_{_{TM}}Tn_{_{TM}}.$$
(20)
In conclusion, the very simple TM is able to reproduce some of the properties emerged in the analysis of the Lennard-Jones liquid system, indicating that the use of simplified analytical model of the PES is a promising field of research.

## VII Conclusions

The analysis of the PES has been very useful in the investigation of supercooled dynamics of simple model liquids. Besides the two main approaches in the investigation of the PES, the instantaneous normal mode and inherent structures approaches, we have discussed an additional sampling of the PES explored during the dynamic evolution of the system in configuration space: it consists in the analysis of the minima of the square gradient of the potential energy, $`W=|V|^2`$. The characteristics of these sampled points seem to have a very close relationship with the transport properties of the system. The main result we obtained is the characterization of the dynamic transition temperature $`T_c`$ as the temperature at which the behavior of quantities related to the minima of $`W`$ changes on character: more specifically the number of negative eigenvalues of the Hessian at these points extrapolates to zero at $`T_c`$. The minima of $`W`$ sampled at different temperatures appear to contain the relevant information about the diffusive directions and the process of minimization of $`W`$ eliminates the non-diffusive negative directions usually present at instantaneous configurations. Recently, an interesting technique to evaluate the non-diffusive directions as been adopted in keyes\_cgf . In some previous works we called these points saddles, as the absolute minima of the function $`W`$ are true saddle points of $`V`$. However, following Doye and Wales doye\_wales  a closer inspection has revealed that very often the sampled minima during dynamic evolution are local minima of $`W`$, that correspond to inflection points of $`V`$ (the energy profile along some normal mode directions at these points is an inflection direction). The fact that the number of inflection directions at a given local minimum of $`W`$ is small, allows us to think that the main properties of these points are determined by the saddle directions (that give the value of the order, i.e. the number of negative curvatures). A further evidence of the relevance of the minima of $`W`$ is obtained analyzing the diffusivity. The temperature behavior of the number of negative curvatures is related to that of the diffusion coefficient, indicating a close relationship between them. It could be that the main information about the dynamic processes is in the saddle points and the relevant properties of the local minima of $`W`$ are good approximations of the properties of these true saddle points. But there is also the possibility that the true saddles are not so relevant in the description of the dynamics, at least not more relevant than points with inflection directions (local minima of $`W`$).

All these findings lead to the following viewpoint about the supercooled dynamics: the relevant diffusion processes happen in a subspace of the total configuration space with dimension proportional to the number of saddle directions at the sampled minima of $`W`$, and, due to the presence of these “open” directions, the relevant diffusive processes are not activated in energy. Below $`T_c`$ things are different, as the number of negative directions at sampled minima of $`W`$ is vanishingly small, indicating that the system is with high probability in a basin of an inherent structure (minimum of $`V`$). In order to change basin and diffuse, now the system has to find the “good” directions leading to other minima (an entropic process, considering the whole $`3N`$ dimensional space of the degrees of freedom, that could correspond to an energy activated process in a suitable $`n`$ dimensional subspace - with $`n3N`$).

Another important result obtained refers to the structural organization of the PES, i.e. the relationship among the relevant points of it. A very simple structure of the PES emerges: saddles of order $`n+1`$ lie above saddles of order $`n`$ by a constant quantity $`\mathrm{\Delta }e`$ and adjacent inherent structures are equally spaced in configuration space. The above observed regularities suggest the possible use of simplified models of the PES in order to capture relevant characteristics of the supercooled dynamics of liquids. The model we analyzed here is the Trigonometric Model, introduced by Madan and Keyes madan\_keyes , a model of independent sinusoidal degrees of freedom. The behavior of the calculated diffusivity and saddle order seems to indicate that the model, despite the absence of cooperativity, is able to reproduce important aspects of liquid dynamics close to $`T_c`$. The study of more complex and realistic PES models, with interactions among the degrees of freedom, is a interesting development of this kind of analysis.

###### Acknowledgements.

We acknowledge financial support from INFM - Iniziativa Calcolo Parallelo. We thank S. Ciuchi, C. De Michele and D.J. Wales for usefull discussions.


format: text

text:
# Macroscopic spin tunneling and quantum critical behavior of a condensate in double-well potential

## Abstract

In a previous work pu , we have shown that a spinor condensate confined in a periodic or double-well potential exhibits ferromagnetic behavior due to the magnetic dipole-dipole interactions between different wells, and in the absence of external magnetic field, the ground state has a two-fold degeneracy. In this work, we demonstrate the possibility of observing macroscopic quantum spin tunneling between these two degenerate states and show how the tunneling rate critically depends on the strength of the transverse field.

Tunneling, a process in which a system penetrates into a classically forbidden region (e.g., a potential barrier), is an intrinsically quantum effect with no classical counterpart. Of all tunneling effects, macroscopic quantum tunneling (MQT), the tunneling of a macroscopic variable of a macroscopic system leggett , represents a particularly intriguing and interesting scenario as it touches the boundary between the classical and the quantum world and may help shed light on the quantum-classical interface. As such, it is one of the most striking manifestations of quantum mechanics. The tunneling of large magnetic spins has recently received much attention, both theoretically and experimentally book ; qtm94  in view of its promise as one of the few realistic candidates for an experimental demonstration of MQT, and also because of its connection to quantum computing. Despite considerable efforts, though, there are not many clear and definitive demonstrations of MQT in spin systems so far — although there have been some indications that tunneling might be the underlying reason for some observed results qtm94 ; exp . The reasons for this difficulty are as follows: First, most theories apply to single system, while the conventional magnetic materials used in the experiments contain many domains, each possessing its own set of parameters such as magnetic anisotropy and barrier energy; second, due to the difficulty of cooling the samples down to ultracold temperatures, thermal processes cannot be completely excluded; and third, the spins in solid materials are inevitably imbedded in a crystal matrix and the spin-matrix interaction matrix  complicates the physical picture. For these reasons, macroscopic spin tunneling remains one of the most anticipated, yet elusive, quantum phenomena.

In this paper, we demonstrate the possibility of observing spin tunneling in a spinor atomic condensate trapped in a double-well potential, thereby eliminating most of these difficulties. We remark at the outset that while the inter-well tunneling of condensates has been previously considered, all previous studies focused instead on the tunneling of an external degree of freedom, the condensate center-of-mass motion milburn ; smerzi ; jessen ; zhang .

In an earlier paper pu , we showed that because of the long range magnetic dipole-dipole interaction, a spinor atomic condensate trapped in a one-dimensional periodic lattice potential or a double-well potential behaves as a ferromagnet. In the absence of external magnetic fields, the ground state spins of the “mini-condensates” confined in individual wells all align parallel to each other and along the lattice direction (the $`z`$-axis), giving rise to a spontaneous magnetization along $`z`$. For the sake of simplicity, we restrict the present discussion to a double-well potential, each well containing $`N`$ spin-1 condensate atoms. In the tight-binding approximation milburn , they are described by the zero-temperature spin Hamiltonian pu

$$H=\lambda _a^{}(𝐒_1^2+𝐒_2^2)3\lambda S_1^zS_2^z+\lambda 𝐒_1𝐒_2h(S_1^x+S_2^x),$$
(1)
where $`𝐒_i`$ is the total spin of the condensate in the $`i^{\mathrm{th}}`$ well, and $`S_i^{x,y,z}`$ its cartesian components. The first term in (1) represents the on-site Hamiltonian of the spinor condensate. It includes short-range nonlinear spin-exchange interactions law , where the parameter $`\lambda _a^{}`$ is related to the $`s`$-wave scattering lengths ho  and needs to be negative. The second and third terms in (1) arise from the site-to-site dipole-dipole interaction pu , where $`\lambda \gamma _B^2\mu _0/(4\pi r^3)`$ for pure magnetic dipolar interaction, with $`\gamma _B`$ being the gyromagnetic ratio, $`\mu _0`$ the vacuum permeability and $`r`$ the distance between the two wells. The value of $`\lambda `$ can be greatly enhanced by the light-induced optical dipolar interaction if one chooses appropriate laser fields to form the potential well spinwave ; optical . The last term describes the effect of an external transverse magnetic field, taken to be along the $`x`$-axis without the loss of generality. Here $`h=\gamma _BB`$, with $`B`$ being the strength of the applied field. In the following, we assume that the nonlinear short-range atom-atom interaction is strong enough that the first term in Hamiltonian (1) dominates over the magnetic dipolar interaction note . As a result, the total spin quantum number for each mini-condensate is fixed to its maximum value $`N`$ — the number of particles in each well. We can therefore neglect the first term in the Hamiltonian (1), since it is a constant of motion and commutes with the remaining terms.

Before discussing quantum mechanical spin tunneling, let us first investigate the classical situation. The Hamiltonian in that limit is still given by (1), except that the spins are now $`c`$-numbers that can be represented by vectors of fixed length $`N`$ in spin space. For zero field, $`h=0`$, it is easy to see that the classical ground state is two-fold degenerate with $`𝐒=𝐒_1=𝐒_2`$ pointing along either the $`z`$\- or $`(z)`$-axis. Under the influence of a weak transverse field along the $`x`$-direction, the two ground states move away from the $`z`$-axis and towards the $`x`$-axis while remaining in the $`xz`$-plane, as shown in Fig. 1. For $`0h<h_c^{(c)}=3N\lambda `$, the two minima are located at

$$\theta =\pi /2\pm \mathrm{cos}^1\left(h/h_c^{(c)}\right),$$

where $`\theta `$ is the angle between $`𝐒`$ and the $`z`$-axis. For $`hh_c^{(c)}`$, the two minima merge along the $`x`$-axis, and the degeneracy is removed. Hence $`h_c^{(c)}`$ can be regarded as the classical critical field strength, beyond which the system is completely polarized by the external field.

Let us now turn to a quantum mechanical description of the system. Our goal is to investigate whether or not tunneling is present in the classically degenerate regime (i.e., when $`0h<h_c^{(c)}`$). We will present both a full numerical calculation and analytical results using the instanton technique.

A well-known consequence of the tunneling between two degenerate states is the lifting of their degeneracy: The two new eigenstates are a symmetric and an antisymmetric superposition of the original states characterized by an energy difference (or tunneling splitting) $`\mathrm{\Delta }\epsilon `$ inversely proportional to the tunneling rate. The quantity of interest to determine the occurrence of tunneling is therefore this energy difference between the two lowest eigenstates of the Hamiltonian. We determine it by expanding the Hamiltonian (1) onto the basis spanned by $`S_1^zS_2^z`$, and evaluate numerically the eigenvalues of the resultant $`(2N+1)^2\times (2N+1)^2`$ matrix. Fig. 2(a) summarizes the result of this analysis. It shows that $`\mathrm{\Delta }\epsilon `$ is essentially zero for small values of $`h`$, but becomes finite when $`h`$ exceeds a threshold value $`h_c^{(q)}`$ — the quantum critical field strength. This means that for $`0hh_c^{(q)}`$, quantum mechanics agrees with classical mechanics in that the system is degenerate. However, for $`h_c^{(q)}<h<h_c^{(c)}`$, even though the system is still degenerate in the classical picture, the presence of tunneling removes the degeneracy in the quantum treatment. Fig. 2(b) displays $`h_c^{(q)}`$ as a function of $`N`$, from which one sees that $`h_c^{(q)}`$ increases with $`N`$ and approaches $`h_c^{(c)}`$ as $`N`$ tends to infinity. In other words, as $`N`$ increases, the system, as expected, behaves more and more classically.

To gain some analytical insight, we first notice that Hamiltonian (1) can be rewritten as

$$H=\frac{3\lambda }{4}[(S^z)^2(S^z)^2]+\frac{\lambda }{2}𝐒^2hS^x$$

where $`𝐒=𝐒_1+𝐒_2`$ and $`𝐒^{}=𝐒_1𝐒_2`$. For $`h<h_c^{(c)}`$, $`𝐒_1`$ and $`𝐒_2`$ are tightly bound together, such that $`S^z0`$ and $`𝐒`$ is approximately a constant of motion with quantum number $`S=2N`$. Hence, neglecting the constant terms, the effective Hamiltonian of the system reads

$$H_{\mathrm{eff}}=\frac{3\lambda }{4}(S^z)^2hS^x.$$
(2)
Hamiltonian (2) describes a quantum spin with the easy-axis anisotropy in a transverse field and it has been extensively studied in the context of spin tunneling book . Using the instanton technique, which has been proved to be quite accurate for $`S>5`$, the tunneling splitting between the two classically degenerated ground states can be expressed as garg

$$\mathrm{\Delta }\epsilon =p\omega \sqrt{S_c/(2\pi )}e^{S_c},$$
(3)
where $`p`$ is a prefactor on the order of unity \[by fitting the numerical results with Eq. (3), we find $`p2.75`$\], $`\omega =h_c^{(c)}x`$ is the typical instanton frequency with $`x=\sqrt{1(h/h_c^{(c)})^2}`$, and $`S_c=2N\mathrm{ln}[(1+x)(1x)]4Nx`$ is the classical action. The solid curves in Fig. 2(a) represent the values calculated using Eq. (3), which agree well with the numerical results.

A more intuitive understanding of this threshold behavior in quantum tunneling can be obtained from a quantum perturbation viewpoint. Tunneling results from the transverse magnetic field. If we regard the last term in Hamiltonian (1) as a perturbation to the rest of the Hamiltonian, which possesses two degenerated ground states $`|N,N`$ and $`|N,N`$, where we have labelled the states with the eigenvalues of $`S_1^z`$ and $`S_2^z`$, the tunneling level splitting is then, in the high order perturbation theory qm , given by the shortest chain of matrix elements and energy denominators connecting these two states garanin . For the situation at hand, this chain can be represented as

$`|N,N`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|N1,N+|N,N1)|N1,N1{\displaystyle \frac{1}{\sqrt{2}}}(|N2,N1+|N1,N2)`$
$`\mathrm{}{\displaystyle \frac{1}{\sqrt{2}}}(|N,N+1+|N+1,N)|N,N.`$
Since the spin has to travel a distance of $`4N`$ to reach from one ground state to the other, the tunneling appears minimally in the $`(4N)^{\mathrm{th}}`$ order of perturbation theory. The energy splitting associated with this chain can be calculated as

$`\mathrm{\Delta }\epsilon `$ $`=`$ $`16\lambda \left({\displaystyle \frac{h}{2\sqrt{2}\lambda }}\right)^{4N}\left[{\displaystyle \frac{\mathrm{\Gamma }\left(3/2+\sqrt{N^2+1/4}N\right)}{\mathrm{\Gamma }\left(N+\sqrt{N^2+1/4}1/2\right)}}\right]^2`$
$``$ $`8\lambda N^2\left({\displaystyle \frac{eh}{4\sqrt{2}N\lambda }}\right)^{4N},`$
where the last line is obtained under the limit $`N1`$. The power law dependence on the field strength indicates that $`\mathrm{\Delta }\epsilon `$ is vanishingly small for large $`N`$ unless the term inside the bracket exceeds one, resulting in a critical field strength at $`(4\sqrt{2}/e)N\lambda 0.7h_c^{(c)}`$. Although this treatment explains the threshold behavior, it fails to predict accurately the critical field strength. In particular, it fails to predict the $`N`$-dependence of the ratio $`h_c^{(q)}/h_c^{(c)}`$ \[see Fig. 2(b)\] . This failure can be easily understood since for $`h`$ close to the critical point, it is no longer valid to regard the transverse field as a perturbation.

We now turn to the detection of tunneling. It is not practical to directly measure the energy difference between the two lowest energy states, since $`\mathrm{\Delta }\epsilon `$ is relatively small compared to the total ground state energy. Instead, we use a method which is perfectly suited for the situation at hand, where thermal effects are negligible, but may not be appropriate for more traditional solid-state systems. Our proposed detection scheme starts with the system prepared in one of the degenerated ground state, say, $`|N,N`$. The external field $`h`$ is then slowly ramped up from zero to a final value, $`h_f<h_c^{(c)}`$. Classically, the macroscopic spin would simply adiabatically follow the instantaneous ground state closest to the initial state. Quantum mechanically, though, this is only true as long as $`h_f<h_c^{(q)}`$. But for $`h_f>h_c^{(q)}`$, the system is in the quantum tunneling regime and does not reach an equilibrium. Rather, it oscillates back and forth between the two classically degenerate states, a phenomenon termed as macroscopic quantum coherence. The observation of such oscillation will be a direct signature of quantum tunneling, but it requires the effect of dissipation to be small. Otherwise, spin relaxation will result, which is normally the case for the experiment on solid magnetic materials. The long coherence time associated with the atomic condensate, however, makes it ideal for this purpose. In fact, macroscopic center-of-mass oscillation has already been observed in ultracold atoms trapped in optical lattices jessen . Fig. 3 depicts the numerical simulation of the measurement scheme in both the classical and quantum treatments. The quantum mechanical results are obtained by time evolving the Schrödinger equation in spin space with Hamiltonian (1), using the Crank-Nicolson method, while the classical results are obtained by solving the dynamical equation

$$\frac{d𝐒_i}{dt}=\gamma _B𝐒_i\times \frac{\delta H}{\delta 𝐒_i},$$
(4)
where $`H`$ is the energy functional and has the same form as Hamiltonian (1) but with the spin operators treated as $`c`$-numbers. Fig. 3(a) is for $`h_f<h_c^{(q)}`$. There is no qualitative difference between the quantum and classical results, a steady state is reached in the end (small oscillations around this steady state persists due to nonadiabiticity). Fig. 3(b) is for $`h_f>h_c^{(q)}`$. The classical result is quite similar to that of Fig. 3(a). However, the quantum calculation clearly shows the oscillations of the system between the two macroscopically distinct states. These two states differ by a minus sign in the expectation value of $`S^z`$, which can be easily measured experimentally using, for example, Stern-Gerlach technique jessen .

In conclusion, we have shown that the quantum macroscopic tunneling of spin is possible in a spinor condensate trapped in double-well potential. Compared to the more conventional solid state magnetic materials, our system possesses several decisive advantages. First and perhaps foremost, this is an exceedingly clean system characterized by a few simple parameters, without the complication of domain separations, and with well-understood microscopic physics; It is amenable to exquisite experimental control; A single parameter — the transverse field strength — is capable of switching on and off the tunneling. Typical temperatures of the atomic condensate is in the nanoKelvin regime, hence the thermal activation normally presented in the solid materials can be safely neglected notespin . Due to these reasons, we believe that this system is indeed ideal for studying quantum magnetism in general pu ; spinwave , and macrosopic quantum tunneling in particular. Finally we remark that although we have considered a double-well potential in this work, we expect similar behavior for a condensate confined in a one-dimensional periodic lattice potential.

We would like to thank Dr. Hui Hu for many valuable discussions and correspondence, and Dr. Christian Tanguy for providing us the final expression of $`\mathrm{\Delta }\epsilon `$ in the perturbation treatment. This work is supported in part by the US Office of Naval Research under Contract No. 14-91-J1205, by the National Science Foundation under Grant No. PHY00-98129, by the US Army Research Office, by NASA, and by the Joint Services Optics Program.


text:
# Radiation Driven Implosion of Molecular Cloud Cores

## 1 Introduction

Massive stars have a strong impact on the evolution of their parental molecular clouds (MCs). Their fast winds and their ionizing radiation cause an important feedback of energy and momentum into the circumstellar environment.

One of the most interesting effects massive stars could have on MCs is “sequential star formation” \[Elmegreen & Lada 1977\]. When massive stars heat their environment to equilibrium temperatures $``$ 10 000 K via their ionizing UV radiation, the ionized region is divided from the surrounding molecular gas by a very thin transition layer called ’ionization front’. Due to the high pressure in the hot, ionized gas, a shock front is driven into the cooler molecular gas which sweeps up a thin, dense layer of shocked material between it and the ionization front. This layer is likely to be gravitationally unstable and thus may form a new generation of stars on timescales of several million years.

Molecular clouds are highly irregular, with observed structures from large scales down to even the smallest observable scales. There is ongoing debate whether the morphology of molecular clouds can be described by a fractal. Broad linewidths in molecular clouds also suggest irregular supersonic motion, which indicates the presence of turbulence. In the presence of ionizing radiation, theoretical models predict that density enhancements will likely be excavated by the ionization front \[Elmegreen, Kimura & Tosa 1995\]. They will subsequently protrude into the ionized gas, and are directly exposed to the ionizing radiation.

The dynamical behaviour (Radiation Driven Implosion, RDI) in such exposed globules was studied numerically by Sandford, Whitaker & Klein \[Sandford, Whitaker & Klein 1982, Sandford, Whitaker & Klein 1984\], who pointed out that the mass distribution of stars formed in a perturbed molecular cloud affected by RDI may sensitively depend on the spatial scales of the perturbations. A similar approach was followed by Lefloch & Lazareff \[Lefloch & Lazareff 1994\]. In their models, the time evolution can be subdivided into two stages: the implosion phase, lasting for about 100 kyrs, and, if gravitationally stable, a subsequent quiet, quasi-static cometary phase, in which the globule is slowly eroded by ionization on timescales of several Myrs. Through semi-analytical models, Bertoldi \[Bertoldi 1989\] and Bertoldi & McKee \[Bertoldi & McKee 1990\] constrained the ranges for mass and ionizing flux in which gravitational instability in the compressed globules should be expected. More recently, Hollenbach & Bally (1998) and López-Martin et al. (2001) investigated the photoevaporation of dense clumps of gas by an external source. The hydrodynamics of photoionized columns in the Eagle Nebula has been modeled by Williams et al. (2001).

Sugitani, Tamura & Ogura \[Sugitani et al. 1999, Sugitani et al. 2000\] observed a sample of 89 bright-rimmed clouds (BRCs), apparently objects with properties as mentioned above, and found an excess in the luminosities and luminosity-to-cloud ratios of embedded IRAS sources when compared to sources in isolated globules, which is an indication for enhanced star formation in their BRCs. Near-IR imaging of 44 BRCs revealed that young stellar objects seem to be aligned along the axis towards the ionizing cluster. There is even evidence for an age gradient with older stellar objects closer to the OB cluster and the youngest objects well inside the globules, aligned with the IRAS sources. This result is consistent with sequential star formation while the shock front advances further and further into the molecular cloud.

The interplay between dynamics and self-gravity in exposed globules has not yet been treated in detail. Using the new, fast SPHI method described in Kessel & Burkert \[Kessel-Deynet & Burkert 2000\], we present 3D calculations of imploding clumps including self-gravity. We find evidence that initial small-scale density perturbations can lead to a substantial stabilization against gravitational collapse during the implosion phase and the subsequent cometary phase, even if the parameters of the globules suggest them to be highly unstable according to the analytical estimations of Bertoldi & McKee \[Bertoldi & McKee 1990\].

## 2 Initial conditions

### 2.1 Density distribution

The initial conditions are static solutions to the equation of Emden \[Emden 1907\], called “Bonnor-Ebert spheres” \[Bonnor 1956\]. In these spheres, which are supported by an external pressure at the truncation radius $`r_\mathrm{t}`$, the gravitational force is everywhere balanced by the pressure gradient, which is vanishing at the center. The scale-free form of the differential equation describing this situation is

$$\frac{1}{\xi ^2}\frac{\mathrm{d}}{\mathrm{d}\xi }\left(\xi ^2\frac{\mathrm{d}\psi }{\mathrm{d}\xi }\right)=\mathrm{exp}(\psi ).$$
(1)
The relation between the parameters $`\xi `$ and $`\psi `$ on the one hand and physical quantities density $`\rho `$ and radius $`r`$ on the other hand is given by

$$\rho =\lambda \mathrm{exp}(\psi ),r=\beta ^{1/2}\lambda ^{1/2}\xi ,$$

where $`\lambda `$ is the maximum density for $`\psi =0`$ and

$$\beta =\frac{kT}{4\pi G\mu }$$

is the temperature of the isothermal gas. Here, $`k`$ is Boltzmann’s Constant, $`G`$ the Gravitational Constant and $`\mu `$ the mean molecular weight.

Stability analysis shows that there is a critical value $`\psi _\mathrm{c}`$ for the truncation radius of the sphere. Below $`\psi _\mathrm{c}`$, the sphere is stable against radial density perturbations. Above, it commences gravitational collapse on a few gravitational timescales if perturbed. At $`\psi _\mathrm{c}`$, the sphere contains about one Jeans mass of matter, and the density contrast from the center to the pressure boundary is $`14`$. Since we don’t expect such supercritical cores to survive long enough in molecular clouds to be influenced by ionizing radiation in relevant number, we instead restrict ourselves to the critical case with the truncation radius $`\psi _\mathrm{t}=\psi _\mathrm{c}`$. This enables us to study whether compression by ionizing radiation can force one initial Jeans mass into gravitational collapse. If not, we can conclude that induced collapse for cases with masses less than a Jeans mass should also be highly improbable. Additional simulations of initially Jeans-stable clumps with $`\psi _\mathrm{t}<\psi _\mathrm{c}`$ will be studied in a subsequent paper.

The maximum particle density of H atoms in the cores is assumed to be $`\lambda =1000\mathrm{cm}^3`$ and the equilibrium temperature in the neutral cloud material $`T_\mathrm{c}=10\mathrm{K}`$ to mimic the typical conditions observed in molecular clouds.

The integration volume has the shape of a cube, which contains the Bonnor-Ebert sphere, embedded in a surrounding gas of constant density $`n_\mathrm{s}`$. This gas provides the boundary pressure for the sphere. The diameter of the sphere is $`0.8`$ times the size of the box. Note that this configuration alone would be unstable against gravitational collapse due to the additional mass of the surrounding gas, in which pressure gradient and gravitational force are not in balance. However, the timescale of this collapse is much longer than the dynamical timescale which is imposed on the system by the large sound velocity of about $`13`$ km s<sup>-1</sup> in the hot, ionized medium. The setup chosen is not thought to provide a stable starting configuration. Our intention is rather to start with a reasonable density distribution, in which the initial core contains about one Jeans mass.

We allow material to stream freely away from the initial integration volume, i.e. no boundary conditions are imposed.

### 2.2 Initial density perturbations

To investigate the effects of density irregularities, we impose in some calculations Gaussian random perturbations on the initial density distribution as described by Klessen & Burkert \[Klessen & Burkert 2000\]. The density distribution is determined by the power spectrum $`P(k)k^\alpha `$. This is achieved by applying the Zel’Dovich \[Zel’Dovich 1970\] approximation. Besides the unperturbed case without and with self-gravity (in the following cases A and B), we ran perturbed cases with $`\alpha =3`$ (more power on large scales, case C) and $`\alpha =1`$ (more power on small scales, case D). The amplitude of the perturbations was chosen to be highly non-linear, with a mean relative variance of $`\sigma _\mathrm{r}=\sigma /\rho =0.5`$.

### 2.3 Parameters of the radiation field

The ionizing flux falls perpendicular onto one of the cube surfaces with a rate of photons per time and area, $`F_0`$. We specify the ”strength” $`C`$ of the radiation field as the ratio between the cube size and the path length $`l`$ which the radiation can penetrate into the surrounding gas via an r-type ionization front. This length is given by the balance between $`F_0`$ and the recombinations along $`l`$:

$$F_0=_0^ln(s)^2\alpha _\mathrm{B}𝑑s=n_\mathrm{s}^2\alpha _\mathrm{B}l,$$
(2)
where $`\alpha _\mathrm{B}`$ denotes the recombination coefficient adopting the ”on the spot” approximation \[Baker & Menzel 1962\]. In other words, the parameter $`C`$ determines how much of the density enhancement is excavated by the ionization front before it turns from weak r-type into weak d-type. From that moment on the medium dynamically reacts to the pressure difference across the ionization front. For the calculations in this paper, we chose $`C=.97`$.

## 3 Results

### 3.1 Time evolution without self-gravity: case A

The time evolution in case A (Figs. 1,2) is analogous to the models of Lefloch & Lazareff \[Lefloch & Lazareff 1994\]. One can distinguish between four different phases:

1. Immediately after the ionizing flux is switched on, a weak R-type ionization front is driven supersonically into the medium. After $`1.5`$ kyrs, the denser part of the cloud is excavated and surrounded by a region of ionized gas. The ionization front stalls and becomes nearly sonic. It converts into a weak D-type front. The gas can now react dynamically to the pressure imbalance between ionized and unionized gas (Fig. 1, first panels).
2. In the second phase, this imbalance drives a shock front into the cloud with velocities $`5\mathrm{k}\mathrm{m}/\mathrm{s}`$. The shock accumulates a thin layer of dense shocked material (Fig. 1, second panels). It thereby converges towards the axis of the finger-like structure (“globule”) due to the curvature of the ionization front.
3. After 130 Myrs, the globule reaches a stage of maximum compression (Fig. 1, third panels). Caused by the overshoot in pressure, the cloud starts reexpanding and
4. enters a fourth phase, in which gas is being ionized continuously away from the surface facing the ionizing source (Fig. 2). This leads to a mass loss of $`10^4\mathrm{M}_{}/\mathrm{yr}`$. Gas is leaving the ionization front with velocities of $`13\mathrm{km}/\mathrm{s}`$, causing a rocket effect which accelerates the gas in the tip into the direction of photon flow. During this acceleration, gas is collected inside the globule in a lid-shaped feature facing the ionizing source. It is subject to dynamical instabilities which lead to the disintegration of the cloud after $`600`$ kyrs. The remaining fragments are then ionized completely after a few 10 kyrs.

### 3.2 Time evolution with self-gravity, no initial density perturbations: case B

Placing our models with $`R_{\mathrm{pc}}^2/S_{49}4.71`$ and $`m_{}40`$ into the Figure 14 in Bertoldi & McKee \[Bertoldi & McKee 1990\] (reprinted here as Fig. 5), they lie far in the regime where gravitational instability is expected. Simulations with smooth initial conditions including self-gravity verify this conclusion. The overall evolution is indistinguishable to the one in case A, until the point of maximum compression is reached. At that moment gravitational collapse sets in, which can be seen in the sudden rise in the maximum density of the system (see Fig. 6).

### 3.3 Time evolution with self-gravity, large-scale and small-scale initial density perturbations: cases C and D

In these models, the evolution fundamentally differs from the one in model B. During the compression phase, density enhancements are amplified in the shocked layer, similar to what is observed in the simulations by Elmegreen et al. \[Elmegreen, Kimura & Tosa 1995\] (Figs. 3, 4) and those of Williams et al. (2001) (Fig. 6f). They are reminiscent of the observed small-scale structures seen in M16 (Hester et al. 1996; White et al. 1999). Contrary to case B, no collapse occurs at the moment of maximum compression. We rather observe an evolution analogous to the phases 3 and 4 of case A. After a delay of 300 kyrs, the layer starts fragmenting by gravitational instability.

## 4 Discussion

### 4.1 Jeans Criterion

As already pointed out in Chapter 2.1, the initial conditions were chosen such that the integration volume contains $`1`$ Jeans mass. Looking at the Jeans mass $`M_\mathrm{J}`$ and Jeans length $`L_\mathrm{J}`$,

$$M_\mathrm{J}=\frac{4}{3}\pi ^{5/2}\left(\frac{kT}{Gm}\right)^{3/2}\rho ^{1/2};L_\mathrm{J}=\sqrt{\frac{15kT}{2\pi \rho Gm}},$$
(3)
and the freefall timescale,

$$T_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32G\rho }},$$
(4)
it follows that $`M_\mathrm{J}`$ as well as $`T_{\mathrm{ff}}`$ drop with increasing density. At the moment of highest compression in model A, regions with densities higher than a threshold density $`\rho _\mathrm{m}=510^{17}\mathrm{g}\mathrm{cm}^3`$ contain $`5.3M_{}`$. For this $`\rho _\mathrm{m}`$, one finds $`M_\mathrm{J}=2.7\mathrm{M}_{}`$, $`L_\mathrm{J}=710^3\mathrm{pc}`$ and $`T_{\mathrm{ff}}=10\mathrm{kyrs}`$. Compared to these values, the compressed region in the tip of the cloud contains 2 Jeans masses and is large enough, and $`T_{\mathrm{ff}}`$ is small compared to other timescales in the system. One can thus expect that the globule is unstable to gravitational collapse at this moment. This explains why we observe collapse when including gravity in model B.

One should however use the Jeans masses criterion with care, as the cases C and D illustrate. It only considers the balance between internal energy and gravity and is thus a static approach, neglecting the effects of internal motion. It is as such only applicable in very few cases.

### 4.2 Reason for different behaviour of case B with respective to cases C and D

Since the initial conditions are chosen such that they contain $`1`$ Jeans mass for all cases, the stabilization of cases C and D cannot be caused by gas pressure alone. The only mechanism which differs in the simulations and which can be responsible is undirected gas motion, i.e. disordered kinetic energy.

Fig. 7 shows the total amount of kinetic energy $`E_{\mathrm{kin}}`$ relative to the center of mass compared to the gravitational energy $`E_{\mathrm{grav}}`$ in the system for case C. During the compression phase, more and more material is collected by the post-shock layer and accelerated to velocities $`5`$ km/s. As a consequence, $`E_{\mathrm{kin}}`$ grows steadily until the shock converges along the axis of the configuration. Then follows a short phase of strong dissipation caused by the colliding shocks. During reexpansion and acceleration of the globule, there is still a residual amount of $`E_{\mathrm{kin}}`$ left in the globule.

The difference between cases C and D with respect to model B is that in the latter case the gas flow is strongly converging in the tip of the globule, whereas in the former the shock front is refracted by the density fluctuations and thus a substantial amount of $`E_{\mathrm{kin}}`$ is transformed into undirected motion on scales small enough to support the globule. Additionally, the knots formed in the compression layer are relatively stable entities, which behave like a swarm of bullets. They tend to pass by each other without merging. Both effects inhibit gravitational collapse during maximum compression for cases C and D. In 2D calculations using cylindrical symmetry, these bullets would be represented by concentric rings and their cross section for collision would be artificially enhanced. It is thus doubtful whether the delay could be observed at all in 2D calculations.

In the late phases in models C and D, all models show gravitational instability in the dense layer which was built up in the tip of the cloud facing the source, which leads to fragmentation and disruption of the cloud (Fig. 4, lowest panel). This happens at different evolutionary times for several models we calculated, varying between 420 kyrs and 600 kyrs. We don’t show these late stages in Fig. 6, since then over 90 per cent of the gas is already ionized. This means that numerical resolution is poor due to low particle numbers in the neutral gas. A second reason why we doubt the reliability of the time when instability occurs is that it will strongly depend on the dissipation rate of the undirected kinetic energy. Since our numerical scheme uses artificial viscosity, we expect numerical dissipation which cannot represent the real value. However, the important result we trust in is the delay in gravitational collapse after maximum compression which we observe in all models with perturbed initial conditions.

The delay in the collapse for models C and D was neither observed nor anticipated in earlier work. Bertoldi & McKee \[Bertoldi & McKee 1990\] derived the conditions for which gravitational collapse can be expected by analysing the balance between the gas pressure in stratified clouds and the external pressure and recoil. In their Figure 14, reprinted here as Fig. 5, all the models presented here fall far into the unstable region. Note that our globules are far from stratified; the timescale for stratification is roughly given by the sound crossing time,

$$t_{\mathrm{sc}}=2r/a_\mathrm{c},$$
(5)
where $`r`$ is the radius of the globule and $`a_\mathrm{c}`$ the sound velocity in the neutral gas. In our case,

$$t_{\mathrm{sc}}\frac{0.5\mathrm{pc}}{0.3\mathrm{pc}/\mathrm{Myrs}}1\mathrm{M}\mathrm{y}\mathrm{r},$$
(6)
which is substantially longer than the timescales considered here.

### 4.3 Properties of the “bullets” in case C

In Fig. 8, the distribution of thermal and kinetic virial coefficients for the six most massive “bullets” in case D are shown after 100 kyrs evolutionary time, shortly before maximum compression. The coefficients were estimated by taking the ratio of internal or kinetic energy, respectively, compared to the gravitational energy at a radius, where the bullet is by its higher density still well separated from the embedding medium. The virial coefficients are all well above 10, which indicates that the bullets are no gravitationally bound entities, but are pressure-confined. This assumption is even affirmed due to the fact that all of the bullets are situated in the shocked, dense layer which is formed by the overpressure in the neighbouring ionized medium and the recoil caused by gas traversing the ionization front with high velocity.

The hydrodynamic lifetime of the bullets, estimated from their extent and the sound velocity in the neutral medium, are found to lie between 100 and 400 kyrs, which is comparable to the time the shock front needs to penetrate the whole globule.

### 4.4 Simplifications and caveats of the models

#### 4.4.1 Equation of state

In the calculations presented above, we use a very simplified equation of state. We assume that the heating and cooling timescales are short compared to the other important timescales, and that the equilibrium temperature is constant at $`T_\mathrm{c}=10\mathrm{K}`$ in the neutral and $`T_\mathrm{i}=\mathrm{10\hspace{0.17em}000}\mathrm{K}`$ in the ionized regions. For material in the ionization front, we couple the temperature to the ionization fraction by

$$T=T_\mathrm{c}(1x)+T_\mathrm{i}x.$$
(7)
This ansatz is unsatisfactory in several respects. First, the temperature of recombination zones is severely underestimated, since the temperatures of recently recombined gas tend to lie in the regime of several thousand degrees Kelvin. We expect that this leads to overestimated densities in the shadow region “behind” the main clump. We thus expect the mass reservoir to be overestimated as well, leading to longer lifetimes and increased probability for gravitational collapse in the late stages after rebounce.

It would be desirable to implement a detailed heating and cooling model, which takes into account the temperatures of dust and gas separately. This implies the application of a more general method for determining the radiation field generated by the massive stars. One possibility is application of flux-limited diffusion, which was proposed by Levermore & Pomraning \[Levermore & Pomraning 1981\] and already successfully implemented into grid codes \[Yorke & Kaisig 1995\]. This issue is currently being addressed.

#### 4.4.2 Diffuse recombination field

The model used for solving the radiation transfer of ionizing photons uses the simplification called “on the spot assumption”. Lyman continuum photons produced by electrons falling into the ground state of protons are themselves capable of ionizing H atoms. We assume that this happens very close to the site of their production and that following the exact equations of radiation transfer for these photons can be neglected. It was shown by Yorke et al. \[Yorke & Welz 1996\] that this assumption is not valid in special cases where the diffuse recombination radiation leads to the penetration of ionization fronts into regions which are shielded from the direct illumination by the ionizing stars. Those regions can thus be heated to temperatures around $`3000\mathrm{K}`$. As in the paragraph above, this again affects the regions where ionized gas streams into the shadows and starts recombining (see for example Fig. 3, lower panels).

We expect that the recombination of this material will be delayed, if not prevented, by the Lyman continuum photons. This can also be addressed as soon as an implementation of the Flux-Limited Diffusion Approximation for SPH is at hand.

## 5 Summary

In this paper we demonstrate that the physics of Radiation Driven Implosion strongly depends on the properties of the initial cloud material. In the presence of initial perturbations, the evolution is highly nonlinear and the fate of the imploding cloud is not easy to foresee. We find that gravitational collapse need not necessarily occur in those objects, depending on the presence, size and strength of initial density perturbations. The redistribution of kinetic energy from the implosion into undirected motion can stabilize the globules for a significant amount of time. Unfortunately, it is difficult to estimate this delay from the models presented here, since it depends on the decay rate of the kinetic energy, and the SPH method used here can reproduce this rate only insufficiently.

The delay could be relevant for star formation if the globule is stabilized for such a long time that the ionizing radiation can totally erode the globule. If gravitational collapse occurs before, we expect the collapsing regions to decouple from the gas around them. The remaining gas will then be further accelerated by the rocket effect from the ionization front, whereas the decoupled cores will lag behind and eventually be overtaken by the ionization front. This process could then repeat for several times, as long as there is still mass in the reservoir. We expect the outcome to be very similar to the observations of Sugitani et al. \[Sugitani et al. 1999, Sugitani et al. 2000\], where young stellar objects are aligned to the axis of the globules with an age gradient from the oldest objects far away from the globules and the youngest ones still sitting inside. By using the numerical technique called “sink particles” \[Bate 1995\], we are confident to be able to reproduce this scenario. One of the next steps should be to take into account the measured internal motion in molecular clouds. We will present calculations with clouds initially supported by turbulent motion in a subsequent paper.

## Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (DFG), grants Bu 842/4-1 and Bu 842/4-2. The 3D projections in Figs. 14 were computed with the help of starsplatter by Joel Welling and John Dubinski (http://www.psc.edu/Packages/StarSplatter\_Home/).