text:
# References

On P-wave meson decay constants
in the heavy quark limit of QCD

A. Le Yaouanc, L. Oliver, O. Pène and J.-C. Raynal
Laboratoire de Physique Théorique***Unité Mixte de Recherche UMR 8627 - CNRS
Université de Paris XI, Bâtiment 210, 91405 Orsay Cedex, France

V. Morénas

Laboratoire de Physique Corpusculaire
Université Blaise Pascal - CNRS/IN2P3, F-63000 Aubière Cedex, France

## Abstract

In previous work it has been shown that, either from a sum rule for the subleading Isgur-Wise function $`\xi _3(1)`$ or from a combination of Uraltsev and Bjorken SR, one infers for $`P`$-wave states $`|\tau _{1/2}(1)||\tau _{3/2}(1)|`$. This implies, in the heavy quark limit of QCD, a hierarchy for the production rates of $`P`$-states $`\mathrm{\Gamma }(\overline{B}_dD\left(\frac{1}{2}\right)\mathrm{}\nu )\mathrm{\Gamma }(\overline{B}_dD\left(\frac{3}{2}\right)\mathrm{}\nu )`$ that seems at present to be contradicted by experiment. It was also shown that the decay constants of $`j=\frac{3}{2}`$ $`P`$-states vanish in the heavy quark limit of QCD, $`f_{3/2}^{(n)}=0`$. Assuming the model of factorization in the decays $`\overline{B}_d\overline{D}_s^{}D`$, one expects the opposite hierarchy for the emission rates $`\mathrm{\Gamma }(\overline{B}_d\overline{D}_s\left(\frac{3}{2}\right)D)\mathrm{\Gamma }(\overline{B}_d\overline{D}_s\left(\frac{1}{2}\right)D)`$, since $`j=\frac{1}{2}`$ $`P`$-states are coupled to vacuum. Moreover, using Bjorken SR and previously discovered SR involving heavy-light meson decay constants and IW functions, one can prove that the sums $`\underset{n}{}\left(\frac{f^{(n)}}{f^{(0)}}\right)^2`$, $`\underset{n}{}\left(\frac{f_{1/2}^{(n)}}{f^{(0)}}\right)^2`$ (where $`f^{(n)}`$ and $`f_{1/2}^{(n)}`$ are the decay constants of $`S`$-states and $`j=\frac{1}{2}`$ $`P`$-states) are divergent. This situation seems to be realized in the relativistic quark models à la Bakamjian and Thomas, that satisfy HQET and predict decays constants $`f^{(n)}`$ and $`f_{1/2}^{(n)}`$ that do not decrease with the radial quantum number $`n`$.

LPT Orsay 01-36

June 2001

In previous work ,  we have pointed out a problem between experiment and the heavy quark limit of QCD for the semi-leptonic decays of $`B`$ mesons to $`L=1`$ excited states $`D_{0,1}\left(\frac{1}{2}\right)`$, $`D_{1,2}\left(\frac{3}{2}\right)`$. In a few words, the argument is as follows.

From Bjorken SR

$$\rho ^2=\frac{1}{4}+\underset{n}{}|\tau _{1/2}^{(n)}(1)|^2+2\underset{n}{}|\tau _{3/2}^{(n)}(1)|^2$$
(1)
and the recently demonstrated Uraltsev SR

$$\underset{n}{}|\tau _{3/2}^{(n)}(1)|^2\underset{n}{}|\tau _{1/2}^{(n)}(1)|^2=\frac{1}{4}.$$
(2)
one obtains,

$`{\displaystyle \underset{n}{}}|\tau _{3/2}^{(n)}(1)|^2={\displaystyle \frac{\rho ^2}{3}}`$ (3)
$`{\displaystyle \underset{n}{}}|\tau _{1/2}^{(n)}(1)|^2={\displaystyle \frac{1}{3}}\left(\rho ^2{\displaystyle \frac{3}{4}}\right).`$ (4)
One can see that $`\underset{n}{}|\tau _{3/2}^{(n)}(1)|^2`$ is proportional to $`\rho ^2`$ and that $`\underset{n}{}|\tau _{1/2}^{(n)}(1)|^2`$ is proportional to the deviation of $`\rho ^2`$ from the lower bound $`\frac{3}{4}`$. Then, there is little room left for $`\underset{n}{}|\tau _{1/2}^{(n)}(1)|^2`$, as it has been pointed out recently from a SR obtained for the subleading function $`\xi _3(1)`$ . The expected hierarchy for the form factors $`|\tau _{3/2}(1)|\tau _{1/2}(1)`$, that can be inferred from the precedent equations, implies that $`\overline{B}_dD^{}\mathrm{}\nu `$ and $`\overline{B}_dD^{}\pi `$ (assuming factorization, a reasonable assumption in view of the recent papers on non-leptonic decays with emission of a light meson) are dominated by the narrow resonances

$`\mathrm{\Gamma }\left(\overline{B}_dD_{1,2}\left({\scriptscriptstyle \frac{3}{2}}\right)\mathrm{}\nu \right)\mathrm{\Gamma }\left(\overline{B}_dD_{0,1}\left({\scriptscriptstyle \frac{1}{2}}\right)\mathrm{}\nu \right)`$ (5)
$`\mathrm{\Gamma }\left(\overline{B}_dD_{1,2}\left({\scriptscriptstyle \frac{3}{2}}\right)\pi \right)\mathrm{\Gamma }\left(\overline{B}_dD_{0,1}\left({\scriptscriptstyle \frac{1}{2}}\right)\pi \right).`$ (6)
On the other hand, it was demonstrated that the decay constants of $`j=\frac{3}{2}`$ $`P`$-wave mesons vanish, $`f_{3/2}^{(n)}=0`$, in the heavy quark limit of QCD , . This result can be summarized in the statement that we expect the “broad” $`D^{}`$ resonances ($`j=\frac{1}{2}`$) to have a much larger decay constant than the “narrow” ones ($`j=\frac{3}{2}`$). This to be contrasted to the opposite expected hierarchy for form factors, $`|\tau _{3/2}||\tau _{1/2}|`$ , that can be inferred from equations (3), (4). This hierarchy in the production is expected to be opposite to the one due to the selection rule $`f_{3/2}^{(n)}=0`$,

$$\mathrm{\Gamma }\left(\overline{B}_d\overline{D}_{s1,2}\left({\scriptscriptstyle \frac{3}{2}}\right)D\right)\mathrm{\Gamma }\left(\overline{B}_d\overline{D}_{s0,1}\left({\scriptscriptstyle \frac{1}{2}}\right)D\right)$$
(7)
where there is emission of $`\overline{D}_{sJ}(j)`$. Of course, this is only a qualitative statement, because factorization in the decays (7) is just a model and is not in a firm status as in the light meson emission case (6) . To summarize, “broad” resonances ($`j=\frac{1}{2}`$) are expected to be suppressed in the production, while “narrow” resonances ($`j=\frac{3}{2}`$) are expected to be suppressed in the emission. The BaBar experiment has started looking at $`B_d(\overline{D}K)D`$  where such decays with emitted excited states might be seen, but the statistics has to be improved.

Let us now make some further remarks on decay constants of excited states. Imposing duality to $`\mathrm{\Delta }\mathrm{\Gamma }`$ in the $`B_s^0`$-$`\overline{B}_s^0`$ system in the heavy $`b`$, $`c`$ quark limit, the following sum rules have been demonstrated, in the heavy quark limit of QCD  :

$`{\displaystyle \underset{n}{}}{\displaystyle \frac{f^{(n)}}{f^{(0)}}}\xi ^{(n)}(w)=1`$ (9)
$`{\displaystyle \underset{n}{}}{\displaystyle \frac{f_{1/2}^{(n)}}{f^{(0)}}}\tau _{1/2}^{(n)}(w)={\displaystyle \frac{1}{2}}`$
The decay constants for $`S`$-states $`f^{(n)}`$ and for $`j=\frac{1}{2}`$ $`P`$-states $`f_{1/2}^{(n)}`$ are properly defined in ref. . The sum rules (9)-(9) are strong constraints, as they hold for any value of $`w`$.

The $`w`$-dependence in eqns. (9), (9) was obtained by considering the two-body intermediate states of the type $`D_s\overline{D}_s`$ (ground state and excited states) for the diagrams of the spectator quark type. As explained in ref. , in this kind of diagrams the transition $`\overline{B}_sD_s\overline{D}_sB_s`$ occurs by $`\overline{D}_s`$ emission and $`D_s`$ production by the $`\overline{B}_s`$ ($`\overline{s}`$ quark being spectator), followed by $`D_s`$ emission and $`\overline{D}_s`$ production by the $`B_s`$. In the heavy quark limit, the expression for the contribution of one intermediate state is proportional to a quantity of the type $`f_{D_s}^2[\xi (w_c)]^2`$ where $`f_{D_s}`$ is a generic ground state or excited $`D_s`$ meson decay constant and $`\xi (w_c)`$ is a generic Isgur-Wise function $`\xi ^{(n)}(w)`$, $`\tau _{1/2}^{(n)}(w)`$ or $`\tau _{3/2}^{(n)}(w)`$ , taken at the fixed value of $`q^2=m_c^2`$ :

$$w_c=\frac{m_b^2+m_c^2q^2}{2m_bm_c}=\frac{m_b}{2m_c}.$$
(10)
Varying the ratio $`m_b/m_c`$, the $`w`$-dependence in the IW functions appears. Identifying the sum over the exclusive modes with the contribution to the quark box diagram having the same topology , one obtains, considering the vector or the axial weak current, the sum rules (9) and (9).

In an earlier paper , we had demonstrated that indeed duality for $`\mathrm{\Delta }\mathrm{\Gamma }_{B_s}`$ occurs for $`N_c=3`$ in the heavy quark limit in the particular Shifman-Voloshin conditions $`\mathrm{\Lambda }_{QCD}m_bm_cm_b,m_c`$ .

On the other hand, the sum rules (9), (9) were studied  within the relativistic quark models of the Bakamjian-Thomas (BT) type , that satisfy HQET relations for form factors and decay constants. The sum rules are satisfied for different values of $`w`$, although the numerical convergence becomes worst as $`w`$ increases.

From (9), (9), using Schwartz inequality, one can obtain the lower bounds :

$`\left[{\displaystyle \underset{n}{}}\left[\xi ^{(n)}(w)\right]^2\right]\left[{\displaystyle \underset{n}{}}\left({\displaystyle \frac{f^{(n)}}{f^{(0)}}}\right)^2\right]\left({\displaystyle \underset{n}{}}{\displaystyle \frac{f^{(n)}}{f^{(0)}}}\xi ^{(n)}(w)\right)^2=1`$ (11)
$`\left[{\displaystyle \underset{n}{}}\left[\tau _{1/2}^{(n)}(w)\right]^2\right]\left[{\displaystyle \underset{n}{}}\left({\displaystyle \frac{f_{1/2}^{(n)}}{f^{(0)}}}\right)^2\right]\left({\displaystyle \underset{n}{}}{\displaystyle \frac{f_{1/2}^{(n)}}{f^{(0)}}}\tau _{1/2}^{(n)}(w)\right)^2={\displaystyle \frac{1}{4}}`$ (12)
Considering first these inequalities at $`w=1`$, one can see, from $`\xi ^{(n)}(1)=0`$ for $`n0`$, that (11) does not provide any useful constraint. However, from (12) at $`w=1`$ and (4) one obtains the bound for the $`j=\frac{1}{2}`$ $`P`$-wave decay constants

