text:
# Continuous demixing at liquid-vapor coexistence in a symmetrical binary fluid mixture

## I Introduction and background

Symmetrical binary fluid mixtures are two-component fluid models in which the configurational energy is invariant with respect to interchange of the two particle species. One example is the symmetrical Lennard-Jones (LJ) mixture in which the interactions between particles of species 1 and 2 are controlled by a LJ potential with scale parameters $`\sigma _{11}=\sigma _{22}=\sigma _{12}=\sigma `$ and interaction strengths $`ϵ_{11}=ϵ_{22}=ϵϵ_{12}`$. The phase diagram of the mixture is uniquely specified by the ratio of interaction strengths between the unlike and like species: $`\delta ϵ_{12}/ϵ`$.

Symmetrical mixtures have been the subject of considerable recent attention on account of their surprisingly rich phase behaviour PAN88 ; RECHT ; DEMIGUEL ; WILDING97 ; WILDING98 . The phase diagram of the mixture is typically represented in terms of the temperature $`T`$, the overall number density $`\rho =\rho _1+\rho _2`$, and the concentration $`c=\rho _2/(\rho _1+\rho _2)`$. Most attention has focused on the case of equal species concentration, $`c=0.5`$. Within this particular symmetry plane of the full phase diagram, the $`\rho `$-$`T`$ dependence of the phase behaviour on the value of $`\delta `$ has been systematically studied using simulation and mean field theory for the case of a symmetrical square-well mixture WILDING98 . The findings of that work are summarized in fig 1. of ref. WILDING98 . Three distinct classes of phase diagram were found depending on the choice of the parameter $`\delta `$. Specifically, for large $`\delta <1`$ (fig. 1(a) of ref. WILDING98 ) there occurs a “$`\lambda `$-line” of critical demixing transitions which intersects (and is truncated by) the liquid-vapor line at a critical end point (CEP). For small $`\delta >0`$, on the other hand, the line of demixing transitions intersects the liquid-vapor line at the liquid-vapor critical critical point, forming a tricritical point (fig. 1(c) of ref. WILDING98 ). Intermediate between these two regimes, one observes both a liquid-vapor critical point and a tricritical point at higher densities (fig. 1(b) of ref. WILDING98 ). In the simulations this latter regime was observed to occur for $`0.65\delta 0.68`$, while within the particular mean field theory employed it occurred for $`0.605\delta 0.708`$.

It has been suggested WILDING98  that the scenario of phase behaviour outlined above for the square well mixture, might in fact be generic to all symmetrical fluids. On the theoretical side, support for this proposal has come from studies of variety of model systems, including mean field studies of a lattice fluid model SCHOEN  and integral equation theories KAHL01 ; PASCHINGER1 ; PASCHINGER2 ; CACCAMO  of the hard core Yukawa fluid mixture (HCYFM). The same overall scenario was initially reported for the hard core $`+`$ Lennard-Jones mixture by Pini et al PINI00  in a study employing the Hierarchical Reference Theory (HRT)– a powerful fluctuation-inclusive approach based on renormalization group techniques PAROLA91 . By contrast, the available simulation data for models other than the square well fluid is less comprehensive and has, in the main part, concentrated on the regime of large $`\delta `$. Specifically, MC simulations of the symmetrical LJ fluid WILDING97  found evidence for a CEP (and associated coexistence curve singularities) for the particular case $`\delta =0.7`$. Similarly a CEP has been reported in simulation studies of the HCYFM for $`\delta =0.9`$ CACCAMO .

Very recently, however, the situation with specific regard to the CEP regime has become less clear following two theoretical studies which question its existence in the expected range of $`\delta `$. In the first of these, Antonevych et al. ANTON02  applied the modified hypernetted chain (MHNC) integral equation theory to the symmetrical LJ fluid, but found no evidence for a CEP at $`c=0.5`$. Instead they reported that the intermediate regime (fig. 1(b) of ref. WILDING98 ) persists right up to $`\delta =0.81`$, (the limit of their study). This value of $`\delta `$ is considerably greater than that ($`\delta =0.7`$) at which CEP behaviour was originally found in the simulations of the same model. On this basis, Antonevych et al. suggested that the demixing transition could always be first order at liquid-vapor coexistence, although they speculated that it might be only weakly so. Furthermore, they questioned the ability of simulation to distinguish this possibility from a true CEP.

In a separate study, Pini et al. PINI02  have applied the HRT approach to the HCYFM. They too observed no CEP, instead reporting (in accord with the finding of Antonevych et al. ANTON02 ) that the intermediate regime persists (albeit weakly) to $`\delta =0.8`$, the limit of their study. Moreover, these authors report NOTE  that in their previous study of the LJ mixture PINI00 , the apparent observation of a CEP was an artifact stemming from the low resolution of the grid in the $`\rho `$-$`c`$ plane on which the HRT equations were solved.

Clearly, therefore, the evidence emerging from the various theoretical and simulation studies is contradictory with regard to the existence of a CEP and the range of values of $`\delta `$ for which it occurs. Most significant in our view is the failure of HRT to observe a CEP for $`\delta <0.8`$. Given the fluctuation-inclusive nature of this approach (and its success in other contexts, see e.g. WILDING01 ), one cannot simply ascribe discrepancies between the simulations and the MHNC study of the LJ fluid ANTON02  to the neglect of critical fluctuations by the latter.

Questions have been raised too regarding the reliability of simulation in distinguishing a CEP from a weak first order transition ANTON02 . While, in our opinion, the original work WILDING97  does provide compelling evidence for CEP behaviour, it is certainly the case that owing to the prevailing computational constraints, the range of system sizes studied was smaller than one might have wished. This in turn precluded an investigation of corrections to scaling and the approach to the limiting critical behaviour.

In the light of the above considerations, it would seem worthwhile to revisit the original system in which a CEP was observed, with a view to performing a more comprehensive determination of the nature of the demixing transition at liquid-vapor coexistence. To this end we have carried out a detailed Monte Carlo (MC) finite-size scaling (FSS) simulation study of the symmetrical LJ fluid with $`\delta =0.7`$, for which it was previously claimed that a CEP occurs WILDING97 . The new study has been executed along similar lines to the original one, but employed a considerably larger range of systems sizes, thus permitting a more sophisticated finite-size scaling analysis. The results demonstrate (unambiguously, we believe) that the model does indeed exhibit a CEP of the expected Ising type.

## II Method and results

The simulation methodology employed here is broadly similar to that described in ref. WILDING97  and we refer the reader to that paper for a detailed description. Briefly, we have performed a MC simulation study of the liquid-vapor coexistence line of the symmetrical LJ fluid mixture with $`\delta =0.7`$ in the neighborhood of the demixing transition. The simulations were carried out within the grand canonical ensemble and the chemical potentials of the two species were constrained to be equal ($`\mu _1=\mu _2=\mu `$), implying on symmetry grounds that $`c=0.5`$. The system was confined to a cubic box of linear dimension $`L`$ having periodic boundary conditions. In all, seven system sizes were studied having $`L=10\sigma ,12.5\sigma ,15\sigma ,17.5\sigma ,20\sigma ,25\sigma `$ and $`30\sigma `$. The largest of these systems contains some $`1.6\times 10^4`$ particles in the liquid phase at the demixing point. The range of system sizes studied here is to be compared with that of the original study WILDING97 , for which the largest size attained was $`L=12.5\sigma `$. The greater range of sizes studied here was partly facilitated by the use of high performance parallel computers.

The liquid-vapor coexistence curve $`\mu _{\mathrm{cx}}(T)`$ was obtained using multicanonical extended sampling techniques BERG . These allow both the liquid and vapor phases to be sampled (and thus connected) in a single simulation run. A suitable form for the requisite preweighting function was obtained using a variant of the recently proposed Wang-Landau method WANGLAN . To locate the coexistence line, the equal peak area criterion was applied to the measured form of the number density distribution function $`p_L(\rho )`$ EW ; WILDING95 . Histogram reweighting techniques were employed to fine-tune this procedure and to facilitate the mapping of a portion of the liquid-vapor coexistence curve in the neighborhood of the demixing transition.

