text:
# Shapes of tree representations of spin-glass landscapes

## 1 Introduction

The notion of energy (fitness) landscapes has played a crucial role in the development of many areas of physics and biology such as disordered systems, neural networks, combinatorial optimization problems , RNA folding , and evolutionary change , to mention only a few. In particular, considerable effort has been devoted to the study of the interplay between the geometry of the landscape and the nature of the relaxation dynamics, searching heuristic or evolutionary process unfolding on the landscape. However, the inherent high-dimensionality of these landscapes poses a serious hindrance to the characterization of their topology. In fact, most of the studies have focused on the statistical characterization of a few local properties of the landscape by looking at, e.g., the auto-correlation function of unbiased walks over the configuration space  or the energy distribution of local minima , while a satisfactory description of a landscape should address also the (relative) energy differences of the local minima, the height of the barriers between these minima, as well as the distribution of saddle-points .

The idea of condensing all the landscape information into a tree structure, termed barrier tree, was introduced in the context of RNA and protein folding , and spin-glass models . The advantage of barrier trees, whose leaves represent the local minima and the internal nodes the lowest-energy saddles connecting those minima, is that they are both visually appealing (much information can be obtained from just looking at them, see figure 1), as well as mathematically well-defined, lending themselves to rigorous analysis . However, a general quantitative measure to characterize unambiguously different kinds of barrier trees (and hence energy landscapes) remains to be obtained. For instance, the size-frequency distribution of low energy saddles $`\psi (w)w^D`$, where $`w=w(s)`$ is the fraction of minima that can be connected through saddle $`s`$, does not provide a good measure because $`D2`$ regardless of whether the barrier tree results from a spin-glass landscape or whether it is generated randomly .

The situation seems to be similar in the analysis of phylogenetic trees, where it is also believed that the shape of the tree contains valuable clues about the evolutionary process , but so far no *single* satisfactory measure of tree-shape has been constructed . In this contribution we apply five measures of tree shape that were originally used to study phylogenetic trees (see, e.g., ) to barrier trees resulting from the Ising $`p`$-spin model. These measures provide statistical information about the shape of the barrier tree, mainly its symmetry or balance, but ignore the lengths of the branches, i.e., the height of the barriers between minima. While the extreme statistics of these heights provide useful information on the performance of local search algorithms such as simulated annealing , a measure based solely on the shape of the barrier tree seems more adequate to classification purposes since the shape is probably insensitive to variations in minor details of the underlying energy landscape. Tree shapes are important in the biological context also because many methods of phylogeny estimation, including parsimony, do not produce branch lengths. We find that all five measures can be used to tell random from spin-glass trees, but only one measure can distinguish between different spin-glass trees.

A few cautionary remarks regarding the relevance of the concept of energy landscape to the understanding of the thermodynamic properties of spin glasses are in order. The characterization of an energy landscape is usually based on a few key elements – local minima, energy barriers and saddles – which, in turn, are rigorously defined in terms of the spin configuration space, the energy assignment to each spin configuration, and the neighborhood relation in the configuration space (see, e.g., Sec. 2). However, there seems to be no simple relation between those elements and the thermodynamic phases (equilibrium states) of the spin-glass model, since these phases are clearly independent of the neighborhood relation between spin configurations or the relaxation dynamics. To illustrate this point, we note that though the landscape of the Ising spin glass with short-ranged interactions in three dimensions is as complex as its counterpart of infinite-range interactions , its thermodynamics can be successfully studied under the assumption of a trivial ergodicity breaking at low temperature, in which there is a single thermodynamic phase. (We refer the reader to Ref.  for a lucid discussion of these controversial issues.) Nonetheless, the study of the organization of the metastable states (local minima) of spin glass models has a long tradition in the physics literature of disordered systems, beginning with the work of Bray and Moore more than two decades ago . It is from this perspective that the statistical studies of energy landscapes in general, and the present work in particular, should be considered.

In the next section, a brief overview of the Ising $`p`$-spin model and energy landscapes is given. Section 3 then reviews the notion of barrier trees. In section 4, phylogenetic trees and several tree shape statistics used to study them are discussed. In particular, we point out that the subtree that connects two leaves in a barrier tree corresponds to the evolutionary path of minimum fitness cost, in contrast with the traditional phylogenetic trees for which the subtrees are determined by the similarity between the leaves, regardless of the barrier height between them. The results of applying these measures to barrier trees of $`p`$-spin energy landscapes are presented in section 5. The main conclusions are then summarized in the final section of the paper.

## 2 Energy landscapes of $`p`$-spin models

Consider a system of $`N`$ Ising spins $`s=(s_1,\mathrm{},s_N)`$ where $`s\{1,+1\}`$, with the following energy function

$$H_p(s)=\underset{1i_1<i_2\mathrm{}<i_pN}{}J_{i_1i_2\mathrm{}i_p}s_{i_1}s_{i_2}\mathrm{}s_{i_p}.$$
(1)
Here, $`1pN`$ and the $`J_{i_1i_2\mathrm{}i_p}`$ are i.i.d. random variables from a Gaussian distribution with mean 0 and variance $`p!/(2N^{p1})`$. Thus, each of the $`2^N`$ possible spin configurations is assigned an energy value that is completely specified by the (magnetic) interactions within all possible subsets of $`p`$ spins. This spin glass is known as the $`p`$-spin model . In the regime of large $`p`$ and $`N`$ with $`Np`$ the energies of any two distinct spin configurations become independent random variables, distributed by a Gaussian of mean $`0`$ and variance $`N`$, so that the random energy model (REM) is recovered in this limit . Here the scaling of the variance with $`N`$ guarantees the extensivity of the free-energy. Another important limiting case is $`p=2`$, which corresponds to the Sherrington-Kirkpatrick (SK) model .

We call two spin configurations $`s`$ and $`t`$ *neighbors* if they differ in only one of the $`N`$ spins, i.e., when $`_{i=1}^N\left|s_it_i\right|=2`$. In other words, neighboring configurations can be turned into each other by a single spin flip. When $`p`$ is small, neighboring spin configurations will have highly correlated energy values, since the one spin $`s_i`$ that is different between the two configurations only influences a small number of the possible subsets of $`p`$ spins. When $`p`$ increases, the energy values of neighboring spin configurations will become less correlated, becoming completely uncorrelated in the limit $`p\mathrm{}`$.

The energy landscape of a $`p`$-spin model consists of the configuration space $`V`$ of the $`2^N`$ possible spin configurations, with the single spin flip neighborhood relation imposed on it, and where the energy value $`H_p(s)`$ of each spin configuration $`s`$ is considered to be its “height”. This gives rise to the intuitive image of a more or less mountainous landscape with peaks, valleys, and saddle points. A “walk” on this landscape consists of moving from one neighboring spin configuration to another, climbing up a peak or going down a valley, or perhaps just moving around randomly. A local search algorithm such as simulated annealing can be seen as performing such a walk, in search of the lowest valley. In addition, $`p`$-spin landscapes have been used to model evolutionary processes (see, e.g., ) since they form a class of tunably rugged landscapes similar to Kauffman’s Nk-model . In this context, evolution is described as an “adaptive” walk on the energy landscape.

## 3 Barrier trees

A local minimum in a $`p`$-spin energy landscape is simply a spin configuration $`s`$ that has a lower energy than all of its neighbors. A *path* $`\stackrel{}{p}_{st}`$ between two configurations $`s`$ and $`t`$ is a sequence of neighboring configurations, starting at $`s`$ and ending at $`t`$. In other words, it represents a series of single spin flips that transforms configuration $`s`$ into $`t`$. Note that there exist multiple paths between any pair of configurations $`s`$ and $`t`$. A *saddle point* between two local minima $`s`$ and $`t`$ is then defined as the minimum from the set of maximum energy values along each possible path $`\stackrel{}{p}_{st}`$ between $`s`$ and $`t`$; see e.g. . So, the energy value $`E[s,t]`$ of this saddle point is

$$E[s,t]=\underset{\stackrel{}{p}_{st}}{\mathrm{min}}\left\{\underset{z\stackrel{}{p}_{st}}{\mathrm{max}}\left\{H_p(z)\right\}\right\}.$$
(2)
The barrier $`B(s)`$ of a local minimum $`s`$ is defined as the height of the lowest saddle point that connects $`s`$ with a local minimum $`t`$ of lower energy,

$$B(s)=\underset{t}{\mathrm{min}}\left\{E[s,t]H_p(s)|H_p(t)<H_p(s)\right\}.$$
(3)
The information about the energy values of a landscape’s local minima and the barriers that connect these local minima can be represented by a barrier tree. In such a tree, the leaves of the tree represent the local minima, and the internal nodes represent the saddles, with the barrier sizes given by the length of the branches connecting the local minima to their corresponding saddles. Figure 1 shows an example of a barrier tree for an $`N=10`$ and $`p=2`$ $`p`$-spin landscape. There are 10 local minima in this landscape (labeled 1 to 10 in the tree), with 9 saddle points (the internal nodes). The length of each branch in the tree indicates the height of the corresponding barrier (this value is shown along the branch).

The algorithm for constructing these barrier trees is presented in . It is implemented in the barriers program111The source code is available at http://www.tbi.univie.ac.at/$``$ivo/RNA/Barriers/., which constructs the tree from a sorted list of energy values of all spin configurations in the landscape. The program barriers is used here to generate barrier trees of $`p`$-spin landscapes, see also  for applications to spin glass problems.

## 4 Phylogenetic trees

Phylogenetic trees are often used to study the historical relations within or between groups of biological species; see e.g. . The currently existing species (or subspecies) form the leaves of the tree, and two related species are linked through their last common ancestor, which forms an internal node in the tree. The length of the branch between a species and its ancestor indicates how long ago the speciation event occurred that led to the current species. “Dead-ends” in this tree represent extinction events. So, the shape of a phylogenetic tree contains information about patterns of speciation and extinction (and possibly about the rates of these occurrences), and thus tells us something about the evolution of different species.