$$\underset{n}{}\left(\frac{f_{1/2}^{(n)}}{f^{(0)}}\right)^2\frac{3}{4\rho ^23}$$
(13)
that contains the IW slope $`\rho ^2`$ in the right-hand side. Although the bound (13) is very weak, it deserves a few comments. Its r.h.s., and also its l.h.s. diverge as $`\rho ^2\frac{3}{4}`$ and it reflects the fact that the excited $`P`$-wave mesons $`D_{0,1}(j=\frac{1}{2})`$ do indeed couple to vacuum, respectively through the vector and axial currents for $`J=0,1`$, as was already proved from the sum rule (9). This is to be contrasted with the selection rule $`f_{3/2}^{(n)}=0`$   that applies to the excited mesons $`D_{1,2}(j=\frac{3}{2})`$. In the example of the non-relativistic quark model, also given as an illustration in ref. , both sides of the inequality (13) are of $`O(v^2/c^2)`$, since $`f_{1/2}^{(n)}`$ is of $`O(v/c)`$ (formula (17) of ) and $`\rho ^2`$ is of $`O(c^2/v^2)`$ (formula (52) of ). In this case, as can be easily seen using completeness, the l.h.s. of (13) is infinite, proportional to the divergent integral $`_0^{\mathrm{}}p^4𝑑p`$, two powers of $`p`$ coming from the current and $`p^2`$ from the measure d$`\stackrel{}{p}`$.

Actually, one can demonstrate this divergence also in field theory, in the heavy quark limit of QCD by considering arbitrary large $`w`$. From Bjorken SR for any value of $`w`$ ,

$$\frac{w+1}{2}\underset{n}{}|\xi ^{(n)}(w)|^2+(w1)\left[2\underset{n}{}|\tau _{1/2}^{(n)}(w)|^2+(w+1)^2\underset{n}{}|\tau _{3/2}^{(n)}(w)|^2\right]+\mathrm{}=1$$
(14)
one obtains, since all terms in this sum are definite positive

$$\underset{n}{}|\xi ^{(n)}(w)|^2\frac{2}{w+1}\underset{n}{}|\tau _{1/2}^{(n)}(w)|^2\frac{1}{2(w1)}$$
(15)
From (11) and (12) for any $`w`$ one gets

$$\underset{n}{}\left(\frac{f^{(n)}}{f^{(0)}}\right)^2\frac{w+1}{2}\underset{n}{}\left(\frac{f_{1/2}^{(n)}}{f^{(0)}}\right)^2\frac{w1}{2}$$
(16)
Since the l.h.s. of these inequalities is independent of $`w`$, that can be made arbitrarily large, one concludes that the sums $`\underset{n}{}\left(\frac{f^{(n)}}{f^{(0)}}\right)^2`$, $`\underset{n}{}\left(\frac{f_{1/2}^{(n)}}{f^{(0)}}\right)^2`$ must diverge.

This situation seems to be realized in quark models à la Bakamjian and Thomas. The decay constants $`f^{(n)}`$ and $`f_{1/2}^{(n)}`$ were computed within the BT quark models for different quark-antiquark potentials up to $`n=10`$, and the convergence of the SR (9) and (9) was studied for different values of $`w`$ . The error on the decay constants induced by a truncation procedure in the calculation increases strongly with $`n`$. The decay constants $`f^{(n)}`$ and $`f_{1/2}^{(n)}`$ do not decrease with increasing $`n`$. For the most sophisticated Godfrey-Isgur potential, that describes all $`q\overline{q}`$, $`q\overline{Q}`$ and $`Q\overline{Q}`$ quarkonia, including spin dependent effects , one obtains the decay constants of the Table for the $`S`$-states and the $`j=\frac{1}{2}`$ $`P`$-states .

| Radial excitation | $`\sqrt{M}f^{(n)}`$ (GeV3/2) | $`\sqrt{M}f_{1/2}^{(n)}`$ (GeV3/2) |
| --- | --- | --- |
| $`n=0`$ | 0.67(2) | 0.64(2) |
| $`n=1`$ | 0.73(4) | 0.66(4) |
| $`n=2`$ | 0.76(5) | 0.71(5) |
| $`n=3`$ | 0.78(9) | 0.73(8) |
| $`n=4`$ | 0.80(10) | 0.76(11) |

Table :  Decay constants $`f^{(n)}`$, $`f_{1/2}^{(n)}`$ for radial excitations $`S`$-states and $`j=\frac{1}{2}`$ $`P`$-states in the GI spectroscopic model. The estimated error is in parentheses , and $`M`$ is the bound state mass.

In view of the values of the decay constants in the Table, besides the qualitative hierarchy (7), one expects for the $`n=0`$ states, assuming factorization, neglecting Penguin diagrams and also the mass differences between the $`P`$ states and the ground state :

$`{\displaystyle \frac{\mathrm{\Gamma }(\overline{B}_dD_{s0}^{}\left(\frac{1}{2}\right)D^+)}{\mathrm{\Gamma }(\overline{B}_dD_s^{}D^+)}}{\displaystyle \frac{\mathrm{\Gamma }(\overline{B}_dD_{s0}^{}\left(\frac{1}{2}\right)D^+)}{\mathrm{\Gamma }(\overline{B}_dD_s^{}D^+)}}1`$
$`{\displaystyle \frac{\mathrm{\Gamma }(\overline{B}_dD_{s1}^{}\left(\frac{1}{2}\right)D^+)}{\mathrm{\Gamma }(\overline{B}_dD_s^{}D^+)}}{\displaystyle \frac{\mathrm{\Gamma }(\overline{B}_dD_{s1}^{}\left(\frac{1}{2}\right)D^+)}{\mathrm{\Gamma }(\overline{B}_dD_s^{}D^+)}}1.`$ (17)
These relations follow from the approximate numerical equality $`f^{(0)}f_{1/2}^{(0)}`$ in the Table and the heavy quark relations among, respectively, the decay constants of $`D`$, $`D^{}`$ and $`D_0\left(\frac{1}{2}\right)`$, $`D_1\left(\frac{1}{2}\right)`$ mesons .

In conclusion, we have underlined a disymmetry in the prediction of the rates of production and emission of $`P`$-wave heavy-light mesons, and we have undertaken a theoretical discussion of decay constants of excited heavy mesons.


format: text

text:
# Monte Carlo Methods for Small Molecule High-Throughput Experimentation

## 1 Introduction

High-throughput synthesis is now established as one of the methods for the discovery of new drugs, materials, and catalysts. High-throughput, or combinatorial, methods allow for simultaneous creation of a large number of structurally diverse and complex compounds, generalizing the traditional techniques of single compound synthesis. Parallel synthesis and split/pool synthesis on solid phase, for example, are two commonly used methods for combinatorial synthesis. After combinatorial synthesis of the desired number of compounds, high-throughput screening is used to identify the few molecules optimally possessing the property of interest. Among the high-throughput methodologies, small molecule combinatorial chemistry is the most developed and has been applied successfully in areas such as transition metal complexation, chemical genetic screening, catalysis, and drug discovery.

High-throughput chemistry can be viewed as a search over a multi-dimensional space of composition variables for molecules possessing a high degree of the desired function, or figure of merit. Indeed, the parallel synthesis and split/pool synthesis methods search the composition space in a regular, grid-like fashion. As the complexity of the molecular library grows, the number of dimensions in the composition variable space grows, and with a grid-like method, the number of compounds that must be synthesized to search the space grows exponentially. Synthesis and screening of mixtures of compounds can partially alleviate the dimensional curse. However, a mixture approach raises the question of how to deconvolute and interpret the results. The greater the degree of mixing, the stronger the synergistic effects can be in the mixture, and the more difficult it is to identify individual compounds responsible for the activity.

The challenge of searching the composition space in an efficient way has led to extensive efforts in the rational design of combinatorial, or high-throughput, libraries. A basic assumption in library design is that structurally similar compounds tend to display similar activity profiles. By designing libraries with maximum structural diversity, the potential for finding active compounds in the high-throughput screenings can be enhanced. This design approach requires a quantitative account of the structural and functional diversity of the library, and many descriptors have been developed. Optimization of a library to maximum diversity is then driven by a reliable statistical method. Several structurally diverse libraries have been successfully designed along these lines. For example, strategies have been presented to optimize the structural diversity of libraries of potential substituents or entire molecules by using stochastic optimization of diversity functions and a point mutation Monte Carlo technique. Peptide libraries have been designed by using topological descriptors and quantitative structure-activity relationships combined with a genetic algorithm and simulated annealing. Diverse libraries of synthetic biodegradable polymers have been designed by using molecular topology descriptors and a genetic algorithm. Similarly, peptoid libraries have been designed by using multivariate quantitative structure-activity relationships and statistical experimental design.

The question of how an initial library should be redesigned for subsequent rounds of high-throughput experimentation in light of the results of the first round of screening remains unanswered. In this paper, we suggest that Monte Carlo methods provide a natural means for library redesign in high-throughput experimentation. The Monte Carlo method is a well-known statistical method for sampling large spaces efficiently and ergodically. There is a striking analogy between searching configuration space for regions of low free energy in a Monte Carlo simulation and searching composition space for regions of high figure of merit in a combinatorial chemistry experiment. Importantly, Monte Carlo methods do not suffer the curse of dimensionality. A Monte Carlo approach should, therefore, be exponentially more efficient than a regular, grid-like method for libraries of complex molecules. Indeed, a Monte Carlo approach to materials discovery proves to be dramatically more efficient than does a grid-based approach. The application of Monte Carlo methods to small molecule high-throughput experimentation differs from the conventional computer simulation technique in several aspects. First, the variables in molecular high-throughput experimentation are discrete, and no continuous moves are available. Second, multiple simultaneous searches of the variable space are performed in high-throughput experimentation when screening the large libraries. Finally, temperature has a natural meaning in Monte Carlo computer simulation, whereas “temperature” is simply a control parameter in a Monte Carlo protocol for high-throughput experimentation. In principle, the temperature in the protocol serves to specify how strong is the differentiation between compounds with low and high figures of merit. In practice, temperatures that are too low may cause the method to become trapped in local optima unless a sufficiently powerful Monte Carlo scheme is used.

Despite the many successes of high-throughput experimentation, the method has been criticized as a simple machinery, lacking incorporation of a priori knowledge when compared with the traditional synthetic approach. A priori knowledge, such as chemical intuition, previous database or experimental information, well-known theory, patentability, or other specific constraints, are indispensable to an efficient library design and are the traditional province of the synthetic chemist. Fascinatingly, the Monte Carlo approach to high-throughput experimentation can naturally incorporate such knowledge in the experimental design through the technique of biased Monte Carlo.

Genetic algorithms are the computational analog of Darwinian evolution. Typically, a genetic algorithm consists of three basic processes: crossover, mutation, and selection. In the crossover step, new compounds are generated by mixing the compositions of parent compounds. In the mutation step, individual molecules are changed at random. In the selection step, the best molecules are identified for the next round. The application of genetic algorithms to combinatorial synthesis and library design has achieved considerable success. Nonetheless, unlike Monte Carlo algorithms, genetic algorithms do not satisfy detailed balance. Because of this, genetic algorithms can not be guaranteed to sample properly the variable space or to locate optimal molecules. Furthermore, in most experiments, one wants to identify several initially promising molecules in the hope that, among them, a few can survive further stringent screenings, such as patentability or lack of side effects. In the genetic approach, however, all the molecules in the library tend to become similar to each other due to the crossover step. While diversity can be encouraged in a genetic approach, diversity can never be guaranteed. The Monte Carlo approach, on the other hand, can maintain or even increase the diversity of a molecular library, due to the satisfaction of detailed balance.