Fig. 1 shows the measured form of the liquid-vapor coexistence line $`\mu _{\mathrm{cx}}(T)`$ NOTE0 . To elucidate the nature of the demixing transition along the liquid branch of this tightly determined coexistence line, we have measured the probability distribution function $`p_L(m)`$ of the demixing order parameter $`m`$, where $`m=(N_1N_2)/(N_1+N_2)`$, with $`N_1`$ and $`N_2`$ the instantaneous counts of the respective particle species. The form of $`p_L(m)`$, for the $`L=25\sigma `$ system at temperatures spanning the demixing transition, is shown in fig. 2. We find that $`p_L(m)`$ evolves smoothly between the strongly double peaked form (corresponding to a demixed liquid) at low temperature and a singly peaked form (corresponding to a mixed liquid) at higher temperatures.

To quantify this evolution more precisely, we have measured the temperature dependence of the fourth order cumulant ratio: $`U_L=1m^4/3m^2^2`$, which provides a dimensionless measure of the shape of a distribution. The results (fig 3) demonstrate (for all system sizes) a smooth evolution of the concentration distribution from the ordered to the disordered phase. The rate of change increases with the system size, but there is no evidence of a jump discontinuity. More significantly, there is a fairly well defined crossing point of the curves for $`U_L0.46`$ close to the reduced temperature $`T^{}1/ϵ=0.959`$. Such an intersection reflects scale invariance in the concentration distribution and does not occur at a first order phase transition. We note further that the magnitude of the cumulant ratio at the intersection point corresponds closely to the known Ising universal value, $`U_L^{\mathrm{Ising}}=0.465`$ TSYPIN00 .

Although the results of fig. 3 constitute strong evidence for a critical demixing transition, a close inspection reveals that the cumulant crossings do not occur precisely at a unique temperature for all system sizes. Such behaviour mirrors that observed in other contexts (see eg. WILDING95 ; WILDING01 ) where it was found to be attributable to corrections to finite-size scaling. Since these corrections have a known universal scaling form, one can attempt to fit them, thereby permitting an extrapolation to the thermodynamic limit. The procedure involves determining, for each system size, the coexistence state point for which the form of $`p_L(m)`$ on the liquid branch best matches the (independently known TSYPIN00 ) fixed point order parameter distribution function $`p^{}(m)`$ appropriate to the Ising universality class. The matching temperature as a function of system size is expected to scale like $`T_c(L)L^{(1+\theta )/\nu }`$ WILDING95 , where $`\theta 0.50(2)`$ is the Ising correction-to-scaling exponent JUSTIN  and $`\nu =0.6294(4)`$ is the correlation length exponent FERRENBERG . Fig. 4 shows the data expressed in this form, which indeed exhibits the anticipated scaling behaviour. Extrapolation of a linear fit to $`L=\mathrm{}`$ yields an estimate for the critical temperature $`T_{\mathrm{cep}}=0.9595(3)`$, which occurs for $`\mu _{\mathrm{cep}}=3.704(2)`$.

It is instructive to plot the finite-size forms of $`p_L(m)`$ at the estimated critical point in order to expose the scale and character of corrections to finite-size scaling. Figure 5 compares the distributions for the $`L=12.5\sigma `$ and $`L=25\sigma `$ system sizes with that of $`p^{}(m)`$ TSYPIN00 . One sees that the results for the larger system size closely match the limiting form.

Finally, we have considered the behaviour of the number density along the liquid branch at coexistence. Were the demixing transition first order along the coexistence line, one would expect that coupling between the concentration and the number density would engender a first order phase transition in the liquid branch number density (cf. fig. 1(b) of ref. WILDING98 ). This would be manifest as a discontinuity in the density of the coexisting liquid phase and by a splitting of the liquid-phase peak in the density distribution $`p_L(\rho )`$. Fig. 6 shows the measured form of the liquid-phase peak for the $`L=25\sigma `$ system size, for several coexistence state points spanning the demixing temperature. Also shown (inset) is the average liquid phase density. The distribution of the density in the liquid phase is singly peaked and we find no evidence of a discontinuity in the average density for any of the system sizes studied NOTE1 .

## III Discussion and conclusions

In summary, we have performed a finite-size scaling analysis of the demixing transition of a symmetrical LJ fluid along the liquid branch of the liquid-vapor coexistence line. The results demonstrate that for the case $`\delta =0.7`$, this transition is continuous and of the expected Ising type DIEHL . Accordingly we can conclude that this system does indeed exhibit a critical end point, as originally proposed in ref WILDING97 . This clear finding should provide a benchmark for testing the reliability of theoretical approaches to the phase behaviour of symmetrical mixtures.

With regard to the original motivation for this study, namely the failure of both the MHNC integral equation study and HRT to observe a CEP in the expected range of $`\delta `$ ANTON02 ; PINI02 , it is difficult to comment definitively at this stage on the source of the discrepancies. We note, however, that problems with the accuracy of the MHNC approach for studying the liquid branch of the HCYFM have previously been reported CACCAMO . With regard to the HRT, it would certainly be useful if this method could be applied to other symmetrical fluid models in order to determine whether there is any systematic model dependence on the qualitative features of the phase behaviour.


format: text

text:
# Asymptotic equivalence of the jackknife and infinitesimal jackknife variance estimators for some smooth statistics

## 1 Introduction

This note concerns the asymptotic behavior of the jackknife variance estimator $`v_{jack}`$, especially regarding its relationship to the infinitesimal jackknife variance estimator $`v_{ijack}`$. We consider, in particular, the variance estimates $`v_{jack}`$ and $`v_{ijack}`$ for smoothly trimmed L-statistics and for smooth functions of the sample mean. We prove that $`v_{jack}`$ and $`v_{ijack}`$ are often asymptotically equivalent to one another in the sense that

(1.1)
$$v_{jack}v_{ijack}=\mathrm{O}_\mathrm{p}\left(n^h\right)$$

for some $`h>0`$. The equivalence of $`v_{jack}`$ and $`v_{ijack}`$ can sometimes be used to prove that $`v_{jack}`$ is asymptotically normal, but — remarkably — it also holds even when these estimators are not asymptotically normal.

These issues are relevant to the following scenario of statistical practice: One wishes to estimate some functional $`T(p)`$ of an unknown population distribution $`p`$, and, to this end, one draws $`n`$ samples from the population and uses $`T(ϵ_n)`$ to estimate $`T(p)`$, where $`ϵ_n`$ is the empirical distribution of the $`n`$ samples. One would then like an estimate of the sampling variance of $`T(ϵ_n)`$. Two widely-used nonparametric estimates of this variance are $`v_{jack}`$ and $`v_{boot}`$, the jackknife and bootstrap variance estimates. Having obtained one of these estimates, one naturally desires to know how accurate it is. To assess the accuracy of the usual Monte Carlo approximation of $`v_{boot}`$, one may use the jackknife-after-bootstrap technique of Efron (1992). This note concerns the asymptotic behavior of $`v_{jack}`$ and a closely related estimator $`v_{ijack}`$, the infinitesimal jackknife.

Beran(1984) showed that $`v_{jack}`$, $`v_{ijack}`$, and $`v_{boot}`$ are asymptotically equivalent, and asymptotically normal, if the functional $`T(p)`$ has a well-behaved second-order functional derivative. While the proof of Beran (1984) requires a strong statement of the Dvoretsky-Kiefer-Wolfowitz inequality to handle the asymptotics of $`v_{boot}`$, the asymptotic equivalence of $`v_{jack}`$ and $`v_{ijack}`$ is easier to prove, and does not require second-order differentiability of $`T`$. Indeed, we shall see that $`v_{jack}`$ and $`v_{ijack}`$ are asymptotically equivalent if $`T`$ has a well-behaved first-order derivative and $`v_{ijack}`$ is consistent as an estimator of the variance of $`T(ϵ_n)`$.

The equivalence of $`v_{jack}`$ to $`v_{ijack}`$ can help one to determine the asymptotic variance of the former. For example, Gardiner and Sen (1979) have carefully studied the asymptotic normality of $`v_{ijack}`$ in the context of L-statistics. We shall see in Section 4 that their work establishes the asymptotic normality of $`v_{jack}`$, too, thanks to the equivalence of $`v_{jack}`$ and $`v_{ijack}`$ for variance estimation of L-statistics.

Sometimes $`v_{jack}`$ and $`v_{ijack}`$ are asymptotically equivalent even when they are not asymptotically normal. Consider, for example, the estimation of the variance of a function $`g`$ of the sample mean. If $`g`$ is once, but not twice, continuously differentiable, then $`v_{jack}`$ and $`v_{ijack}`$ may not be asymptotically normal, yet they will satisfy (1.1) as long as $`g^{}`$ is Hölder continuous of order $`h`$.

