text:
# Separation quality of a geometric ratchet

## I Introduction

During the last decade, ratchets have been the subject of intense research efforts (for a recent comprehensive review see ). Ratchets are able to produce a directed current of particles although no net average force is applied. Besides the fundamental interest in such a somewhat counterintuitive physical phenomenon, their analysis is important both for the description of natural nonequilibrium transport processes (like “Brownian motors” in cells ) and for concrete technical applications as rectifiers and separation devices. The various types of ratchets considered so far include rocking, flashing and correlation ratchets, where a temporally periodic force, periodic switching of a potential and colored non-thermal noise, respectively, induce directed transport in an asymmetric potential . Apart from these, there are “geometric” ratchets, which do not necessarily require any time-dependent forcing but consist, instead, of a two-dimensional periodic array of asymmetric obstacles , see Fig. 1. Particles are driven by a constant external force through the array while they are undergoing diffusive motion. Because of the asymmetry of the obstacles, the particles’ average drift velocity acquires a component perpendicular to the direction of the external force, which constitutes the ratchet effect. Since this is easier to realize experimentally than the time-dependent ratchets, it has been proposed and already demonstrated  as a device for the separation of charged biomolecules which were subject to an external electric field and underwent diffusive motion in an array of micrometer-sized obstacles produced by a lithographic process.

In this article, we want to take up the latter example and analyze specifically the quality of the separation effect and its dependence on various parameters. A similar numerical analysis has been carried out in Ref. , where the ratchet effect was investigated for a smooth periodic potential. However, we will emphasize that optimizing the ratchet effect alone is not equivalent to optimization of the separation quality. Apart from that, we will discuss the criterion for assessment of the separation quality, point out several “trivial” possibilities for optimization and analyze the effective change in the diffusion tensor brought about by the presence of the obstacles. The discussion will center around a model of particles undergoing drift-diffusive motion on a discrete lattice. This model is analyzed numerically using both a Monte-Carlo scheme as well as the numerical solutions of a master equation.

The remainder of the article is organized as follows: First, we will review briefly the geometric ratchet used for the separation of particles, point out the distinction between microscopic and macroscopic drift velocities and diffusion tensors and establish a measure of separation quality. Then, we will set up our lattice model. Using the results of the numerical simulations, we will present the dependence of the macroscopic diffusion tensor and drift velocity on the parameter characterizing the force applied to the particles, as well as on the shape of the obstacles. Finally, we will optimize the separation quality for a situation of two particle species that are subject to different forces, for a restricted set of obstacle shapes.

## II Separation in a geometric ratchet

In a “geometric ratchet”, particles drift and diffuse in a periodic potential, where the potential inside each elementary cell is asymmetric, leading to a ratchet effect. In the present work, we will only consider the special case of a two-dimensional array of impenetrable obstacles. The drift velocity and diffusion constant of the particles are assumed to be the same everywhere outside the obstacles. While the diffusion constant is fixed, the drift velocity depends not only on the mobility but also on the force applied to the particles. If the latter derives from an electric field, it will be proportional to the charge of the particles, which is important for effecting a separation of different particle species, in the way it has been demonstrated in Ref. . In this work, differently charged biomolecules have been injected at a corner of a periodic array of obstacles. An electric field is applied to the setup so that the particles are subject to a force pointing downward. Due to the asymmetry of the obstacles, the average drift velocity has a horizontal component. Therefore, after the particles have traversed several rows of obstacles, the center of the particle distribution is deflected by a certain amount. When the magnitude of this deflection is sufficiently different for two particle species, they may be separated by collecting the particles arriving at a certain row below the injection point. This permits a continuous separation of particles, in contrast to electrophoresis. The quality of separation does not depend on the strength of the ratchet effect alone but rather on the difference in deflections for two given species.

### A Macroscopic drift and diffusion

If the center of the particle distribution moves at a “macroscopic” average drift velocity $`\overline{\stackrel{}{v}}`$, then the slope of the line it traces (starting from the point of injection) is given by the ratio $`\overline{v}_x/\overline{v}_y`$. The average deflection in the final row is obtained by multiplying this ratio with the height of the array. Obviously, if the average drift velocity vector for one of the species is just proportional to that of the other one, no separation can result, regardless of the strength of the ratchet effect (i.e., the magnitude of the slope) itself.

Apart from the average drift, the particle distribution will also undergo diffusion, with a “macroscopic diffusion tensor” $`\overline{D}`$ that will be different from the microscopic one, due to the presence of the obstacles (which may hinder the expansion of the particle cloud in some directions, for example). It is important to know about the diffusive spreading of the distribution, since it affects the separation quality. A large difference in deflections for the two species will be useless if it is bought at the price of a large width of the respective distributions, which will overlap so that they cannot be separated unambiguously.

The general functional dependence of $`\overline{\stackrel{}{v}}`$ and $`\overline{D}`$ can be obtained from dimensional analysis. The microscopic parameters entering are the microscopic drift velocity $`v`$ (pointing along the $`y`$-direction), the (isotropic) diffusion constant $`D`$, the height of the obstacles $`h`$, and a collection of parameters describing their shape $`S`$ (including the aspect ratio). The only possibility of forming a dimensionless parameter out of a combination of $`D,v`$ and $`h`$ is given by

$$\xi \frac{vh}{D}.$$
(1)
It is proportional to the microscopic drift velocity $`v`$, and therefore to the microscopic mobility $`\mu `$ multiplied by the force $`F`$. The parameter $`\xi `$ will be used to present the results of our numerical simulations of a lattice model and to compare them with the real physical parameters. The most general form of macroscopic drift velocity and diffusion tensor is given by:

$`\overline{\stackrel{}{v}}`$ $`=`$ $`v\stackrel{}{v}_0(\xi ,S)`$ (2)
$`\overline{D}`$ $`=`$ $`DD_0(\xi ,S).`$ (3)
Here, $`\stackrel{}{v}_0`$ and $`D_0`$ are dimensionless vector and tensor functions.

From these considerations, one can already conclude that two species differing only in their mobilities (but not in the forces acting on them) will not become separated, in spite of the ratchet effect. This is due to the Einstein relation $`D=\mu k_BT`$, which ensures that the ratio $`\mu /D`$ will be the same for both species. Since the forces are assumed to be the same, $`v/D`$ is also unchanged, such that the respective values of $`\xi `$ are equal. Therefore, the average drift velocity vector only gets scaled when one passes from one species to the other, so that the drift slope remains the same, as has been discussed above.

For small values of the external force (i.e., $`\xi `$), the macroscopic drift velocity depends linearly on the force and, therefore, on $`v`$ (as long as the linear component is not suppressed due to symmetry). Therefore, a “macroscopic mobility tensor” $`\overline{\mu }`$ relating $`\overline{\stackrel{}{v}}`$ to the force $`\stackrel{}{F}`$ may be defined, such that $`\overline{\stackrel{}{v}}=\left(\overline{\mu }/\mu \right)v\stackrel{}{n}`$, with $`\stackrel{}{n}`$ the direction of the microscopic drift. The numerical calculations (see remarks in section III) will confirm that it fulfills the Onsager symmetry relations: $`\overline{\mu }`$ is found to be symmetric. Likewise the macroscopic diffusion and mobility tensors are connected by the Einstein relation for small values of the external force: $`\overline{D}=\overline{\mu }k_BT`$, where the temperature $`T`$ is obtained from the ratio $`D/\mu `$ of the corresponding microscopic quantities. Physically, this derives from the fact that the equilibrium distribution in a setup with a wall at the bottom is given by the Boltzmann distribution, which carries over from the microscopic density to the coarse-grained density, whose evolution is governed by $`\overline{D}`$ and $`\overline{\mu }`$. Apart from the fact that one has to be still in the regime of linear variation of $`\overline{\stackrel{}{v}}`$ for the Einstein relation to make sense, it can only hold as long as the force is not so strong as to make the density fall off rapidly over the scale of a single obstacle, because then the coarse-graining procedure is no longer justified. This happens approximately at $`\xi 1`$, which is also the condition that has to be reached to see an appreciable separation effect.

### B Analytical estimate

The magnitude of the ratchet effect can be estimated analytically  for high enough external force (large $`\xi `$). Then, one can treat the motion in the direction of the external force as deterministic (neglecting diffusion), so that diffusive spreading takes place only in the perpendicular direction. In this simplified picture, the geometric ratchet becomes analogous to a time-dependent one-dimensional “flashing” ratchet, with the $`y`$-coordinate playing the role of time. For a typical obstacle (as it has been used, e.g., in the experiment of van Oudenaarden and Boxer ), the diffusing particle distribution (in the shape of half a Gaussian curve) will be split in two parts by the “top part” of the obstacle, see Fig. 2. Here and in the following, we assume that there is one connected obstacle. The left part proceeds downward further on, while the right part moves one cell to the right. Therefore, the slope $`\overline{v}_x/\overline{v}_y`$ is given essentially by the percentage of particles that have moved to the right in such an “elementary step”, i.e., by an integral over the respective part of the Gaussian distribution. For large $`\xi `$ it can be approximated by an exponential (which becomes a good approximation if the magnitude of the following exponent exceeds $`2`$):

$$P(\xi )\frac{\overline{v}_x}{\overline{v}_y}\frac{2b}{\sqrt{\pi }w}\sqrt{\frac{h^{}}{h}}\frac{1}{\sqrt{\xi }}\mathrm{exp}\left(\frac{w^2}{4hh^{}}\xi \right).$$
(4)
Here, the dependence on $`\xi `$ has been made explicit, while all the other dimensionless parameters are ratios of lengths determining the shape $`S`$ of the obstacle (see Fig. 2). Obviously, if the force becomes very large, the particles will only move down inside “channels”, since they do not have time to spread to the left or right. Then the slope $`\overline{v}_x/\overline{v}_y`$ becomes very small, as is expressed by this formula. On the other hand, for very small forces neglecting the possibility of diffusing backwards in $`y`$-direction (or more than one cell in $`x`$-direction) renders this estimate invalid. Qualitatively, however, it is correct that the slope tends to a constant in the limit of vanishing force $`\xi 0`$. Therefore, separation is ineffective both at very small forces (since the slopes of two species are the same) and at very large forces (since the slopes differ but are both very small). Since the prefactor multiplying $`\xi `$ in the exponent is of the order of one, the optimum separation quality for the obstacle shape discussed here will also be reached when $`\xi _1`$ and $`\xi _2`$ are of the order of one. This is confirmed by the numerical analysis below.

### C Quality of separation

In the long-time limit, the ensemble of particles of a given species (having started at the injection point) assumes the form of a two-dimensional Gaussian distribution, which drifts at a velocity $`\overline{\stackrel{}{v}}`$. Therefore, the distribution of particles along the $`x`$-coordinate of the final row is also a Gaussian, which is centered around some value $`x`$ and has a width $`\sigma `$. We will assume that separation is performed by collecting all the particles up to some point $`x_0`$ in one bin and the rest in another bin. Ideally, the two bins would only receive particles of a single species ($`1`$ or $`2`$). Due to the overlap of the two distributions, this is not possible and there is a certain percentage of particles that are attributed to the wrong bin. Qualitatively, the optimal choice of $`x_0`$ is one where this “error” is minimized. However, since there are two different types of errors (percentage of particles $`1`$ attributed to bin $`2`$ and vice versa), no unambiguous definition of the optimal choice of $`x_0`$ and the corresponding optimal separation quality exists. We suggest to take the separation of $`x_1`$ and $`x_2`$ and divide it by the maximum of the widths $`\sigma _1`$ and $`\sigma _2`$ to arrive at a measure of separation quality which is easily evaluated:

$$Q\left|x_1x_2\right|/\mathrm{max}(\sigma _1,\sigma _2).$$
(5)
For the parameters considered here it seems to be appropriate.

Another possible definition consists in replacing the maximum by the geometric mean of the two widths:

$$Q^{}\left|x_1x_2\right|/\sqrt{\sigma _1\sigma _2}.$$
(6)
Note, however, that there are situations when this may be a misleading measure, particularly when one of the widths $`\sigma _1`$ or $`\sigma _2`$ is much larger than the other. In these cases, $`Q`$ will probably be better suited.

Other measures have been used in the literature. For example, in Refs. , the authors essentially asked how large the relative difference in diffusion constants of two species should be in order to have a separation that exceeds the spread of one of the two distributions.

In any case, any reasonably defined optimal separation quality can only be a function of two dimensionless parameters: the measure $`Q`$ (or $`Q^{}`$) used here and the ratio $`\sigma _1/\sigma _2`$.

The goal is to optimize the separation quality for two given species of particles by varying the applied force and the shape of the obstacles (including their size and aspect ratio). The following parameters are naturally assumed to be fixed: the microscopic mobilities and diffusion constants of the species (connected by the Einstein relation for a given temperature), the ratio $`\lambda `$ of the forces (equal to the ratio of charges in an electric field) and the total height $`H`$ of the periodic array in $`y`$-direction. The latter will be dictated by practical considerations (if it could be made arbitrarily large, one would do this to get an ideal separation effect).

We want to argue when and why a nontrivial optimum is to be expected. It has been pointed out above that choosing the force too large or too small results in a decrease of separation quality. Likewise, one could make the ratchet effect itself arbitrarily strong, by choosing a slanted line of large slope $`dx/dy`$ (as one big obstacle). However, both species would drift along that line and could not be separated either.

In order to obtain a more quantitative understanding, we use the general functional forms given for $`\overline{\stackrel{}{v}}`$ and $`\overline{D}`$, and insert them into the measure $`Q`$ introduced above. More precisely, we use the scaling expressions $`x=H\overline{v}_x/\overline{v}_y=Hx_0(\xi ,S)`$ and $`\sigma =\sqrt{DH/v}\sigma _0(\xi ,S)`$ for both species, which are assumed to have the same microscopic values of $`\mu `$ and $`D`$ but different forces acting on them, such that $`v_2=\lambda v_1\lambda v`$ and $`\xi _2=\lambda \xi _1\lambda \xi `$. Thus, we obtain:

$`Q`$ $`=`$ $`{\displaystyle \frac{\left|x_1x_2\right|}{\mathrm{max}(\sigma _1,\sigma _2)}}`$ (7)
$`=`$ $`\sqrt{{\displaystyle \frac{H}{h}}}{\displaystyle \frac{\sqrt{\xi }\left|x_0(\xi )x_0(\lambda \xi )\right|}{\mathrm{max}(\sigma _0(\lambda \xi )/\sqrt{\lambda },\sigma _0(\xi ))}}.`$ (8)
After the fraction in the last line has been optimized by varying both $`\xi `$ and the shape $`S`$ of the obstacles, $`Q`$ could be made arbitrarily large by having the height $`h`$ of an individual obstacle go to zero (at fixed aspect ratio), such that the array (of fixed height $`H`$) contains more rows. Since $`\xi =vh/D`$ must remain constant, this means that the microscopic drift velocity $`v`$ has to go to infinity. Physically, the enhancement of separation quality can be understood in the following way: Although the slopes will remain the same, the relative size of the diffusive spread decreases like $`\sqrt{h/H}`$. Note as well that, due to the same reason, the maximum of the fraction does not exist, strictly speaking. It is possible to keep $`h`$ fixed but still have an effective increase in the number of rows by placing an array of miniaturized obstacles inside an “elementary” cell. However, in practice there are obvious restrictions on the force that can be applied to the particles, as well as on the minimum size of the obstacles. Therefore, it is only possible to choose the force as large as possible and the corresponding value of $`h`$ in such a way that $`\xi `$ takes on the optimal value under these restrictions.

It is also helpful to consider the behavior of $`Q`$ for the analytical estimate given in Eq. (4), which is valid for large values of $`\xi `$. After $`H/h`$ rows (and a time $`t=H/\overline{v}_y`$), the average deflection is $`\overline{v}_xt=\overline{v}_yPt`$ and the variance of the resulting Poisson distribution is $`2\overline{D}_{xx}t=bHP`$. This gives the relation $`\overline{D}_{xx}/\overline{v}_x=b/2`$, which is confirmed by the numerical results below (for obstacles of type “A” and “B” in Fig. 6). We obtain:

$$Q\sqrt{\frac{H}{b}}\frac{|P(\xi )P(\lambda \xi )|}{\mathrm{max}(\sqrt{P(\xi )},\sqrt{P(\lambda \xi )})}.$$
(9)
By setting $`P(\xi )\mathrm{exp}(\gamma \xi )`$, with some exponent $`\gamma `$, we can find the approximate behavior of $`Q`$: It drops like $`\mathrm{exp}(\xi \gamma /2)`$ at large $`\xi `$, regardless of $`\lambda (>1)`$. On the other hand, the alternative definition $`Q^{}`$ given above would only decrease like $`\mathrm{exp}(\xi \gamma (3\lambda )/4)`$. For $`\lambda >3`$, this may even rise at large $`\xi `$, which is probably an indication that, for this regime, the expression for $`Q^{}`$ is not any longer a good measure of separation quality: The difference in spread of the two species grows too fast.

## III The lattice model

We have set up a model to study numerically the diffusion of particles under the influence of an external force in a periodic array of obstacles: A single particle is positioned on a point of a two-dimensional lattice which consists of square fields. During each time-step, the particle either changes its position to one of the neighboring squares (with a certain probability) or it remains on its original square.

The probabilities to move to the right or to the left are both equal to $`\mathrm{\Gamma }`$. In the absence of obstacles, there is consequently no net flow of particles in the horizontal direction. The probability to move downward (i.e. in positive $`y`$-direction) is $`\mathrm{\Gamma }+\alpha `$, the probability to move upward is $`\mathrm{\Gamma }\alpha `$ (see inset of Fig. 1). The probability to remain on the original square is therefore equal to $`14\mathrm{\Gamma }`$. Because of the different probabilities to move up- and downward, a net flow results in the vertical direction.

The following relations hold for the microscopic drift velocity $`v`$ in y-direction and the microscopic diffusion constant $`D`$:

$`v`$ $`=`$ $`2\alpha `$ (10)
$`D`$ $`=`$ $`\mathrm{\Gamma }.`$ (11)
We assume lengths to be measured in lattice constants and time in units of the elementary time-step. For the interpretation of the results, only the value of the dimensionless parameter $`\xi `$ is needed:

$$\xi =vh/D=2\alpha h/\mathrm{\Gamma },$$
(12)
where $`h`$ is the height of an obstacle (elementary cell of the array) measured in lattice constants. From this expression, it is clear that increasing the spatial resolution ($`h`$) means decreasing $`\alpha /\mathrm{\Gamma }`$, such that this ratio vanishes in the physically relevant continuum limit. Actually, we have already assumed $`\alpha \mathrm{\Gamma }`$ in writing down Eq. (11). Otherwise, one would have to take into account that the diffusion resulting from the model stated above is anisotropic, with $`D_{xx}=\mathrm{\Gamma }`$ and $`D_{yy}=\mathrm{\Gamma }2\alpha ^2`$. In our Monte-Carlo simulation, we have usually chosen $`\mathrm{\Gamma }`$ to be of the order of $`0.1`$ and $`\alpha /\mathrm{\Gamma }`$ to be less than about $`0.1`$.

Obstacles are represented by “forbidden” squares. If the particle would have to move onto such a square, it does not move and remains on its original square. All obstacles are arranged periodically.

Particles start in the top left square of the array. They continue moving until they have reached the final row, where their $`x`$-coordinate is saved. Using the results of many runs, the average deflection $`x`$ and the standard deviation $`\sigma `$ of the resulting distribution can be calculated, see Fig. 1. Alternatively, a particle runs for a given number of time steps and its final coordinates are registered. In this way, the macroscopic drift velocity $`\overline{\stackrel{}{v}}`$ and diffusion tensor $`\overline{D}`$ can be obtained.

In order to test the Monte-Carlo simulation (which has been implemented in C++), we have verified that the distribution of deflections in the final row is Gaussian (when coarse-grained over several obstacles) and becomes independent of the precise starting position after a sufficiently large time (number of rows of the array), and all quantities show the correct scaling behavior described by Eq. (3), provided $`\alpha /\mathrm{\Gamma }`$ is chosen small enough. The lattice resolution (connected with $`\alpha /\mathrm{\Gamma }`$) has been chosen such that the results do not depend on it appreciably any more.

Furthermore, at small values of the external force (i.e., of $`v`$ and $`\xi `$), the macroscopic drift velocity depends linearly on this force and fulfills the Onsager symmetry relations. Likewise the Einstein relation holds: $`\overline{D}=\overline{\mu }k_BT`$.

## IV Master equation

While every quantity of interest (macroscopic drift velocity and diffusion tensor, average deflection and spread of distribution) can be calculated using the Monte-Carlo simulation, it is nevertheless useful to consider a master equation solution as well. This is both because reaching a high statistical accuracy requires a large number of Monte-Carlo runs and because a discussion of the master equation yields additional physical insights.

The particle distribution $`p(x,y)`$ is defined on the lattice, with integer coordinates $`x`$ and $`y`$. In a single time-step, the distribution changes by

$`\delta p(x,y)`$ $`=`$ $`\mathrm{\Gamma }(p(x+1,y)+p(x1,y))`$ (15)
$`+(\mathrm{\Gamma }\alpha )p(x,y+1)+(\mathrm{\Gamma }+\alpha )p(x,y1)`$
$`4\mathrm{\Gamma }p(x,y),`$
where the quantities on the right-hand side are to be evaluated at time $`t`$ and $`\delta pp(t+1)p(t)`$. The equation shown here holds for every site $`(x,y)`$ which has no neighboring obstacle sites. For each “forbidden” site, the corresponding incoming and outgoing rates have to be left out. The temporal continuum limit is obtained by letting $`\mathrm{\Gamma }`$ and $`\alpha `$ tend to zero, with their ratio kept fixed.

At large times, an ensemble of particles that has started at the injection point will be spread over many obstacle cells. Viewed on the scale of only a few cells (much less than the total spread), the distribution $`p`$ is periodic. Therefore, we can calculate the average drift velocity by solving for the stationary distribution $`p`$ defined inside the cell of a single obstacle, imposing periodic boundary conditions. $`\overline{\stackrel{}{v}}`$ is then given by the average velocity inside the cell, i.e., by

$$\overline{\stackrel{}{v}}=\underset{\text{links}}{}\stackrel{}{j}.$$
(16)
Here, the sum runs over all links connecting adjacent sites and the “current density” $`\stackrel{}{j}`$ along each link is obtained by multiplying the values of $`p(x,y)`$ on the two connected sites by the transition probabilities of a jump along the link, taking into account the direction of the link and leaving out blocked sites. The distribution $`p`$ itself is assumed to be normalized: $`p(x,y)=1`$.

Actually, every site with no neighboring obstacle sites contributes just $`p(x,y)\stackrel{}{v}`$ to this sum, where $`\stackrel{}{v}`$ is the microscopic drift velocity. In this sense, the deviation of the macroscopic drift $`\overline{\stackrel{}{v}}`$ from $`\stackrel{}{v}`$ is seen to arise only from the boundaries of the obstacle. This can be understood most easily in the continuum limit, where $`\stackrel{}{j}=D\stackrel{}{}\rho +\stackrel{}{v}\rho `$ ($`\rho `$ being the continuous distribution). An integration by parts yields:

$$\overline{\stackrel{}{v}}=_{cell}\stackrel{}{j}d^2\stackrel{}{r}=\stackrel{}{v}+D_{obst.}\rho \stackrel{}{n}𝑑s.$$
(17)
Here, the second integral extends over the boundary of the obstacle and $`\widehat{n}`$ is the normal vector of that boundary (pointing away from the obstacle). Therefore, those boundaries where the density $`\rho `$ “piles up” (see Fig. 3 below) contribute most to the change in drift velocity. The second integral can therefore be interpreted as being proportional to the “force” due to the obstacle, impeding and deflecting the free flow of particles.