The traditional approaches to reconstruct or infer a phylogenetic tree given the extant species (i.e., the leaves) are based on sequence similarity, i.e., configurational overlap only, which can be justified by the (usually implicit) assumptions of a flat fitness landscape and a diffusive behavior in sequence space . However, it seems intuitive that, regardless of the similarity between two sequences, if the fitness costs of all possible evolutionary paths connecting them are high, they must be put far apart in the phylogenetic tree. In the barrier tree approach both similarity and fitness cost are taken into account to yield the evolutionary path of minimum fitness cost connecting two species. Thus it may be viewed as a generalization of the maximum parsimony principle to rugged fitness landscapes. We recall that maximum parsimony chooses the tree (or trees) that require the fewest evolutionary changes (spin flips, in the present context). In figure 2 we give an example in which the distance-based tree and barrier tree have different topologies.

Various methods have been proposed to analyze phylogenetic trees. Here, we consider five statistics that were used in  to measure tree symmetry and balance. The trees are assumed to be binary trees with $`n`$ leaves (or species) and thus $`n1`$ internal nodes, with the root being the last common ancestor of all $`n`$ species. Let $`d(i,j)`$ be the graph-theoretical distance between two nodes of the tree, i.e., the number of edges along the path that connects them. Furthermore, we denote the root of the tree by $`\mathrm{}`$. The *height* of a leaf $`l`$ is $`h_l=d(\mathrm{},l)`$. Equivalently, $`h_l`$ is the number of internal nodes between leaf $`l`$ and the root $`\mathrm{}`$ (inclusive). For each interior node $`i`$ we have two subtrees with $`r_i`$ and $`s_i`$ leaves, respectively. We assume $`r_is_i`$. The *subtree-height* of an interior node $`i`$ is $`m_i=\mathrm{max}_{lT_i}d(i,l)`$ where the maximum is taken over all leaves $`l`$ in the subtree $`T_i`$ below $`i`$, i.e., the subtree of which $`i`$ is the root.

With this notation we may define the following five characteristic values for the shape of a binary rooted tree:

1. $`H=\frac{1}{n}_{l=1}^nh_l`$ is the *average height* of a leaf in the tree222In  this quantity is denoted by $`\overline{N}`$..
2. $`\sigma _H=\sqrt{\frac{1}{n}_{l=1}^n\left(h_lH\right)^2}`$, is the standard deviation of the leaf height.
3. $`C=\frac{2}{n(n3)+2}_{i=1}^{n1}\left(r_is_i\right)`$ is a measure for the *imbalance* of the tree. Up to normalization it is the same as Colless’s imbalance measure . A closely related measure of this type, which we will not use here, essentially amounts to averaging $`(r_is_i)/(r_i+s_i)`$ instead, see e.g. .
4. $`B_1=_i\mathrm{}1/m_i`$ is the average inverse subtree height, where the sum is taken over all $`n2`$ internal nodes $`i`$ excluding the root $`\mathrm{}`$.
5. $`B_2=_{l=1}^n2^{h_l}h_l`$ is an alternatively weighted average leaf-height.

In  the variance $`\sigma _H^2`$ was considered, but here we will use the standard deviation. Both $`H`$ and $`\sigma _H`$ have larger values for more asymmetric trees. In  it is shown that the expected value of $`H=2_{k=2}^n1/k`$ for random trees with $`n`$ leaves. Values of $`H`$ larger than this indicate trees more asymmetric than a random tree. For a completely symmetric tree, $`\sigma _H=0`$, while it has a maximum value for a completely asymmetric tree.

The imbalance measure $`C`$ examines the internal nodes of a tree. It “weighs” the subtrees branching out from each internal node by counting and comparing the number of leaves in each subtree. These weight differences are then averaged and normalized over all internal nodes of the tree. The value of $`C`$ increases from 0 for a completely symmetric tree to 1 for a completely asymmetric tree. The quantity $`B_1`$ looks at the longest possible path $`m_i`$ from each internal node $`i`$ to any of the leaves in its subtree. The statistic $`B_2`$ is based on an index of information content. For highly asymmetric trees $`B_2`$ will quickly converge to a value of 2. For a completely symmetric tree, it will be equal to $`\mathrm{log}_2(n)`$, where $`n`$ is the number of leaves in the tree. Both $`B_1`$ and $`B_2`$ have smaller values for increasingly asymmetric trees.

## 5 Results

Here, we apply the statistics presented in the previous section to barrier trees that result from $`p`$-spin energy landscapes. The parameter values used are $`N=10`$, $`12`$, $`15`$, $`18`$, $`20`$, and $`22`$, and $`p=2`$, $`3`$, and $`\mathrm{}`$ (REM). For each combination of $`N`$ and $`p`$, $`100`$ independent landscapes were generated randomly, and the barrier trees for each of these landscapes were constructed. The tree statistics reported here are the averages of the 100 trees for each parameter combination.

Three of the statistics ($`H`$, $`\sigma _H`$, and $`B_1`$) are exponential in $`N`$, the number of spins. For example, the relation between $`H`$ and $`N`$ is $`H=0.43\mathrm{exp}(0.20N)`$ for $`p=2`$. The data for $`H`$ (symbols) and the exponential fits (solid lines) are shown in a semi-log plot in figure 3a for the different values of $`p`$. The data for $`\sigma _H`$ and $`B_1`$ (not shown) are similar, but with different slopes. It turns out that the average number of leaves, $`n`$, in the barrier trees is also exponential in the number of spins $`N`$, with the same slope as for $`H`$. Analytical values for this exponent, defined as $`lim_N\mathrm{}\mathrm{ln}n/N`$, for the $`p`$-spin models are derived e.g. in  (see also ). The numerical values from the simulations reported here are compiled in table 1.

The imbalance $`C`$ is sub-linear in $`N`$, as shown in figure 3b, and eventually converges to $`1`$ for large $`N`$ (so the logarithm fitting cannot remain valid for all $`N`$). However, the value of $`B_2`$ very quickly converges to $`2`$ with increasing $`N`$, as shown in Figure 3c, and thus it does not seem to be a very useful measure to distinguish the trees. It should be noted at this point that the number $`n`$ of leaves, i.e., local minima, depends not only on $`N`$ but also on the type of the landscape. We will therefore consider tree measures as a function of $`n`$ rather than $`N`$ below. In fact, $`n`$ is the natural parameter in the analysis of phylogenetic trees, obtained by simply counting the number of leaves in the tree.

The data presented so far explicitly assumed that we already know the values of $`N`$ and $`p`$ of the energy landscapes. Suppose we do not know these values, but all we have are several instances (say, $`100`$) of barrier trees that were derived from some $`p`$-spin model. Is it possible to say something about the underlying landscape from just the barrier tree? (The problem faced by the phylogeneticists is even harder since they have access to a single tree of arbitrary size only .)

In figure 4a, the values of $`H`$ of all landscapes are lumped together and plotted against $`n`$, the (average) number of leaves in the tree, in a log-log plot. This results in a straight line, indicating a power-law behavior. In fact, the slope in this case is 1, since the slopes of $`H`$ and $`n`$ when plotted against $`N`$ are equal (see table 1). The plots for $`\sigma _H`$ and $`B_1`$ (not shown) are similar, but with slopes slightly larger and slightly smaller than 1, respectively.

In this case, it is not possible to distinguish the trees that result from, e.g., the $`p=2`$ and the $`p=3`$ landscapes: they all fall on the same line. Since $`B_2`$ converges to a fixed value for larger trees, it is not very useful for distinguishing different trees either. However, when the imbalance statistic $`C`$ is plotted against $`n`$, there is a clear distinction, as can be seen in figure 4b. In particular, if one plots $`C`$ against $`\mathrm{ln}\mathrm{ln}n`$ (i.e., if $`n`$ is on a double logarithmic scale), the data is fitted by straight lines with different slopes for the different values of $`p`$ (data not shown). So, $`C`$ clearly can be used as a statistic to distinguish and classify different barrier trees, and thus their underlying landscapes.

Next we address the question: How different are barrier trees of $`p`$-spin landscapes from random trees? To answer it, we generated random binary trees with $`n=10^i`$ leaves for $`i=1`$, $`2`$, $`3`$, $`4`$, and $`5`$. For each value of $`i`$, 100 random trees were generated and the same five statistics were calculated. The random trees were generated as follows. First, create $`n`$ nodes (the leaves) and put them in a set $`A`$. Next, remove two random nodes $`x`$ and $`y`$ from $`A`$, create a new node $`z`$ and make $`x`$ and $`y`$ its two children, and put $`z`$ in the set $`A`$. Repeat this procedure until there is only one node left in $`A`$, which will be the root of the tree.

Figure 5a shows the data for $`H`$ against $`n`$ (the number of leaves) in a semi-log plot for random trees. Clearly, $`H`$ depends logarithmically on $`n`$. The results for $`\sigma _H`$ and $`B_2`$ (not shown) are similar, but with different slopes. Figure 5b shows the data for $`C`$ against $`n`$ on a log-log plot. The fit to the data is a power law with exponent $`0.81`$. The results for $`B_1`$ (not shown) are similar but with a (positive) exponent very close to 1.

The results for random trees are quite different from those for the barriers trees of $`p`$-spin landscapes. In the $`p`$-spin case, $`H`$ and $`\sigma _H`$ depend exponentially on $`n`$, whereas for random trees this dependence is logarithmic. For $`C`$, the reverse is true. Moreover, for random trees, $`C`$ decreases with increasing $`n`$, while it increases for $`p`$-spin barrier trees. For $`B_1`$, in both cases the dependence is exponential, but with different exponents. Finally, for $`B_2`$ both cases are completely different, with a (downward) convergence to the value $`B_2=2`$ for $`p`$-spin trees, and a logarithmic increase for random trees. Clearly, barrier trees from $`p`$-spin landscapes are much more asymmetric than random trees. A similar, though qualitative, conclusion was reached by considering the size-distribution of minima connected through a high-energy (i.e., closer to the root) saddle-point .