In this paper, we propose several strategies for small molecule high-throughput experimentation derived by analogy with Monte Carlo methods. We compare these Monte Carlo protocols to the genetic algorithm approach. In order to make this comparison and to demonstrate the effectiveness of the Monte Carlo approach, we perform simulated high-throughput experiments. In section 2, we introduce a random energy model that we use as a surrogate for experimental measurement of the figure of merit. The random energy model is not fundamental to the protocols; it is introduced as a simple way to test, parameterize, and validate the various searching methods. In an experimental implementation, the random energy model would be replaced by the value returned by the screen. In section 3, we introduce the Monte Carlo protocols and provide a means to calculate the diversity of a library. In section 4, we compare the various protocols. In section 5, we discuss some implications of these results. We conclude in section 6.

## 2 Space of Variables and A Random Energy Model for Small Molecule High-Throughput Experimentation

To quantitatively describe the molecules in a high-throughput library, we uniquely characterize each molecule by its composition, such as the identity of the core and substituents. For specificity, we will consider the figure of merit of interest to be a binding constant, but our results should be generically valid. A schematic view of our model is presented in Figure 1. For simplicity, we consider the small molecule to consist of one core, drawn from a library of cores, and six binding substituents, each drawn from a single library of substituents. Numerous energetic interactions could exist between this molecule and the substrate. It is commonly believed that descriptors can be directly related to compound performance. A large class of descriptors, such as one-dimensional, two-dimensional, three-dimensional, and BCUT descriptors, has been used to measure the diversity between substituents, cores, and molecules in the literature. To simplify, we will limit ourselves to a set of six weakly correlated descriptors for each substituents and core. For example, the descriptors could be hydrogen bond donors, hydrogen bond acceptors, flexibility, an electro-topological calculation, clog$`P`$, and aromatic density.

To carry out the simulated experiments, we need a figure of merit function that mimics the experimental step of measuring the figure of merit. Once constructed or “synthesized,” the molecules are scored by the model, which takes the composition or molecular descriptors as input. A random energy model can mimic the generic features of an experimental figure of merit. For example, the NK model is used to model combinatorial chemistry experiments on peptides, the block NK and generalized NK models are used to model protein molecular evolution experiments, and the random phase volume model is used to model materials discovery.

The basic building block for our random energy model is a random polynomial of $`n`$ descriptors, $`x_1,\mathrm{},x_n`$:

$$F(x_1,\mathrm{},x_n,\{G\})=\underset{k=0}{\overset{q}{}}\underset{\stackrel{i_1+\mathrm{}+i_n=k}{i_1,\mathrm{},i_n0}}{}f_{i_1\mathrm{}i_n;k}^{\frac{1}{2}}\xi ^kG_{i_1\mathrm{}i_n}x_1^{i_1}\mathrm{}x_n^{i_n}$$
(1)
Here $`q`$ is the degree of the polynomial, $`G_{i_1\mathrm{}i_n}`$ are the fixed coefficients of the polynomial, and $`f_{i_1\mathrm{}i_n;k}`$ are constant symmetry factors,

$$f_{i_1\mathrm{}i_n;k}=\frac{k!}{i_1!\mathrm{}i_n!}$$
(2)
We choose to take the square root of $`f`$ here because we consider each term of the unsymmetrized polynomial to be random. The scale factor $`\xi `$ is used to equalize roughly the contributions from each term of the polynomial. Since $`x_i`$ will be drawn from a Gaussian random distribution of zero mean and unit variance, we set

$$\xi =\left(\frac{x^q}{x^2}\right)^{\frac{1}{q2}}=\left(\frac{q!}{(q/2)!2^{q/2}}\right)^{\frac{1}{q2}}$$
(3)
We use $`q=6`$ and $`n=6`$ in our random energy model.

The random energy model accounts for contributions to substrate binding arising from interactions between the substrate and core and from interactions between the substrate and each of the substituents. In addition, synergistic effects between the substituents and core are incorporated. Consider, for example, a molecule made from core number $`m`$ from the core library and substituents number $`s_1,\mathrm{},s_6`$ from the substituent library. The core is characterized by six descriptors, $`D_1^{(m)},\mathrm{},D_6^{(m)}`$. Similarly, each substituent is characterized by six descriptors, $`d_1^{(s_i)},\mathrm{},d_6^{(s_i)}`$. We denote the core contribution to binding by $`E_\mathrm{C}`$ and the substituent contributions by $`E_\mathrm{S}`$. We denote the contribution due to synergistic substituent-substituent interactions by $`E_{\mathrm{SS}}`$ and the contribution due to synergistic core-substituent interactions by $`E_{\mathrm{CS}}`$. The total contribution to the figure of merit is, then,

$$E=E_\mathrm{S}+E_\mathrm{C}+E_{\mathrm{SS}}+E_{\mathrm{CS}}$$
(4)
Each of these factors is given in terms of the random polynomial:

$`E_\mathrm{S}`$ $`=`$ $`\alpha _1{\displaystyle \underset{i=1}{\overset{6}{}}}F(d_1^{(s_i)},\mathrm{},d_6^{(s_i)},\{G_\mathrm{S}\})`$ (5)
$`E_\mathrm{C}`$ $`=`$ $`\alpha _2F(D_1^{(m)},\mathrm{},D_6^{(m)},\{G_\mathrm{C}\})`$ (6)
$`E_{\mathrm{SS}}`$ $`=`$ $`\alpha _3{\displaystyle \underset{i=1}{\overset{6}{}}}h_iF(d_{j_1}^{(s_i)},d_{j_2}^{(s_i)},d_{j_3}^{(s_i)},d_{j_4}^{(s_{i+1})},d_{j_5}^{(s_{i+1})},d_{j_6}^{(s_{i+1})},\{G_{\mathrm{SS}}\})`$ (7)
$`E_{\mathrm{CS}}`$ $`=`$ $`\alpha _4{\displaystyle \underset{i=1}{\overset{6}{}}}h_iF(d_{k_1}^{(s_i)},d_{k_2}^{(s_i)},d_{k_3}^{(s_i)},D_{k_4}^{(m)},D_{k_5}^{(m)},D_{k_6}^{(m)},\{G_{\mathrm{CS}}\})`$ (8)
where the $`\{G_\mathrm{S}\},\{G_\mathrm{C}\},\{G_{\mathrm{SS}}\}`$, and $`\{G_{\mathrm{CS}}\}`$ are four sets of fixed random Gaussian variables with zero mean and unit variance. The $`\alpha _i`$ are constants to be adjusted so that the synergistic terms will contribute in desired percentages, and $`h_i`$ is a structural constant indicating the strength of the interaction at binding site $`i`$. The interaction strengths $`h_i`$ are chosen from a Gaussian distribution of zero mean and unit variance for each site on each core. Only synergistic interactions between neighboring substituents are considered in $`E_{\mathrm{SS}}`$, and it is understood that $`s_7`$ refers to $`s_1`$ in eq 7. In principle, the polynomial in eq 7 could be a function of all 12 descriptors of both substituents. We assume, however, that important contributions come from interactions among three randomly chosen distinct descriptors of substituent $`s_i,d_{j_1}^{(s_i)},d_{j_2}^{(s_i)},`$ and $`d_{j_3}^{(s_i)},`$ and another three randomly chosen distinct descriptors of substituent $`s_{i+1},d_{j_4}^{(s_{i+1})},d_{j_5}^{(s_{i+1})},`$ and $`d_{j_6}^{(s_{i+1})}`$. Similarly, we assume that core-substituent contributions come from interactions between three randomly chosen distinct descriptors of the substituent, $`d_{k_1}^{(s_i)},d_{k_2}^{(s_i)},`$ and $`d_{k_3}^{(s_i)},`$ and another three randomly chosen distinct descriptors of the core, $`D_{k_4}^{(m)},D_{k_5}^{(m)},`$ and $`D_{k_6}^{(m)}`$. Both $`j_i`$ and $`k_i`$ are descriptor indices ranging from 1 to 6. Assuming that we have integrated out the degrees of freedom of the substrate, these indices depend only on the core.

The parameters in the random energy model are chosen to mimic the complicated interactions between a small molecule and a substrate. We choose to focus on the case where these interactions are unpredictable, which is typical. That is, in a typical experiment, it would not be possible to predict the value of the screen in terms of molecular descriptors. Indeed, when rational design fails, an intelligent use of high-throughput experimentation is called for. The task of library design and redesign, rather than single molecule design, is the one we address in the next section.

## 3 Monte Carlo Strategies and Diversity Measurement

Before initiating the Monte Carlo protocol, we first build the core and substituent libraries. We denote the size of the core library by $`N_\mathrm{C}`$ and the size of the substituent library by $`N_\mathrm{S}`$. In a real experiment, the six descriptors would then be calculated for each core and substituent. In the simulated experiment, the values of the six descriptors of each substituent and core are extracted from a Gaussian random distribution with zero mean and unit variance. In the simulated experiment, we also associate two sets of random interaction descriptor indices to each core for the interaction terms in eqs 7 and 8.

To give a baseline for comparison, we first design the library using a random construction. New molecules are constructed by random selection of one core and six substituents from the libraries. Since the properties of each substituent and core are assigned randomly, this first library should be reasonably diverse and comparable to examples in literature.

For the Monte Carlo schemes, the initial molecular configurations are assigned randomly as before. The library is modified by the Monte Carlo protocol in subsequent rounds of high-throughput experimentation. Two kinds of move are possible for each molecule in the library, core changes and substituent changes. Either the core is changed with probability $`p_{\mathrm{core}}`$, or one of the six substituents is picked randomly to change. We denote the probability of changing from core $`m`$ to $`m^{}`$ by $`T(mm^{})`$ and from substituent $`i`$ to $`i^{}`$ by $`t(ii^{})`$. The new configurations are updated according to the acceptance rule at $`\beta `$, the inverse of the protocol temperature. All the samples are sequentially updated in one Monte Carlo round.

For the simple Metropolis method, the transition matrices are

$`T(mm^{})`$ $`=`$ $`1/N_\mathrm{C}`$ (9)
$`t(ii^{})`$ $`=`$ $`1/N_\mathrm{S}`$ (10)
and the acceptance rule is

$$\mathrm{acc}(on)=\mathrm{min}[1,\mathrm{exp}(\beta \mathrm{\Delta }E)]$$
(11)
To make use of the idea that smaller moves are accepted more often, we could try to choose a modified substituent or core that is similar to current one, that is, we could use a transition matrix weighted towards those substituents or cores close to the current one in the six-dimensional descriptor space. Interestingly, this refinement turns out not to work any better than does the simple random move. It seems that even a small move in the descriptor space is already much larger than the typical distance between peaks on the figure of merit landscape.

Biased Monte Carlo methods have been shown to improve the sampling of complex molecular systems by many orders of magnitude. In contrast to conventional Metropolis Monte Carlo, trial moves in biased schemes are no longer chosen completely at random. By generating trial configurations with a probability that depends on a priori knowledge, the moves are more likely to be favorable and more likely to be accepted. As we are dealing with a discrete configurational space, the implementation of biased Monte Carlo in this case is relatively simple. First, we need a biasing term for both substituents and cores. Since the form of this term is not unique, we can proceed in several different ways. One strategy is to bias our choice of core and substituent on the individual contributions of the cores and substituents to the figure of merit. We might know, or be able to estimate, these contributions from theory. For the random energy model, for example, we know