For the numerical solution, we first set up the (sparse) matrix corresponding to the transition rates on the lattice with periodic boundary conditions and then solve for the stationary solution $`p`$. This is done by setting $`p(x,y)=1`$ on an arbitrary non-blocked site and striking out the respective column and row in the homogeneous linear system of equations, such that it becomes inhomogeneous and nonsingular. We have used a sparse matrix bi-conjugate-gradient solver from the LAPACK package of linear algebra routines . For a typical cell size of $`128\times 128`$, the number of nonvanishing matrix entries is about $`10^5`$. In the end, $`p`$ is normalized and the current density and macroscopic drift velocity are calculated. This type of calculation has already been performed by the authors of Ref. , for a model with a smoothly varying periodic potential.

It is also possible to calculate the macroscopic diffusion tensor $`\overline{D}`$ using the master equation. This involves the solution of an inhomogeneous linear equation, with the linear operator defined by Eq. (15) and the inhomogeneity derived from the solution $`p`$ of the homogeneous equation. The derivation of the equation and the formula for $`\overline{D}`$ can be found in the Appendix. To the best of our knowledge, such an analysis has not been carried out before.

## V Numerical results

In this section, we present the results of both the Monte-Carlo simulations and solutions of the master equation.

In the Monte-Carlo simulation, the relative statistical accuracy of the macroscopic drift velocity ($`\delta \overline{v}_x/\overline{v}_x`$ etc.) is given approximately by $`\sqrt{D/\left(v^2Nt\right)}`$, while that of $`\overline{D}`$ is estimated to be $`1/\sqrt{N}`$, where $`N`$ is the number of samples and $`t`$ the number of time-steps. Independence of the detailed initial conditions is reached when the diffusive spread $`\sqrt{Dt}`$ becomes much larger than the width of an obstacle cell (typically $`160`$). We have chosen values of $`t10^7`$ and $`N10^3`$ in order to fulfill these criteria.

### A Density distribution and flow field

The particle distribution $`p`$ and corresponding current density $`\stackrel{}{j}`$ resulting from the solution of the master equation for a typical obstacle shape are depicted in Fig. 3, both for a low and a high value of $`\xi `$. This corresponds to the optimal situation for the separation of two species, where one of them is almost not deflected at all (high force, i.e., high $`\xi `$), while the other one has an appreciable probability to go to a neighboring cell to the right, due to diffusion around the top part of the obstacle.

### B Macroscopic drift velocity and diffusion tensor

The results for the components of the macroscopic drift velocity are displayed in Fig. 4. These data have been obtained for the obstacle of Fig. 3, using both the master equation with periodic boundary conditions and the Monte-Carlo simulation. As has been explained in a previous section, the onset of a nonlinear dependence of $`\overline{v}_y`$ and $`\overline{v}_x`$ on $`\xi `$ at rather low $`\xi `$ implies that the obstacle is well-suited for obtaining a separation effect (rather than merely a ratchet effect).

In Fig. 5, the components of the (symmetrized) macroscopic diffusion tensor $`\overline{D}`$ have been plotted vs. $`\xi `$, for the same obstacle. At low values of the driving force, all components of $`\overline{D}`$ are generally reduced compared with the microscopic diffusion constant $`D`$, since the obstacles hinder the free diffusion of particles. At higher forces, the $`y`$-component is increased, while the $`x`$-component $`\overline{D}_{xx}`$ is further reduced. This can be understood in the following way: The motion proceeds inside vertical “channels”, such that the particle cannot move easily in horizontal ($`x`$) direction, while the diffusion in $`y`$-direction is more or less free. Comparing the data for $`\overline{D}_{xx}`$ and $`\overline{v}_x`$ shows that $`\overline{D}_{xx}/\overline{v}_xb/2`$ at larger values of $`\xi `$ (to the right of the maximum of $`\overline{v}_x(\xi )`$), as expected: both quantities decrease exponentially. Note that the off-diagonal component $`\overline{D}_{xy}`$ of the macroscopic diffusion tensor changes sign at about $`\xi 8`$, which seems to approximately coincide with the sign-change in $`\overline{v}_x/\overline{v}_y`$ (see Fig. 4). We have not come up, however, with an explanation for the approximate correlation between $`\overline{D}_{xy}/D`$ and $`\overline{v}_x/\overline{v}_y`$ yet.

The statistical accuracy of the Monte-Carlo results for $`\overline{D}`$ is worse than that of the $`\overline{v}`$ results, as expected. The deviation between the results of the master equation and Monte-Carlo simulation at larger values of $`\xi `$ is reduced when the grid resolution is enhanced (i.e., when $`\alpha /\mathrm{\Gamma }`$ gets smaller for fixed $`\xi `$). In the examples shown here, $`\alpha /\mathrm{\Gamma }`$ takes on a maximum value of about $`0.1`$.

### C Optimization of separation quality

For the obstacle discussed above, the slope $`\overline{v}_x/\overline{v}_y`$ (proportional to the average deflection $`x`$ in the last row) is shown in Fig. 6, together with those of other obstacle shapes discussed further below. The exponential decay is consistent with the analytical estimate derived in a previous section, see Eq. (4). From the slope of the logarithmic plot, a value of about $`0.17`$ has been obtained for the prefactor $`w^2/(4hh^{})`$ in the exponent of Eq. (4), which is roughly consistent with the geometrical parameters of the obstacle. Given the slope and the spread $`\sigma `$ (derived from the components of $`\overline{D}`$), one can obtain the separation quality $`Q`$ defined in Eq. (6), if one assumes some ratio $`\lambda =\xi _2/\xi _1`$ of the forces acting on the two species (see Fig. 7).

### D Influence of the obstacle shape

Both the magnitude of the ratchet effect and the separation quality depend very much on the shape of the obstacle. We have not performed a systematic search over obstacle shapes for a kind of “global” optimization of the separation quality because of the numerical effort involved. However, there are a few general properties resulting from certain geometrical features. These are illustrated by the numerical results for the slope $`\overline{v}_x/\overline{v}_y`$ and the separation quality $`Q`$ plotted in Figs. 6 and 7. They can be summarized as follows:

The vertical “wall” at one side of the obstacle “A” depicted in Fig. 3 acts to prevent particles from diffusing back to the left, thereby increasing $`\overline{v}_x`$ and leading to a better ratchet and separation effect. This can be seen by comparing against a version with a shorter wall (“B”). The triangular “roof” of the obstacle splits the particle distribution in two halves as it drifts downwards. If the external force is high, the particles do not have time to diffuse a sufficient distance to the right and will be deflected back by the left side of this “roof”, therefore streaming downward, with no net deflection to the right. The horizontal position of the upper tip determines the strength of the force where this transition takes place: If it is moved to the left (“C”), a much higher force is necessary. In that case, the slope $`\overline{v}_x/\overline{v}_y`$ falls off more slowly with increasing force $`\xi `$. At the same time, the value of the slope is generally increased at small forces, since more particles are deflected one cell to the right.

In order to illustrate the difference between having a strong ratchet effect and a good separation effect, we have tried a triangular obstacle (“E”), which yields a comparatively large slope that, however, does not change very strongly with applied force. By flipping the triangle along the horizontal axis, an obstacle with reflection symmetry is created (“D”). This has the peculiar feature  that, for symmetry reasons, $`\overline{v}_x`$ is an even function of the microscopic velocity $`v_y`$ (driving force $`\xi `$), so the linear mobility at low driving force vanishes. In principle, this nonlinear dependence of $`\overline{v}_x`$ on $`\xi `$ is well-suited for achieving a separation effect. However, it must be kept in mind that at low values of the external force the diffusive motion is dominant, so the spread and the overlap of the particle distributions of the two species in the final row is significant. Therefore, the separation quality decreases towards $`\xi =0`$.

The separation quality can become zero for a special value of the external force whenever the slope shown in Fig. 6 has an extremum as a function of microscopic velocity $`\xi `$, such that two different $`\xi `$ can produce the same slope. This occurs for two of the obstacles (D and E in Fig. 7).

For most of the obstacle shapes considered here, there is at most a local maximum of the separation quality at low forces. The global maximum is expected to occur at much higher values of the force, which may be unattainable in the experiment (and are difficult to reach in a numerical simulation with a finite lattice resolution). However, for a well-suited obstacle like the one depicted in Fig. 3, the quality peaks at moderate values of the force.

## VI Conclusions

In this work, we have analyzed a geometric ratchet consisting of a two-dimensional array of obstacles, where particles perform drift-diffusive motion under the action of a constant external force. We have carried out numerical calculations using both a Monte-Carlo simulation and a master equation solution in order to obtain the dependence of the “macroscopic” drift velocity and diffusion tensor on external force and obstacle shape. Using these results, we have quantified the quality of the separation effect that can be achieved when two species of particles with differing microscopic mobilities are injected into the array. Our results show the strong dependence on several features of the shape of the obstacles and demonstrate the distinction between a strong ratchet effect and a good separation effect.

## VII Acknowledgments

We thank Hanno Gassmann for helpful discussions.

## VIII Appendix: Calculating the macroscopic diffusion tensor using the master equation

Solving the master equation for the probability density inside a cell containing a single obstacle allows one to obtain easily a numerically exact result for the macroscopic drift velocity $`\overline{\stackrel{}{v}}`$. It is given by the average flow velocity of particles inside the cell, i.e., the integral (or sum) over the current density (see Eqs. (16) and (17)). However, to assess the quality of separation of a given geometric ratchet, it is equally important to know the macroscopic diffusion tensor, which governs the spreading of the macroscopic particle density (i.e., the density averaged over many obstacles). Its evaluation using the Monte-Carlo simulation requires a large number of samples, to obtain a good statistical accuracy. Therefore, it is desirable to calculate $`\overline{D}`$ using the master equation as well. The steps involved in the derivation of $`\overline{D}`$ are slightly more involved than the straightforward calculation of $`\overline{\stackrel{}{v}}`$.

Our strategy is to derive the equation of motion for the macroscopic density $`\overline{\rho }`$,

$$_t\overline{\rho }=\overline{D}_{ij}_i_j\overline{\rho }\overline{\stackrel{}{v}}\stackrel{}{}\overline{\rho },$$
(18)
from the analogous equation for the microscopic density $`\rho `$, making use of the slow variation of $`\overline{\rho }`$.

If the particle density $`\rho `$ is spread over many obstacle cells (as is the case after waiting for a sufficiently long time), it may, in a first approximation, be described by

$$\rho =\rho _0\overline{\rho },$$
(19)
where $`\rho _0`$ refers to the detailed density that varies on the scale of a single obstacle but is periodic throughout the array. $`\rho _0`$ has been obtained before by solving the master equation for a single cell, using periodic boundary conditions at the borders of the cell and the restriction for the current density $`\stackrel{}{j}_0`$ to run parallel to the walls of the obstacle. If $`\overline{\rho }`$ were constant, this would constitute a (not normalizable) stationary periodic solution to the Fokker-Planck equation. However, $`\overline{\rho }`$ is assumed to vary (very slowly), such that this is not a stationary solution and does not fulfill the boundary conditions exactly. Therefore, it has to be supplemented by “correction terms”, which depend (necessarily linearly) on the spatial derivatives of $`\overline{\rho }`$. Consequently, we take $`\rho `$ to be given by:

$$\rho =\rho _0\overline{\rho }+g^i_i\overline{\rho }+K^{ij}_i_j\overline{\rho }+\mathrm{},$$
(20)
Here, $`g^i`$ and $`K^{ij}`$ are as yet unknown periodic functions (that vary on the scale of a single obstacle). We emphasize that, of course, many different microscopic densities $`\rho `$ yield the same macroscopic density $`\overline{\rho }`$. Therefore, $`\rho `$ is not uniquely specified if $`\overline{\rho }`$ is given. However, in the long-time limit assumed here, $`\rho `$ has “equilibrated” and the deviations from $`\rho _0`$ are in a one-to-one correspondence with the lowest order spatial derivatives of $`\overline{\rho }`$.