Finally, we investigate the structure of subtrees of barrier trees. Instead of calculating the tree statistics on the entire tree, they are calculated on subtrees starting at some internal node of the tree. This way, it can be determined whether the tree has a self-similar structure or not (i.e., whether subtrees look similar to the tree as a whole). Three different instances of an $`N=20,p=3`$ $`p`$-spin landscape were taken and their barrier trees were constructed. Next, the five statistics were calculated on the subtrees starting from each of the $`n1`$ internal nodes of these trees (a binary tree with $`n`$ leaves has $`n1`$ internal nodes). Each internal node in a barrier tree represents a saddle point in the underlying energy landscape, and the internal nodes are characterized by the energy values of their corresponding saddle points. Obviously, nodes higher up in the tree (closer to the root) will have higher energy values than nodes lower down (or deeper) in the tree (closer to the leaves).

Figure 6a shows the results for $`H`$. The statistics for the subtrees are plotted against the energy value of the internal node that forms the root of the subtree. The data for the three different landscapes can be distinguished from the slightly different range in energy values, but the overall shape is the same for all three. Clearly, the values for $`H`$ depend strongly on the energy value of the root of the subtree. The results for $`\sigma _H`$ and $`B_1`$ (not shown) are similar. Figure 6b shows the results for the statistic $`C`$, where we have discarded subtrees with $`n=2`$ since $`C`$ is not defined for them. We note that $`C=1`$ for all subtrees with $`n=3`$ leaves. We find that subtrees with a high-energy root are extremely unbalanced. The same is true for the subtrees with very few leaves that we find near the global optimum. In the intermediate regime we find nodes with very balanced subtrees. For $`B_2`$, the values vary much more widely, but with most of the points falling on or near the $`B_2=2`$ line. The overall results for $`p=2`$ and for the REM are very similar.

¿From these plots it is clear that the structure of subtrees does not reflect the structure of the tree as a whole. Instead, there is a rather well-defined dependence on the structure of a subtree and its *depth* in the tree, i.e., the energy value of the internal node that forms its root. The landscape structure around local minima and saddle points with a relatively high energy value is therefore significantly different from the structure around local minima and saddle points with a low energy. Since this structure is correlated with the energy value of the local minima and saddle points, this information could possibly be used to guide local search algorithms.

## 6 Conclusions

One of the main outcomes of the mean-field (replica) theory of spin glasses was the prediction that, in the low-temperature phase, the phase space is broken into infinitely many pure states or valleys . Finding universal, in the sense of model independent, features of the distribution of valleys is one of the main goals of the theory . Unfortunately, the concept of valley in the replica theory is not easily related to the more tangible concepts of local minima and saddles commonly used to characterize complex landscapes. The reason being that the replica valleys are weighted by their Boltzmann factors which, at least for Gaussian distributed couplings, results in a single valley (up to an overall spin flip) at zero temperature, thus contrasting with the exponentially large number of local minima of the landscape. In that sense, the replica theory is of little use in the characterization of energy landscapes of spin-glass models. To investigate the organization of the local minima one has then to resort either to the annealed estimates of the correlations between local minima in the thermodynamic limit  or to the exact numerical calculation of the barrier trees for relatively small system sizes . As pointed out before, the latter seems a more convenient approach for the purpose of classifying families of landscapes or spin-glass models. Furthermore, the barrier tree also contains information about the structure of subspaces in the landscape around local minima and saddle points, since the structure and symmetry of the subtrees are clearly correlated with their depth in the tree (or the energy level of the node that forms the root of the subtree). This information can probably be used to guide local search algorithms such as simulated annealing in their search for the lowest (or global) minimum.

In searching for efficient measures to characterize landscapes, it was recently shown that the size-frequency distribution of the number of leaves $`w`$ connected by a saddle-point $`s`$ is too robust to be a useful measure, as it yields the same power-law $`\psi (w)w^2`$ for both $`p`$-spin barrier trees and random trees . In this contribution, we improve considerably the landscape systematics by considering five measures proposed originally to characterize the shape of phylogenetic trees . Three of the measures used, namely, $`H`$, $`\sigma _H`$ and $`B_1`$ were proved independent of the underlying spin-glass landscape when plotted against the number of leaves $`n`$, as evidenced by the “data collapse” illustrated in figure 4a. Only one of the measures, the imbalance $`C`$, can be used to differentiate between, e.g., $`p=2`$-type and $`p=3`$-type energy landscapes. It is remarkable, however, that all five measures yield completely different results for random trees, owing to a different scaling with $`n`$. In particular, barrier trees generated from energy landscapes of $`p`$-spin models appear to be much more asymmetric than random trees, and the asymmetry increases with increasing $`N`$ (the number of spins) or $`n`$ (the number of leaves). It remains a challenge to find (if it exists) a disordered spin system whose associated barrier trees exhibit balance and symmetry properties akin to those of the random trees or, perhaps an easier task, that violate the scaling law shown in figure 4a.

## Acknowledgments

Thanks to Christoph Flamm at the University of Vienna, Austria, for his help with the barriers program. This research was supported by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), project 99/09644-9. The work of J.F.F. is supported in part by CNPq and WH is supported by FAPESP.

## References


format: text

text:
# Crystal and magnetic structure of LaTiO3: evidence for non-degenerate 𝑡_{2𝑔}-orbitals

(, DRAFT)

## Abstract

The crystal and magnetic structure of LaTiO3 has been studied by x-ray and neutron diffraction techniques using nearly stoichiometric samples. We find a strong structural anomaly near the antiferromagnetic ordering, TN=146 K. In addition, the octahedra in LaTiO3 exhibit an intrinsic distortion which implies a splitting of the t2g-levels. Our results indicate that LaTiO3 should be considered as a Jahn-Teller system where the structural distortion and the resulting level splitting are enhanced by the magnetic ordering.

LaTiO3 has been studied already in the seventies and was thought to be a text book example of a Mott-insulator with antiferromagnetic order imada . Ti is in its three-valent state with a single electron in the t2g-orbitals of the 3d-shell. The titanate is hence an electron analog to the cuprates with a single hole in the 3d-shell. However, the t2g-orbitals in the LaTiO3 are less Jahn-Teller active and, therefore, the orbital moment may not be fully quenched in the titanate. The physics of the orbital degree of freedom has recently reattracted attention to this material keimer ; khaliullin .

The ordered moment in LaTiO3 amounts to 0.46 $`\mu _B`$ which is much smaller than the value of 1$`\mu _B`$ expected for a single electron with quenched orbital moment meijer . Quantum fluctuations can explain only about 15% reduction in the 3D-case. A straight-forward explanation could be given in terms of spin-orbit coupling, as an unquenched orbital moment would align antiparallel to the spin-moment in the titanate. However, in a recent neutron scattering experiment the magnon spin gap was observed at 3.3 meV, and it was argued that the strong interaction of an orbital moment with the crystal lattice implies a much larger value for the spin gap keimer . An orbital contribution to the ordered moment in LaTiO3 was hence excluded. On the basis of standard theories, however, even the G-type antiferromagnetic ordering in LaTiO3 may not be explained without a spin-orbit coupling. Instead one expects ferromagnetism 5 ; 6  related with the orbital degeneracy. Under the assumption of a specific structural distortion, Moshizuki and Imada recently presented a successful model for the antiferromagnetic order in LaTiO3 moshizuki . However, there is no experimental evidence for such a distortion. The puzzling magnetic properties of LaTiO3 led Khaliullin and Maekawa to suggest a novel theoretical description for RETiO3 based on the idea of an orbital liquid. They were able to explain many of the magnetic characteristics of LaTiO3 khaliullin , but the presumed orbital fluctuations have not been observed fritsch . Therefore, magnetism in LaTiO3 still remains an open issue.

We have reanalyzed the crystal and magnetic structure of LaTiO3 by x-ray and by neutron diffraction samples with almost perfect stoichiometry. First, we find a clear structural anomaly at the Néel-ordering and, second, the shape of the octahedra in this compound is not ideal but distorted. From these observations we conclude that LaTiO3 has to be considered as a soft Jahn-Teller system thereby explaining many of its unusual magnetic properties.

Samples of LaTiO3 were prepared following the standard solid state reaction; however the last process was performed in a floating zone image furnace. This procedure gave large single crystals with small mosaic spread. Since twinning would induce severe problems in any single crystal diffraction experiment, we crushed parts of these single crystals in order to perform powder diffraction studies. Special care was attracted to the stoichiometry of these samples. LaTiO3 has the tendency to incorporate La and Ti vacancies which are frequently denoted by an excess of oxygen in the formula LaTiO3+δ. By adjusting the ratio of TiO and TiO2 in the starting materials lichtenberg  one may vary the final oxygen content of the samples. Stoichiometry was verified by thermogravimetric analyzes and by determining the Néel-temperature in a SQUID magnetometer. For the stoichiometric sample studied here we obtain TN=146 K, amongst the highest values reported so far for LaTiO3. The temperature dependence of the lattice parameters was determined on a Siemens D5000 x-ray diffractometer using Cu-Kα radiation. Powder neutron diffraction experiments were performed at the Orphée-reactor in Saclay, using the high-resolution diffractometer 3T.2 ($`\lambda `$=1.2Å) and the high-flux diffractometer G4.1 ($`\lambda `$=2.4Å).

At room temperature, LaTiO3 exhibits the crystal structure of GdFeO3-type (space group $`Pbnm`$), see Fig. 1. This structure arises from the ideal perovskite (space group $`Pm\overline{3}m`$) by tilting the TiO6-octahedron around a c-axis (the subscript $`c`$ denotes the notation with respect to the cubic perovskite lattice) followed by a rotation around the c-axis (or c axis). One may note that while tilt distortions are opposite for any neighboring pair of octahedra, octahedra neighboring along the c-axis rotate in the same sense.

Fig. 2 shows the temperature dependence of the orthorhombic lattice parameters and that of the crystal volume. There are clear anisotropic anomalies near the antiferromagnetic ordering. Upon passing into the ordered state the lattice elongates along $`a`$ and shrinks along $`b`$, whereas for the c-direction and for the lattice volume no anomalous changes can be resolved. Anomalies in the thermal expansion have also been observed independently in a capacitance dilatometer experiment bruns .