After the necessary definitions are presented in the next section, we prove in Section 3 that $`v_{jack}`$ and $`v_{ijack}`$ are asymptotically equivalent as estimators of the variance of smooth functions of the sample mean. In Section 4 we discuss the asymptotic normality of $`v_{jack}`$ as an estimator of the variance of trimmed L-statistics.

## 2 Background and definitions

Let $`p`$ be a probability measure on a sample space $`𝒳`$. Given $`n`$ samples from $`𝒳`$, sampled independently under the probability law $`p`$, one desires to estimate the value $`T(p)`$ of some real functional $`T`$ on the space $`𝒫(𝒳)`$ of all probability measures on $`𝒳`$. Denote by $`ϵ_n`$ the map that converts $`n`$ data points $`x_1,x_2,\mathrm{},x_n`$ into the empirical measure

(2.1)
$$ϵ_n(x_1,x_2,\mathrm{},x_n)=\frac{1}{n}\underset{i=1}{\overset{n}{}}\delta (x_i)$$

where $`\delta (x_i)`$ denotes a point-mass at $`x_i`$. The plug-in estimate of $`T(p)`$ given the data $`𝐱=(x_1,\mathrm{},x_n)`$ is

(2.2)
$$T_n=T(ϵ_n(𝐱)).$$

Suppose $`T_n`$ is an asymptotically normal estimator of $`T(p)`$, i.e., suppose the distribution of $`n^{1/2}(T_nT(p))`$ tends to $`𝒩(0,\sigma ^2)`$. The jackknife is a computational technique for estimating $`\sigma ^2`$: one transforms the $`n`$ original data points into $`n`$ pseudovalues and computes the sample variance of those pseudovalues.

Given the data $`𝐱=x_1,x_2,\mathrm{},x_n`$, the jackknife pseudovalues are

$$Q_{ni}=nT_n(ϵ_n)(n1)T(ϵ_{ni})i=1,2,\mathrm{},n$$

with $`ϵ_n`$ as in (2.1) and

(2.3)
$$ϵ_{ni}=\frac{1}{n1}\underset{ji}{}\delta (x_j).$$

The jackknife variance estimator is

(2.4)
$$v_{jack}(x_1,x_2,\mathrm{},x_n)=\frac{1}{n1}\underset{i=1}{\overset{n}{}}\left(Q_{ni}\overline{Q_n}\right)^2$$

where $`\overline{Q_n}=\frac{1}{n}Q_{nj}`$. The variance estimator $`v_{jack}`$ is said to be consistent if $`v_{jack}\sigma ^2`$ almost surely as $`n\mathrm{}`$. Sufficient conditions for the consistency of $`v_{jack}`$ are given in terms of the functional differentiability of $`T`$. An early result of this kind states that $`v_{jack}`$ is consistent if $`T`$ is strongly Fréchet differentiable (Parr(1985)), and it is now known that $`v_{jack}`$ is consistent even if $`T`$ is only continuously Gâteaux differentiable as defined in Shao(1993).

A functional derivative of $`T`$ at $`p`$, denoted $`T_p`$, is a linear functional that best approximates the behavior of $`T`$ near $`p`$ in some sense. For instance, a functional $`T`$ on the space of bounded signed measures $`(𝒳)`$ is Gâteaux differentiable at $`p`$ if there exists a continuous linear functional $`T_p`$ on $`(X)`$ such that

$$\underset{t0}{lim}\left|t^1\left(T(p+tm)T(p)\right)T_p(m)\right|=0$$

for all $`m(𝒳)`$. The concept of Hadamard differentiability is more relevant to statistical asymptotics, for the fluctuations of $`T(ϵ_n)`$ about $`T(p)`$ are asymptotically normal if $`T`$ is Hadamard differentiable at $`p`$. A functional $`T:𝒫()`$ is Hadamard differentiable at $`p`$ if there exists a continuous linear functional $`T_p`$ on $`()`$ such that

$$\underset{t0}{lim}\left|t^1\left(T(p+tm_t)T(p)\right)T_p(m)\right|=0$$

whenever $`\{m_t\}_t`$ is such that $`\underset{t0}{lim}m_t=m`$ and $`m_t()=0`$ for all $`t`$, the topology on $`()`$ being the one induced by the norm $`m=\underset{t}{sup}\left\{\left|m((\mathrm{},t])\right|\right\}`$. If $`T`$ is Hadamard differentiable at $`p`$, the variance of $`n^{1/2}T(ϵ_n)`$ tends to

(2.5)
$$\sigma ^2=𝔼_p\varphi _p^2$$

as $`n\mathrm{}`$, where $`\varphi _p(x)`$ is the influence function

(2.6)
$$\varphi _p(x)=T_p(\delta (x)p)$$

(this can be shown via the Delta Method using Donsker’s theorem (van der Waart(1998))). The infinitesimal jackknife estimator (Jaeckel(1972)) of $`\sigma ^2`$ is obtained by substituting the empirical measure $`ϵ_n`$ for $`p`$ in (2.5):

(2.7)
$$v_{ijack}=𝔼_{ϵ_n}\varphi _{ϵ_n}^2.$$

## 3 Functions of the mean

When $`q`$ is a measure, we denote $`xq(dx)`$ by $`\overline{q}`$ if the integral is defined. Let $`gC^1()`$ and let

$$T(m)=g\left(\overline{m}\right)$$

be defined for all finite signed measures $`m`$ with finite first moment. The functional derivative at $`m`$ of $`T`$, evaluated at $`q`$, is $`T_m(q)=g^{}\left(\overline{m}\right)\overline{q}`$; the influence function (2.6) is $`\varphi _m(x)=g^{}\left(\overline{m}\right)\left(x\overline{m}\right)`$. Suppose that $`x_1,x_2,\mathrm{}`$ are iid $`p`$, and $`p`$ has a finite second moment. Let $`T_n`$ denote the plug-in estimator defined in (2.2). Then the asymptotic variance of $`n^{1/2}\left(T_nT(p)\right)`$ is

(3.1)
$$\sigma ^2=g^{}(\overline{p})^2\left\{x^2p(dx)\overline{p}^2\right\}.$$

Let $`v_{jack}`$ and $`v_{ijack}`$ denote the jackknife and infinitesimal jackknife variance estimates of $`\sigma ^2`$.

###### Proposition 3.1

If $`g^{}`$ is Hölder continuous of order $`h>1/2`$ (with global Hölder constant) and $`p`$ has a finite moment of order $`2+2h`$ then

$$v_{jack}v_{ijack}=\mathrm{O}_\mathrm{p}\left(n^h\right).$$

Proof: Setting $`\mathrm{\Delta }_{ni}=\left(Q_{nj}\overline{Q_n}\right)\varphi _{ϵ_n}(x_i)`$, formula (2.4) for $`v_{jack}`$ yields

(3.2)
$$v_{jack}=𝔼_{ϵ_n}\varphi _{ϵ_n}^2+\frac{1}{n1}𝔼_{ϵ_n}\varphi _{ϵ_n}^2+\frac{2}{n1}\underset{i=1}{\overset{n}{}}\varphi _{ϵ_n}(x_i)\mathrm{\Delta }_{ni}+\frac{1}{n1}\underset{i=1}{\overset{n}{}}\mathrm{\Delta }_{ni}^2.$$

The second term on the right hand side of (3.2) is $`\mathrm{O}_\mathrm{p}\left(1/n\right)`$ since

$$𝔼_{ϵ_n}\varphi _{ϵ_n}^2=\frac{1}{n}\underset{i=1}{\overset{n}{}}\varphi _{ϵ_n}^2(x_i)=\frac{1}{n}\underset{i=1}{\overset{n}{}}g^{}\left(\overline{ϵ_n}\right)^2(x_i\overline{ϵ_n})^2$$

converges almost surely to $`\sigma ^2`$.

To control the last two terms on the right hand side of (3.2) we need a bound on $`\mathrm{\Delta }_{ni}`$. Recall the notation $`ϵ_{ni}`$ of (2.3). Since $`g`$ is differentiable, $`g\left(\overline{ϵ_{nj}}\right)g\left(\overline{ϵ_{ni}}\right)=g^{}\left(\eta _{ji}\right)\left(\overline{ϵ_{nj}}\overline{ϵ_{ni}}\right)`$ for some $`\eta _{ji}`$ between $`\overline{ϵ_{ni}}`$ and $`\overline{ϵ_{nj}}`$, so that