Our further strategy is as follows: We will first rederive the known result for the macroscopic drift velocity $`\overline{\stackrel{}{v}}`$, which is the constant coefficient appearing in the part of the equation of motion for $`\overline{\rho }`$ that contains the spatial derivatives of first order: $`_t\overline{\rho }=\overline{\stackrel{}{v}}\stackrel{}{}\overline{\rho }`$. To this end, we will insert the Ansatz (20) into the equation of motion for $`\rho `$, keeping only the terms containing up to first derivatives of $`\overline{\rho }`$ and eliminating $`g^i`$ using the boundary condition. Using the result for $`\overline{\stackrel{}{v}}`$, we will obtain an inhomogeneous linear equation for $`g^i`$ itself, which must be solved numerically. The relation between $`g^i`$ and the macroscopic diffusion tensor $`\overline{D}`$ will be obtained by going to the second order in the spatial derivatives of $`\overline{\rho }`$. Although this involves the unknown function $`K^{ij}`$, we will be able to eliminate it in the same way that $`g^i`$ had been eliminated in the first step.

Let us first derive a boundary condition for $`g^i`$. The current densities for the first and second term on the right hand side of Eq. (20) are given by

$`\stackrel{}{j}`$ $`=`$ $`\stackrel{}{j}_0\overline{\rho }D\rho _0\left(\stackrel{}{}\overline{\rho }\right)`$ (21)
$`(\stackrel{}{j}_1)_l`$ $`=`$ $`D(_lg^i)(_i\overline{\rho })+\stackrel{}{v}_lg^i(_i\overline{\rho })Dg^i_l_i\overline{\rho },`$ (22)
Here, $`\stackrel{}{j}_0=(D\stackrel{}{}+\stackrel{}{v})\rho _0`$ is the current density of $`\rho _0`$ alone. We demand $`\stackrel{}{j}+\stackrel{}{j}_1`$ to be parallel to the obstacle wall, keeping only terms including first order derivatives of $`\overline{\rho }`$ and then canceling these terms. This leads to the following boundary condition for $`g^i`$:

$$\widehat{n}\stackrel{}{j}_{g^i}=D\widehat{n}_i\rho _0,$$
(24)
where $`\stackrel{}{j}_{g^i}=(D\stackrel{}{}+\stackrel{}{v})g^i`$ is the current density related to $`g^i`$ and $`\widehat{n}`$ is the outer normal vector of the obstacle wall.

In the next step, we calculate $`_t\rho =\stackrel{}{}(\stackrel{}{j}+\stackrel{}{j}_1+\mathrm{})`$ up to first order in the spatial derivatives of $`\overline{\rho }`$ and demand it to equal $`\rho _0_t\overline{\rho }\rho _0\overline{\stackrel{}{v}}\stackrel{}{}\overline{\rho }`$, which is essentially the drift term for $`\overline{\rho }`$. Physically, this equation means that a nonvanishing slope of $`\overline{\rho }`$ will lead to an overall increase (or decrease) of the microscopic density $`\rho `$ inside an obstacle cell. The detailed shape of the distribution within that cell is not changed, only its magnitude. After dropping the overall factor $`_i\overline{\rho }`$, we arrive at

$$\stackrel{}{}\stackrel{}{j}_{g^i}=2(\stackrel{}{j}_0)_i+(\stackrel{}{v}_i+\overline{\stackrel{}{v}}_i)\rho _0.$$
(25)
It is not necessary to know $`g^i`$ in order to obtain $`\overline{\stackrel{}{v}}`$. We integrate both sides of this equation over the whole cell, assuming periodic boundary conditions for $`g^i`$ (as well as for $`\rho _0`$). The boundary term resulting from the walls of the obstacle contains $`g^i`$, but it can be transformed using Eq. (24), such that we end up with an equation where $`g^i`$ has been eliminated:

$$D\widehat{n}_i\rho _0𝑑s=2(\stackrel{}{j}_0)_id^2\stackrel{}{r}+(\stackrel{}{v}_i+\overline{\stackrel{}{v}}_i).$$
(26)
The integral on the left-hand side runs along the obstacle wall, while that on the right hand side extends over the whole cell. Since

$$(\stackrel{}{j}_0)_id^2\stackrel{}{r}=\stackrel{}{v}_i+D\widehat{n}_i\rho _0𝑑s,$$
(27)
we have

$$\overline{\stackrel{}{v}}_i=(\stackrel{}{j}_0)_id^2\stackrel{}{r},$$
(28)
as had been assumed in the main text already (see Eqs. (16) and (17)).

Inserting $`\overline{\stackrel{}{v}}`$ into Eq. (25) yields an inhomogeneous linear partial differential equation for $`g^i`$ which has to be solved numerically (assuming periodicity and the boundary condition Eq. (24)). Note that $`g^i`$ is only determined up to a constant multiple of $`\rho _0`$, since $`\rho _0`$ solves the homogeneous equation. However, as we will see, this does not affect the result for the diffusion tensor $`\overline{D}`$ to be derived from $`g^i`$. Further remarks concerning the numerical solution of the master equation on the discrete lattice can be found at the end of this appendix.

The current density related to the part of $`\rho `$ that involves second derivatives of $`\overline{\rho }`$ (see Eq. (20)) is equal to

$$(\stackrel{}{j}_2)_l=(_i_j\overline{\rho })\left[(D\stackrel{}{}_l+\stackrel{}{v}_l)K^{ij}\right]+\mathrm{},$$
(29)
where we have neglected higher derivatives of $`\overline{\rho }`$. We arrive at a boundary condition for $`K^{ij}`$ at the walls of the obstacle in the same way as for $`g^i`$, by demanding $`\stackrel{}{j}_{tot}\stackrel{}{j}+\stackrel{}{j}_1+\stackrel{}{j}_2`$ to be parallel to the wall. This time, we keep only the terms including second derivatives of $`\overline{\rho }`$. This leads to

$$\left[\widehat{n}\left(D\stackrel{}{}+\stackrel{}{v}\right)\right]K^{ij}=D\widehat{n}_jg^i.$$
(30)
To this order, the time-derivative of $`\rho `$, $`_t\rho =\stackrel{}{}\stackrel{}{j}_{tot}`$, includes both the diffusion of $`\overline{\rho }`$ and the drift of the term $`g^i_i\overline{\rho }`$:

$$\stackrel{}{}\stackrel{}{j}_{tot}=\overline{\stackrel{}{v}}_lg^i_l_i\overline{\rho }\rho _0\overline{D}_{li}_l_i\overline{\rho }+\mathrm{},$$
(31)
keeping only the second order with respect to the spatial derivatives of $`\overline{\rho }`$ on both sides of the equation.

As before, we integrate this equation over the cell and use the boundary condition Eq. (30) at the obstacle walls to eliminate $`K^{ij}`$. The resulting expression for $`\overline{D}`$ then is given by:

$$\overline{D}_{ji}=D\left(\delta _{ji}\widehat{n}_jg^i𝑑s\right)+(\overline{\stackrel{}{v}}_j\stackrel{}{v}_j)g^id^2\stackrel{}{r}.$$
(32)
Note that adding $`\lambda \rho _0`$ to $`g^i`$ (with an arbitrary constant $`\lambda `$) does not affect the result for $`\overline{D}`$, due to Eqs. (27) and (28).

For the numerical solution, it is, in principle, possible to discretize the equations (24) and (25) for $`g^i`$ as well as the expression (32) for $`\overline{D}`$. However, this is guaranteed to coincide with the results of the Monte-Carlo simulation only in the continuum limit (where, e.g., $`\alpha /\mathrm{\Gamma }0`$). In order to have a better agreement even when one is not yet in the continuum limit, it is advisable to start directly from the discretized master equation and redo the steps of the derivation shown here for the discrete lattice.

The equation which has actually been solved numerically to arrive at $`g^i`$ is given by:

$$Lg^i=p\overline{\stackrel{}{v}}_i\mathrm{\Gamma }(p^{i+}p^i)\delta _{i2}\alpha (p^{i+}+p^i).$$
(33)
Here, $`Lg^i`$ corresponds to $`\stackrel{}{}\stackrel{}{j}_{g^i}`$. Therefore, $`L`$ is the matrix kernel which is also used for solving the homogeneous equation, $`Lp=0`$, including the same treatment of obstacle walls and periodic boundary conditions (see the right-hand side of Eq. (15) for the definition of $`L`$). $`p`$ is evaluated at the “current site” (the site which the left hand side refers to), while $`p^{i+}`$ and $`p^i`$ are evaluated at the neighboring sites, in positive or negative direction $`i`$ (=$`1,2`$, corresponding to $`x,y`$), respectively. At obstacle walls, these neighboring sites may turn out to be “forbidden”, in which case $`p^{i+}`$ or $`p^i`$ vanishes. This implements the discrete version of the boundary condition discussed above. (Note that according to the convention used here, the microscopic drift velocity is assumed to point in negative $`y`$-direction if $`\alpha `$ is positive).

In order to evaluate $`\overline{D}`$, we must carry out a sum over all sites at the wall of the obstacle, i.e., those which have forbidden sites as neighbors. This sum is denoted by $`_W`$. The sum extending over all allowed sites in the cell is denoted by $``$:

$$\overline{D}_{ji}=\mathrm{\Gamma }\left(\delta _{ji}\underset{W}{}\widehat{n}_jg^i\right)+(\overline{\stackrel{}{v}}_j\stackrel{}{v}_j)g^i\alpha \delta _{j2}\underset{W}{}g^i.$$
(34)
The last term vanishes in the continuum limit but is important to ensure that $`\overline{D}`$ does not change on adding a homogeneous solution $`\lambda p`$ to $`g^i`$.


format: text

text:
# Comment on “Strain effect and the phase diagram of La1-xBaxMnO3 thin films”

## Abstract

Recent experiments in La1-xBaxMnO3 thin films by Zhang et al. \[Phys. Rev. B 64, 184404 (2001)\] showed that the ferromagnetic Curie temperature $`T_c`$ is enhanced (reduced) by tensile (compressive) strain. The results are regarded to be anomalous because it is generally understood by the authors that the tensile strain should lead to a reduction of the electron hopping $`t`$ (and thus also $`T_c`$) due to an elongation of the (in-plane) Mn-O bond length. In this comment such a general understanding of the strain effects is disapproved. It is argued that the tensile strain leads primarily to a larger Mn-O-Mn bond angle (with the Mn-O bond length nearly unchanged), thus the pure geometric consideration on $`t`$ is already in the correct direction towards the experimental findings by Zhang et al. At the same time, an indirect strain effect under consideration of the dynamic electron-phonon coupling is proposed, which may improve the quantitative explanation.

Very recently, the magnetic and transport properties for La1-xBaxMnO3 $`(x=0.050.33)`$ thin films grown on SrTiO3 substrate were reported by Zhang et al.  Depending on the doping level $`x`$ the films may suffer tensile $`(x0.2)`$ or compressive $`(x>0.2)`$ strain. It was found that the ferromagnetic transition temperature $`T_c`$ is enhanced (reduced) by the tensile (compressive) strain. The findings are regarded to be anomalous because a general understanding by the authors is “Tensile strain elongates the in-plane Mn-O bond length $`d`$, reducing $`t_0`$ (due to $`t_0d^{3.5}`$) and thus $`T_c`$; in contrast, compressive strain raises $`T_c`$” (see the second paragraph in Sec. III.E). In order to explain their unusual results, Zhang et al. proposed a mechanism involving the orbital degree of freedom, see Ref.  for details. The idea is that an elongation of the in-plane Mn-O distance due to the tensile strain helps increase the $`d_{x^2y^2}`$ component of the orbital state, thus in average the electron hopping would be enhanced because the hopping integral between $`d_{x^2y^2}`$ orbitals is larger than that for $`d_{3z^2r^2}`$ orbitals. As consistent results with theirs, Zhang et al. mentioned several other experimental findings: in Pr0.5Ca0.5MnO3/SrTiO3 thin films the charge ordering (CO) becomes less stable under tensile strain; and in La0.67Ca0.33MnO3/LaAlO3 thin films it is favored by compressive strain. Also, the same results as found by Zhang et al. in La1-xBaxMnO3 films were obtained by others.