The anomalous temperature dependence of the $`a`$ and $`b`$ lattice parameters is not restricted to a temperature range close to TN, see Fig. 3. In particular the $`a`$-parameter which is almost constant between 200-500 K deviates from a normal thermal expansion over a large temperature interval. Here, it is instructive to analyze the orthorhombic splitting. The lattice in $`Pbnm`$ results from that of an ideal cubic perovskite (denoted by subscript $`c`$) by setting $`a=\sqrt{2}a_c`$, $`b=\sqrt{2}a_c`$ and $`c=2a_c`$. Considering a rigid tilt of the ideal octahedron one obtains the non-equalities $`c\frac{1}{\sqrt{2}}(a+b)`$ due to the rotation and the tilt and $`ab`$ due to the tilt distortion only. However, tilting of an ideal octahedron around the $`b`$-axis would lead to $`b>a`$, whereas LaTiO3 exhibits just the opposite behavior below about 650 K. The sign of the orthorhombic distortion, $`ϵ=\frac{ba}{a+b}`$, already indicates that the octahedron shape in LaTiO3 is not ideal.

We have determined the crystal structure on our almost stoichiometric sample using high resolution neutron diffraction. The results are given in table 1. The splitting in the Ti-O-bond distances is small though not negligible. There is, however, a large splitting in the O2-O2 edge lengths of the octahedron basal plane, which has escaped previous attention. At 8 K the edges along $`a`$ and $`b`$ differ by about 4% . This distortion is comparable in its strength to that observed for the Ti-O-bond distances in YTiO3 which is considered as being orbitally ordered imada . In LaTiO3, the octahedron elongation perpendicular to the tilt axis over-compensates the shortening of $`a`$ through a rigid tilt. The negative sign of the orthorhombic strain in LaTiO3 is hence caused by the elongation of the octahedron along $`a`$. Interestingly, this effect is found to be strengthened at the antiferromagnetic transition; the behavior of the O2-O2 bond distances accounts for the anisotropic anomalies : an elongation along $`a`$ and a shrinking along $`b`$. Comparing our new results with those of ref. eitel ; maclean  one recognizes that the splitting in Ti-O-bond lengths, the splitting in the O2-O2 edges as well as the anomalies near $`T_N`$ are less pronounced in the samples of ref. eitel ; maclean . Non-stoichiometry seems to significantly reduce the long-range distortion of the TiO6-octahedron. This was confirmed by our own studies on samples with excess oxygen which revealed a continuous reduction of the orthorhombic strain and the anomalies around $`T_N`$.

Meijer et al. performed neutron diffraction studies on single crystals in order to determine the magnetic structure meijer . According to the symmetry of the GdFeO3 structure bertaut ; meijer  a G-type antiferromagnetic order may be associated with a weak ferromagnetic moment. However, Meijer et al. could not distinguish between the following two configurations due to twinning : either a G-type moment along $`a`$ and ferromagnetism along $`c`$ or a G-type moment along $`c`$ and ferromagnetism along $`a`$. Using the high flux diffractometer data we get evidence that the ordered antiferromagnetic moment points along the $`a`$-direction, see Fig. 4, although some underlying weak nuclear reflection complicates the analysis. Constraining the crystal structure to the high resolution results given in Table I, we find a low temperature ordered moment of 0.57(5) $`\mu _B`$ slightly higher than previous studies.

The crystal structure within the RETiO3 series is determined through the ionic radius of the RE. A smaller radius implies increasing tilt and rotation distortions. In Fig. 5a) we plot the orthorhombic splitting against the square of the tilt angle. There is a clear deviation from the proportionality expected in Landau-theory bruce , which we attribute to the anomalous stretching of the octahedron basal plane. Without this anomalous elongation one would expect an orthorhombic splitting in LaTiO3 three times larger and of positive sign. The anomalous octahedron elongation at the magnetic ordering is, therefore, about an order of magnitude smaller than the total deformation in LaTiO3, see Fig. 3 and table I. In Fig. 5b) we compare the octahedron deformation in LaTiO3 with that seen in the RETiO3-series using the data of MacLean et al. maclean . In LaTiO3 we find a small variation in the TiO-bond distances, described by $`r_{TiO}=\frac{TiO_{long}}{TiO_{short}}`$, but a strong variation in the basal plane edge lengths described by $`r_{O2O2}=\frac{O2O2_{long}}{O2O2_{short}}`$. In contrast, in YTiO3 the variation in the Ti-O-distances is strong and that in the O2-O2-lengths is small, as it was recently confirmed heester . Fig. 5b) shows that these distortions change continuously in the RETiO3 series. The two distinct deformations of the octahedra imply different orbital ordering schemes. The variation of the Ti-O-distances arises from the shift of an O2-site towards a Ti-site; in consequence, the distance to the opposite Ti is elongated and antiferroorbital-type ordering results which has been directly observed in YTiO3 akimitsu . In contrast the elongation of the basal planes along $`a`$ in LaTiO3 implies ferro-orbital-type ordering. Following Goodenough-Kanamori rules these orderings imply ferromagnetic and antiferromagnetic coupling, respectively. Most interestingly, the cross-over between the two types of structural deformation note-struc  coincides with the cross-over from antiferromagnetic to ferromagnetic ordering greedan ; amow1 ; amow2  in the RETiO3-compounds. With decreasing RE ionic radius NdTiO3 and SmTiO3 exhibit still antiferromagnetic order with small $`T_N`$ amow1 ; amow2 , whereas GdTiO3 is ferromagnetically ordered. A similar magnetic crossover is seen in the La1-xYxTiO3-series okimoto . This coincidence between the structural and magnetic crossovers clearly documents the relevance of the octahedron deformation for the magnetism in RETiO3.

In the following we analyze the effect of the octahedron basal plane elongation on the $`t_{2g}`$-energy levels. Qualitatively one may argue, that the octahedron elongation splits the degenerate $`xz`$ and $`yz`$-levels into $`\frac{1}{\sqrt{2}}(xz+yz)`$ and $`\frac{1}{\sqrt{2}}(xzyz)`$. In the notation used here the orbital coordinates refer to the cubic perovskite lattice, i.e. $`x,y,z`$ are parallel to the Ti-O-bonds. We have calculated the $`t_{2g}`$-level splitting within a full Madelung-sum point charge model taking into account the second order covalent contribution. Using the experimentally determined crystal structure, we obtain a sizeable splitting: one level is about 0.24eV below the two others which remain almost degenerate. The crystal field splitting is about one order of magnitude stronger than the spin orbit coupling, and we have to conclude that the $`t_{2g}`$-orbitals are not degenerate in LaTiO3 ; but the splitting is not sufficiently strong to completely quench the orbital moment. The occupied orbital can be constructed as a linear combination of the three standard $`t_{2g}`$-orbitals : $`0.77(\frac{1}{\sqrt{2}}|xz+yz>)\pm 0.636|xy>`$. The admixture of the $`xy`$-orbital to $`\frac{1}{\sqrt{2}}|xz+yz>`$ alternates in its sign along the $`c`$-axis and is constant in the $`a,b`$-plane. The occupied orbital is oriented in the orthorhombic $`b,c`$-plane and has almost the shape of the usual $`3z^2r^2`$ orbital pointing near the cubic c direction comment  (56o angle with the $`c`$-axis compared to 54.7o for the c-direction), see Fig. 5c). The lowered orbital points to the center of one of the octahedron triangles-faces which are formed by two O2’s and one O1. It is interesting to note that the O2-O1 distances of this triangle are elongated compared to the average. This means that the TiO6-octahedron is compressed along the direction of the occupied orbital. Since the occupied orbital is not oriented along one of the bonds the resulting antiferromagnetic exchange is not expected to be very anisotropic. Assuming the calculated orbital splitting and quantum fluctuations we obtain an ordered moment of 0.72 $`\mu _B`$. The ordered moment will be further reduced as covalence will diminish the level splitting, admix orbital contributions and transfer a part of the spin to the oxygen-sites which are not active in the G-type antiferromagnetic ordering. Similar to the proposal by Moshizuki and Imada moshizuki  we may hence explain most of the magnetic properties of LaTiO3 basing on the structural distortion observed in the almost stoichiometric compound. The distortion assumed in reference moshizuki  has not the same geometry as the one we find but is comparable in magnitude.

Since LaTiO3 exhibits intermediately strong t2g-level splitting, the orbital occupancy, the admixture of an orbital moment and the resulting antiferromagnetic exchange interaction depend sensitively on the crystal structure. On one side the system may gain magnetic energy in the ordered phase either by an enhancement of the magnetic interaction parameter or via spin-orbit coupling. On the other side there will be a loss of elastic energy associated with the additional structural distortion. In the ordered state there should be a trend to increase the orbital ordering in agreement with experiment; according to our calculation the level splitting is enhanced by 6% between 155 and 8 K.

In summary we have found a sizeable deformation of the TiO6-octahedron in LaTiO3 which causes a relevant splitting in the t2g-levels (about 0.24eV).

We gratefully acknowledge discussions with M. Grüninger, G. Khaliullin and M. Laad. This work was supported by the Deutsche Forschungsgemeinschaft through SFB 608.


text:
# Driven depinning of strongly disordered media and anisotropic mean-field limits

## Abstract

Extended systems driven through strong disorder are modeled generically using coarse-grained degrees of freedom that interact elastically in the directions parallel to the driving force and that slip along at least one of the directions transverse to the motion. A realization of such a model is a collection of elastic channels with transverse viscous couplings. In the infinite range limit this model has a tricritical point separating a region where the depinning is continuous, in the universality class of elastic depinning, from a region where depinning is hysteretic. Many of the collective transport models discussed in the literature are special cases of the generic model.

Nonequilibrium transitions from stuck to moving states underlie the physics of a wide range of phenomena, from fracture and earthquake rupture to flux flow in type-II superconductors Fisher98 . Various dynamical models have been proposed in different contexts. One class of models, overdamped elastic media pulled by an applied force $`F`$, has been studied extensively. These predict a nonequilibrium phase transition from a pinned state to a sliding state at a critical value $`F_T`$ of the driving force. Generically, this depinning transition displays critical behavior as in equilibrium continuous transitions footcontinuous , with the medium’s mean velocity $`v`$ acting as the order parameter. In overdamped elastic media, the sliding state is unique and no hysteresis can occur AAMunique . Universality classes have been identified, which are distinguished, for example, by the range of interactions or by the periodicity (or nonperiodicity) of the pinning forces.