$`e^{(i)}`$ $`=`$ $`\alpha _1F(d_1^{(i)},\mathrm{},d_6^{(i)},\{G_\mathrm{S}\})`$ (12)
$`E^{(m)}`$ $`=`$ $`\alpha _2F(D_1^{(m)},\mathrm{},D_6^{(m)},\{G_\mathrm{C}\})`$ (13)
where $`e^{(i)}`$ is the bias energy to the substituent $`i`$ in the library, and $`E^{(m)}`$ the bias energy of core $`m`$ in the library. Alternatively, we can estimate the contribution of each substituent or core to the figure of merit experimentally. A electrospray ionization source coupled to a mass spectrometer, for example, can serve this purpose. To measure the contributions, we do a pre-experiment on 10000 randomly constructed molecules. This number of compounds will give on average each substituent 60 hits and each core 667 hits. By averaging the figure of merit of the molecules containing a particular substituent or core over the total number of hits, we can obtain experimental estimates of $`e^{(i)}`$ and $`E^{(m)}`$. Using these two methods of bias, we construct three different types of biased Monte Carlo schemes: theoretical biased move, experimental biased move, and mixed biased move. In theoretical bias, both $`e^{(i)}`$ and $`E^{(m)}`$ are from the random energy model. In experimental bias, both $`e^{(i)}`$ and $`E^{(m)}`$ are calculated from the pre-experiment. In mixed bias, $`e^{(i)}`$ comes from the random energy model, while $`E^{(m)}`$ comes from the pre-experiment.

These biases tend to exhibit a large gap between a few dominant cores and substituents and the rest. To ensure the participation of more substituents and cores in the strategy, we introduce cutoff energies for the substituent and core, $`e_\mathrm{c}`$ and $`E_\mathrm{c}`$. We arbitrarily choose $`e_\mathrm{c}`$ to be the 21st lowest substituent energy and $`E_\mathrm{c}`$ to be the 4th lowest core energy. The biased energy, $`e_\mathrm{b}^{(i)}`$, for the $`i^{th}`$ substituent is