Although the orbital-relevant mechanism proposed by Zhang et al. sounds good, we disagree to their general understanding of the strain effects, and further question the reality of the proposed mechanism. In this comment, we point out that the so called general understanding, which is extensively adopted in the literature,  is purely a misunderstanding. In fact, a simple correction to it, although not sufficient for a quantitative explanation, is already in the right direction toward the experimental findings by Zhang et al. At the same time, another mechanism in consideration of the electron-phonon (e-p) coupling is proposed to improve the quantitative explanation. In the following mainly the tensile strain is discussed and the case for the compressive strain is the reverse.

In thin films, the tensile strain leads to an elongation of the in-plane lattice parameter $`a`$. But it does not simply signify an elongation of the (in-plane) Mn-O bond length $`d`$. A possible change of the Mn-O-Mn bond angle $`\theta `$ was ignored by Zhang et al. and many others.  There is a simple geometric relation between $`a`$ and $`d`$: $`a=2d\mathrm{sin}(\theta /2)`$. Moreover, the hopping $`t`$ between two nearest neighbor Mn sites is $`\theta `$-dependent itself, which reaches the maximum when $`\theta =180^o`$. Totally, the direct dependence of $`t`$ on $`d`$ and $`\theta `$ is given by the following empirical formula:

$$td^{3.5}\mathrm{sin}(\theta /2).$$
(1)
From Eq. (1), if the angle $`\theta `$ is assumed independent of strain, the elongation of the lattice parameter $`a`$ (and the corresponding $`d`$) will lead to a reduction of $`t`$ as told by Zhang et al. Once it is noted that the angle $`\theta `$ may be modified, however, whether the hopping $`t`$ (and further $`T_c`$) is reduced or not by the tensile strain is not obvious even if the length $`d`$ is assumed elongated. Thus the above general understanding is not general at all.

Actually, there seems no direct experimental evidence that the bond length $`d`$ is elongated by the tensile strain in manganite thin films. On the contrary, recent measurements in La1-xCaxMnO3 $`(x=0.31,0.39)`$ thin films grown on SrTiO3 and LaAlO3 substrates,  with extended x-ray absorption fine structure, showed that the bond length $`d`$ is fixed at the value for the bulk materials, no matter whether the films are under tensile or compressive strain and the thickness of them. What is really modified is just the bond angle $`\theta `$. It becomes larger (smaller) under the tensile (compressive) strain, being consistent with the elongation (contraction) of $`a`$. A similar result was also observed in La0.67Sr0.33MnO3/LaAlO3 films.  In addition, it was found in bulk materials R1-xAxMnO3 (R=La, Pr, Nd, Y; A=Ca, Sr, Ba) with $`x=0.3,0.33`$, that changes in the average radius of the ions at (R, A) site (or equivalently the tolerance factor) only affect the mean Mn-O-Mn bond angle, with the mean Mn-O distance nearly unchanged. Based on all these findings, it is plausible to believe, for the current La1-xBaxMnO3 thin films, that the primary factor modified by strain is the angle $`\theta `$. Then a conclusion, which is exactly contrary to the so called general understanding, that is, the hopping $`t`$ is enhanced by the tensile strain due to an enlarged $`\theta `$, may be naturally obtained. In this way, the experimental findings by Zhang et al. are not so strange as claimed. Simultaneously, the orbital-relevant mechanism proposed by them, which is based on the hypothetical elongation of $`d`$, is not valid.

It is interesting to further make a quantitative calculation based on the above argument. When the bond length $`d`$ is invariant from strain, Eq. (1) will reduce to $`ta`$. If we assume the Curie temperature $`T_ct`$, we have $`T_ca`$. This means that the relative change of $`T_c`$, i.e., $`\mathrm{\Delta }T_c/T_c`$ is in the same order as that of the in-plane lattice parameter $`a`$, i.e., $`\mathrm{\Delta }a/a`$. In Ref. , for La1-xBaxMnO3 $`(x=0.05)`$ thin films which suffer tensile strain, the out-of-plane lattice parameter with different film thickness $`h`$ is given: $`3.872`$Å for $`h=85`$nm and $`3.862`$Å for $`h=20`$nm. The relative change is calculated by $`(3.8623.872)/3.8720.3\%`$. Then the corresponding relative change of $`a`$ should have roughly the same magnitude but opposite sign, i.e., $`\mathrm{\Delta }a/a+0.3\%`$. On the other hand, the Curie temperature $`T_c`$ increases from $`130`$K for $`h=85`$nm to $`180`$K for $`h=20`$nm, with a relative change $`\mathrm{\Delta }T_c/T_c=(180130)/130+38\%`$. It is clear that the change of $`\theta `$ (or $`a`$) is too small to account for the significant change of $`T_c`$, although this consideration has already given the correct tendency of the $`T_c`$ variation.

From the above paragraph, it has turned out that the pure geometric consideration on $`t`$ as shown by Eq. (1) is not enough for a quantitative explanation of the experimental results in Ref. , thus other cooperative mechanisms are still required. Due to the intrinsic complexity of the concrete films, there are various factors to influence $`T_c`$ such as nonstoichiometry, disorder, phase segregation, etc. Here we focus on discussing the $`T_c`$ variation caused by strain since the understanding of it is indispensable in thin films. Above we have addressed the direct effect due to strain on the electron hopping (and $`T_c`$). In the following we propose an indirect strain effect under consideration of the dynamic e-p coupling as a potential way to improve the quantitative explanation. Actually, since the e-p coupling is generally important in manganites, it is natural to ask the possible relevance of it to the change of $`T_c`$, as recently discussed by Chen et al. in La0.9Sr0.1MnO3 films. Explicitly, we consider that an electron at site $`i`$ is coupled to a local stretching mode which involves the stretching motions at least from the in-plane oxygen ions, for example, the mode $`B_{2g}(1)`$ with notations in Ref. . The coupling is of the form:

$$fc_i^{}c_iQ_i,$$
(2)
or in boson representation ($`\mathrm{}=1`$):

$$gc_i^{}c_i(b_i+b_i^{})$$
(3)
with $`g=f/\sqrt{2m\omega }`$. Here $`c_i,c_i^{}`$ are fermionic operators and $`b_i,b_i^{}`$ are bosonic operators, $`Q_i`$ represents the normal coordinate describing the stretching mode considered, $`f,g`$ are coupling constants, $`m`$ is the oxygen mass, and $`\omega `$ is the normal vibration frequency. In consideration of the e-p coupling (3), the bare hopping $`t`$ will be renormalized into the polaron hopping

$$\stackrel{~}{t}t\mathrm{exp}[g^2/\omega ^2]=t\mathrm{exp}[f^2/(2m\omega ^3)]$$
(4)
by an extra factor. The relative change of $`T_c`$ now results from that of $`\stackrel{~}{t}`$, which may become significant as shown below. To analyze the possible change of $`\omega `$ due to the modification of $`\theta `$, we introduce the central idea proposed recently by Egami and co-workers in a different context. If the Mn-O-Mn bond is straight, the vibration of the middle oxygen, when compressing one Mn-O bond, will inevitably stretch the other one. On the other hand, if the Mn-O-Mn bond is sufficiently bent, the contraction of one Mn-O bond can be accomodated by unbending the Mn-O-Mn triangle, without stretching the other Mn-O bond. Therefore, in the former case, the oxygen feels stronger restoring force when it vibrates, i.e., the former case corresponds to a larger effective elastic constant $`K`$ or vibration frequency $`\omega `$. In our picture, the tensile strain tends to straighten the in-plane Mn-O-Mn bond, and thus will lead to a larger $`\omega `$ based on the above analysis. At the same time, the coupling constant $`f`$ is considered to be unaffected in view of the fixed $`d`$. Then the renormalization factor will be raised. In particular, due to the exponential dependence shown in Eq. (4), a tiny hardening of $`\omega `$ will substantially increase the hopping $`\stackrel{~}{t}`$. Therefore, it is expected that the quantitative explanation addressed above can be largely improved. A full quantitative explanation will be left to the future when more experimental data are available, e.g., the change of $`\omega `$ with strain. As the first step, the experimental refinement to the local atomic configuration is needed to clarify the change of $`\theta `$ and/or $`d`$ with strain for the specific La1-xBaxMnO3 films discussed currently.

Finally, we point out that all the above ideas indicate the potential ways in which the electron itineracy can be modified by strain, thus they are of general significance for the understanding of the strain effects in manganite thin films. As mentioned above, the strain effects might not always be the dominant for experimental explanations in concrete films, nevertheless, they must be seriously studied. We mention that the (in)stability of the CO in the manganite films refered to in the first paragraph can be self-consistently explained by our proposal.

In summary, the so called general understanding of strain effects in manganite thin films is disapproved. It is argued that strain primarily leads to a change of the in-plane Mn-O-Mn bond angle, thus the pure geometric consideration on the electron hopping can already give the correct tendency of the $`T_c`$ variation found by Zhang et al. in La1-xBaxMnO3 thin films. Furthermore, in order to improve the quantitative explanation, an indirect strain effect under consideration of the dynamic electron-phonon coupling is proposed.

The author would like to acknowledge J. Xu, C. Schneider, T. Kopp, and R. Y. Gu for valuable discussions and J. Zhang for communications. This work was financially supported by the Deutsche Forschungsgemeinschaft through SFB 484, the National Natural Science Foundation of China (Grant No. 19904007) and the Texas Center for Superconductivity and Advanced Materials at the University of Houston.


text:
# Dimensional cross-over and charge order in half-doped Manganites and Cobaltites

## Abstract

We propose a generic model for understanding the effect of quenched disorder on charge ordering in half-doped Manganese- and Cobalt-Oxides with different crystal structures. Current experimental observations (Table 1) are explained in the light of the global phase diagram of the model (Fig. 2).

In the past decade much theoretical work was devoted to elucidating the nature and origin of charge ordering (CO) in the doped transition metal oxides, particularly in manganites, which are of great experimental and practical interest Millis-ManganiteOrder ; Daggoto01-ManganiteTheoryReview ; 01-ManganiteExperimentalReview . Much attention was focused on studying the interplay of charge correlations with orbital and spin degrees of freedom which are often important for the CO TokuraNagaosa2000 . Surprisingly, however, no coherent approach to understanding the effect of quenched disorder, such as introduced by a random distribution of dopant ions, on charge order (distinguished from orbital order) has been proposed so far.

Dopant disorder is a source of an unavoidable random electrostatic potential which couples linearly to charge density fluctuations 01-ManganiteExperimentalReview . In a strongly correlated electron system close to criticality such perturbation is of crucial importance. We propose that essential account for this random potential can be formulated in terms of a generic model (1), mapping the CO interaction on an effective Ising problem. We construct the phase diagram of this model, Fig. 2, and find that it adequately describes the essential features of the CO structure in half-doped manganese and cobalt oxides. In particular, a dimensional cross-over in our model gives a natural explanation to the puzzling experimental observation that CO in the layered perovskites is short-range Sternlieb1996 ; Zaliznyak2000 ; Wakabayashi2001\_LaSr1.5MnO , while it is long-range in the pseudo-cubic materials Radaelli1997\_LaCaMnO3 ; Zimmermann99\_PrCa0.5MnO ; Nakamura1999\_NdSrMnO3 .