The elastic medium model is often inadequate to describe many systems which exhibit plasticity (e.g., due to topological defects in the medium) or inertial effects (underdamping.) The dynamics in plastic systems can be inhomogeneous, with coexisting pinned and moving regions. The depinning transition may be discontinuous (first order), possibly with macroscopic hysteresis. Several mean-field models of driven extended systems with locally underdamped relaxation or phase slip have been proposed in the literature strogatz88&89 ; levy92 ; levy94 ; NV97 ; Fisher98 ; MMP00 ; SF01 .

In this paper, we present a model of driven disordered systems that encompasses several of the models discussed in the literature. This model incorporates the anisotropy of the sliding state in the plastic flow region that results either from flow along coupled channels oriented in the direction of drive (e.g., as in the moving smectic phase BMCMR97 ) or in layered materials such as the high-$`T_c`$ cuprate superconductors. We restrict ourselves to systems with a periodic structure along the direction of motion, such as charge density waves (CDWs), 2D colloidal arrays, and vortex lattices in type-II superconductors and consider only the dynamics of a scalar displacement field. This model uses coarse-grained degrees of freedom that are solid-like regions. These volumes slip relative to each other in the transverse dimensions, due to the presence of small scale defects (phase slips, dislocations, grain boundaries) at their boundaries, but remain elastically coupled in the longitudinal dimensions. Our current results are for the mean-field limits, in the transverse dimensions at least, with transverse viscous interactions. The studies carried out so far of these types of models for finite range interactions MCMKD02 ; SF02  suggest that the mean field approximation described here may give the correct topology for the phase diagram, although there will be corrections to the critical behavior in finite dimensions.

Anisotropic slip model. Elastic models have their limits in describing real materials. Coppersmith argued that elastic models with weak disorder have arbitrarily large strains, as large rare regions with atypically low pinning result in large displacement gradients coppersmith90&91 ; sncajm91 . This leads to the breakdown of the elastic model, though the low density of such regions may preclude its easy observation. Simulations, imaging and noise experiments for moderate and strong disorder also indicate the breakdown of the elastic model faleski96 ; nori96 ; myers99 ; olson\_hyst01 ; tonomura99 ; troyanovski99 ; Maeda02 . These calculations and observations often show a sliding state that consists of coherent regions moving at different average velocities, with extended defect structures along their boundaries. Furthermore, the dynamics is strongly anisotropic, with slip occurring preferentially where the shear deformations are largest, at the borders of channels that are on average aligned with the direction of motion. Guided by this work, we propose an anisotropic description of the inhomogeneous dynamics in the presence of defects.

We consider a $`d=d_{}+d_{}`$-dimensional medium composed of $`d_{}`$-dimensional elastic channels coupled via interactions that allow for slip of the channels in the remaining $`d_{}`$ directions. The system is driven by a uniform force $`F`$ along one of the $`d_{}`$ directions. A cartoon of a $`1+1`$-dimensional system is shown in Fig. 1. We choose to discretize space in both transverse and longitudinal directions, using integer vectors $`i`$ for the $`d_{}`$-dimensional intra-layer index and $`\mathrm{}`$ for the $`d_{}`$-dimensional layer index. The local displacement along the direction of motion is $`\varphi _{\mathrm{}}^i(t)`$. Assuming overdamped dynamics, the equation of motion in the laboratory frame convective  is

$`\dot{\varphi }_{\mathrm{}}^i(t)=K{\displaystyle \underset{j}{}}(\varphi _{\mathrm{}}^j\varphi _{\mathrm{}}^i)+F+h_{\mathrm{}}^iY(\varphi _{\mathrm{}}^i\beta _{\mathrm{}}^i)+\sigma _{\mathrm{},i}^\alpha ,`$ (1)
where the dot denotes a time derivative (we have chosen to scale time so that the damping constant is unity) and $`j`$ ranges over sites $`j`$ that are nearest neighbor to $`i`$. The first term on the right hand side represents an elastic intra-channel coupling of strength $`K`$. The third term is the pinning force. The function $`Y(x)`$ has period $`1`$, and the $`\beta _{\mathrm{}}^i`$ are random phases chosen independently and uniformly in $`[0,1)`$. The random pinning strengths $`h_{\mathrm{}}^i`$ are chosen from a probability distribution $`\rho (h)`$. Finally, $`\sigma _{\mathrm{},i}^\alpha `$ represents an interaction (of type $`\alpha `$) that allows for local slips of neighboring channels. In this paper, we assume a linear viscoelastic ($`V`$) stress-strain relation,

$$\sigma _{\mathrm{},i}^V=_{\mathrm{}}^t𝑑t^{}\underset{m}{}J_m\mathrm{}^\alpha (tt^{})[\varphi _m^i(t^{})\varphi _{\mathrm{}}^i(t)],$$
(2)
where $`m`$ indexes the layers neighboring $`\mathrm{}`$. The stress transfer function Fisher98 ; MMP00  $`J_m\mathrm{}^\alpha `$ is generally nonlocal in space and/or time. We focus on a simplified version of the viscoelastic model where the local slip force is a purely viscous coupling of strength $`\eta `$: $`\sigma _{\mathrm{},i}^\mathrm{V}=\eta _m(\dot{\varphi }_m^i\dot{\varphi }_{\mathrm{}}^i)`$. With this coupling, Eq. (1) is a simplified “viscous slip” form of the hydrodynamics of a driven viscoelastic medium that flows in response to large-scale shear, but responds elastically within the layers to long-wavelength compressions boonyip80 . Gradients of the displacement transverse to the channels or layers correspond to shear deformations arising from defect structures between the channels. Some justification for modeling the slip that results from moving defects as viscous couplings comes from work mcmks02  showing that that the hydrodynamics of a two-dimensional crystal with free dislocations is identical to that of a viscoelastic fluid. In contrast, gradients of displacement longitudinal to the drive mainly yield compressional deformations (exactly so when $`d_{}=1`$.) The elastic response to long wavelength compressions in fluids is intimately related to the fact that compressions are associated with fluctuations in the conserved density, suggesting that the inclusion of compressional forces of strength $`K`$ is necessary to describe the driven dynamics of systems with a conserved number of particles, such as vortex lattices. The coarse-grained model described by Eq. (1) allows us to investigate the competition between elastic interaction and plastic flow in controlling the topology of the nonequilibrium phase diagram.

Mean-field treatment. One mean-field approximation for the viscous slip model is obtained by taking the limit of infinite range interactions in both the transverse and longitudinal directions. Each discrete displacement then couples to all others only through the mean velocity, $`v=N^1\dot{\varphi }_{\mathrm{}}^i`$, and the mean displacement, $`\overline{\varphi }=N^1\varphi _{\mathrm{}}^i`$. We look for solutions moving with stationary velocity: $`\overline{\varphi }=vt`$. Since all displacements $`\varphi `$ are coupled, they can now be indexed by their disorder parameters $`\beta `$ and $`h`$, rather than the spatial indices $`i`$ and $`\mathrm{}`$. The mean-field dynamics is governed by the equation

$$(1+\eta )\dot{\varphi }(\beta ,h)=K\left(vt\varphi \right)+F+\eta v+hY(\varphi \beta ).$$
(3)
When $`K=0`$, the mean field velocity is determined by the self-consistency condition $`\dot{\varphi }(h)_h=v`$, where the subscript $`h`$ indicates an average over the distribution of pinning strengths $`\rho (h)`$. When $`K0`$, the mean field velocity is found by imposing $`\varphi (h,\beta )vt_{h,\beta }=0`$.

It is useful to review the case where $`\eta =0`$ and $`K0`$. In this limit, Eq. (3) reduces to the mean field theory of a driven elastic medium worked out by Fisher and collaborators NF92 . For the piecewise harmonic pinning $`Y(\varphi )=1/2\varphi `$ for $`0\varphi 1`$, no moving solution exists for $`F<F_T=\frac{h^2}{2(K+h)}_h`$. Just above threshold the mean velocity has a universal dependence on the driving force, with $`v(FF_T)^\beta `$. The critical exponent $`\beta `$ depends on the shape of the pinning force: $`\beta =1`$ for the discontinuous piecewise harmonic force and $`\beta =3/2`$ for generic smooth forces. Using a functional RG expansion in $`4ϵ`$ dimensions, Narayan and Fisher NF92  showed that the discontinuous force captures a crucial intrinsic discontinuity of the large scale, low-frequency dynamics, giving the general result $`\beta =1ϵ/6+𝒪(ϵ^2)`$, in reasonable agreement with numerical simulations in two and three dimensions myersSIM ; aamSIM . For simplicity and to reflect the “jerkiness” of the motion in finite-dimensional systems at low velocities, we use piecewise harmonic pinning.

When $`\eta >0`$, the nature of the depinning differs qualitatively from the $`\eta =0`$ case in that hysteresis in the dynamics can take place. This can be shown starting from the $`\eta =0`$ solutions $`v_{\eta =0}(F)`$, with $`v_{\eta =0}`$ a well-defined function of $`F`$. The solution for general $`\eta `$ can be found by substituting the effective drive force $`G=F+\eta v`$ for $`F`$ in the $`v_{\eta =0}(F)`$ relation and scaling velocity down by $`1+\eta `$. The linear transformation $`F=G\eta v`$ then gives the general $`v(F)`$ curve. One result is that when $`1+\eta ^1<\mathrm{max}_F(dv_{\eta =0}/dF)`$, the $`v`$ vs. $`F`$ curve is multivalued. Two cases can now be considered. The viscous-elastic case $`K0`$ has quite different behavior from the purely viscous case $`K=0`$, due to the distinct self-consistency conditions imposed by elasticity (or conservation in the flow direction.)