$$Q_{ni}\overline{Q_n}=\frac{n1}{n}\underset{j=1}{\overset{n}{}}\left(g\left(\overline{ϵ_{nj}}\right)g\left(\overline{ϵ_{ni}}\right)\right)=\frac{n1}{n}\underset{j=1}{\overset{n}{}}g^{}\left(\eta _{ji}\right)\left(\overline{ϵ_{nj}}\overline{ϵ_{ni}}\right).$$

Therefore, since $`\varphi _{ϵ_n}(x_i)=g^{}\left(\overline{ϵ_n}\right)(x_i\overline{ϵ_n})=\frac{1}{n}_jg^{}\left(\overline{ϵ_n}\right)(x_ix_j)`$,

$`\mathrm{\Delta }_{ni}=\left(Q_{ni}\overline{Q_n}\right)\varphi _{ϵ_n}(x_i)`$ $`=`$ $`{\displaystyle \frac{n1}{n}}{\displaystyle \underset{j=1}{\overset{n}{}}}g^{}\left(\eta _{ji}\right)\left(\overline{ϵ_{nj}}\overline{ϵ_{ni}}\right){\displaystyle \frac{1}{n}}{\displaystyle \underset{j=1}{\overset{n}{}}}g^{}\left(\overline{ϵ_n}\right)(x_ix_j)`$
$`=`$ $`{\displaystyle \frac{1}{n}}{\displaystyle \underset{j=1}{\overset{n}{}}}\left(g^{}\left(\eta _{ji}\right)g^{}\left(\overline{ϵ_n}\right)\right)(x_ix_j).`$
But $`g^{}`$ is Hölder continuous of order $`h`$ and $`|\eta _{ji}\overline{ϵ_n}|<\mathrm{max}\{|\overline{ϵ_{nj}}\overline{ϵ_n}|,|\overline{ϵ_{ni}}\overline{ϵ_n}|\}`$, so the identity $`(n1)\left(ϵ_nϵ_{ni}\right)=\delta _{x_i}ϵ_n`$ implies that

$$\left|g^{}\left(\eta _{ji}\right)g^{}\left(\overline{ϵ_n}\right)\right|C\left(|\overline{ϵ_{nj}}\overline{ϵ_n}|^h+|\overline{ϵ_{ni}}\overline{ϵ_n}|^h\right)C(n1)^h\left(|\overline{ϵ_n}x_j|^h+|\overline{ϵ_n}x_i|^h\right),$$

where $`C`$ is a global Hölder constant for $`g^{}`$. It follows that

$$\left|\mathrm{\Delta }_{ni}\right|=C(n1)^h\frac{1}{n}\underset{j=1}{\overset{n}{}}\left(|\overline{ϵ_n}x_j|^h+|\overline{ϵ_n}x_i|^h\right)\left(|\overline{ϵ_n}x_j|+|\overline{ϵ_n}x_i|\right).$$

With this bound on $`\mathrm{\Delta }_{ni}`$, and assuming that $`p`$ has a finite moment of order $`2(1+h)`$, it may be shown that

$$\frac{1}{n}\underset{i=1}{\overset{n}{}}\mathrm{\Delta }_{ni}^2=\mathrm{O}_\mathrm{p}(n^{2h}),$$

and then, by the Cauchy-Schwartz inequality, that

$$\left|\frac{1}{n}\underset{i=1}{\overset{n}{}}\varphi _{ϵ_n}(x_i)\mathrm{\Delta }_{ni}\right|=\mathrm{O}_\mathrm{p}\left(n^h\right).$$

Substituting the preceding estimates in (3.2) completes the proof. $`\mathrm{}`$

Consider $`g(x)=x\text{sgn}(x)x^2`$. This function has a Lipschitz continuous derivative but does not have a second-order derivative at $`0`$. By Proposition 3.1,

$$v_{jack}v_{ijack}=\mathrm{O}_\mathrm{p}\left(1/n\right).$$

However, for some population distributions having mean $`0`$, one can prove that $`v_{ijack}`$ is not asymptotically normal, and simulations suggest that $`v_{jack}v_{boot}`$ is $`\mathrm{O}_\mathrm{p}\left(1/\sqrt{n}\right)`$ rather than $`\mathrm{O}_\mathrm{p}\left(1/n\right)`$. This example shows that $`v_{jack}`$ is more closely related to $`v_{ijack}`$ than it is to $`v_{boot}`$.

## 4 Trimmed L-statistics

Suppose that $`\mathrm{}:(0,1)`$ is supported on $`[\alpha ,1\alpha ]`$ for some $`0<\alpha <1/2`$, and let

(4.1)
$$L(p)=_0^1P^1(s)\mathrm{}(s)𝑑s.$$

Here $`P^1`$ denotes the quantile function for $`p`$, i.e., $`P^1(s)=\mathrm{min}\{x:P(x)s\}`$ for $`0<s<1`$ where $`P`$ denotes the cdf of $`p`$. A plug-in estimate for $`L`$ is called a trimmed L-statistic, or a trimmed linear combination of quantiles.

If the weight function $`\mathrm{}`$ is continuous then $`L`$ is Hadamard differentiable at all $`p𝒫()`$ (see, e.g., Lemma 22.10 of van der Waart (1998)), and so the L-statistics are asymptotically normal (an original reference is Stigler (1974)). The asymptotic variance $`\sigma ^2`$ of the L-statistics may be estimated by $`v_{jack}`$, which converges almost surely to $`\sigma ^2`$ if $`\mathrm{}`$ is continuous (Parr(1985), Shao and Tu (1995)). The jackknife and infinitesimal jackknife would seem to be the only nonparametric methods of consistent variance estimation for L-statistics, aside from the bootstrap (Shucany and Parr(1982)).

We turn now to the question of the asymptotic normality of $`v_{jack}`$. In this regard, a variant of the L-functional (4.1) has been treated in the literature, namely

(4.2)
$$(p)=x\mathrm{}(P(x))p(dx).$$

If $`P`$ is continuous and strictly increasing then $``$ of (4.2) equals $`L`$ of (4.1). Beran (1984) proves that $`v_{jack}`$ for $``$ is asymptotically normal — and so is $`v_{boot}`$ — if $`\mathrm{}`$ is continuously differentiable and $`p`$ has bounded support. Section 2.2.3 of Shao and Tu (1995) incorrectly claims that $`v_{jack}`$ for $``$ is asymptotically normal if $`\mathrm{}`$ is Hölder continuous of order greater than $`1/2`$, and it also wrongly claims that the asymptotic variance equals $`\text{Var}\left(\varphi _p^2\right)`$, where $`\varphi _p`$ is the influence function of $``$. A detailed discussion of those errors is given in an unpublished technical report (Gottlieb (2001)). Nevertheless, a reworking of Definition 2.6 and Theorem 2.7 in Shao and Tu (1995) leads us to the following general proposition, which will presently be applied to the case where the $`T(ϵ_n)`$ are L-statistics:

###### Proposition 4.1

Let $`ϵ_n`$ denote the empirical distribution of $`n`$ iid samples from $`p`$, and let $`v_{jack}`$ and $`v_{ijack}`$ denote the jackknife and infinitesimal jackknife estimates of the variance of $`T(ϵ_n)`$.

Let $`q^{}q`$ denote the supremum of the absolute value of the difference between the cdf’s of $`q^{}`$ and $`q`$. Suppose that there exist positive constants $`C`$, $`\delta `$, and $`h`$ such that

(4.3)
$$T(q^{})=T(q)+_qT(q^{}q)+(q^{},q)$$

for all $`q^{},q`$ with $`q^{}p,qp<\delta `$, where the remainder $`|(q^{},q)|Cq^{}q^{1+h}`$. Then

$$v_{jack}v_{ijack}=\mathrm{O}_\mathrm{p}\left(n^h\right)$$

if $`v_{ijack}`$ is bounded in probability.

The straightforward proof of this proposition proceeds like the proof of Proposition 3.1 above, except that (4.3) is used to bound $`\mathrm{\Delta }_{ni}`$ in (3.2).