At half-doping, the ground state of Co oxides and of many isostructural Mn oxides manifests a checkerboard-type planar charge modulation accompanied with a Jahn-Teller (JT) distortion of the MO6 (M = Mn, Co, etc.) octahedra shown in Fig. 1, Sternlieb1996 ; Zaliznyak2000 ; Wakabayashi2001\_LaSr1.5MnO ; Radaelli1997\_LaCaMnO3 ; Zimmermann99\_PrCa0.5MnO ; Nakamura1999\_NdSrMnO3 ; comment1 . The electrons are highly localized in a charge-ordered state with alternating Co2+ and Co3+ valence (Mn4+ and Mn3+ in manganites), and the in-plane O2- ions move towards the higher-valence M ions. To be explicit, we present the model construction in terms of a La1.5Sr0.5CoO4 Zaliznyak2000  which has no superimposed further symmetry breaking due to orbital order.

As a general framework, we argue that hopping dynamics, orbital, magnetic, local Coulomb and Jahn-Teller interactions can all be integrated out into an effective charge-ordering interaction which then competes with the random charge potential introduced by the dopant ions. A perfectly ordered state has a two-fold degeneracy associated with the choice of sublattice occupied by the higher/lower valence ions. We define a local Ising variable, $`\tau \left(\text{r}\right)=\pm 1`$, to denote the two alternative ionic configurations, Fig. 1(c). A single-domain CO state corresponds to a ferromagnetic order of $`\tau \left(\text{r}\right)`$. Topological disorder is introduced by domain walls as shown in Fig. 1(d). In principle, the model can be enriched by introducing the ‘vacancies’, $`\tau \left(\text{r}\right)=0`$, which represent cells with even number of electrons, but we shall exclude such possibilities in the analysis and discussion of this paper.

Because of the Jahn-Teller distortion, there is an elastic strain energy cost for nearest neighbor cells with opposite $`\tau `$ values, Fig. 1(d). This strain energy can be described by an effective short-range ferromagnetic interaction $`J\left(𝐫𝐫^{}\right)\tau \left(\text{r}\right)\tau \left(𝐫^{}\right)`$ JahnTeller-Review75 . Without a loss of generality we may approximate $`J\left(𝐫𝐫^{}\right)`$ by a nearest-neighbor coupling. The potential introduced by dopant ions randomly favors one or the other configuration of a two-ion cell, and therefore effectively acts as a random field $`h\left(\text{r}\right)\tau \left(\text{r}\right)`$. Importantly, depending on the crystal structure the effective coupling may be anisotropic, $`J_cJ_{ab}`$. Hence, for the charge order, the effect of elastic crystal strain and random dopant distribution is captured by an effective *anisotropic* random-field-Ising-model (RFIM),

$$H=J_{ab}\underset{\mathrm{𝐫𝐫}^{}_{ab}}{}\tau _𝐫\tau _𝐫^{}J_c\underset{\mathrm{𝐫𝐫}^{}_c}{}\tau _𝐫\tau _𝐫^{}+\underset{𝐫}{}h_𝐫\tau _𝐫,$$
(1)
where $`\mathrm{𝐫𝐫}^{}_{ab}`$ and $`\mathrm{𝐫𝐫}^{}_c`$ denote nearest-neighbor cells in the $`ab`$-plane and along the $`c`$-axis respectively, see Fig. 1 (e), (f). To our knowledge, all previous analytical and numerical studies focused on the isotropic $`3D`$ RFIM. Here we make the first attempt aimed at a qualitative and quantitative understanding of the anisotropic RFIM (1).

Random field alters both lower critical dimension Imry-Ma1975\_Disorder  and critical exponents RFI-NewExponents  of an Ising system. In particular, for weak random field, $`\mathrm{\Delta }/J1`$, the isotropic $`3D`$ RFIM orders at a finite temperature $`T_c>0`$, while the $`2D`$ RFIM has no long range order (LRO) even at $`T=0`$ ($`J`$ and $`\mathrm{\Delta }=\left\{\text{rms of }h_𝐫\right\}`$ parameterize the Ising coupling and disorder potential, respectively). This prompts an interesting question: what is the fate of a $`3D`$ anisotropic random field Ising model (ARFIM) of (1)? The simplest way of connecting the renormalization group (RG) flows leads us to the phase diagram shown in Fig. 2.

To connect our model analysis with the experiment, we summarize the charge order correlation lengths in the $`ab`$-plane, $`\xi _{ab}`$, and along the $`c`$-axis, $`\xi _c`$, measured in several representative perovskite oxides in Table 1. Correlations in the pseudo-cubic systems I-III (Fig. 1 (e)) are resolution-limited (typically this means $`\xi 2000`$ Å Radaelli1997\_LaCaMnO3 ; Zimmermann99\_PrCa0.5MnO ), and CO is apparently long-range. On the other hand, in the layered systems IV-V (Fig. 1 (f)) the CO is finite-range, with highly anisotropic $`3D`$ correlations $`\xi _{ab}\xi _c`$, which implies $`J_{ab}J_c`$. Large in-plane correlation length $`\xi _{ab}1`$ indicates weak planar disorder, $`\mathrm{\Delta }J_{ab}`$. Therefore, these materials probe the most interesting and complicated region of the phase diagram.

Because the isotropic $`3D`$ RFIM is long-range ordered for $`\mathrm{\Delta }<\mathrm{\Delta }_c^{3D}2.27J`$ Fisher01-RFIM3D-numerical , it is clear that for $`\mathrm{\Delta }<J_c<J_{ab}`$ the system is in the regime of weak $`3D`$ disorder and has a $`3D`$ LRO. Thus, we examine an intriguing case,

$$J_c\mathrm{\Delta }<J_{ab},$$
(2)
where the $`c`$-axis coupling is in the regime of *strong* disorder, while the $`ab`$-planes are in the *weak* disorder limit.

For $`J_c=0`$ the system decouples into independent $`2D`$ RFIM planes, and the ground state is disordered. This form of disordered ground state is perturbatively stable. For an infinitesimal inter-plane coupling $`J_c/\mathrm{\Delta }1`$, the $`ab`$-planes remain uncorrelated on all length scales ($`c`$-axis correlation length $`\xi _c<1`$), and each plane may still be treated as an effectively independent $`2D`$ RFIM. Therefore, a critical line of phase transitions exists in the $`(J_{ab}/\mathrm{\Delta },J_c/J_{ab})`$ phase diagram, separating the disordered and $`3D`$-LRO phases (solid line in Fig. 2). With increasing $`J_c/J_{ab}`$ on the disordered side we enter a critical regime, where the ground state is still disordered but with significant and highly anisotropic correlations in all spacial directions. In the isotropic $`2D`$ and $`3D`$ RFIM Villain-RFIM\_review  the ground state correlation length scales with $`\kappa =\mathrm{\Delta }/J`$ as $`\xi ^{2D}\mathrm{exp}\left(\frac{1}{\sigma }\kappa ^2\right)`$, and $`\xi ^{3D}\left(\kappa \kappa _c\right)^\nu `$ (for $`\kappa \kappa _c`$), respectively. What are the correlation lengths in a highly anisotropic $`3D`$ RFIM (1), with $`J_c\mathrm{\Delta }<J_{ab}`$?

The phase diagram, Fig. 2, instructs that upon a RG transformation model (1) scales either to a $`3D`$ LRO fixed point, or towards a $`3D`$ strong disorder fixed point. Currently, there is no generally satisfying analytical derivation of the RG equations for the RFIM for $`d<4`$. However, because $`d=2`$ is the lower critical dimension of the RFIM, 3D random field critical fixed point should be accessible perturbatively in $`2+\epsilon `$ expansion, Cardy-RGbook . Therefore, we derive perturbative RG equations for the anisotropic $`3D`$ RFIM (1) starting from a semi-phenomenological scaling ”ansatz” for the $`2D`$ RFIM Villain-RFIM\_review ; Cardy-RGbook . These equations, which naturally account for the RG flow of the coupling anisotropy, allow us to derive and analyze all qualitative features of the phase diagram of Fig. 2.

For the isotropic $`d`$-dimensional RFIM a simple scaling analysis of Imry-Ma, $`\left\{J(l)Jl^{d1};\mathrm{\Delta }\left(l\right)\sqrt{l^d}\mathrm{\Delta }\right\}`$, Imry-Ma1975\_Disorder ; Cardy-RGbook , predicts the relative disorder $`\kappa =\frac{\mathrm{\Delta }}{J}`$ to scale as

$$\kappa \left(l\right)=\frac{\mathrm{\Delta }\left(l\right)}{J\left(l\right)}\frac{\sqrt{l^d}}{l^{d1}}\kappa .$$
(3)
The marginality of $`\kappa \left(l\right)`$ in 2D is associated with $`d=2`$ being the lower critical dimension. It is broken by subdominant contributions from the interface roughening at the domain boundaries which lead to a disordered ground state Cardy-RGbook ; Villain-RFIM\_review . For weak disorder, $`\kappa 1`$, the correlation length is $`\xi =A\mathrm{exp}\left\{+\frac{1}{\sigma }\kappa ^2\right\}`$, with $`\sigma 1`$ (clearly, this expression fails for strong disorder, where $`\xi 1`$). By demanding that $`\kappa \left(l\right)`$ scales consistently with the correlation length, $`\xi \left(l\right)=\xi /l`$, and with the expression $`\xi \left(l\right)=A\mathrm{exp}\left\{+\frac{1}{\sigma }\kappa ^2\left(l\right)\right\}`$, we obtain a scaling equation

$$\kappa ^2\left(l\right)=\kappa ^2\left[1\kappa ^2\sigma \mathrm{ln}\left(l\right)\right],$$
(4)
i.e. scaling of the relative disorder in 2D is governed by the logarithmic corrections to Eq. (3). The disordered ground state is manifested by the flow of $`\kappa ^1\left(l\right)=\frac{J\left(l\right)}{\mathrm{\Delta }\left(l\right)}0`$. In the context of the model (1), equation (4) describes the RG flow for $`J_c=0`$.

To proceed with deriving the qualitative features of the RG flows shown in Fig. 2 for $`0<J_cJ_{ab}`$ we note that the anisotropy of interactions is, in fact, an irrelevant perturbation at the 3D critical point. This is reflected by the fact that all stable fixed points are located on the isotropic line in the phase diagram, at $`J_c/J_{ab}=1`$. At the mean-field level (and without random fields) the anisotropy of the interactions can be removed by appropriate anisotropic re-scaling of the coordinates Cardy-RGbook . Here we implement this idea by devising an *anisotropic* real-space RG transformation, effectively coarse-graining to anisotropic blocks $`l_{ab}\times l_{ab}\times l_c`$. We derive the RG equations along the lines of $`2+\epsilon `$ expansion, and demonstrate explicitly the flow towards the isotropic fixed points.

We parametrize the anisotropy as $`\alpha =\left(\frac{J_c}{J_{ab}}\right)^x`$, $`x>0`$, and define the RG transformation through an infinitesimal anisotropic re-scaling, $`\{d(\mathrm{ln}l_{ab})=d(\mathrm{ln}l);d(\mathrm{ln}l_c)=\alpha (l)d(\mathrm{ln}l)\}`$. This determines the shape of the coarse-grained blocks to be gradually varying, flowing from the anisotropic towards isotropic. We use the Imry-Ma approximation Cardy-RGbook  to derive the lowest order contributions to the flow equations for the model parameters, and supplement it by the appropriate higher order correction to reproduce the scaling equation (4) for the 2D RFIM case, $`\alpha =0`$. The resulting perturbative RG equations are