$$e_\mathrm{b}^{(i)}=\{\begin{array}{cc}e^{(i)}\hfill & \text{if }e^{(i)}>e_\mathrm{c}\hfill \\ e_\mathrm{c}\hfill & \text{otherwise}\hfill \end{array}$$
(14)
And the biased energy, $`E_\mathrm{b}^{(m)}`$, for the $`m^{th}`$ core is

$$E_\mathrm{b}^{(m)}=\{\begin{array}{cc}E^{(m)}\hfill & \text{if }E^{(m)}>E_\mathrm{c}\hfill \\ E_\mathrm{c}\hfill & \text{otherwise}\hfill \end{array}$$
(15)
To correct for this bias, we introduce Rosenbluth factors. Since the transition probabilities are the same at each Monte Carlo step, we have a constant Rosenbluth factor for the substituent

$$w(n)=w(o)=\underset{i=1}{\overset{N_\mathrm{S}}{}}\mathrm{exp}(\beta e_\mathrm{b}^{(i)})$$
(16)
The probability of transition from substituent $`i`$ to $`i^{}`$ is

$$t(ii^{})=\frac{\mathrm{exp}(\beta e_\mathrm{b}^{(i^{})})}{w(n)}$$
(17)
In the same way, the Rosenbluth factor for the core is

$$W(n)=W(o)=\underset{m=1}{\overset{N_\mathrm{C}}{}}\mathrm{exp}(\beta E_\mathrm{b}^{(m)})$$
(18)
The probability of transition from core $`m`$ to $`m^{}`$ is

$$T(mm^{})=\frac{\mathrm{exp}(\beta E_\mathrm{b}^{(m^{})})}{W(n)}$$
(19)
Finally, we define the remaining, non-biased part of the figure of merit to be

$$E_\mathrm{b}=EE_\mathrm{b}^{(m)}\underset{i=1}{\overset{6}{}}e_\mathrm{b}^{(s_i)}$$
(20)
To satisfy the detail balance, the acceptance rule becomes

$$\mathrm{acc}(on)=\mathrm{min}[1,\mathrm{exp}(\beta \mathrm{\Delta }E_\mathrm{b})]$$
(21)
We add a swap move that attempts to exchange fragments between two molecules to the set of Monte Carlo moves. We denote the probability of attempting a swap instead of a single-molecule move as $`p_{\mathrm{swap}}`$. In a swap move, the cores or a pair of substituents may be swapped between two randomly selected molecules. The probability of switching the core or substituent at the same position is given by $`p_{\mathrm{swap}_\mathrm{C}}`$ and $`p_{\mathrm{swap}_\mathrm{S}}`$, respectively. The crossover event from genetic algorithms could also be introduced in the swap moves, but this additional move did not improve the results. The acceptance rule for swapping is $`\mathrm{acc}(on)=\mathrm{min}[1,\mathrm{exp}(\beta \mathrm{\Delta }E)]`$.

Parallel tempering is known to be a powerful tool for searching rugged energy landscapes. In parallel tempering, the samples are divided into $`k`$ groups. The first group of samples is simulated at $`\beta _1`$, the second group is at $`\beta _2`$, and so on, with $`\beta _1<\beta _2<\mathrm{}<\beta _k`$. At the end of each round, samples in group $`i<k`$ are allowed to exchange configurations with samples in group $`i+1`$ with probability $`p_{\mathrm{ex}}`$. The corresponding acceptance rule for a parallel tempering exchange is

$$\mathrm{acc}(on)=\mathrm{min}[1,\mathrm{exp}(\mathrm{\Delta }\beta \mathrm{\Delta }E)]$$
(22)
where $`\mathrm{\Delta }\beta =\beta _i\beta _{i+1}`$ and $`\mathrm{\Delta }E`$ is the difference in energy between the sample in group $`i`$ and the sample in group $`i+1`$. It is important to notice that this exchange step is experimentally cost-free. Nonetheless, this step can be dramatically effective at facilitating the protocol to escape from local minimum. The number of groups, the number of samples in each group, the value of $`\beta _i`$, and the exchange probability, $`p_{\mathrm{ex}}`$, are experimental parameters to be tuned.

For comparison, we compare these Monte Carlo protocols to a standard genetic algorithm approach. In the genetic algorithm, as in the Monte Carlo strategy, we perform multiple rounds of experimentation on a large set of compounds. The difference between the Monte Carlo and the genetic algorithm lies in how the library is redesigned, that is, how the compounds are modified in each round. In the genetic algorithm, first we randomly select two parents. Then we list the explicit composition of each molecule, i.e. core, substituent 1, …, substituent 6. After aligning the sequences from the two parents, we make a random cut and exchange part of the sequences before the cut. We also allow for random changes, or mutation, in the cores or substituents of the offsprings. Finally, since the population is doubled by crossover, we select the better half of the molecules to survive this procedure and continue on to the next round.

The diversity of the library as it passes through the rounds of high-throughput experimentation is an important quantity. We calculate the diversity, $`𝒟`$, as the standard deviation of the library in the 42-dimensional descriptor space.

$$𝒟^2=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left[\underset{j=1}{\overset{6}{}}(D_j^{(m(i))}D_j)^2+\underset{j=1}{\overset{6}{}}\underset{k=1}{\overset{6}{}}(d_j^{(s_k(i))}d_j^{(s_k)})^2\right]$$
(23)
where $`m(i)`$ is the index of the core of molecule $`i`$, $`s_k(i)`$ is the index of substituent $`k`$ of molecule $`i`$, and $`j`$ is the index for the descriptor. The average value in each descriptor dimension is given by $`D_j=N^1_{i=1}^ND_j^{(m(i))}`$ and $`d_j^{(s_k)}=N^1_{i=1}^Nd_j^{(s_k(i))}`$. The diversity of the library will change as the library changes. A larger library will generally possess a higher absolute diversity simply due to the increased number of compounds. This important, but trivial, contribution to the diversity is scaled out by the factor of $`1/N`$ in eq 23.

## 4 Results

To gauge how the synergistic terms in the figure of merit affect the efficiency of the Monte Carlo protocols, we consider three models with increasingly important synergistic effects. We do this by adjusting the $`\alpha _i`$ in eqs 58, so that the absolute values of the terms are on average in the ratio $`E_\mathrm{S}:E_\mathrm{C}:E_{\mathrm{SS}}:E_{\mathrm{CS}}=1:1:0.5:0.3`$ in model I, $`1:1:1:0.6`$ in model II, and $`1:1:2:1.2`$ in model III. Finally, we set $`\alpha _1=0.01`$ arbitrarily in model I. To maintain the same statistical magnitude of total energy, we set $`\alpha _1=0.00778`$ in model II and $`\alpha _1=0.00538`$ in model III.

The size of the library is fixed at $`N_\mathrm{C}=15`$ and $`N_\mathrm{S}=1000`$. The compositional space of this model has $`15\times 1000^6`$ distinct molecules. Clearly, it is impossible to search exhaustively even this modestly complex space. We fix the total number of molecules to be synthesized at 100000, that is, all protocols will have roughly the same experimental cost. Specifically, 100000 molecules will be made in the random library design protocol, while in the case of the Monte Carlo or genetic protocols, the number of molecules times the number of simulation rounds is kept fixed at 100000.

To locate optimal parameters for the protocols, we perform a few short pre-experiments. We first fix the energy coefficients in the energy function and the descriptors of the substituent and core libraries. For simple Metropolis, we find that it is optimal to use 10 samples with 10000 rounds, suggesting that the system is still far from equilibrium at the random initial configuration. With the biased Monte Carlo method, we find that 100 samples and 1000 rounds is optimal. We focus on systems with 1000 or 100 rounds, since fewer rounds are typically preferred in experiments. It is more difficult to achieve effective sampling in the system with 100 rounds, and so we use this system when setting optimal parameter values. For parallel tempering, it was optimal to have the samples divided into three subsets, with 30% of the population at $`\beta _1`$, 40% at $`\beta _2`$, and 30% at $`\beta _3`$. The optimal parameters are listed in Table I for each model. Determination of these parameter values corresponds experimentally to gaining familiarity with the protocol on a new system.

The various Monte Carlo schemes are compared with the random selection method and the genetic algorithm. Once the optimal parameters are chosen, the coefficients of the energy function and the descriptor values of the substituent and core libraries are generated differently in each simulated experiment. The simulation results for the three models are shown in Figures 24. Each data point in the figure is an average over 20 independent runs. This averaging is intended to give representative performance of the protocols on various figures of merit of experimental interest. Since there is much randomness in the results, the standard deviation of the average is shown as well.

## 5 Discussion

Although the average absolute values of the figures of merit in the three models are adjusted to be equal, the stronger the interaction terms, the greater the figure of merit we can find. For instance, the biased schemes find values of $`E`$ in the range 30–40 in model I, but values in the range 40–50 in model II and values in the range 50–60 in model III. This suggests that the figure of merit landscape has changed in detail as the synergistic effects in the model are adjusted. It is clear that for all systems, the Metropolis methods perform better than does random selection. The system with 1000 molecules and 100 rounds is not well-equilibrated by the Metropolis schemes, and an experiment with 100 molecules and 1000 rounds significantly improves the optimal compounds identified. However, by incorporating a priori knowledge, the biased Monte Carlo schemes are able to equilibrate the experiment with either 1000 or 100 rounds. Interestingly, the theoretical bias and experimental bias methods yield similar results. This strongly suggests that a minimal number of pre-experiments can be very useful, both for the understanding of the structure of the figure of merit landscape and for improving the performance in future rounds.

The results produced with the composite moves including swap and parallel tempering are slightly improved relative to those from the plain Monte Carlo schemes. Typically, however, these composite moves significantly improve the sampling of a rough landscape. Indeed, swapping and crossover moves are very effective in protein molecular evolution, where the variable space is extremely large. Perhaps the variable space is not so large in small molecule high-throughput experimentation that these composite moves are required. Alternatively, the random energy model may underestimate the ruggedness of the landscape. The landscape for RNA substituents, for example, is estimated to be extremely rough, and composite moves may prove more important in this case.

The genetic algorithm is relatively easy to use. It does not satisfy detailed balance, however, so there is no theoretical guarantee of the outcome. The optimal figures of merit identified are, nonetheless, comparable to those from the better Monte Carlo methods for all three models. However, due to the crossover and selection steps in the genetic algorithm, the molecules in the library tend to become similar to each other, which prevents this scheme from sampling the whole variable space. To help elucidate this point, diversity measurements for model I are shown in Figure 5. It is clear that the genetic algorithm has reduced the diversity of the library by 60% relative to the biased schemes. Interestingly, the Monte Carlo simulations actually increase the diversity from the initial random configurations. The biased schemes tend to bring the system to equilibrium relatively quickly, and the diversity measurements are similar for the 100 and 1000 round experiments. For the Metropolis method, on the other hand, an experiment with 100 rounds is less diverse than an experiment with 1000 rounds. The genetic algorithm approach finds less favorable figure of merit values in the 100 compound 1000 round experiment, presumably due to a greater sensitivity to the $`\sqrt{10}`$ reduction in the absolute diversity relative to the already small absolute diversity in the 1000 compound 100 round experiment.

The greater the number of potentially favorable molecules in the library space, the greater the diversity of the experimental library will be for the Monte Carlo methods. The genetic algorithm, on the other hand, will tend to produce a library that contains many copies of a single favorable molecule. A key distinction, then, is that a Monte Carlo strategy will sample many compounds from the figure of merit landscape, whereas a genetic algorithm will tend to produce a single molecule with a favorable figure of merit value. How strongly the compounds with high figures of merit are favored in the Monte Carlo strategy is determined by the protocol temperature, since the probability of observing a compound with figure of merit $`E`$ is proportional to $`\mathrm{exp}(\beta E)`$. The sampling achieved by the Monte Carlo methods is important not only because it assures that the composition space is thoroughly sampled, but also because it assures that the library of final hits will be as diverse as possible.

The Monte Carlo methods perform equally well on all three models. The three models were introduced to gauge the impact of unpredictable, synergistic effects in the experimental figure of merit. It might be expected that the a priori bias methods would perform less well as the synergistic effects become more pronounced. That the biased methods perform well even in model III suggests that the Monte Carlo approach may be rather robust. In other words, even a limited amount of *a priori* information is useful in the Monte Carlo approach to library redesign.

## 6 Conclusion

Monte Carlo appears to be a fruitful paradigm for experimental design of multi-round combinatorial chemistry, or high-throughput, experiments. A criticism of high-throughput experimentation has been its mechanical structure and lack of incorporation of a priori knowledge. As shown here, a biased Monte Carlo approach handily allows the incorporation of a priori knowledge. Indeed, our simulation results reveal that biased Monte Carlo schemes greatly improve the chances of locating optimal compounds. For the moderately complex libraries considered here, the bias can be determined equally well by experimental or theoretical means. Although the compounds identified from a traditional genetic algorithm are comparable to those from the better Monte Carlo schemes, the diversity of identified molecules is dramatically decreased in the genetic approach. Genetic algorithms, therefore, are less suitable when the list of good molecules is further winnowed by a secondary screen, a tertiary screen, patentability considerations, lack of side effects, or other concerns. Interestingly, composite Monte Carlo moves such as swap or parallel tempering bring only a slight improvement to the plain biased Monte Carlo protocols, possibly due to the relatively small size of the composition space in small molecule high-throughput experimentation. Presumably, as the complexity of the library is increased, these composite moves will prove more useful for the more challenging figures of merit. Although we have here chosen the initial library configurations at random, the sophisticated initial library design strategies available in the literature can be used, and they would complement the multi-round library redesign strategies presented here.

## Acknowledgments

This research was supported by National Science Foundation through grant number CTS–9702403.


text:
# Influence of oxygen isotope exchange on the ground state of manganites

## I INTRODUCTION

Manganese oxides of general formula R1-xAxMnO3 (where R is a rare-earth ion and A is an alkaline-earth ion) have attracted considerable interest due to their various fascinating properties. They exhibit a rich phase diagram, which includes ferromagnetic (FM) metallic, antiferromagnetic (AFM) insulating, charge and orbitally ordered phases depending on doping $`x`$, average ionic radius $`<r_\mathrm{A}>`$ on R and A sites, magnetic field, and temperature. In manganites, charge, spin, orbital, and lattice degrees of freedom are interconnected, and their balance determines the electronic and magnetic state of these materials. In some compositions a competition of opposing interactions associated with these different degrees of freedom induces an electronic phase separation, which is considered to be an intrinsic properties of these manganites . Another prominent manifestation of the spin-charge-lattice interplay in manganites is the unusual and large oxygen isotope effect . Substitution of 16O by 18O lowers significantly the Curie temperature $`T_\mathrm{C}`$ of the ferromagnetic compositions of La1-xCaxMnO3 , increases the charge ordering (CO) transition temperature $`T_{\mathrm{CO}}`$ for the transitions from FM metallic to AFM charge-ordered state  and even induces a metal-insulator transition . Many models have been proposed to explain the isotope effect in manganites: a small polaron model , a bipolaron model , change of interatomic distance by lattice vibrations , nonadiabatic behavior of the oxygen ions , isotope dependence of the nonstoichiometry in manganites . In spite of significant experimental and theoretical effort, a physical picture that would consistently explain the influence of oxygen isotope substitution on different states (FM and CO) is presently lacking. To better understand the nature of the oxygen isotope effect in manganites, it is important to know how it affects the ground state. Studies of low temperature specific heat can address this issue, since the specific heat carries information about principal excitations.

In this paper we report a systematic study of specific heat, magnetization and resistivity in FM metallic composition of La1-xCaxMnO3 ($`x=0.2`$, 0.3, and 0.375), charge-ordered La0.5Ca0.5MnO3, as well as (La1-yPry)0.67Ca0.33MnO3 ($`0<y<1`$). The latter compound changes from FM metal for $`y=0`$ to CO AFM insulator for $`y=1`$ . For the compositions $`y<0.75`$ the (La1-yPry)0.67Ca0.33MnO3 exhibits a percolative transition from a higher temperature CO insulator to a metal with decreasing temperature . Analysis of low temperature magnetization, resistivity , and synchrotron x-ray diffraction data  indicates metallic and CO phase coexistence in the $`y<0.75`$ compositions of these materials. One goal of this paper is to compare the low temperature behavior of these various samples in order to investigate the influence of oxygen isotope exchange on the properties of both single phase and phase separated manganites. We have found that the most prominent isotope effect is observed in the phase separated systems, which might be the key to understanding the origin of the giant oxygen isotope effect in manganites.

## II EXPERIMENTAL DETAILS

Ceramic samples of La1-xCaxMnO3 ($`x=0.2`$, 0.3, 0.375, and 0.5) and (La1-yPry)0.67Ca0.33MnO3 ($`0<y<1`$) were prepared by a standard solid state reaction technique. X-ray powder diffraction showed that all samples are single phase and good quality. For oxygen isotope exchange we used a standard procedure : the 16O and 18O samples were prepared from the same pellet of the starting material, and were simultaneously treated at T = 1000 C and pressure 1 atm in different closed quartz tubes, one filled with 16O and the other with 18O. We estimate the extent of the exchange of 16O by 18O to be at least 85 %, as determined by weight change.

The specific heat was measured in the temperature range 2-20 K and magnetic field range 0-8.5 T by relaxation calorimetry. In the temperature range 20-300 K the specific heat was measured by a Quantum Design PPMS. The specific heat measurements have an absolute accuracy of $`\pm 3`$%. The magnetization was measured with a commercial SQUID magnetometer in magnetic field range 0-5.5 T and with a Quantum Design PPMS in fields up to 9 T. Resistivity was measured by a standard four-probe technique.

## III RESULTS

### A Low temperature specific heat of La1-xCaxMnO3 with 16O and 18O

In our samples, substitution of 16O by 18O decreases the Curie temperature (determined from the temperature dependence of the magnetization in a magnetic field of 50 G) of La1-xCaxMnO3 significantly: TC(16O) - TC(18O) = 17.7 $`\pm 0.5`$ K, 11 $`\pm 0.5`$ K, 6 $`\pm 0.5`$ K for $`x=0.2`$, 0.3, and 0.375 respectively. The $`x`$ dependence of the oxygen isotope effect and its values are consistent with prior work . Our low temperature specific heat measurement of the 16O and 18O samples of these compositions is the same within the experimental error ($`\pm 2`$%) in the temperature range from 2 K to 10 K, as shown in Fig. 1a,b for $`x=0.2`$ and 0.3. Since these materials are in the FM metallic state at low temperatures we fit the $`x=0.2`$ and 0.3 data to a form $`C=\gamma T+\beta T^3+\delta T^{3/2}`$ for $`4.2<T<8`$ K, where $`\gamma T`$ is the charge carrier contribution, $`\beta T^3`$ the phonon contribution, and $`\delta T^{3/2}`$ the ferromagnetic spin-wave contribution. This fitting range was chosen, because for $`T<4`$ K the hyperfine contribution from the nuclear magnetic levels of Mn ions may appear, and for $`T>\mathrm{\Theta }_\mathrm{D}/50`$, where $`\mathrm{\Theta }_\mathrm{D}`$ is the Debye temperature, the higher terms of the lattice expansion ($`\beta _5T^5`$, $`\beta _7T^7`$ etc.) could be present. We found that values of the respective fitting parameters for 16O and 18O samples ($`\gamma _{16}=5.7\pm 0.1`$ mJ/mole K2, $`\gamma _{18}=5.6\pm 0.1`$ mJ/mole K2, $`\beta _{16}=0.178\pm 0.003`$ mJ/mole K4, and $`\beta _{18}=0.177\pm 0.003`$ mJ/mole K4 for $`x=0.2`$ samples and $`\gamma _{16}=5.8\pm 0.1`$ mJ/mole K2, $`\gamma _{18}=5.8\pm 0.1`$ mJ/mole K2, $`\beta _{16}=0.180\pm 0.003`$ mJ/mole K4, and $`\beta _{18}=0.177\pm 0.003`$ mJ/mole K2 for $`x=0.3`$ samples) coincide within experimental error. The best fit requires $`\delta =0`$. As was noted in previous work , it is difficult to resolve the ferromagnetic spin-wave contribution to the specific heat in FM metallic manganites due to its small value and the presence of the $`\gamma T`$ contribution.

Although we cannot resolve the spin-wave contribution from our data, it should not affect our conclusions, since inelastic neutron scattering  shows that the spin dynamics remains the same for 16O and 18O La1-xCaxMnO3 in this composition range at low temperatures. The values of the charge carrier ($`\gamma T`$) and lattice contribution ($`\beta T^3`$) are close to that found in other FM metallic manganites .

Although we did not observe a change in the lattice contribution to the specific heat in the temperature range $`4.2<T<8`$ K, the 18O samples should have a larger lattice contribution, since the frequencies of the lattice vibrations should be smaller for the heavier oxygen. We do not expect any measurable contribution from optical modes in this temperature range. For estimates of the changes in the lattice contribution with oxygen isotope substitution we consider $`\mathrm{\Theta }_\mathrm{D}M^{1/2}`$, and therefore $`\beta M^{3/2}`$. This change in $`\beta `$ corresponds to only a 2% (below our resolution limit) increase of total specific heat for $`T=6`$ K, since the specific heat at low temperatures is dominated by the charge carrier contribution. However, at higher temperatures, where the specific heat is dominated by the lattice contribution, the specific heat of the 18O sample is larger, than that of the 16O sample (Fig. 1a, inset), which corresponds to a larger lattice contribution. The magnitude of this increase is consistent with estimates for the change in the lattice specific heat described above. After oxygen back-exchange (16O $``$ 18O; 18O $``$ 16O) the 16O -18O pairs of samples showed the same results: a significant change in $`T_\mathrm{C}`$ and no change in the low temperature specific heat within experimental error.

In the charge-ordered compound La0.5Ca0.5MnO3 oxygen isotope substitution strongly increases the charge ordering temperature , which coincides for this compound with the temperature of the transition from FM to AFM state: TCO(16O) - TCO(18O) = 11.7 $`\pm 0.5`$ K (determined from the temperature dependence of the magnetization measured in $`B=50`$ G and the temperature dependence of the resistivity), but the low temperature specific heat remains unchanged (Fig. 1c). Thus, for single phase metallic and CO compositions the specific heat does not show any change in principal excitations (except for the expected increase of the lattice contribution for 18O samples which we observed for $`T>11`$ K, as mentioned above), which indicates that the ground state of these materials remains unchanged after oxygen isotope substitution.

### B Low temperature specific heat of (La1-yPry)0.67Ca0.33MnO3 ($`0<y<1`$)

To study the influence of oxygen isotope substitution on the materials with metallic and CO phase coexistence we chose (La1-yPry)0.67Ca0.33MnO3 ($`0<y<1`$). Before presenting the experimental results for 16O - 18O samples, it is important to understand how compositional ($`y`$) variations affect the low temperature properties. Fig. 2 shows the temperature dependence of the resistivity, $`\rho (T)`$, of our (La1-yPry)0.67Ca0.33MnO3 ($`0<y<1`$) samples. The magnetization values at $`\mu _0H=1`$ T and $`T=5`$ K (after zero field cooling (ZFC)) for different $`y`$ is shown in the inset of Fig. 2. For the $`y=0.4`$, 0.5, and 0.6 compositions the temperature dependence of the resistivity changes from insulating at higher temperatures to metallic at low temperatures. The magnetization values at $`\mu _0H=1`$ T (slightly above saturation field, but below the field where the CO is affected by the magnetic field) decreases with increase of $`x`$ indicating that the fraction of the FM phase decreases for compositions richer in Pr (Fig. 2, inset). This behavior of the resistivity and magnetization is consistent with previous work .

The low temperature specific heat of (La1-yPry)0.67Ca0.33MnO3 ($`0<y<1`$) is shown in Fig. 3. The specific heat of La0.7Ca0.3MnO3 in the temperature range from 4 K to 19 K is given by

$$C=\gamma T+\beta T^3+\beta _5T^5$$
(1)
where the higher term of the lattice contribution $`\beta _5T^5`$ ($`\beta _5=0.08\pm 0.001`$ $`\mu `$J/mole-K6) is required to fit the data in this temperature range. We express the specific heat of Pr0.65Ca0.35MnO3 in the temperature range from 4 K to 12 K in the following form :

$$C=\gamma T+\beta T^3+C^{}(T)$$
(2)
The first term in this electrically insulating sample originates from spin and charge disorder . The second term represents the lattice and the AFM spin-wave contribution to the specific heat. The third term, $`C^{}(T)`$, is an anomalous contribution present only in the CO state , which manifests itself as upward curvature in the $`C/T`$ vs. $`T^2`$ plot. All the samples at intermediate compositions ($`y=0.4`$, 0.5, and 0.6) exhibit the $`C^{}(T)`$ term even though the low temperature resistivity of $`y=0.4`$

and 0.5 samples is similar to metallic La0.7Ca0.3MnO3. This indicates a presence of CO for metallic (at low T) compositions of (La1-yPry)0.67Ca0.33MnO3 ($`y=0.4`$, 0.5, and 0.6) and, therefore, coexistence of CO insulating and

metallic phases. The presence of charge ordering below the metal-insulator transition was also observed by synchrotron x-ray diffraction in $`y0.5`$ .

The specific heat of (La1-yPry)0.67Ca0.33MnO3 increases consistently with increase of the Pr content $`y`$, suggesting the increase of the volume fraction of the CO phase. To test this assumption we consider the fraction of the metallic phase in (La1-yPry)0.67Ca0.33MnO3 to be proportional to the saturation magnetization value . Since the magnetization is almost temperature independent in the temperature range from 2 K to 20 K, the fractions of the FM metallic phase $`f_{\mathrm{met}}=M(y)/M(0)`$ in this temperature region are $`0.80\pm 0.02`$, $`0.62\pm 0.01`$, and $`0.45\pm 0.01`$ for $`y=0.4`$, 0.5, and 0.6 respectively (Fig. 2, inset). If we suppose that the FM fraction has the specific heat of La0.7Ca0.3MnO3 ($`C_{\mathrm{met}}`$) and the CO fraction has the specific heat of Pr0.65Ca0.35MnO3 ($`C_{\mathrm{CO}}`$), then the specific heat of such a two phase system can be expressed as $`C_{\mathrm{tp}}(T)=f_{\mathrm{met}}C_{\mathrm{met}}+(1f_{\mathrm{met}})C_{\mathrm{CO}}`$. These calculated $`C_{\mathrm{tp}}(T)`$ curves for $`y=0.4`$, 0.5, and 0.6 samples are shown in Fig. 3a as dashed lines. These curves show a reasonable quantitative agreement with the experimental data, which argues strongly for the coexistence of CO and metallic phases at low temperatures in these materials. However, this approach does not give an exact agreement with the experimentally observed specific heat. A fit of the data to the form of Eq. (2) gives much better agreement with the data as shown by the solid lines in Fig. 3b. The values of the fitting parameters are shown in Table I, where we consider $`C^{}(T)`$ to be a contribution from nonmagnetic excitations with dispersion relation $`ϵ=\mathrm{\Delta }+Bq^2`$ ($`\mathrm{\Delta }`$ is an energy gap and $`q`$ is a wave vector) . These results show that the specific heat of $`y=0.4`$, 0.5, and 0.6 samples is not just the linear combination of $`C_{\mathrm{met}}`$ and $`C_{\mathrm{CO}}`$. Since the phase separation in (La1-yPry)0.67Ca0.33MnO3 is electronic (not chemical), the presence of both La and Pr ions could modify the metallic and CO phase in these materials, and CO and metallic phases in phase separated regions of (La1-yPry)0.67Ca0.33MnO3 differ from CO and metallic phases of Pr0.65Ca0.35MnO3 and La0.7Ca0.3MnO3 respectively. The following analysis of the specific heat of (La1-yPry)0.67Ca0.33MnO3 samples will be done using Eq.(2), since the volume fraction approach gives less close agreement with the data.

### C Specific heat, resistivity, and magnetization of (La1-yPry)0.67Ca0.33MnO3 with 16O and 18O

#### 1 Low temperature specific heat and resistivity

The low temperature specific heat and resistivity of $`y=0.4`$ and 0.5 samples with 16O and 18O is shown in Fig. 4. The resistivity of these compositions is affected significantly by the oxygen isotope substitution (Fig. 4, insets): the onset of metallic behavior (the peak temperature for cooling run) decreases by 50 K for the $`y=0.5`$ sample and by 28 K for the $`y=0.4`$ sample; the low temperature value of the resistivity increases approximately by one order in magnitude for both compositions. The

low temperature resistivity of metallic La1-xCaxMnO3 did not show such dramatic effects . This behavior is consistent with previous work: substitution of 16O by 18O favors the insulating state .

Unlike in the single phase metallic and CO compositions of La1-xCaxMnO3 (Fig. 1), in phase separated (La1-yPry)0.67Ca0.33MnO3 the substitution of 16O by 18O induces substantial changes in the low temperature specific heat (Fig. 4). The specific heat, dominated by the $`C^{}(T)`$ contribution in most of the temperature range, is larger for the 18O samples. This increase of the specific heat for the 18O sample is similar to the increase of the specific heat with an increase of the Pr content $`y`$. The specific heat of the $`y=0.6`$ sample is shown in Fig. 4a for comparison. A result of the fits of our data to Eq. (2) is given in Table I. All fitting parameters are close to those of the higher $`y`$ value samples. These results indicate that the amount of the CO phase is larger in 18O samples.

#### 2 High temperature specific heat and magnetization

The specific heat of $`y=0.5`$ 16O and 18O samples in the temperature range from 2 K to 270 K is shown in Fig. 5a. The room temperature value of the specific heat, 110 J/mole K, is in agreement with the room temperature values reported previously in manganites . For both samples we observed a broad anomaly associated with transition to a magnetically and charge ordered state at $`T220`$ K. Similar behavior was reported for an 16O sample of similar composition , where the anomaly in C(T) was attributed to the CO detected by x-ray diffraction. The temperature dependence of the magnetization measured in $`\mu _0H=50`$ G and the inverse susceptibility is shown in the inset of Fig. 5a. The onset of the magnetic transition at $`220`$ K in both samples is best seen from the temperature dependence of the inverse susceptibility. The high temperature inverse susceptibility of both samples (which extrapolates to intersect the x-axis at positive value) indicates a FM interaction. However, for the 18O sample, the value of the inverse susceptibility increases (the magnetization drops) just below the transition at 220 K, which is indicative of a transition to the AFM state. This suggests a presence of competing FM and AFM interactions leading to a two phase coexistence in these materials. Ferromagnetic ordering in the 18O sample is established below $`100`$ K. The 16O sample exhibits a higher value of the magnetization below 220 K than expected for the AFM state, perhaps due to the presence of a fraction of the FM phase. The magnetization of the 16O sample increases further below 120 K indicating that a larger volume fraction becomes FM. These magnetization data of the 16O sample are consistent with , where the transition to a FM state from the CO AFM state and corresponding insulator to metal transition was shown to be of a percolative nature. The substitution of 16O by 18O does not change the temperature of the magnetic and CO transitions in this material but decreases the volume fraction of the FM phase. The specific heat does not show an anomaly associated with the metal-insulator transition for both samples, which is consistent with percolative nature of the transition occurring gradually over a wide temperature range.

The low temperature specific heat (2 K$`<T<`$19 K) of the 18O sample is higher than that of 16O sample (Fig. 4a). This tendency persists also at higher temperatures (see Fig. 5a,b). This difference in specific heat of 18O and 16O samples is not primarily caused by an enhanced lattice contribution from the heavier atom: the difference in the lattice contributions at $`T=70`$ K, estimated from the Debye function, should be $`1`$%, while the observed difference is $`6`$%. A similar behavior is found for La0.5Ca0.5MnO3, which has an excess specific heat compared with metallic La0.7Ca0.3MnO3 (see Fig. 5c). This implies that the excitations responsible for the excess specific heat in the CO state  are also contributing at higher temperatures. In the case of (La1-yPry)0.67Ca0.33MnO3 ($`x=0.5`$) samples, the larger specific heat of the 18O sample (which comes from the excitations in the CO states) indicates a larger volume fraction of the CO phase over a wide temperature range.

#### 3 Effect of a magnetic field

The field dependence of the magnetization for $`y=0.5`$ and 0.4 16O and 18O samples at $`T=5.5`$ K is shown in Fig. 6. Samples were cooled in zero magnetic field (ZFC) and the hysteresis loops in a magnetic field up to 9 T were taken. After ZFC the samples appear to be in a state with smaller than full spin alignment magnetization value. This indicates the FM and AFM phase coexistence in these materials. The values of the ZFC magnetization at $`\mu _0H=1`$ T (where the saturation for the FM phase is reached) are smaller for the materials with 18O which corresponds to the larger fraction of the AFM phase in the 18O samples. At $`\mu _0H1.5`$ T the magnetic field starts to affect the CO (and associated AFM ordering) which causes an increase of the magnetization. In the 18O samples the transition is completed in magnetic fields higher than in 16O samples. Even when the transition is completed ($`\mu _0H>8`$ T), the magnetization value of the 18O samples is somewhat smaller than that of the 16O samples. These magnetization measurements show that substitution of 16O by 18O not only increases the fraction of the AFM CO phase, as shown by specific heat, but also lowers its energy, since it requires a higher magnetic field to stabilize the FM state than for the 16O samples. As in other CO compositions, once the FM phase is stabilized, the material remains in FM state even though the magnetic field is removed, and only heating the sample above the CO transition temperature can return it to the state with lower magnetization value.

Application of magnetic field of 8.5 T reduces the resistivity of the (La1-yPry)0.67Ca0.33MnO3 $`y=0.5`$ 16O and 18O samples (Fig. 7, insets). The resistivity of both samples in $`\mu _0H=8.5`$ T has a metallic temperature dependence (for $`T<300`$ K) with values essentially the same for 16O and 18O and close to the resistivity of the metallic La0.7Ca0.3MnO3.

In $`\mu _0H=8.5`$ T the specific heat of both samples drops significantly (Fig. 7, open circles), and tends to retain this low value when the magnetic field is reduced to 4 T (Fig.7a, crosses) and 0 T (triangles). This irreversible behavior is clearly demonstrated in the $`C(H)`$ plot (Fig. 8). This is similar to the behavior of the magnetization - only by heating the sample above the CO transition temperature does the specific heat return to its higher value. We reported similar magnetic and thermal history dependencies associated with “melting” of the CO in magnetic field for Pr0.65Ca0.35MnO3 . As in Pr0.65Ca0.35MnO3, the specific heat at $`\mu _0H=8.5`$ T (magnetic field sufficient to induce a transition to the FM state (Fig. 6)) is significantly larger than $`C(T)`$ of metallic La0.7Ca0.3MnO3 (dotted line in Figs. 7a and b) and exhibits the characteristic upward curvature of the CO state in the $`C/T`$ vs. $`T^2`$ plot. The specific heat in 8.5 T is similar to the specific heat of samples with a smaller fraction of the CO: for the 18O sample, $`C(T)`$ is close to the zero field $`C(T)`$ of 16O $`y=0.5`$ (dashed line in Fig. 7a), and for the 16O sample $`C(T)`$ is close to the zero field $`C(T)`$ of $`y=0.4`$ (dashed line in Fig. 7a). The results of fitting the 8.5 T data to Eq. 2 (see Table I) gives values close to the fitting parameters of the samples with the smaller fraction of the CO.

Unlike the magnetization, the $`C(H)`$ (when the field is applied to a ZFC sample) at $`T=5.5`$ K and 10 K shows no saturation in a magnetic field of 8 T (Fig.8). If we attempt to relate the value of the magnetization to the fraction of the metallic phase and the rest to the fraction of the CO phase, as we did above for different $`y`$, we can calculate the field dependence of the specific heat as $`C_{\mathrm{tp}}(H)=f_{\mathrm{met}}C_{\mathrm{met}}+(1f_{\mathrm{met}})C_{\mathrm{CO}}`$. This calculated specific heat of the two phase mixture $`C_{\mathrm{tp}}(H)`$ is not in agreement with the experimental data (Fig. 8): the ZFC calculated $`C_{\mathrm{tp}}(H)`$ (dotted lines in Fig. 8) decreases rapidly in the field interval from 4 T to 6 T and saturates at $`\mu _0H=7`$ T, while the ZFC experimental $`C(H)`$ data have approximately the same slope up to $`\mu _0H=8`$ T without a sign of saturation. As we noted before, the magnetization does not recover its full FM value even in magnetic field of 9 T, which probably means that some fraction of the AFM phase still remains at this field. These results indicate that the CO is not completely destroyed by the magnetic field sufficient to induce the FM transition. This phenomenon appears to be common for all CO manganites.

The magnetic field has a different effect on the samples with different oxygen isotopes. $`C(T)`$ is higher for the 18O sample due to larger fraction of the CO phase remaining in this magnetic field, since $`C(T)`$ (Fig. 7a) is close to to the zero field $`C(T)`$ of the 16O with $`y=0.5`$, while the $`C(T)`$ of the 16O sample in 8.5 T is close to $`C(T)`$ of the $`y=0.4`$ sample which has smaller fraction of the CO phase. This has qualitative agreement with the magnetization: larger magnetic field is needed to be applied to the 18O samples to induce the same effect on the charge ordering as for the 16O samples.

## IV DISCUSSION

Our study shows that the principal excitations of single phase systems such as metallic compositions of La1-xCaxMnO3 and insulating charge-ordered La0.5Ca0.5MnO3 are not affected by oxygen isotope exchange. A recent theoretical model , which suggests that the ground state of metallic manganites is a polaronic Fermi-liquid, predicts an increase of the effective mass (and hence $`\gamma `$) by about 3-6 % for the 18O material. Other authors point out that a large oxygen isotope effect on the effective mass of the charge carriers can only be expected for systems close to a local structural instability . Recent neutron scattering experiments have shown that the spin-wave dynamics and therefore the exchange interaction is the same for ferromagnetic La1-xCaxMnO3 with 16O and 18O . Also, Mossbauer effect measurements have shown that the exchange interaction does not change with oxygen isotope substitution . These results also show that the metallic ground state of these materials is not affected by the oxygen isotope exchange, which is in agreement with our specific heat data because of the following reasons. Within the double exchange model, the exchange interaction is determined by the an effective transfer integral $`t_{\mathrm{eff}}`$ for an electron hopping between Mn ions. In the tight-binding approximation, the density of states at the Fermi level, $`N(E_\mathrm{F})`$, is inversely proportional to $`t_{\mathrm{eff}}`$. Since the $`N(E_\mathrm{F})`$ is directly proportional to $`\gamma `$, materials with the same exchange interaction should have the same $`\gamma `$. In contrast, drastic changes in the low temperature properties are observed for the phase separated (La1-yPry)0.67Ca0.33MnO3, associated, as we have shown, with an increase of the CO phase fraction in the 18O materials.

A large isotope effect was observed for the ferromagnetic transition in metallic compositions of La1-xCaxMnO3, but only a negligible effect was found for La1-xSrxMnO3 . The charge ordering temperature changes significantly with isotope exchange for La0.5Ca0.5MnO3 and Nd0.5Sr0.5MnO3, but not for Pr0.5Ca0.5MnO3 . So what is in common for the large isotope effect materials? We believe it is the fact that the phase transition is first order and a phase separation is observed in the temperature region where the transition occurs. The FM transition in La1-xCaxMnO3 is believed to be the first order because the heat capacity peak associated with FM transition shifts significantly with application of magnetic field  and the spin dynamics observed in the vicinity of the FM transition is not consistent with a second order transition . Moreover, there is strong evidence for phase separation in the vicinity of the Curie temperatures for La1-xCaxMnO3 . In contrast, La1-xSrxMnO3 does not exhibit a phase separation near the FM transition and the phase transition in this material is second order . The charge ordering transition for La0.5Ca0.5MnO3 and Nd0.5Sr0.5MnO3 is a first order transition from a FM metal phase to a AFM CO phase accompanied by a hysteresis, showing the presence of both phases in the transition region. But Pr0.5Ca0.5MnO3 has a transition from a charge disordered to a charge ordered phase without signs of phase separation  and this transition is probably second order.

As was pointed out by some authors , the polaron effect alone  is too small to induce the observed changes associated with oxygen isotope exchange. Perhaps the polaron bandnarrowing effect produced by 18O induces localization effects in the phase separated region, which stabilizes the CO phase and leads to the increase of the volume fraction of the CO phase, as we observe for (La1-yPry)0.67Ca0.33MnO3. This leads to a decrease of the Curie temperature and an increase of the CO temperature in the phase separated system. In this picture, the largest effects should be observed near the percolation threshold, which is consistent with our observation that (La1-yPry)0.67Ca0.33MnO3 with $`y=0.5`$ exhibit larger changes with oxygen isotope exchange than $`y=0.4`$, and with the fact that for $`y=0.75`$ the oxygen isotope exchange induces a metal to insulator transition .

The experimental facts discussed above suggest: 1) the free energy of the CO phase is lower for 18O than for 16O, 2) the free energy of the metallic phase does not change much with oxygen isotope exchange, and 3) the large isotope effect on the transition temperature occurs when the transition is of first order (phase separation occurs). Fig. 9a shows a schematic free energy diagram in the vicinity of the FM - CO transition which occurs for La0.5Ca0.5MnO3 and Nd0.5Sr0.5MnO3 compositions. Since the transition is first order, the free energy curve $`F(T)`$ of the CO state crosses the $`F(T)`$ of the FM state at $`T_{\mathrm{CO}}`$ and the CO state has a lower value of $`F(T)`$ at $`T<T_{\mathrm{CO}}`$. The $`F(T)`$ of the CO phase of 18O material is lower than the 16O material. Therefore, the $`F(T)`$ of the CO phase of 18O material intersects the $`F(T)`$ of the FM phase at higher temperature, $`T_{\mathrm{CO}}`$(18O) $`>`$ $`T_{\mathrm{CO}}`$(16O). Using this picture, the change in the CO temperature $`\mathrm{\Delta }T_{\mathrm{CO}}`$ = $`T_{\mathrm{CO}}`$(18O) - $`T_{\mathrm{CO}}`$(16O) can be approximated as $`\mathrm{\Delta }T_{\mathrm{CO}}`$ = $`\delta /\mathrm{\Delta }S`$, where $`\delta `$ is the difference between the free energy of the 16O and 18O CO phases, and $`\mathrm{\Delta }S`$ is the difference in entropy (d$`F`$/d$`T`$) of the CO and FM phases at the transition temperature.