Now, let $``$ be a trimmed L-functional of the form (4.2) whose weight function $`\mathrm{}`$ is Hölder continuous of order $`h`$. Upon integrating the right hand side of (4.2) by parts, it becomes easy to verify that $``$ admits the expansion (4.3) near any $`p`$. Since $`v_{ijack}`$ converges almost surely, Proposition 4.1 implies that $`v_{jack}`$ and $`v_{ijack}`$ are asymptotically equivalent.

This equivalence allows us to conclude that $`v_{jack}`$ is asymptotically normal for many L-functionals, for Gardiner and Sen (1979) have found hypotheses that guarantee the asymptotic normality of $`v_{ijack}`$ for generalized L-functionals of the form (4.1). They begin by assuming that the cdf $`P`$ of the population distribution is continuous. In this case

$$v_{ijack}=𝔼_{ϵ_n}\varphi _{ϵ_n}^2=\mathrm{}(P_n(y))\left[P_n(yz)P_n(y)P_n(z)\right]\mathrm{}(P_n(z))𝑑y𝑑z.$$

In order to make contact with the work of Gardiner and Sen (1979), let us suppose that $`P`$ is continuous and strictly increasing. Their hypotheses are general enough to apply to non-trimmed L-statistics, but too complicated to be repeated here. Suffice it to say that their theorem applies under our current assumptions that $`\mathrm{}`$ is trimmed and that $`P`$ has no jumps or flats, if it is assumed in addition that $`P`$ does not have very heavy tails and that $`\mathrm{}`$ is piecewise continuously differentiable with Hölder continuity of order greater than $`1/2`$ at the cusps. (Imagine, for example, a piecewise-linear $`\mathrm{}`$ whose graph is shaped like a desert mesa; this is one of the weight functions recommended in Stigler (1973) for smoothly trimmed means.) In these cases $`v_{jack}`$ is asymptotically normal as well, by Proposition 4.1.

## 5 Acknowledgments

The author would like to thank Steve Evans for his advice and encouragement, and also Rudolf Beran. Support from the ESI and WPI in Vienna and the START project Nonlinear Schrödinger and quantum Boltzmann equations (FWF Y-137) is acknowledged.


text:
# Eta Production in the pp Scattering

## 1 Introduction

Zagreb-UCLA-ANL meson nucleon partial wave analysis opened the possibility to explain the number of intermediate energy data. Important result of the PWA are the eta nucleon elastic T matrices. The three-body processes can serve us to lower the error of the two-body T matrices. The idea is to take two-body T matrices from existing PWA’s and apply them to the three-body processes.

The most convenient process for that purpose is the $`pppp\eta `$. It is the isospin 1/2 process, what enables us to test the half of the amplitudes at a time. Dominant contribution to the production process comes through meson-exchange, and decay of nucleon resonances.

## 2 Model

Our model has been described earlier. Because of the large uncertainties in the treatment of final state interaction of the considered theoretical models, we have decided to ignore its influence at the moment. The FSI is significant near threshold, so energies from 100 to 500 MeV would be the most convenient place to test this model. Unfortunately, in this range, experimental data are very hard to get. Preliminary results are given in the Figure 1.

## 3 Discussion

Light and dark lines show constructive and destructive total cross section, respectably. Enormous discrepancy below few MeV are expected because we have completely disregarded final state interaction. Our preliminary results again do not show rho meson dominance - which would be expected from most of the other models. To get better picture about what is going on, we shall calculate further differential cross sections, and total cross sections for pn to $`\eta `$pn. Completed results will be reported elsewhere.


text:
# Dynamics and Scaling of 2D Polymers in a Dilute Solution

## Abstract

The breakdown of dynamical scaling for a dilute polymer solution in 2D has been suggested by Shannon and Choy \[Phys. Rev. Lett. 79, 1455 (1997)\]. However, we show here both numerically and analytically that dynamical scaling holds when the finite-size dependence of the relevant dynamical quantities is properly taken into account. We carry out large-scale simulations in 2D for a polymer chain in a good solvent with full hydrodynamic interactions to verify dynamical scaling. This is achieved by novel mesoscopic simulation techniques.

dynamic scaling, diffusion, polymers, mesoscopic simulation, hydrodynamics  preprint: Preprint submitted to …

The dynamics of polymer chains has attracted attention for decades already. In 3D, polymer dynamics exhibits rich and complex behavior which depends on the solvent conditions and polymer concentration degennes ; doi\_and\_edwards . The 2D case, however, has attracted much less attention. Recently, it has been realized that it has important applications in the field of colloids and biomolecules. Examples include the 2D diffusion of DNA oligonucleotides confined to interfaces maier99  and the lateral diffusion of lipids and proteins along biological interfaces cicuta01  such as cell membranes. Further, the dynamics of polymers in 2D is of major importance in thin films whose thickness is less than the size of the polymer. Wetting, surface adhesion, and flow in confined geometries are examples oron97  of this broad and fundamental field.

An important feature of essentially all the 2D diffusion processes in soft matter is that they take place in a solvated environment, which implies that the role of the hydrodynamic interaction (HI) cannot be disregarded. The HI originates from interactions mediated by the solvent in the presence of momentum conservation. Hydrodynamics plays a major role in a wide range of applications, and hence the understanding of its impact on the dynamics of polymer systems has received a great deal of attention. In 3D the effects of hydrodynamics are well understood: it is well known that the dynamics of polymers in dilute solution is well described by the Zimm model doi\_and\_edwards . In 2D, however, the situation becomes significantly more complicated as will be discussed below.

To understand the dynamics of polymer chains, with or without hydrodynamics, a common technique is to apply the theory of dynamical scaling doi\_and\_edwards . Two key quantities here are the radius of gyration $`R_g`$ and the center-of-mass (CM) diffusion coefficient $`D`$ of the chain. In the dilute limit, they follow the scaling relations $`R_gN^\nu `$ and $`DN^{\nu _D}`$, with corresponding scaling exponents $`\nu `$ and $`\nu _D`$, respectively. Another central quantity is the intermediate scattering function defined as

$$S(\stackrel{}{k},t)=(1/N)\underset{m,n}{}\mathrm{exp}\left\{i\stackrel{}{k}\left[\stackrel{}{r}_m(t)\stackrel{}{r}_n(0)\right]\right\},$$
(1)
where $`N`$ is the degree of polymerization, $`\stackrel{}{k}`$ is a wave vector, and $`\{\stackrel{}{r}_n\}`$’s are the positions of the monomers. This function should then scale as doi\_and\_edwards

$$S(k,t)=k^{1/\nu }F(tk^x),$$
(2)
where $`x`$ is the dynamical scaling exponent related to the other exponents through the relation

$$x=2+\nu _D/\nu .$$
(3)
This is valid for $`k(2\pi /R_g,2\pi /a)`$, where $`a`$ is the size of a monomer. Equations (2) and (3) are the cornerstones of dynamical scaling of polymers.

In the purely dissipative case, the values of the scaling exponents for polymer chains are well understood doi\_and\_edwards . In the dilute limit the simple Rouse model gives $`\nu =1/2`$ and $`\nu _D=1`$. When proper volume exclusion is taken into account, $`\nu =3/4`$ in 2D and approximately $`3/5`$ in 3D, while $`\nu _D=1`$ still holds for dilute 3D systems and for all polymer concentrations in 2D ala-nissila96 .

However, when the HI is taken into account, the situation becomes dramatically different. While in 3D theory and numerical simulations agree with the prediction of the Zimm equations that $`\nu =\nu _D`$ (i.e. $`x=3`$) dunweg91 ; dunweg93b ; pierleoni91 ; ahlrichs99 , in 2D the situation is less clear. What has been established both theoretically degennes ; doi\_and\_edwards  and computationally shannon97 ; vianney90  is that in good solvent conditions $`\nu =3/4`$ in 2D, as in the case of no HI.

The situation with $`\nu _D`$ is more subtle, however. Using lattice-gas simulations Vianney et al. vianney90  found a large positive value of $`\nu _D=0.78\pm 0.05`$. The molecular dynamics (MD) simulations of $`S(k,t)`$ by Shannon and Choy shannon97 , in turn, gave $`x=2`$ which would imply that $`\nu _D=0`$, if Eq. (3) holds. However, from their MD data for $`D`$ vs. $`N`$ they concluded that $`\nu _D>0`$ thus contradicting the scaling law. They also solved the Zimm equations numerically in 2D and verified the result that $`x=2`$, but found that now $`\nu _D<0`$ shannon97 . These results prompted the authors of Ref. shannon97  to suggest that dynamical scaling is broken for 2D polymers. Essentially, the very basis of polymer dynamics is being questioned.