$`{\displaystyle \frac{\mathrm{d}\left(\kappa ^2\right)}{\mathrm{d}\mathrm{ln}\left(l\right)}}`$ $`=\alpha \kappa ^2+\sigma \kappa ^4+𝒪\left(\alpha ^2\right)`$ (5)
$`{\displaystyle \frac{\mathrm{d}\alpha }{\mathrm{d}\mathrm{ln}\left(l\right)}}`$ $`=x\left(1\alpha \right)\alpha `$ (6)
We obtain an unstable 2D fixed point at $`\{\alpha ^{}=0,\kappa ^2=0\}`$, and stable isotropic 3D fixed points at $`\alpha ^{}=1`$. The $`3D`$ critical point is at $`\{\alpha ^{}=1,\kappa ^2=\kappa _c^2=\alpha ^{}/\sigma \}`$. With $`\kappa _c2.27`$ from numerical simulations Fisher01-RFIM3D-numerical , the extrapolation of our perturbative equations to $`\alpha =1`$ would imply $`\sigma 0.2`$, which is too small. This indicates that while qualitatively correct, and valid in the perturbative region of parameters, RG equations (5),(6) are inadequate for quantitative estimates.

It is straightforward to integrate Eq. (6) and obtain $`\alpha (l)^1=1+(\alpha (1)^11)l^x`$, the scaling flow of the anisotropy. In principle, the phase transition line can be determined by a further numerical integration. The important asymptotes can be derived analytically,

$`\alpha (\kappa ^2)C_0e^{\frac{x}{\sigma \kappa ^2}},`$ $`\alpha 1,`$ (7)
$`\alpha (\kappa ^2)1+C_1(\kappa ^2\kappa _c^2),`$ $`1\alpha 1,`$ (8)
and show that neglecting $`𝒪\left(\alpha ^2\right)`$ term in Eq. (5) is indeed justified for $`\alpha 1`$, but not for $`\alpha 1`$.

To obtain quantitative estimates of the anisotropic correlation lengths as a function of $`J_c`$, $`J_{ab}`$ and $`\mathrm{\Delta }`$, we use an alternative approach. In the absence of the exact RG equations for the 3D ARFIM, we resort to a real-space rescaling procedure where we employ the Migdal-Kadanoff bond-moving technique Kadanoff76-BondMoving  combined with the scaling ansatz (4). We adopt conventional approximations in ignoring the generation of longer range interactions Villain-RFIM\_review  and using the “majority rule” Kadanoff76-BondMoving  for block spin variable. Hence, we follow only the transformation of effective nearest neighbor block interactions $`J_{ab}(l)`$, $`J_c(l)`$ and effective disorder parameter $`\mathrm{\Delta }(l)`$.

Our strategy is to map the anisotropic RFIM onto a ”solvable” isotropic $`3D`$ RFIM, whose known properties enable us to derive the estimates for $`\xi _{ab}(J_{ab},J_c,\mathrm{\Delta })`$ and $`\xi _c(J_{ab},J_c,\mathrm{\Delta })`$. Because the correlation length scales as $`\xi \left(l\right)=\xi /l`$, a *bare anisotropic* $`3D`$ model with finite $`\xi _{ab}>\xi _c1`$ is mapped onto an *effective isotropic* $`3D`$ model with $`\xi _{ab}=\xi _c`$ when the coarse-grained block size $`l_{ab}\times l_{ab}\times l_c`$ is chosen with $`l_{ab}=\xi _{ab}/\xi _c`$ and $`l_c=1`$. In other words, the cross-over from quasi-$`2D`$ to isotropic $`3D`$ scaling is intuitively implemented in two steps. First, a $`2D`$ transformation to blocks of size $`l_{ab}\times l_{ab}\times 1`$, and then a usual isotropic $`3D`$ scaling (see Fig. 3).

Performing the Migdal-Kadanoff transformation to $`l_{ab}\times l_{ab}\times 1`$ blocks we first move the $`J_c`$ bonds Kadanoff76-BondMoving  and thus create an effective block interaction $`J_c\left(l_{ab}\right)=J_cl_{ab}^2`$, Fig. 3B. Now that the c-axis coupling is removed, each planar $`l_{ab}\times l_{ab}`$ section can be integrated out in a $`2D`$ fashion using (4). This approximation is equivalent to assuming that fluctuations on length scale $`l_{ab}<\xi _{ab}/\xi _c`$ in $`ab`$-planes are uncorrelated in the c-axis direction. In addition, we keep the Imry-Ma approximation for scaling of the disorder $`\mathrm{\Delta }\left(l_{ab}\right)`$ under $`2D`$ block transformation Cardy-RGbook . Indeed, a perturbative expression for the surface tension per unit length on the scale $`l`$ in $`2D`$ RFIM Villain-RFIM\_review , $`\mathrm{\Sigma }(l)=J\left[1\sigma \kappa ^2\mathrm{ln}(l)\right]`$, indicates that logarithmic corrections to $`\kappa \left(l\right)`$ come from the effective weakening of the inter-block interaction $`J\left(l\right)`$ as domain boundaries adjust to the random field. We thus arrive at the final form of our approximate scaling equations for block $`l_{ab}\times l_{ab}\times 1`$,

$`J_c\left(l_{ab}\right)`$ $`=J_cl_{ab}^2`$ (9)
$`J_{ab}\left(l_{ab}\right)`$ $`=J_{ab}l_{ab}\sqrt{1\sigma \kappa ^2\mathrm{ln}\left(l_{ab}\right)}`$ (10)
$`\mathrm{\Delta }\left(l_{ab}\right)`$ $`=\mathrm{\Delta }l_{ab}.`$ (11)
For $`l_{ab}=\xi _{ab}/\xi _c`$, the correlation lengths in the transformed model are isotropic $`\xi _{ab}\left(l_{ab}=\xi _{ab}/\xi _c\right)=\xi _c`$. Since the disorder remains isotropic, it must be that the rescaled interactions are also isotropic. Thus, by imposing a self-consistency condition, $`J_{ab}\left(l_{ab}=\xi _{ab}/\xi _c\right)=J_c\left(\xi _{ab}/\xi _c\right)`$, we obtain a mapping of the bare *anisotropic* RFIM (1) onto *isotropic* $`3D`$ RFIM with effective interactions

$$J_{ab}\left(l_{ab}=\xi _{ab}/\xi _c\right)=J_c\left(\xi _{ab}/\xi _c\right)=\left(\frac{\xi _{ab}}{\xi _c}\right)^2J_c,$$
(12)
and $`\mathrm{\Delta }\left(\xi _{ab}/\xi _c\right)=\left(\frac{\xi _{ab}}{\xi _c}\right)\mathrm{\Delta }`$. Using (10),(11),(12) and expression for the correlation length of the $`3D`$ RFIM, $`\xi _{3D}=\mathrm{\Lambda }\left(\kappa _c^2\kappa ^2\right)^\nu =\xi _c`$, we obtain

$`\left({\displaystyle \frac{J_{ab}}{\mathrm{\Delta }}}\right)^2`$ $`=\kappa _c^2\left({\displaystyle \frac{\mathrm{\Lambda }}{\xi _c}}\right)^{1/\nu }+\sigma \mathrm{ln}\left({\displaystyle \frac{\xi _{ab}}{\xi _c}}\right)`$ (13)
$`\left({\displaystyle \frac{J_c}{\mathrm{\Delta }}}\right)^2`$ $`=\left({\displaystyle \frac{\xi _c}{\xi _{ab}}}\right)^2\left[\kappa _c^2\left({\displaystyle \frac{\mathrm{\Lambda }}{\xi _c}}\right)^{1/\nu }\right].`$ (14)
Substituting $`\xi _c=1`$, we obtain the crossover line $`\xi _c1`$ in the phase diagram,

$$\frac{J_c}{J_{ab}}=\frac{\mathrm{\Delta }}{J_{ab}}\sqrt{\kappa _c^2\mathrm{\Lambda }^{1/\nu }}e^{\frac{1}{\sigma }\left[\left(\frac{J_{ab}}{\mathrm{\Delta }}\right)^2+\mathrm{\Lambda }^{1/\nu }\kappa _c^2\right]}.$$
(15)
With $`\kappa _c2.27`$, $`\mathrm{\Lambda }0.1`$ and $`\nu 1.37`$ from numerical simulations, Fisher01-RFIM3D-numerical , equations (13),(14) allow us to retrodict the effective model parameters. For La0.5Sr1.5CoO4 (Table 1, V) we find $`\frac{\mathrm{\Delta }}{J}0.85`$, and $`\frac{J_{ab}}{J_c}>53`$, indeed satisfying the relations $`J_c\mathrm{\Delta }<J_{ab}`$ under which our approximations are valid. For La0.5Sr1.5MnO4 (Table 1, IV) $`\frac{\mathrm{\Delta }}{J_{ab}}0.64`$ and $`\frac{J_{ab}}{J_c}49`$. In both cases the anisotropy of the bare coupling parameters is remarkably high! Yet, we think it is quite realistic for the layered systems. Because neighbor $`ab`$-planes are shifted so that MO6 octahedra in one plane fit in-between those in the other (see Fig. 1(f)), and the inter-plane spacing is rather large, breathing-type distortions of the octahedra accompanying the CO are weakly coupled between the planes. On the other hand, octahedra within each plane form a corner-sharing network, so the coupling of in-plane distortions is strong. In the pseudocubic perovskite structure, the MO6 octahedra share apical oxygens along the $`c`$-axis as well (see Fig. 1(e)), and crystal fields are more isotropic. Note that, for a system to have a LRO groundstate, interactions need not be isotropic. In fact, the phase diagram instructs us that only very strongly anisotropic systems will be disordered if the disorder potential is weak (as exemplified by Eq. (7), and the above estimates of $`J_{ab}/J_c`$ in the layered compounds).

In conclusion, we argued on very general grounds that charge order observed in half-doped manganites and their isostructural relatives is described by the *anisotropic* ($`J_cJ_{ab}`$) $`3D`$ random field Ising model (1). We constructed the schematic phase diagram of this model, Fig. 2, and supported its intuitive structure by perturbative RG equations (5),(6). We also obtained the quantitative estimates of the effective model parameters as a function of the measured correlation lengths using a Migdal-Kadanoff block transformation scheme combined with a phenomenological scaling equation (4) for a 2D RFIM.

Long-range charge order is observed in the pseudo-cubic materials at low temperatures. In contrast, a disordered ground state with anisotropic correlations is found in the layered compounds. This disparity is naturally explained in our model: because of the strong anisotropy of the layered materials these systems reside in the different regions of the phase diagram. The anisotrpic correlation lengths measured in La0.5Sr1.5CoO4 and La0.5Sr1.5MnO4 corroborate our conclusion. We note that a disordered $`3D`$ RFIM has well known features which may be verified experimentally. Furthermore, our analysis may be applied to thin films, where even in pseudocubic compounds correlations in one direction are limited by the film thickness. The latter, together with the model parameters, will then determine the CO *finite* correlation length *within the film*. Finally, we note that charge disorder in the form of dislocations, as in Fig. 1(d), has inevitable effects on the superimposed orbital and spin orders which have thus far not been elaborated.

###### Acknowledgements.

We thank S. A. Kivelson, A. Tsvelik, F. Essler and J.P. Hill for fruitful discussions, and acknowledge the financial support of DOE#DE-AC02-98CH10886 (O.Z. and I.Z.) and DOE#DE-FG03-00ER45798 (O.Z.).


text:
# On the large escape of ionizing radiation from GEHRs

## 1 INTRODUCTION