For the purely viscous case $`K=0`$ (the $`\tau =0`$ limit of the model in Ref. MMP00, ), periodic solutions of period $`T(h,G)`$ are found for a single particle moving in a pinning force of strength $`h`$ under the effective drive $`G`$. Self-consistency implies that $`v=[T(h,G)]^1_h`$, giving

$$v=\frac{1}{1+\eta }_0^{2G}𝑑h\rho (h)h\left[\mathrm{ln}\left(\frac{G+h/2}{Gh/2}\right)\right]^1.$$
(4)
The behavior now depends on the shape of the disorder distribution $`\rho (h)`$. For narrow distributions, that vanish below a value $`h_{\mathrm{min}}>0`$, no moving solutions exist for $`F<h_{\mathrm{min}}/2`$. In the extreme case of uniform pinning, $`\rho (h)=\delta (hh_0)`$, the depinning is hysteretic for all nonzero values of $`\eta `$, with a jump from $`v=0`$ to finite $`v`$ at $`F=h_0/2`$. For the more general broad distribution with nonvanishing support at $`h=0`$, no pinned solution exists for $`F>0`$: while the bulk of the degrees of freedom are pinned at small $`F`$, weakly pinned $`\varphi `$ can respond to small drives. Typical $`v`$ vs. $`F`$ curves for $`\rho (h)=e^h`$ are shown in Fig. 2(a). The lines in Fig. 2(b) represent the critical forces where there is a macroscopic jump in the average velocity, indicating dynamic hysteresis. For $`\eta >\eta _c`$ and $`F_{}>F>F_{}`$ “fast” and “slow” sliding states can coexist. The critical point $`(\eta _c,F_c)`$ in Fig. 2(b) is in the universality class of the liquid-gas transition (and of the field-driven random field Ising model MCMKD02 .) Near the critical point the mean velocity has universal scaling, $`vv_c(\eta \eta _c)^{1/\delta }(FF_c)^\beta `$, with $`v_c=v(F_c,\eta _c)`$, and $`\beta _{MF}=1/2`$, $`\delta _{MF}=3`$. The point $`(\eta _c,F_c)`$ survives in finite dimensions MCMKD02 , with exponents distinct from their mean field values.

For finite long-time elasticity, i.e., when $`K0`$, the behavior changes dramatically. The elastic forces or particle conservation enforce a uniform time-averaged velocity for all degrees of freedom. In this case the long-time uniform-$`v`$ solution to Eq. (3) is $`\varphi (h,\beta ,t)=vt+\stackrel{~}{\varphi }(h,\beta ,t)`$ with $`\stackrel{~}{\varphi }(h,\beta ,t)`$ of period $`1/v`$ in $`t`$. The self-consistency condition is $`\stackrel{~}{\varphi }(\beta ,t)_{\beta ,h}=0`$. The explicit solution can be obtained for the piecewise linear pinning force. (In the mean-field limit the nature of the phase diagram depends some on the shape of the pinning potential.) There is now a pinned phase that is stable stabilitynote  for $`F<F_T=\frac{h^2}{2(K+h)}_h`$. One then solves

$$\delta F=\frac{\eta _c\eta }{1+\eta _c}v+\frac{h^2}{(K+h)}\frac{1}{e^{(K+h)/[(1+\eta )v]}1}_h,$$
(5)
to find the mean velocity in the sliding state ($`0<\delta FFF_T`$), where $`(\eta _c+1)^1=\frac{h^2}{(K+h)^2}_h`$. The velocity-force curves and a phase diagram are shown in Fig. 2(c,d) for $`\rho (h)=e^h`$. There is now a tricritical point at $`(\eta _c,F_c=F_T)`$. For $`\eta <\eta _c`$, a continuous depinning transition at $`F_T`$ separates a pinned state from a sliding state with unique velocity $`v(1+\eta _c)\frac{FF_T}{\eta _c\eta }`$, giving $`\beta =1`$ in MFT. In finite dimensions, this transition is likely to remain in the same universality class as the depinning of an elastic medium ($`\eta =0`$): numerical studies and analysis by Schwarz and Fisher SF02  of a model with local viscous-like terms show that there is no hysteresis for small slip coupling and that the depinning transition exponents are the same as without slip coupling. In our mean-field example, the linear response diverges at $`\eta _c`$, $`v(\eta =\eta _c)1/\mathrm{ln}(FF_T)`$. For $`\eta >\eta _c`$ there is hysteresis with coexistence of stuck and sliding states.

Relations to other models. In the absence of viscous coupling ($`\eta =0`$), Eq. (1) reduces to the conventional $`d_{}`$-dimensional phase-only model of driven CDWs. When $`K=0`$ and $`\eta `$ is finite, the model describes shear deformations in a viscous fluid in a disordered background MMP00 .

The mean-field limit of the anisotropic slip model for finite $`K`$ maps onto a periodic-pinning version of the infinite-range version of the stress overshoot model studied by Schwarz and Fisher SF01 ; SF02 . The dynamics of crack fronts in brittle materials can be dominated by local inertial effects, in which the motion of a crack segment creates a transient stress on other segments, the stress overshoot, resulting in a net forward motion beyond the local minimum of the pinning. In the infinite-range limit the stress transfer yields a global coupling of strength $`M`$ between the mean displacements at two different times, analogue to the effective force per unit velocity $`\eta `$ in our model. The behavior we conjecture for finite dimensions is identical to that obtained by Schwarz and Fisher.

The viscous model ($`\eta 0`$) with finite $`K`$ is also closely related to a model of sliding CDWs that incorporates the coupling of the CDW to normal carriers by adding a global velocity coupling littlewood88 ; levy92 ; levy94  to the conventional Fukuyama-Lee-Rice model. The precise relationship is obtained by considering a mean-field limit $`d_{}\mathrm{}`$ with fixed $`d_{}`$. In this limit, each degree of freedom is coupled equally by slip to all other layers, but elastically coupled only to its nearest neighbors in the longitudinal direction. The equation of motion is then

$$(1+\eta )\dot{\varphi }_{\mathrm{}}^i=K\underset{j}{}(\varphi _{\mathrm{}}^j\varphi _{\mathrm{}}^i)+F+\eta v+h_{\mathrm{}}^iY(\varphi _{\mathrm{}}^i\beta _{\mathrm{}}^i).$$

Eq. (Driven depinning of strongly disordered media and anisotropic mean-field limits) is identical to the equation of motion of a purely elastic CDW with friction $`1+\eta `$ and an effective driving force $`F+\eta v`$. One can then obtain the velocity-force characteristics consistent with Ref. levy94,  simply by translating the analytic results NF92  for CDWs. Near threshold, we obtain $`v(F+\eta vF_T)^\beta `$, with $`\beta =1(4d_{})/6+𝒪[(4d_{})^2]`$. The transition is always hysteretic: the medium depins at $`F^{}=F_T`$ when the force is ramped up and repins at a lower value, $`F^{}F_T\left(\frac{1}{\beta }1\right)\left(\frac{\beta \eta }{1+\eta }\right)^{1/(1\beta )}`$, when the force is ramped down from the sliding state. This hysteresis for all $`\eta >0`$ appears to be a consequence of the global nature of the coupling, in contrast with the results of Schwarz and Fisher for local slip coupling SF02 .

One can include phase-slip couplings $`\sigma ^{\mathrm{PS}}`$ periodic in $`\varphi `$. The anisotropic model can describe the driven dynamics of intrinsically layered systems, such as vortex lattices in cuprate materials. Hysteresis depends on the strength of the correlations in the (weak-coupling) layered direction, as seen in simulations olson\_hyst01 . An anisotropic CDW model was studied in this context by Nattermann and Vinokur NV97 . These authors modeled the transverse slip as a nonlinear coupling that is periodic in the local CDW phase differences; this allows for phase slips, i.e., the collapse of the CDW amplitude between coherent regions strogatz88&89 .

Summary. We have proposed that the plastic flow of extended systems in a disordered background can be modeled using a general anisotropic approach that includes many specific physical systems as special cases. We have discussed the behavior of this model in a mean-field limit. When the degrees of freedom are allowed to move at distinct velocities, a pinned phase is generally not present, though there is hysteresis in the current-drive relation. Even when fluid-like shear takes place, particle conservation gives a sharp depinning transition in flow that takes place along channels. One of our future goals is to establish the connection between the parameters of the model and physical parameters of a given experimental system. This work has been supported in part by the National Science Foundation via grants DMR-9730678 and DMR-0109164.


text:
# New generations of stellar model atmospheres

## 1. Introduction

Determining stellar element abundances play an integral role in most endeavours to improve our understanding of stellar, galactic and cosmic evolution. The term observed abundances is somewhat of a misnomer, however, since the chemical composition of course can not be inferred directly from an observed spectrum. The obtained abundances are therefore never more trustworthy than the models employed to analyse the observations. The derivation of accurate abundances require realistic models of both the spectrum formation region and the spectrum formation process, an obvious fact which unfortunately is often overlooked in the analysis and in assessing the systematic errors. Costly new observational facilities like large telescopes and sophisticated instrumentation will only reach their full potential once combined with the appropriate analysis tools. Today, the uncertainties are no longer dominated by observational errors but by shortcomings in the modelling of stellar atmospheres and line formation. Nevertheless, much progress has been made lately in this respect and here I will highlight a few of these. Unfortunately due to page limitations I am unable to discuss many other impressive works for which I sincerely apologize. I will concentrate on one area each for hot and cool stars: line-blanketed non-LTE modelling and 3D hydrodynamical treatment of convection, respectively. I will conclude the review by making a call to arms to develop non-LTE inversion methods which are urgently needed to enable analysis of very large stellar samples.

## 2. Hot stars

### 2.1. Line-blanketed non-LTE model atmospheres

With the high temperatures encountered in hot star atmospheres, radiative rates mostly overwhelm collisional ones, rendering the assumption of local thermodynamic equilibrium (LTE) unacceptable. As a consequence the rate equations for all atomic level populations (assuming time-independence: statistical equilibrium) must be solved simultaneously with the radiative transfer equation, as well as the hydrodynamical equations of conservation of mass, momentum and energy (which in the absence of a wind simplify to hydrostatic and radiative equilibria), in a plane-parallel or spherical 1D geometry. Since the populations depend on the radiation field, which in turn is determined by the opacities and hence populations, this system of equations is extremely non-local and non-linear with everything in principle depending on everything else, everywhere else. Complicating factors also arise from the fact that the intense radiation field often drive massive stellar winds, which among other things cause severe problems for the treatment of radiative transfer by shifting spectral lines in and out of the shadows of other lines formed in deeper layers.