The entropy change associated with this transition was found to be 5 J/mole K for La0.5Ca0.5MnO3 . To estimate $`\delta `$, we use the experimental fact that the magnetic field sufficient to destroy the CO state and induce a transition to the FM state is higher for 18O than for 16O (see, for example, Fig. 6). The gain in magnetic energy associated with the transition to the FM state is $`MH_{\mathrm{FM}}`$, where $`M`$ is the magnetization of the FM state. The difference between the $`MH`$ values for 18O and 16O materials should correspond to the difference in free energy of the CO state of the 18O and 16O materials. The transition to the FM state for La0.5Ca0.5MnO3 occurs at magnetic fields higher than 9 T , beyond our experimental field limit. Therefore, we approximate the difference in free energy of the CO state of the 18O and 16O materials from our $`M(H)`$ data for (La1-yPry)0.67Ca0.33MnO3 (Fig. 6). Although the transition to the FM state is very broad for both 18O and 16O samples, a saturation occurs for the 18O sample at a magnetic field 2 T higher than for the 16O sample. Thus, $`\delta MH_{\mathrm{FM}}`$(18O) - $`MH_{\mathrm{FM}}`$(16O) $`42`$ J/mole, and estimated increase in the CO transition temperature $`\mathrm{\Delta }T_{\mathrm{CO}}=\delta /\mathrm{\Delta }S8.6`$ K, which is the same order as the observed $`\mathrm{\Delta }T_{\mathrm{CO}}=11.7`$ K. This agreement suggests that the proposed free energy diagram (Fig. 9a) is valid.