The issue of whether or not ionizing star clusters are the sources of radiation responsible for the diffuse ionized gas (DIG) and, moreover, for the background radiation at low redshift (Heckman et al 2001), is still controversial. UV observations of four nearby starburst galaxies (Leitherer et al. 1995; Hurwitz et al. 1997) showed that more than 10% of the ionizing photons could escape from these galaxies. On the other hand, Dove et al. (2000) obtained from theoretical arguments that 7% of photons produced by OB associations would have to escape into the DIG to sustain it and could also be consistent with the estimated flux required to photoionize the Magellanic Stream (Weiner & Williams 1996; Bland-Hawthorn & Maloney 1999). All these aspects at low redshift are undoubtedly important to establish whether massive stars in starburst galaxies can rival QSO’s as sources of the metagalactic background radiation at high redshifts (z $`>`$ 3; see Ricotti & Shull 2000). Furthermore, an estimate of the escape fraction at low redshift is to provide an upper limit to the expected escape fraction at high redshift, when the universe was much denser. Beckman et al. (2000) presented evidence, based on the analysis of the H$`\alpha `$ luminosity function of HII region populations in nearby spiral galaxies, that luminous HII regions may be matter bounded. This would imply that these regions are an important source of photons to the DIG. Here we report strong evidence that points at Giant Extragalactic HII Regions (GEHRs) to be matter bounded, and thus as an important source of ionization at large distances from their exciting stars. Our conclusions are based on evolutionary synthesis models (Leitherer et al. 1999; hereafter ST99) applied to previously analyzed GEHRs with observed Wolf-Rayet (WR) features (Díaz et al. 2000; Castellanos et al. 2002; hereafter referred to as Paper I). The observed parameters inside these regions, i.e. the fluxes and equivalent widths of the WR features, the emission line intensities of \[OII\]($`\lambda `$3727), \[OIII\]($`\lambda `$5007), \[NII\]($`\lambda `$6584), \[SII\]($`\lambda `$6717), and \[SIII\]($`\lambda `$9069) relative to H$`\beta `$, the observed H$`\beta `$ luminosity and equivalent width, cannot be fitted simultaneously unless an important fraction of the ionizing radiation escapes from these regions.

## 2 MODELS AND RESULTS

The analyzed GEHRs are regions H13 in NGC 628, CDT3 in NGC 1232 and 74C in NGC 4258. The three of them, particularly 74C and H13, show spectral signatures of prominent Wolf-Rayet features. The presence of these stars allows to constrain the age of the ionizing clusters to between 3Myr and 6Myr. These values depend on the stellar mass loss rates and, consequently, on the stellar metallicity that ultimately controls when massive O stars enter the Wolf-Rayet phase (Meynet 1995). This metallicity has been assumed to be similar to that found through the analysis of the emission line spectrum (see Paper I). Once the age and metallicity are known, a unique Spectral Energy Distribution (SED) can be provided by evolutionary synthesis models.

To estimate the age of the ionizing stellar clusters from the measured strength and equivalent widths of the observed WR features we used the models of Schaerer & Vacca (1998). Ages between 4 and 4.5 Myr, 3 and 3.5 Myr and 3.5 and 4 Myr are found for regions H13, CDT3 and 74C with oxygen abundances of 0.2, 0.4 and 0.5 solar, respectively. We have then used the ST99 models to predict the SEDs of the ionizing clusters. These models assume the same stellar tracks as those used by Schaerer & Vacca (1998) and, at the metallicities here considered, represent a good approach to the ionizing spectra (see Stasińska et al. 2001).

The predicted SEDs were then used as input to the latest version of CLOUDY (Ferland 1999). The models also require of an electron density (assumed to be constant throughout the nebula), the derived gas metallicity (see paper I) and an ionization parameter, U = Q(H)/(4$`\pi `$R2nec); where Q(H) is the total flux of ionizing photons emitted by the cluster per second, R the distance from the ionizing source to the illuminated face of the cloud, ne the electron density and c the speed of light. In spherical geometry, the average U is proportional to (Q(H)n$`{}_{e}{}^{}ϵ_{}^{2}`$)1/3 where $`ϵ`$ is the gas filling factor. In the case of ionization bounded regions, Q(H) must equal, in the absence of dust, the number of ionizing photons derived from H$`\alpha `$ recombinations. Our models have considered values of the ionization parameter between -4.0 $`<`$ log U $`<`$ -2.0, through the span of the WR phase which, combined with the observed numbers of ionizing photons (Q(H)) imply R values of the order of 1020 cm - 1021 cm ($``$ 30 - 300 pc). The thickness of the ionized gas shell is of the order of 1019, less than 10 % of the total dimensions of the region, which results in a plane-paralel geometry.

Figure 1 shows the run of the emission line intensities relative to H$`\beta `$, and the H$`\beta `$ equivalent width and luminosity as a function of the ionized gas shell thickness. The figure also compares the model predictions in each case with the observed values (horizontal bars).

Figure 1 clearly shows that a consistent fit to all observed quantities can only be found if the regions are considered to be matter bounded, i.e. only if there is an important escape of ionizing radiation from the nebulae. Models in which all photons are absorbed (most right-hand points in Figure 1) lead to large departures from all observed values. However, if one assumes the nebulae to be matter bounded, a satisfactory agreement between all predicted quantities and the observations is found at a common distance to the ionizing cluster for each case. Clearly, for these models we have also demanded for consistency with the observed H$`\alpha `$ luminosity which has been found through an iterative process. Successful models have hydrogen column densities of the order of 1019 cm-2 while ionization bounded models reach column densities close to 1020 cm-2.

The fraction of escaping photons depends both on their energy and the gas column density. Table 1 gives the main physical conditions of the observed regions. Columns 1, 2 and 3 identify the GEHRs and list the observed metallicity and resultant age of the best fitting models for the ionizing clusters. Columns 4 indicates the range in ionization parameter used, and columns 5-7 indicate both, the absolute and relative values of the photon flux escaping in three different energy bins (13.6 - 24.5 eV, 24.5 - 54.4 eV, and larger than 54.4 eV). The log of the total values of photons escaping the nebulae range from 49.80 $`<`$ log Qesc(H) $`<`$ 51.51, which imply an escape fraction between 10% (for region H13) and 73% (for region 74C).

## 3 DISCUSSION AND CONCLUSIONS

The observed emission lines of GEHR have been widely used to derive the properties of the ionizing radiation and link them to those of the stellar populations generating it, but always under the assumption that the ionized region is ionization bounded. In the absence of an independent constraint on the age of this stellar population, such as the presence of WR stars, successful fits to the emission line spectra of most GEHR have been found within a narrow range of ages, between 2 and 2.5 Myr (Bresolin et al. 1999; Dopita et al. 2000). In fact, if similar methods of analysis were applied to the regions modelled here, we would obtain a similar result. For these regions, however, we have accurate determinations of the gas metallicity and their age, provided by the analysis of the WR spectral features, and therefore we can synthesize the corresponding SEDs. However, the derived SEDs are unable to reproduce the observed emission line spectra unless the GEHRs are density bounded. In such a case, we have shown here that a satisfactory fit can be found in the three analyzed cases.

Our results suggest that the leaking of ionizing photons do depend on the evolutionary stage of the regions. Regions 74C and CDT3 show the highest escape fractions, consistent with the prominent observed WR features in the former spectrum, and the high excitation lines (\[FeII\],, \[FeIII\]) in the latter one. However, region H13, which also shows a prominent WR feature, shows a far lower escape fraction which can be understood from its lower metallicity (0.2 solar versus 0.5 solar in 74C and CDT3), affecting the WR phase strength. It is plausible that large photon escaping fractions might be expected in those regions that experience a hard WR phase, this one depending on both the age and metallicity.

An important implication from the escape of ionizing photons concerns the EW(H$`\beta `$) - age relation, that does not hold for matter bounded regions. It is a well known fact that few HII regions show EW(H$`\beta `$) as large as those predicted from evolutionary synthesis models (e.g. Melnick, Terlevich & Eggleton, 1985; Mas-Hesse & Kunth 1998). The leak of Lyman continuum radiation provides a natural explanation in the correct direction.

A close analysis of the energy distribution of the escaping radiation (see Table 1) shows that the largest fraction has energies between 13.6 and 24.4 eV (up to 80% in the case of 74C). However, if escaping to incident fractions are considered, helium ionizing photons are clearly dominant (e.g. up to 84% of the incident photons with energies between 24.4 and 54.4 eV escape from region 74C). This would explain why the observed emission line spectrum is so similar to that produced by a star cluster without WR stars.

Furthermore, the derived large amounts of escaping ionizing photons indicate that GEHRs and stellar OB associations may significantly contribute to the ionization of the diffuse gas layers observed above the disks of spiral galaxies. Though, to our knowledge, no measurements of the diffuse interstellar gas (DIG) are available for NGC 628, NGC 1232 and NGC 4258, recent studies by Ferguson et al. (1996), Oey & Kennicutt (1997) and Zurita et al. (2000) stress the fact that the integrated escaping flux from disk galaxies can attain up to 1041 erg s-1, enough to account for the diffuse H$`\alpha `$ flux and fully compatible with our derived values for single GEHRs.

In our models the lower energy photons emitted by the central star cluster are more effectively absorbed by the HII region gas while a large fraction of the photons with energy between 1.8 to 4 Ry escape the nebula. Therefore, the resulting spectral energy distribution of the escaping photons is slightly harder than that of the ionizing star cluster. For the three studied regions the Q(He)/Q(H) ratio of the escaping photons is: 0.22, 0.20 and 0.19 which should be compared to 0.10, 0.14 and 0.15 for the respective ionizing clusters.

One still unsolved problem for the interpretation of observations of the DIG above the galactic disk is the apparently low value of the HeI$`\lambda `$5876 $`\AA `$ /H$`\beta `$ (Reynolds & Tufte 1995) which would be in contradiction with our derived Q(He)/Q(H) ratios. Our studied regions, however, are not representative of galactic HII regions ionized by a single star, but of high luminosity Giant Extragalactic HII Regions on the disks of spiral galaxies like, for example, NGC891. For the DIG in this galaxy the average value of HeI$`\lambda `$5876 $`\AA `$ /H$`\beta `$ is 0.1 (Rand 1997). According to Bresolin, Kennicutt & Garnett (1999) the value of this ratio for stellar effective temperatures between 40000 and 50000 K is highly dependent of metallicity and reaches the value of 0.1 for a metallicity of about 1/4 solar. In fact, most HII galaxies show values of HeI$`\lambda `$5876 $`\AA `$ /H$`\beta `$ around 0.10 (Izotov & Thuan 1998 and references therein).

In summary, by applying evolutionary synthesis models to well-studied GEHRs with observed WR features, we find that these regions need to be matter bounded in order to reconcile model predictions with observations. The amount of escaping ionizing photons per unit time from these regions range from 49.80 $`<`$ log Qesc(H) $`<`$ 51.51, which implies an escape fraction between 10% and 73% of the available incident photons. These fractions seem to increase with both the strength of the WR phase and the metallicity of the ISM. The implication is thus that the mass of the ionizing clusters is larger than the value derived in a straight manner from the observed H$`\alpha `$ luminosity. It is also clear that the ISM is highly non-uniform and that the matter swept up in the shells produced by the mechanical energy deposited by the WR sources and other massive stars is not sufficient to trap the ionization fronts. This is also the result obtained from numerical calculations that consider both the mechanical energy of massive stellar clusters as well as their ionizing luminosity. Tenorio-Tagle et al. (1999) have shown how once the shell of swept up matter becomes Rayleigh - Taylor unstable and fragments, as it evolves out of a galaxy disk into the halo (the blowout phenomenon), it then allows not only for the venting of the hot (wind and SN) matter into the halo, but also for the leakage of a large fraction of the ionizing radiation. The latter soon establishes a giant conical HII region in the low density halo. The low densities in the halo lead to a long recombination time and thus to a large leakage of photons into the IGM. The situation changes later, once the shock is able to sweep enough halo matter, enhancing locally the number of recombinations in the expanding shell. This leads eventually to the trapping of the ionization front. Consecuently, depending on the stage of the evolution, a large fraction of the ionizing flux escapes the nebulae and is freely available to impact upon the gas at large distances from the host galaxy plane (Collins & Rand 2001) and it is also likely to escape the galaxy and cause an important ionization of the IGM. Clearly, further analysis of this kind must be done to infer whether or not our derived escape fractions are typical of other GEHRs.

## Acknowledgements

We thank the anonymous referee for many comments and suggestions which significantly improved the content of the paper.

G.T.-T. is grateful to an IBERDROLA Visiting Professorship to UAM during which part of this work was completed. This work has been partially supported by DGICYT project AYA-2000-0973.