The last couple of decades has seen very impressive progress in constructing non-LTE model atmospheres (e.g. Hillier & Miller 1999; Pauldrach et al. 2001; Werner 2003 and references therein). Much of this is due to the development of efficient numerical algorithms such as the approximate lambda iteration technique (e.g. Scharmer 1981). The mean intensity $`J_\nu `$ is obtained from the source function $`S_\nu `$ using an iterative procedure $`J_\nu ^n=\mathrm{\Lambda }^{}S_\nu ^n+(\mathrm{\Lambda }\mathrm{\Lambda }^{})S_\nu ^{n1}`$ through an approximate (e.g. local) lambda operator $`\mathrm{\Lambda }^{}`$; $`\mathrm{\Lambda }S_\nu ^{n1}`$ simply denotes a formal radiative transfer solution with the current estimate of $`S_\nu `$ (see e.g. Werner 2003 for a detailed description). Improved convergence properties can also be achieved through preconditioning of the equations. Equally important for the success has been the advent of appropriate atomic data and methods for the statistical representation of them. Data from Opacity Project, Iron Project and the like has allowed line-blanketing from millions of spectral lines involving thousands of atomic levels from different species to be taken into account with $`10^5`$ frequency points using the superlevel concept (Anderson 1989).

Some work in this arena still remains however. Often the radiatively-driven wind is not self-consistently calculated from the theoretical radiative acceleration (the so-called wind-momentum problem of WR and similar stars) and instead a parametrized mass loss rate and velocity law are prescribed (e.g. Castor et al. 1975). Alternatively the calculations are restricted to either the quasi-hydrostatic photosphere or the supersonic wind parts. Observationally there is increasing evidence for both time-dependence and spatial inhomogeneities in the winds of hot stars in the form of shocks and clumps being accelerated outwards (e.g. Eversberg et al. 1998). Needless to say, time-dependent 3D non-LTE calculations is a very daunting task indeed but given the enormous progress achieved lately it is not inconceivable that also these issues will be addressed with the next generation of models in the coming years.

### 2.2. Implications for CNO abundances

The introduction of massive line-blanketing in the new generation of model atmospheres of hot stars directly affects the atmospheric structure, the ionization balance and the predicted emergent flux distribution. Indirectly it also modifies the estimated stellar parameters, in particular the effective temperature $`T_{\mathrm{eff}}`$ and the mass loss rate. For main sequence O stars the inclusion of line-blanketing can lead to a substantial decrease in $`T_{\mathrm{eff}}`$ of 1000-4000 K (Martins et al. 2002) at solar metallicity, which obviously also has an impact on the derived element abundances. Even larger differences are expected for giants and supergiants.

C, N and O are among the relatively few elements for which abundances can be determined in OB stars. With the improved atmospheric and ionization structures, the overall agreement between observations and predictions is in general quite impressive, with the exception of some lines formed in the extended wind (Pauldrach et al. 2001; Crowther et al. 2002). A common denominator in most analyses of Galactic and SMC/LMC O supergiants is the existence of non-solar CNO abundance ratios with a large N enrichment and weak C and O depletion. This indicates mixing of unprocessed and CNO-processed material, possibly related to rotational mixing (Meynet & Maeder 2000). Recently the first results for other external galaxies have appeared based on 8-10m telescopes, again suggesting He and N enrichments (Bresolin et al. 2002).

## 3. Cool stars

### 3.1. 3D time-dependent hydrodynamical model atmospheres

For late-type stars such as the Sun, the most pressing uncertainty stems from our poor understanding of convection. The visible surface from which the stellar spectrum originates represents the transition from convective to radiative energy transport. Traditionally convection has been treated by the mixing length theory (MLT) or some close relative thereof. In addition to being a very rudimentary description, MLT is a local theory, which ignores convective overshoot and radiative transfer effects. A far more realistic approach is to solve the hydrodynamical equations coupled to the equation of radiative transfer. Recently much progress has been made in performing 3D radiative-hydrodynamical simulations of stellar surface convection (e.g. Nordlund & Dravins 1990; Freytag et al. 1996, 2002; Stein & Nordlund 1998; Asplund et al. 1999, 2000a; Asplund & García Pérez 2001; Ludwig et al. 2002), which subsequently can be applied to studies of spectral line formation. These investigations have clearly shown that in many cases standard 1D analyses are very misleading in terms of derived element abundances.

The 3D model atmospheres developed by our group have been computed with a 3D, time-dependent, compressible, explicit, radiative-hydrodynamics code (Stein & Nordlund 1998). The hydrodynamical equations for conservation of mass, momentum and energy are solved on a Eulerian mesh with gridsizes of $`100^3`$. The physical dimensions of the grids are sufficiently large to cover many granules simultaneously and extended enough in the vertical direction to reach almost adiabatic conditions in the bottom. In terms of continuum optical depth the simulations extend at least up to log $`\tau _{\mathrm{Ross}}5`$, which for most purposes are sufficient to avoid numerical artifacts of the open upper boundary on spectral line formation. Periodic horizontal boundary conditions are employed, allowing the simulations to describe a small and representative region of the stellar surface (Fig. 1). The temporal evolution of the simulations cover several convective turn-over time-scales to allow thermal relaxation to be established and to obtain statistically significant average atmospheric structures and spectral line profiles following disk-integration. The input parameters discriminating different models are the surface gravity log $`g`$, metallicity \[Fe/H\] and the entropy of the inflowing material at the bottom boundary. The effective temperature of the simulation is therefore a property which depends on the entropy structure and evolves with time following changes in the granulation pattern.

In order to obtain a realistic atmospheric structure, it is crucial to have the best possible input physics, and properly account for the energy exchange between the radiation field and the gas. The adopted equation-of-state is that of Mihalas et al. (1988), which includes the effects of ionization, excitation and dissociation. The continuous opacities come from the Uppsala package (Gustafsson et al. 1975 and subsequent updates) while the line opacities are from Kurucz (1993). The 3D radiative transfer is solved at each time-step under the assumptions of local thermodynamic equilibrium (LTE, $`S_\lambda =B_\lambda `$) and opacity binning. The opacity binning includes the effects of line-blanketing in a manner reminiscent of opacity distribution functions.

It should be stressed that the here described 3D hydrodynamical model atmospheres do not include any free parameters which are tuned to improve the agreement with observations. Yet the simulations are highly successful in reproducing key observational constraints, such as helioseismology, detailed spectral line shapes, and granulation characteristics (e.g. Stein & Nordlund 1998; Asplund et al. 2000a,b). Thus the usual free parameters (mixing length parameters, micro- and macroturbulence) hampering classical 1D analyses have finally become obsolete with the advent of the new generation of 3D model atmospheres. The major drawbacks with such 3D modelling is their computationally demanding nature and that sofar only a relatively limited stellar parameter space has been chartered. Through the continuous improvement in computer hardware and software the former problem will be steadily relaxed in the coming years, while already work is being undertaken to extend the modelling to other types of stars, most notably A stars (Freytag et al. 1996), M dwarfs (Ludwig et al. 2002), supergiants (Freytag et al. 2002) and red giants.

### 3.2. Implications for CNO abundances

Spectrum formation is a highly non-local and non-linear process. It is therefore not surprising that 3D analyses in general yield element abundances which are not the same as with traditional 1D investigations, given both the different photospheric mean structures and the presence of spatial inhomogeneities. In many cases the differences can be profound with far-reaching consequences. Here I will limit the discussion on the impact of 3D model atmospheres on C, N and O for the Sun and metal-poor stars.

Sun: The solar spectrum offers a multitude of atomic and molecular lines to determine the abundances of C, N and O. Primary diagnostics are: \[C i\] 872.7 nm, high-excitation C i lines, and CH vibration-rotation transitions in the IR for C; permitted N i lines and NH vibration-rotation lines in the IR for N; \[O i\] 630.0 nm, high-excitation O i lines, and OH vibration-rotation and rotation-rotation transitions in the IR for O. Secondary, supportive roles are offered by electronic CH, C2, CN and OH lines. Astonishingly, with 1D model atmospheres the discrepancies between different indicators amount to as much as 0.2-0.3 dex, casting doubts on the use of such models for stars where not all types of lines can be employed. In particular for O the dissonance is striking with the often preferred OH lines (e.g. Grevesse & Sauval 1998) indicating about 0.3 dex higher abundance than the O i 777 nm triplet with the \[O i\] line in between.

We have recently re-analysed all of the above-mentioned solar diagnostics using our 3D hydrodynamical solar model atmosphere, the best possible atomic and molecular data, and high-quality solar atlases. Reassuringly, we find excellent agreement between different lines for C, N and O, in sharp contrast to the 1D case. In particular, the presence of cool pockets in the photosphere due to convection produces stronger molecular lines and hence lower derived abundances compared with standard 1D model atmospheres. Our preliminary results indicate log $`ϵ(\mathrm{C})=8.41\pm 0.03`$, log $`ϵ(\mathrm{N})=7.80\pm 0.04`$, and log $`ϵ(\mathrm{C})=8.66\pm 0.03`$, where the here quoted uncertainties are only the standard deviations between the different abundance indicators (Allende Prieto et al. 2001, 2002; Asplund et al., in preparation). It is noteworthy that this represents a substantial decrease from previously adopted values (e.g. Fig. 2) but is in good agreement with results for nearby B stars, local interstellar medium and solar corona/wind measurements. The resolution of these long-standing discrepancies gives further credibility to the new generation of 3D hydrodynamical model atmospheres.