In this letter, our objective is to determine the validity of dynamical scaling for 2D polymers. To this end, we first present analytic arguments which show that when finite-size effects are properly taken into account, the scaling of $`D`$ with respect to $`N`$ is truly logarithmic, leading to $`\nu _D=0`$ and thus to $`x=2`$. Following this, we carefully extract the exponents $`\nu `$, $`\nu _D`$, and $`x`$ through extensive mesoscopic simulations of a 2D polymer in a good solvent with the full HI included. Our results verify both that $`x=2`$ and the predicted logarithmic scaling of $`D`$, and thus we conclude that dynamical scaling is obeyed in 2D.

To overcome the significant difficulties in simulating polymers with full hydrodynamic interactions, we employ a novel mesoscopic simulation method introduced recently by Malevanets and Kapral (MK) malevanets99 ; malevanets00a . The MK method is essentially a hybrid molecular dynamics scheme, where the polymer chain is treated microscopically while the solvent obeys coarse-grained dynamics. In practice this idea is implemented by choosing the monomer-monomer and monomer-solvent interactions as in MD simulations, while the conservative interactions between the solvent particles are absent as in the ideal gas. This description preserves the hydrodynamic modes through so-called collision rules. Further, it allows for a major speedup compared to other simulation techniques such as MD.

To describe the dynamics of the coarse-grained solvent, time is partitioned into segments $`\tau `$ and the simulation box is divided into collision volumes or cells. The effective interactions between the solvent molecules take place at each $`\tau `$: this is called a collision event. In a collision the velocities of the solvent particles are transformed according to $`\stackrel{}{v}_i(t+\tau )=\stackrel{}{V}+\stackrel{}{\omega }\left[\stackrel{}{v}_i(t)\stackrel{}{V}\right]`$. Here $`\stackrel{}{v}_i`$ is the velocity of the particle $`i`$, $`\stackrel{}{V}`$ is the average velocity of all the particles in the cell the particle $`i`$ belongs to, and $`\stackrel{}{\omega }`$ is a random rotation matrix chosen for that particular cell. It can be shown malevanets99  that this multiparticle collision dynamics conserves the momentum and energy in each collision volume, and thus gives a correct description of the hydrodynamics of the velocity field.

Our model system consists of a polymer chain with $`N`$ monomers immersed in a 2D coarse-grained solvent. The mass of a solvent particle is set to $`m`$, and the monomer mass is $`2m`$. The monomer-monomer and monomer-solvent interactions are described by a truncated Lennard-Jones (LJ) potential:

$$U_{LJ}(r)=\{\begin{array}{cc}4ϵ\left[(\sigma /r)^{12}(\sigma /r)^6\right]+ϵ,\hfill & r2^{1/6}\sigma ;\hfill \\ 0,\hfill & r>2^{1/6}\sigma .\hfill \end{array}$$
(4)
Here $`\sigma `$ and $`ϵ`$ together with $`m`$ define the LJ unit system, where the unit of time is defined as $`\tau _{LJ}=\sigma \sqrt{m/ϵ}`$. In addition to the LJ potential, there is an attractive FENE potential between the nearest-neighbor monomers:

$$U_C(r)=(aR_0^2/2)\mathrm{ln}\left(1r^2/R_0^2\right),$$
(5)
where $`a=7ϵ\sigma ^2`$ and $`R_0=2\sigma `$.

The solvent density was set to $`\rho =0.581\sigma ^2`$, and the simulations were carried out at a temperature $`k_BT=1.2ϵ`$, yielding good solvent conditions. The equations of motion were integrated using the velocity Verlet algorithm with a time step $`\delta t=0.005\tau _{LJ}`$. The choice of the parameters that determine the collision dynamics fixes the properties of the coarse-grained solvent, e.g. its viscosity. Here we set the collision time to $`\tau =\tau _{LJ}`$ and the linear size of the collision volume to $`l_c=2\sigma `$. The random rotation angles were chosen from a uniform distribution in $`[0,2\pi )`$. The size of the polymer chain $`N`$ varies from 20 to 80 monomers, and the linear system size $`L`$ ranges from $`40\sigma `$ up to $`420\sigma `$. Periodic boundary conditions were employed for all system sizes. The CM diffusion coefficient $`D`$ was determined using the memory expansion method presented in Ref. ying98 .

First, we checked the scaling of $`R_g`$ with $`N[20,80]`$ footnote . Our estimate for the scaling exponent of the radius of gyration is $`\nu =0.75\pm 0.02`$, in excellent agreement with theory. Next, we computed the dynamic structure factor $`S(k,t)`$ which is depicted in Fig. 1.

Our data show the best collapse with $`x=2.0\pm 0.1`$, and we find that $`x`$ is not particularly sensitive to either $`N`$ or $`L`$ when $`N<L`$ (data not shown). This confirms the MD results of Shannon and Choy shannon97 , and shows that the 3D Zimm result $`x=3`$ is indeed invalid in 2D.

Next, we address the crucial question of the value of the scaling exponent $`\nu _D`$ for the CM diffusion coefficient. In the case of the corresponding system in 3D dunweg93b , the finite-size dependence of $`D`$ is given by $`D1/L`$. Therefore, in principle it is easy to determine $`D`$ for a fixed chain length $`N`$ by running a series of simulations for different values of $`L`$, and then extrapolating to $`L\mathrm{}`$. By repeating this procedure for several values of $`N`$, the exponent $`\nu _D`$ can be determined.

However, in the 2D case the finite-size effects are much more subtle due to the infinite range of the HI. We have carefully calculated $`D`$ analytically for a 2D polymer in a finite system of size $`L`$ using various approximations, including the approaches presented in Refs. kapral76 ; shannon97 ; dunweg93b . In all cases we find falck03a  that $`D`$ follows the scaling relation

$$DA\mathrm{ln}NB\mathrm{ln}\frac{1}{L},$$
(6)
where $`A`$ and $`B`$ are constants whose values depend on the approximations used falck03a ; footnote2 . This shows that extracting $`\nu _D`$ in the “traditional” sense in the thermodynamic limit $`L\mathrm{}`$ is no longer possible.

To numerically study the scaling of $`D`$, we determined $`D`$ for each $`N[20,80]`$ with different values of $`L`$. For instance, for $`N=30`$ we considered the cases $`L\{60,90,120,150,180,210,240\}`$. For every $`N`$, we examined the behavior of $`D`$ as a function of $`\mathrm{ln}(1/L)`$, and found that the behavior indeed is linear. Furthermore, we can estimate the exponent $`\nu _D`$ in terms of effective diffusion coefficients in the following way. We chose cutoff values $`L_{cut}\{10^2,10^3,10^4,10^5,10^6\}`$, and extrapolated a value $`D(N,L_{cut})`$ for each chain size and cutoff. If the data complied with Eq. (6), we should, when plotting $`D(N,L_{cut})`$ vs. $`\mathrm{ln}N`$, obtain a set of equally spaced straight lines. Each line corresponds to a certain cutoff, and the lines should all have the same slope $`A`$.

As can be seen in Fig. 2, this indeed holds within the statistical uncertainties of our data.

Most importantly, Fig. 2 confirms the prediction of logarithmic scaling of $`D`$ with $`N`$, which means that $`\nu _D=0`$. To quantify this, we can extract the exponent $`\nu _D`$ from $`\mathrm{ln}D(N,L_{cut})`$ vs. $`\mathrm{ln}N`$: for large values of $`L`$, we should have $`\mathrm{ln}D\nu _D\mathrm{ln}N`$. The results in Fig. 2 show that $`\nu _D`$ decreases steadily with $`L`$ as it should. For the largest $`L_{cut}`$ studied here, we find $`\nu _D0.05\pm 0.05`$.

The analysis above reveals the reason for the suggested breakdown of scaling in Refs. shannon97 ; vianney90 . While the result $`x=2`$ is correct, as verified here, the results in the previous studies for $`\nu _D`$ are simply incorrect because the exponent has been extracted without proper finite-size scaling analysis. Thus, we can conclude that dynamical scaling holds for 2D polymers with $`x=2`$, $`\nu =3/4`$, and $`\nu _D=0`$.

Finally, we wish to discuss the issue of long-time tails in 2D diffusion. The velocity autocorrelation function for a tagged particle $`i`$ is defined as $`\varphi (t)=\stackrel{}{v}_i(t+t^{})\stackrel{}{v}_i(t^{})`$ and its integral $`_0^{\mathrm{}}\mathrm{dt}\varphi (t)`$ yields $`D`$. In the presence of the HI it has been shown that $`\varphi (t)t^1`$ alder70 , or $`\varphi (t)[t\sqrt{\mathrm{ln}(t)}]^1`$ wainwright71 , which means that rigorously speaking, $`D`$ is not well defined in 2D. This would seem to invalidate the present scaling arguments. However, there are several ways to resolve this problem. A standard method is to view the diffusion coefficients in 2D as time-dependent quantities $`D(t)`$ wainwright71 , from which one can define effective values of $`DD(t_f)`$ at some finite time $`t_f`$. Alternatively, one can normalize $`D(t)`$’s with an appropriately chosen “bare” diffusion coefficient $`D_0(t)`$ (e.g. that of a monomer ala-nissila96 ). Such normalized values $`D(t)/D_0(t)`$ should converge to a finite result at long times.

However, there are many cases where it is in fact very difficult to observe long-time tails in either experiments vdhoef91  or in simulations. Further complications may arise from the fact that confined geometries can strongly influence the asymptotic decay of $`\varphi (t)`$ hagen97 . In the present case, the issue of long-time tails is settled by recognizing that in dynamical scaling, the absolute values of the diffusion coefficients are irrelevant: only the behavior of $`D`$ as a function of the chain or system size matters. Hence, we can here use such effective values provided that they have been determined in a consistent way. To this end, we have simply determined $`D`$’s over a time interval where the coefficients have converged within numerical error. More precisely, assuming $`D\mathrm{}_D^2/t_D`$, where $`\mathrm{}_D`$ is the distance over which the chain diffuses during the time interval $`t_D`$, the diffusion coefficients have been measured at a point where the chain has diffused a scaled distance $`\mathrm{}_D/R_g=24`$.

The convergence of $`D(t)`$ on these time scales is demonstrated in Fig. 3 for several systems with different values of $`N`$ and $`L`$.

In conclusion, we have applied the MK algorithm to a dilute 2D polymer solution. The MK method itself has proved to be an efficient tool for studies of macromolecular systems, especially in the dilute limit where the computational cost is mainly due to the explicit solvent. The technique has enabled us—at a moderate computational cost—to study system sizes that have not been previously amenable to simulations. This approach together with proper finite-size scaling analysis has allowed us to solve the controversy regarding the dynamical scaling of dilute polymer solutions in 2D with full hydrodynamics. We have found, in contrast to previous arguments shannon97 ; vianney90 , that the exponent relation $`x=2+\nu _D/\nu `$ is valid within numerical error. This justifies the scaling hypothesis, and shows that the anomalous exponent $`x=2`$ found in previous studies is due to the logarithmic scaling of $`D`$ as a function of $`N`$.

Acknowledgments –  This work has been supported in part by the Academy of Finland through its Center of Excellence Program and the National Graduate School in Materials Physics (E.F.).


text:
# A Novel Photonic Material for Designing Arbitrarily Shaped Waveguides in Two Dimensions

## Abstract

We investigate numerically optical properties of novel two-dimensional photonic materials where parallel dielectric rods are randomly placed with the restriction that the distance between rods is larger than a certain value. A large complete photonic gap (PG) is found when rods have sufficient density and dielectric contrast. Our result shows that neither long-range nor short-range order is an essential prerequisite to the formation of PGs. A universal principle is proposed for designing arbitrarily shaped waveguides, where waveguides are fenced with side walls of periodic rods and surrounded by the novel photonic materials. We observe highly efficient transmission of light for various waveguides. Due to structural uniformity, the novel photonic materials are best suited for filling up the outer region of waveguides of arbitrary shape and dimension comparable with the wavelength.

preprint: submitted to Phys. Rev. Lett.

Downsizing is an everlasting dream of researchers in engineering field. Researchers in the field of optics are hoping to find a way to fabricate all-optic integrated circuits by using optical elements comparable with the wavelength. In fact, realization of miniature-sized optical waveguides should soon be possible due to the discovery of photonic crystals (PhCs) JoannopoulosNature97 ; LinScience98 . Because of the periodicity of a dielectric constant, PhCs can be designed to have complete photonic gaps (PGs), range of frequencies for which light in any direction cannot propagate within the PhCs. We can steer light through the waveguides made of PhCs having complete PGs. On the other hand, the shapes of waveguides should be structurally commensurate with the periodicity of the host PhC. This severely obstructs the realization of arbitrarily shaped waveguides whose dimension is comparable with the wavelength.

Let us take two-dimensional waveguides of a PhC composed of periodic dielectric rods, for example. The waveguides are usually formed by removing rods along a line. Therefore, they are composed of a set of segments. Waveguides of $`90^{}`$ or $`60^{}`$ bends can be easily obtained from PhCs of square or triangular lattices. For a bend of an arbitrary angle, however, they become zigzag in shape and excess scattering occurs at the junctions of segments. Successive scattering significantly reduces the transmittance of waveguides composed of many branches and bends. Arbitrarily shaped waveguides, therefore, require photonic materials of maximum structurally uniformity in addition to the complete PGs. In this Letter, we propose novel photonic materials in which parallel dielectric rods are randomly placed in a certain region provided that the distance between the centers of rods is larger than a certain value $`D_{min}`$:

$$|𝐑_i𝐑_j|D_{min},$$
(1)
where $`𝐑_i`$ and $`𝐑_j`$ are the positions of $`i`$-th and $`j`$-th rod center. We call this new photonic material as uniformly distributed photonic scatterers (UDPS). It is noted that UDPS have neither long-range nor short-range order. Nevertheless, we show numerically that UDPS can have complete PGs if rods have sufficient density and dielectric contrast. We also propose a new concept to fabricate arbitrarily shaped waveguides, i. e., we fence the waveguides with side walls of periodic rods and fill up the outer region with UDPS. We observe clear propagation of waveguide modes with large transmittance.

Two examples of UDPS are shown in Figs. 1a and 1b composed of $`N=100`$ and $`N=200`$ rods of radius $`a`$, respectively. Here, we put $`D_{min}=4.0a`$. As the figures show, the distribution of rods becomes more uniform with increase in rod density. Transmittance of UDPS is calculated by assuming the incidence of the plane electric field of wavelength $`\lambda `$ from the upper side of Fig. 1a or 1b. Incident light is scattered multiply by each rod. This scattering is treated analytically by solving the Maxwell equationYousif88 . The solution gives the distribution of electric field and energy flow (Poynting vector). From the average energy flow at line L in Fig. 1a or 1b, we calculate the transmittance $`T`$ normalized by that without rods. $`T`$ is a function of normalized frequency $`\mathrm{\Omega }=2\pi a/\lambda `$ (known as size parameter) and becomes very small in the PG region. In some cases, $`T`$ exceeds unity because of the diffraction due to the finite size of UDPS.

Figure 1c shows values of $`T`$ for UDPS of $`N=100`$, 150, and 200 when the electric field $`𝐄`$ is parallel to the rod axis (TM mode). In all the figures of this Letter, we fix the dielectric constant of rods at $`\epsilon =12`$ which corresponds to that of Si at 1.55 $`\mu m`$ used world-wide in optical communications. We found no PG for UDPS of $`N=100`$ though there were two split dips at $`\mathrm{\Omega }=0.35`$. For UDPS of $`N=150`$, we observed a PG of $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }_c=30\%`$, where $`\mathrm{\Omega }_c`$ and $`\mathrm{\Delta }\mathrm{\Omega }`$ are the central frequency and width of the PG, respectively. Here, we define PG to be the frequency range continuously below $`T=0.01`$. This PG, however, is incomplete in that there appear spiky peaks in the gap region. At these spiky peaks we found an intrusion of energy flow through cracks (non-uniform regions) of UDPS. When $`N`$ becomes 200, PG grows up to $`37\%`$ with the suppression of spiky peaks. It was also verified that this PG is isotropic and therefore complete. We also calculated $`T`$ of UDPS for various $`\epsilon `$ with common values of $`D_{min}=4.0a`$ and area fraction of rods $`f=0.138`$, and found PG for $`\epsilon 5`$. The presence of PGs was also confirmed by results of finite difference time domain (FDTD) calculation. Note that the rod radius becomes 0.11 $`\mu m`$ if we correspond $`\mathrm{\Omega }_c`$ of $`N=200`$ to $`\lambda =1.55\mu m`$. This size can be prepared relatively easily using recently developed microfabrication technique.

The actual fabrication process inevitably involves certain fluctuation in rod position $`\mathrm{\Delta }x`$ and radius $`\mathrm{\Delta }a`$. It is natural to expect from the construction rule of UDPS that the PG is unaffected by $`\mathrm{\Delta }x`$. In contrast, the effect of $`\mathrm{\Delta }a`$ should be investigated. We have also plotted in Fig. 1c the transmittance of UDPS with $`\mathrm{\Delta }a/a=\pm 20\%`$. It was confirmed that PG of $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }_c=30\%`$ can survive such large fluctuation. This means that the PG of UDPS is also considerably robust against radius fluctuation.