In the case of the FM transition in La1-xCaxMnO3 ($`0.2<x<0.4`$) in the vicinity of the FM transition the free energy curve $`F(T)`$ of the insulating state crosses the $`F(T)`$ of the FM state at $`T_\mathrm{C}`$ and the FM state has a lower value of $`F(T)`$ at $`T<T_\mathrm{C}`$ (Fig.9b). The insulating state of La1-xCaxMnO3 ($`0.2<x<0.4`$) shows a signature of short range charge ordering  in the vicinity of the FM transition. We argue that the substitution of 18O for 16O increases the volume fraction of the CO phase in the vicinity of the FM transition, which lowers the free energy of the insulating state. This shifts the crossing of $`F(T)`$ of the insulating state with $`F(T)`$ of the FM state to lower temperature (Fig. 9b), in agreement with the lowering of $`T_\mathrm{C}`$ observed when 18O replaces 16O in these materials.

## V CONCLUSIONS

We have studied the influence of oxygen isotope exchange on the resistivity, magnetization, and specific heat of single phase and phase separated manganites. We have found that the principal excitations in single phase metallic and CO materials do not change with oxygen isotope substitution, while isotope exchange induces large changes in the low temperature properties of phase separated (La1-yPry)0.67Ca0.33MnO3. The variation of the specific heat with Pr content $`y`$ in (La1-yPry)0.67Ca0.33MnO3 samples is consistent with the larger volume fraction of the CO phase for higher $`y`$. Oxygen isotope exchange in these materials induced changes similar to the increase of the Pr content $`y`$ which suggests an increase of the volume fraction of the CO phase in 18O samples.