Metal-poor stars: The differences between 1D and 3D model atmosphere predictions are perhaps most pronounced in metal-poor halo stars (Asplund et al. 1999). Due to the weak coupling through spectral lines between the radiation field and the gas at low metallicities, the temperatures in the optically thin photospheric layers can reach very low values with a proper treatment of convection, significantly below the radiative equilibrium temperatures enforced in standard 1D model atmospheres. For spectral lines sensitive to those atmospheric depths, the resulting differences can be profound. In particular, molecular features, minority ionization stages and low-excitation lines are very vulnerable and significant systematic errors can be expected in 1D analyses. This inadequate atmospheric modelling probably lies at the heart of the long-standing oxygen problem in metal-poor stars, i.e. that different diagnostics support either a flat plateau at \[O/Fe\] $`+0.4`$ (not accounting for the new low solar O abundance mentioned above) or linear rising trend reaching \[O/Fe\] $`+1`$ at \[Fe/H\] $`=3`$. In particular, the application of 3D model atmospheres bring down the very high oxygen abundances from OH lines (Asplund & García Pérez 2001) but also those of the \[O i\] line (Nissen et al. 2002); 3D non-LTE calculations of O i yield similar results as in 1D non-LTE (Asplund et al., in preparation). The final verdict is not yet in, partly due other out-standing factors such as non-LTE effects on OH, but it appears that the truth, not surprisingly perhaps, lies in between the hotly contested plateau and linear trend (Nissen et al. 2002). Similarly, many conclusions about C and N abundances in metal-poors probably have to be revised. Normally such determinations are based on molecular lines such as CH, NH and CN, and hence the abundances are likely significantly over-estimated in 1D. Naturally this would bring down the claimed \[C/Fe\] and \[N/Fe\] trends at low metallicity, as well as the incidence of very C-rich halo stars. Work is currently under way to address these issues.

## 4. Large-scale stellar abundance analyses

### 4.1. The need for semi-empirical model atmospheres

The new generations of stellar model atmospheres described above are both theoretical, as are almost all other models. However, the properties of the stellar atmospheres are encoded in the observed stellar spectra which at least in principle can be deciphered to yield the photospheric structure using inversion techniques. This method has often been employed among solar physicists with the Holweger-Müller (1974) and the VAL-3C (Vernazza et al. 1981) solar models being the most well-known examples. Surprisingly, such semi-empirical model atmospheres have rarely been used for stellar spectroscopy. I believe that the time is now ripe for pursuing such modelling, in fact crucial in order to make optimum use of future satellite and ground-based telescope facilities.

Perhaps the biggest problem with using the above-mentioned theoretical model atmospheres for stellar abundance analyses is their computationally intensive nature. Although less sophisticated model atmospheres are less CPU-demanding, also standard abundance analyses are time-consuming in terms of the work involved. Even today the largest studies of stellar abundances involve only a few hundred stars (e.g. Reddy et al. 2003). To fully disentangle the formation and evolution of for example our own Galaxy would require a far more ambitious approach with much larger stellar samples (Freeman & Bland-Hawthorn 2002). Fortunately this is now becoming feasible with new multi-object high-resolution spectrographs like FLAMES/VLT. Several other, even more ambitious instruments are currently in the design phase or being discussed in detail, such as the RAVE project (Steinmetz et al. 2002) and the GAIA satellite (Perryman et al. 2001). The wealth of information these facilities will provide is truly astonishing, as is the magnitude of data needed to be analysed. Without an efficient automatic procedure to derive the inherent information in the spectra, these projects will not come to an optimal execution. Semi-empirical model atmospheres represent one such promising possibility, provided the spectral resolving power is sufficient.

### 4.2. Ingredients and outlook for the future

Semi-empirical model atmospheres are obtained through an iterative procedure aimed at improving the overall agreement between the calculated and observed spectra. Since the temperature structure is modified in this way the model does not in general fulfill flux constancy as in standard theoretical model atmospheres, nor is there any need to invoke uncertain and questionable recipes for convective energy transport, such as the mixing length theory, which can lead to grossly misleading results (e.g. Asplund et al. 1999; Asplund & García Pérez 2001). In fact, this allows at least in principle a possibility to gain insight to convection and non-standard energy transport, such as acoustic and magnetic heating. With the temperature determined from the observations, the gas pressure is normally fixed by the assumption of hydrostatic equilibrium while the electron pressure is the result of the ionization balance using the current estimate of the chemical composition. Thus knowledge of the surface gravity is necessary from for example parallaxes or photometry unless it can be obtained from satisfying the ionization balance of some key elements. Essentially all existing modelling is based on a 1D plane-parallel geometry but some attempts with multi-component atmospheres have been made, corresponding for example to convective up- and downflows (Frutiger et al. 2000). Naturally, inversions must rely on accurate atomic data to yield realistic results.

The main disadvantages with semi-empirical modelling comes with the numerical difficulties often associated with ill-posed inversion problems and the assumption of LTE normally made in the line formation calculations. The former can to a large extent be overcome by including a wide range of spectral lines with different line formation depths, ionization stages and excitation potentials. Of particular importance, at least for late-type stars, is to include elements with both neutral and ionized lines and strong lines with pressure-damped wings. Although not yet utilized, hydrogen lines as well as various molecular lines have a great potential in this respect. The second major shortcoming can only be resolved by relaxing the assumption of LTE, which automatically makes the inversion procedure much more challenging and cumbersome. Very little work with the exception of the pioneering contribution by Soccas-Navarro et al. (1998) has been done on non-LTE inversion. It is, however, of absolute pivotal importance to include such effects, given the fact that LTE is often a poor assumption for the line formation process, leading otherwise to erroneous atmospheric structures and inferred conclusions. It suffices here to mention the case of Fe i which offers the most useful lines for inversion modelling in late-type stars but suffers from over-ionization, in particular at low metallicities.

Recent endeavours like the Opacity Project and Iron Project have come a long way in addressing the great needs for atomic data for non-LTE calculations. Although it would be advantageous to treat all elements in non-LTE for the inversions, as a first step it may be sufficient to do so for the most important elements for determining the actual atmospheric structure, such as Fe, Mg and Ca, while the line formation of other significant species for late-type stars like H and various molecules (e.g. CO) can likely be well approximated by LTE. Other trace elements which are unimportant for setting the atmospheric structure but whose abundance may be of interest can then be handled either in LTE or non-LTE after the initial inversion has converged. This procedure should restrict the computing time to just a few times a normal non-LTE calculation with a fixed 1D atmosphere. Adopting the superlevel concept to reduce the number of atomic levels for complex atoms like Fe should speed up the computations further and bring the non-LTE inversion concept fully within the realms of possibility for the coming few years. This should enable truly grandiose studies of galaxy formation and evolution possible with current and future observational capabilities.

#### Acknowledgments.

I am grateful to the SOC and LOC for organizing an extremely stimulating conference in honour of Prof. André Maeder, for providing me with a travel grant to attend the meeting and for showing great patience in waiting for my contribution to the proceedings. I am already looking forward to the conference celebrating André’s 70th birthday, which will hopefully take place in a similarly wonderful environment and atmosphere. Even further into the future I hope I am not too old to also participate in the conference held in connection with André’s retirement.

## References

Allende Prieto, C., Lambert, D.L., Asplund, M. 2002, ApJ, 556, L63

Allende Prieto, C., Lambert, D.L., Asplund, M. 2001, ApJ, 573, L137

Anderson L.S., 1989, ApJ, 339, 558

Asplund, M., García Pérez, A.E. 2001, A&A, 372, 601

Asplund, M., Nordlund, Å., Trampedach, R., Stein, R.F. 1999, A&A, 346, L17

Asplund, M., Ludwig, H.-G., Nordlund, Å., Stein, R.F. 2000a, A&A, 359, 669

Asplund, M., Nordlund, Å., Trampedach, R., et al., 2000b, A&A, 359, 729

Bresolin F., Kudritzki R.-P., Najarro F., et al., 2002, ApJL, 577, L107

Castor J.I., Abbott D.C., Klein R., 1975, ApJ, 195, 157

Crowther P.A., Hillier D.J., Evans C.J., et al., 2002, ApJ, 579, 774

Eversberg T., Leṕine S., Moffat A.F.J., 1998, ApJ, 494, 799

Freeman K.C., Bland-Hawthorn J., 2002, ARAA, 40, 487

Freytag B., Ludwig H.-G., Steffen M., 1996, A&A, 313, 497

Freytag B., Steffen M., Dorch B., 2002, AN, 323, 213

Frutiger C., Solanki S.K., Fligge M., Bruls J.H.M.J., 2000, A&A, 358, 1109

Grevesse N., Sauval A.J., 1998, Space Science Reviews, 85, 161

Gustafsson, B., Bell, R.A., Eriksson, K., Nordlund, Å. 1975, ApJ, 42, 407

Hillier D.J., Miller D.L., 1999, ApJ, 519, 354

Holweger, H., Müller, E.A. 1974, Solar Physics, 39, 19

Kurucz, R.L. 1993, CD-ROM, private communication

Ludwig H.-G., Allard F., Hauschildt P.H., 2002, A&A, 395, 99

Meynet G., Maeder A., 2000, A&A, 361, 101

Mihalas, D., Däppen, W., Hummer, D.G. 1988, ApJ, 331, 815

Nissen, P.E., Primas, F., Asplund, M., Lambert, D.L. 2002, A&A, 390, 235

Nordlund, Å., Dravins, D. 1990, A&A, 228, 155

Pauldrach A.W.A., Hoffmann T.L., Lennon M., 2001, A&A, 375, 161

Perryman M.A.C., de Boer K.S., Gilmore G., et al., 2001, A&A, 369, 339

Reddy B.E., Tomkin J., Lambert D.L., et al., 2003, in press \[astro-ph/0211551\]

Scharmer G.B., 1981, ApJ, 249, 720

Soccas-Navarro H., Ruiz Cobo B., Trujillo Bueno J., 1998, ApJ, 507, 470

Stein, R.F., Nordlund, Å. 1998, ApJ, 499, 914

Steinmetz M., 2002, in press \[astro-ph/0211417\]

Vernazza J.E., Avrett E.H., Loeser R., 1981, ApJS, 45, 635

Werner K., 2003, in press \[astro-ph/0208294\]