We have observed in Fig. 1 that an increase in rod density enlarges the PG width. We therefore examined a case of much higher density. Figure 2a shows one example of UDPS of $`D_{min}=2.1a`$, which includes no radius fluctuation. Two examples of radial distribution functions $`g(r)`$ are plotted in Fig. 2b for $`D_{min}=2.1a`$, where $`r`$ is the distance between rod centers. For uniform distribution in two dimensions, $`g(r)`$ is proportional to $`r`$ without showing any peak. Distributions in Fig. 2b are very similar to the uniform case. Nevertheless, we find three distinct PGs at $`\mathrm{\Omega }_c`$=0.54, 0.93 and 1.34 for the TM mode in Fig. 2c. It is remarkable that the UDPS have such wide PGs of higher frequencies. If one uses the third PG, rod radius of $`0.33\mu m`$ is required to utilize $`\lambda =1.55\mu m`$. This facilitates the fabrication significantly. We also find a PG of TE mode ($`𝐄`$ rod axis) at $`\mathrm{\Omega }_c=0.68`$.

Before discussing the origin of PGs in UDPS, let us use UDPS for various waveguides. For this purpose we plot in Fig. 3 the average transmittance $`T`$ over five configurations of UDPS for three cases of sample thickness with common $`D_{min}`$. We can see that PG appears even for very thin UDPS containing three or four rods along the direction of light. This indicates that UDPS have wide applicability to build up waveguides of arbitrary shape and size comparable with the wavelength.

To make best use of this property, we first decide the shape of a waveguide. It can be twisty, as shown in Fig. 4, to enable maximum flexibility in designing. Waveguides are separated by side walls from the surrounding medium. To avoid excess scattering from the side walls, they are chosen to be made of periodic rods in a line. Then, we fill up the outer region with UDPS. It is noted that UDPS are the best materials for the surrounding optical medium. PhCs are not suitable for this purpose because their periodicity conflicts with that of side walls. This mismatch causes non-uniformity of rod density, resulting in excess scattering. This is also the case not only for quasi-periodic PhCs Chan98  but also for photonic materials having short-range order Jin01 .

In the waveguides shown below, the density of UDPS is the same as that in Fig. 1b ($`N=200`$) and radius fluctuation of $`\mathrm{\Delta }a/a=\pm 20\%`$ is introduced. Thus, they have a common PG of $`0.366\mathrm{\Omega }0.494`$ shown by the shaded region in each transmittance. We assume that the TM mode is incident from the upper side. Figure 4a shows the distribution of electric field intensity in a waveguide of $`90^{}`$ bend. Rods are shown by open circles. The intensity increases from blue to red. Energy flow is shown by white arrows. Corresponding frequency is indicated by the arrow in the transmittance in Fig. 4b which shows relatively large values over a wide range within PG. We can clearly observe the propagation mode of large transmittance comparable with that of waveguides made of PhCs. We also found no appreciable change in transmittance for a smaller value of $`\mathrm{\Delta }a`$. Therefore, waveguides of UDPS have a wide tolerance for the fabrication process.

UDPS are not limited to waveguides composed of segments. They can also be used for twisty waveguides whose curvatures are comparable with $`\lambda `$. Figure 4c shows such an example composed of two quarter circles. The distribution of electric field intensity and energy flow are shown at the frequency noted by the arrow in the transmittance in Fig. 4d. The rod density and incident light are the same as those in Fig. 4a. As can be seen in the figure, the electric field flows smoothly downward through the sample. While the corresponding value of $`T`$ is not large (0.739), it can be increased by optimization.

Let us discuss the origin of PGs in UDPS. In a study concerning effects of disorder on PG, it was found that there are two kinds of PG, one that is easily smeared out by disorder and the one that is very robust against disorder Lidorikis00 . The former PGs are formed by the coherent interference of scattered waves from periodic rods like Bragg diffraction in X-rays. The latter are formed by the bonding and anti-bonding states of Mie resonance states within each rod, similar to the electronic bandgaps in semiconductors. Since the latter are formed by local interaction, they are not significantly affected by the fluctuations in position and radius. As a matter of fact, an isolated dielectric rod of $`\epsilon =12`$ has Mie resonance at $`\mathrm{\Omega }`$=0.23, 0.66, 1.06, 1.14 and 1.44 for the TM mode and $`\mathrm{\Omega }`$=0.66, 1.05 and 1.40 for the TE mode. It is likely that the PG of UDPS is a result of interaction of these modes. This is also evidenced by the appearance of PG in Fig. 3 for very thin UDPS which is easily understood from the formation of bonding and anti-bonding states by local arrangement of rods. An important difference, however, exists between electrons and photons in that resonance wavefunctions of photons are not localized exponentially. Rather, they decay in inverse power and have a long-range nature. This long-range nature is responsible for the formation of PGs in UDPS that does not require even a short-range order.

One might think that UDPS is deeply related with two-dimensional disordered systems which are used to investigate Anderson localization of lightJohn90 . Let us discuss this point. In three dimensions, Anderson localization takes place only when the disorder is strong enough to satisfy the Ioffe-Regel criterion. In contrast, even a very small amount of disorder is sufficient in one and two dimensions to invoke Anderson localization. We have evaluated the localization length $`\mathrm{}_{loc}`$ in various UDPS samples of common $`D_{min}=4.0a`$ and different sample thickness $`\mathrm{}`$. We assume the form $`T=T_0\mathrm{exp}(\mathrm{}/\mathrm{}_{loc})`$, and found that within the PG region $`\mathrm{}_{loc}`$ is roughly $`2a`$ comparable with the surface distance $`\mathrm{}_{surf}`$ between nearest neighbor rods. Note generally that $`\mathrm{}_{loc}\mathrm{}_{mfp}\mathrm{}_{surf}`$ (usually $`\mathrm{}_{loc}\mathrm{}_{mfp}`$), where $`\mathrm{}_{mfp}`$ is the mean free path. Therefore, it would not be appropriate to conclude that Anderson localization can explain the observed localization length within the PGs. In fact, PGs have not been detected in two-dimensional disordered systems which show localization of light Localization . On the other hand, $`\mathrm{}_{loc}`$ above PGs is estimated to be $`30a40a`$ comparable with the sample thickness. The effect of localization is usually obvious when $`\mathrm{}_{loc}`$ is comparable with or less than the sample thickness $`\mathrm{}`$. Therefore, the thickness dependence of $`T`$ outside PG reflects the effect of Anderson localization.

There are numerous studies concerning the effect of disorder such as randomness of radii, positions or dielectric constants of rods on the PGs DisorderedPC  and waveguides made of PhCs Bayindir01 . PGs are observed when disorder is not so strong, but they are obviously vestiges of PGs of the underlying lattices. In contrast, there is no underlying lattice for UDPS and no peak is observed in radial distribution functions $`g(r)`$ as shown in Fig. 2b. The discovery of UDPS has three important contributions. Firstly, it gives a conceptual breakthrough of the common belief that either periodicity or short-range order is indispensable for the existence of PGs. Secondly, it enables us to fabricate very easily the wavelength-sized optical waveguides. Lastly, it unveils a new and powerful role of well-controlled randomness which can drastically change the optical feature of photonic materials. There has been no study, to the authors’ knowledge, to recognize such an active role of randomness in photonic materials. Future researches of air-hole type UDPS to make use of TE modes and better design policy to increase uniformity are needed for further development of UDPS.

This research was supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science, and Technology. One of the authors (H. M.) expresses his sincere thanks to R. Ohkawa for his continuous encouragement. We are grateful to K. Ohtaka for valuable discussions.