Similar to other CO materials (Pr1-yCayMnO3, La0.5Ca0.5MnO3) all studied composition of (La1-yPry)0.67Ca0.33MnO3 ($`0.4<y<1`$) have shown the presence of the anomalous $`C^{}(T)`$ contribution to the low temperature specific heat due to the presence of the CO phase. The excess specific heat is also present at higher temperatures ($`T<T_{\mathrm{CO}}`$) for the CO materials, or materials containing larger fraction of the CO phase. Specific heat measurements indicate that the charge ordering is not destroyed completely by the application of a magnetic field of 8.5 T in contrast to the resistivity which shows a compleat “melting” of the CO. The nature of the excitations responsible for the $`C^{}(T)`$ term in the specific heat of the CO materials, and the presence of this term in rather strong magnetic field, has yet to be explained.

We suggest that the giant isotope effect on the Curie temperature and the CO temperature in manganites occurs in materials where these transitions are first order and hence, where the CO phase is present in the vicinity of the phase transition. 18O lowers the energy of the CO phase, which leads to an increase in the volume fraction of the CO phase in the phase separated regions. This leads to a decrease of the Curie temperature and increase of the CO temperature. A detailed understanding of the mechanism lowering the energy of the CO phase in the 18O materials still needs to be developed.

Acknowledgment: We thank A. J. Millis for helpful discussions, S-W. Cheong for providing us with as-grown La0.625Ca0.375MnO3 sample and helpful discussion, Z. Li for the sample preparation. The work at Maryland is supported in part by the NSF-MRSEC, DMR #00-80008. The work in Sherbrooke is supported by the Canadian Institute of Advanced Research (CIAR), the Canadian Foundation for Innovation (CFI), the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Fondation FORCE of the Universite de Sherbrooke.


text:
# How to be master of your domains

## Abstract

Abstract: Prepare a two-species BEC in a perfectly phase-mixed state. By applying Rabi drives, one can tune the range of wavelengths of phase-separating excitations that are dynamically unstable. Yada yada yada, this determines the characteristic sizes of the eventual spin domains. The trapping potential is neglected because it makes life hard, but of course this is a terrible approximation, and the results are of only inspirational value. Since this is thus a rather trivial calculation, the very modest increase in kudos that might accrue from journal publication does not seem to outweigh the grief of having to change the title. This will therefore be an ArXiv exclusive.

Consider a two-species condensate held in a homogeneous trap, and subjected to a Rabi drive. Homogeneous stationary states (there are always at least two solutions) satisfy the time-independent Gross-Pitaevskii equation

$`\mu _1\psi _1`$ $`=`$ $`Ae^{i\mathrm{\Omega }t}e^{i\theta _L}\psi _2+\left(g_{11}\left|\psi _1\right|^2+g_{12}\left|\psi _2\right|^2\right)\psi _1`$ (1)
$`\mu _2\psi _2`$ $`=`$ $`Ae^{i\mathrm{\Omega }t}e^{i\theta _L}\psi _1+\left(g_{12}\left|\psi _1\right|^2+g_{22}\left|\psi _2\right|^2\right)\psi _2.`$ (2)
Here $`A`$ is the driving amplitude, which initially we will take to be real and positive; $`\theta _L`$ is the phase of the driving field; and $`\mathrm{\Omega }`$ is its frequency. (Setting $`A`$ positive may obviously be done without loss of generality, but it will be notationally convenient later to redefine it so as to allow it to be negative.)

Our goal is to determine the Bogoliubov spectrum for small perturbations around the homogeneous stationary states, and then to consider which Bogoliubov modes may be dynamically unstable. Such instabilities will indicate the earliest, linear stages of spontaneous phase separation into spin domains. By varying $`A`$, it will turn out to be possible to control the ranges of wavelengths that are unstable. If we assume that the nonlinear evolution of domain formation does not significantly alter the ‘decisions’ made in the linear epoch, about where each species predominates, then control of the instability wavelengths means control of domain sizes . (This is admittedly a yada theorem ; but it is probably not too far wrong.)

From the viewpoint of BEC theory, it is obviously desirable to incorporate realistic features like an inhomogeneous trapping potential, and accurate values for the independent components of the interaction co-efficients $`g_{ij}`$. To get a feel for the basic concept, however, let us consider the tractable case of no potential, as in (2), and also stipulate $`g_{11}=g_{22}g`$. It will then be a nice notation to define $`g_{12}g\mathrm{\Delta }`$. We then have the two solutions

$`\psi _1`$ $`=`$ $`\sqrt{\rho }e^{i(\theta _0+\theta _L/2)}e^{i(\mu +\mathrm{\Omega }/2)t}`$ (3)
$`\psi _2`$ $`=`$ $`\pm \sqrt{\rho }e^{i(\theta _0\theta _L/2)}e^{i(\mu \mathrm{\Omega }/2)t}`$ (4)
$`\mu `$ $`=`$ $`\rho (g+g_{12})\pm A{\displaystyle \frac{\mathrm{}}{2M\xi ^2}}\left[1+\mathrm{\Delta }/4+\alpha /2\right],`$ (5)
where in the last step we introduce another notation that includes the healing length: $`g_{12}g(1+\mathrm{\Delta })`$, $`A\pm \alpha g\rho `$, and $`2g\rho `$ $`\mathrm{}/(2M\xi ^2)`$ for $`M`$ the atomic mass. In these solutions the densities $`\rho `$ of each of the two condensates are equal and constant, but their phases may differ by the drive phase, or by the drive phase plus $`\pi `$. There are thus two stationary states to consider; and at least to a good first approximation, either of them may be prepared at will, with an appropriate initial Rabi pulse. (A caveat: it will ultimately be an important question how much dynamical instabilities may limit the accuracy of the initial state preparation.) So we have two Bogoliubov problems to solve, depending on the branch of the $`\pm `$ in (5). Since the positive parameter $`A`$ already appears in our solutions, however, we have made things more compact by absorbing the $`\pm `$ into $`\alpha `$, which can thus be positive or negative. So from now on we will discuss only one set of Bogoliubov modes, whose frequencies will depend on our arbitrary real control parameter $`\alpha `$.

If for $`j=1,2`$ we write

$`\psi _j\psi _j\left[1+e^{i(𝐤𝐱\omega t)}(r_{j𝐤}+s_{j𝐤})+e^{i(𝐤𝐱\omega _kt)}(r_{j𝐤}^{}s_{j𝐤}^{})\right]`$

and then linearize the time-dependent Gross-Pitaevskii equation in $`r_{j𝐤}`$ and $`s_{j𝐤}`$, we obtain the Bogoliubov equations

$$\frac{\omega _k}{2g\rho }\left(\begin{array}{c}r_{1𝐤}\\ r_{2𝐤}\\ s_{1𝐤}\\ s_{2𝐤}\end{array}\right)=\left(\begin{array}{cccc}0& 0& k^2\xi ^2\alpha /2& \alpha /2\\ 0& 0& \alpha /2& k^2\xi ^2\alpha /2\\ k^2\xi ^2\alpha /2+1& \alpha /2+1+\mathrm{\Delta }& 0& 0\\ \alpha /2+1+\mathrm{\Delta }& k^2\xi ^2\alpha /2+1& 0& 0\end{array}\right)\left(\begin{array}{c}r_{1𝐤}\\ r_{2𝐤}\\ s_{1𝐤}\\ s_{2𝐤}\end{array}\right).$$
(6)
This is trivially equivalent to

$$\frac{\omega _k}{2g\rho }\left(\begin{array}{c}r_{1𝐤}+r_{2𝐤}\\ s_{1𝐤}+s_{2𝐤}\\ r_{1𝐤}r_{2𝐤}\\ s_{1𝐤}s_{2𝐤}\end{array}\right)=\left(\begin{array}{cccc}0& k^2\xi ^2& 0& 0\\ k^2\xi ^2+2+\mathrm{\Delta }& 0& 0& 0\\ 0& 0& 0& k^2\xi ^2\alpha \\ 0& 0& k^2\xi ^2\alpha \mathrm{\Delta }& 0\end{array}\right)\left(\begin{array}{c}r_{1𝐤}+r_{2𝐤}\\ s_{1𝐤}+s_{2𝐤}\\ r_{1𝐤}r_{2𝐤}\\ s_{1𝐤}s_{2𝐤}\end{array}\right).$$
(7)
Hence we have a ‘sonic branch’ of excitations, in which the two species oscillate in phase, and in which the frequencies are always real (and in fact it can be shown, positive) as long as $`\mathrm{\Delta }`$ is not too negative. For all condensates prepared so far, $`\left|\mathrm{\Delta }\right|1`$, so these modes are not very interesting.

On the other hand, we have an ‘optical branch’ in which the two species oscillate out of phase, and here is the physics we want to see. The frequency of an optical mode of wave number $`𝐤`$ is, for $`k\left|𝐤\right|`$,

$$\omega _k=2g\rho \sqrt{\left(k^2\xi ^2\alpha \right)\left(k^2\xi ^2\alpha \mathrm{\Delta }\right)},$$
(8)
where the sign of any real $`\omega _k`$ does not concern us here. (A negative real frequency for a positive norm mode means only energetic instability, and we are interested in dynamical instabilities, for which the Bogoliubov norm is zero.) What matters to us is that $`\omega _k`$ will be imaginary, indicating a dynamical instability, if

$`\alpha +\mathrm{\Delta }<k^2\xi ^2<\alpha `$ , $`\mathrm{\Delta }<0`$ (9)
$`\alpha <k^2\xi ^2<\alpha +\mathrm{\Delta }`$ , $`\mathrm{\Delta }>0.`$ (10)
So with negative $`\mathrm{\Delta }`$, the two species intrinsically prefer to stay mixed; but by applying a positive $`\alpha `$ we can destabilize them toward phase separation. Positive $`\mathrm{\Delta }`$, on the other hand, means that the condensates would intrinsically prefer to separate, with domain sizes no smaller than $`\xi \mathrm{\Delta }^{1/2}`$ (assuming that the minimum domain size is given by the maximum dynamically unstable $`k`$). By applying negative $`\alpha `$ we can enlarge this minimum domain size; if we make it larger than the sample size, then we effectively frustrate the phase separation.

For either sign of $`\mathrm{\Delta }`$, we can turn on phase separation by applying positive $`\alpha `$. We can make the minimum domain size smaller; and for large enough $`\alpha `$, we can even introduce a maximum domain size. Hence we can truly be masters of our domains.