text:
# Entanglement Energetics at Zero Temperature

(14 November, 2003)

## Abstract

We show how many-body ground state entanglement information may be extracted from sub-system energy measurements at zero temperature. Generically, the larger the measured energy fluctuations are, the larger the entanglement is. Examples are given with the two-state system and the harmonic oscillator. Comparisons made with recent qubit experiments show this type of measurement provides another method to quantify entanglement with the environment.

A many-body quantum system is cooled to zero temperature so that it is forced into its overall nondegenerate ground state. We discuss the measurement of a sub-system Hamiltonian and demonstrate that it can be found in an excited state with a probability that depends on the coupling to its environment. This non-intuitive result is a pure quantum phenomenon: it is a consequence of entanglement sch  of the sub-system with the environment. In fact, we demonstrate that knowledge of the probability to find the system in an excited state can be used to determine the degree of entanglement of the sub-system and bath. Consequently, simple systems with well known isolated quantum mechanical properties (such as the two-state system and harmonic oscillator) become “entanglement-meters”.

There is growing interest of ground-state entanglement in condensed matter physics. Theoretical works on ground state entanglement have addressed entropy scaling in harmonic networks ms , spin-spin entanglement in quantum spin chains arnesen  and quantum phase transitions fazio ; nielsen2 . Entanglement properties of the ground state are also essential in the field of adiabatic quantum computing adiabatic . Furthermore it is interesting to link other ground state properties like the persistent current of small mesoscopic rings pc ; pc2 , or of doubly connected Cooper pair boxes cpbox1 ; cpbox2 ; qubit , or the occupation of resonant states resonant  to the zero-temperature entanglement energetics.

We consider a general Hamiltonian $`H=H_s+H_c+H_E`$, that couples $`(c)`$ the system $`(s)`$ we are interested in to a quantum environment $`(E)`$ such as a network of harmonic oscillators rmp . The lowest energy separable state is $`|\mathrm{\Omega }=|0_s|0_E`$, where $`|0_{\{s,E\}}`$ are the lowest uncoupled energy state of both systems. However, if the system Hamiltonian and the total Hamiltonian do not commute (which is the generic situation), then $`|\mathrm{\Omega }`$ is not an energy eigenstate of the total Hamiltonian. Thus, there must be a lower energy eigenstate ($`|0`$) of the total Hamiltonian which is by definition an entangled state. Because time evolution is governed by the full Hamiltonian, the ground state expectation of any operator with no explicit time dependence will have no time evolution, insuring that any measurement is static in time. This situation is in contrast to the usual starting point of assuming that the initial state is a separable state and studying how it becomes entangled. The reduced density operator of the system is given by tracing out the environmental degrees of freedom, $`\rho =\mathrm{Tr}_E|00|`$. Assuming the full state of the whole system is pure, the reduced density matrix contains all accessible system information, including entanglement of the system with its environment. Because repeated measurements of $`H_s`$ will give different energies as the sub-system is not in an energy eigenstate, we are interested in a complete description of the statistical energy fluctuations. These fluctuations may be described in two equivalent ways. The first way is to find the diagonal density matrix elements in the basis where $`H_s`$ is diagonal. These elements represent the probability to measure a particular excited state of $`H_s`$. A second way is to find all energy cumulants. A cumulant of arbitrary order may be calculated from the sub-system energy generating function, $`Z(\chi )=\mathrm{exp}(\chi H_s)`$ (as always, $`𝒪=\mathrm{Tr}\rho 𝒪`$) so that the nth energy cumulant is given by

$$H_s^n=()^n\frac{d^n}{d\chi ^n}\mathrm{ln}Z(\chi )|_{\chi =0}.$$
(1)
These cumulants give information about the measured energy distribution around the average.

Before proceeding to calculate these energy fluctuations, we ask a general question about entanglement. Given the energy distribution function (the diagonal matrix elements of the density matrix only), can anything be said in general about the purity or entropy of the state? Surprisingly, because we are given the additional information that we are at zero temperature, the answer is yes. If we ever measure the sub-system’s energy and find an excited energy, then we know the state is entangled. Though this statement alone links energy fluctuations with entanglement, a further quantitative statement may be made in the weak coupling limit. The reason for this is the following: the assumptions exponentially suppress higher states, so to first order in the coupling constant, we can consider a two-state system where the density matrix has the form $`\rho _{}=1\alpha p`$, $`\rho _{++}=\alpha p`$, $`\rho _+=\rho _+^{}=\alpha c`$. For vanishing coupling constant $`\alpha =0`$, this just gives the density matrix for the separable state. The linear dependence of $`\rho `$ on $`\alpha `$ holds to first order for the model systems considered below and is the entanglement contribution. If one measures the diagonal elements of $`\rho `$, one obtains $`\rho _{}`$ and $`\rho _{++}`$ as the probability to be measured in the ground or excited state (because $`\alpha `$ is small, there is only a small probability to find the sub-system in the upper state). If we now diagonalize $`\rho `$, the eigenvalues are $`\lambda _\pm =\{1p\alpha ,p\alpha \}+𝒪(\alpha ^2)`$. To first order in $`\alpha `$, the eigenvalues are the diagonal matrix elements, so we may (to a good approximation) write the purity or entropy in terms of these probabilities even if the energy difference remains unknown.

The Qubit. Let us now first evaluate the energy fluctuations of a qubit, a two-state system. The most general (trace 1) spin density matrix is $`\rho =(𝟙+\sigma _x\sigma _x+\sigma _y\sigma _y+\sigma _z\sigma _z)/2`$. A simple measure of the entanglement is given by the purity, $`\mathrm{Tr}\rho ^2=(1/2)(1+X^2+Y^2+Z^2)`$, where $`X_i=\sigma _i`$. It is well known that $`(X,Y,Z)`$ form coordinates in the Block sphere. Purity lies at the surface where $`X^2+Y^2+Z^2=1`$, whereas corruption lies deep in the middle.

We take the system Hamiltonian note1  to be $`H_s=(ϵ/2)\sigma _z+(\mathrm{\Delta }/2)\sigma _x`$. Introducing the frequency $`\mathrm{\Omega }=\sqrt{ϵ^2+\mathrm{\Delta }^2}/\mathrm{}`$ and using the identity $`e^{i\frac{\beta }{2}\widehat{n}\widehat{\sigma }}=I\mathrm{cos}\frac{\beta }{2}i\widehat{n}\widehat{\sigma }\mathrm{sin}\frac{\beta }{2}`$ with $`\beta =\chi \mathrm{\Omega }`$, it is straightforward to show

$$Z(i\chi )=\mathrm{cos}(\mathrm{\Omega }\chi /2)i\frac{\mathrm{sin}(\mathrm{\Omega }\chi /2)}{\mathrm{\Omega }}(ϵ\sigma _z+\mathrm{\Delta }\sigma _x).$$
(2)
The energy probability distribution may be easily found by Fourier transforming Eq. (2), or by tracing in the diagonal basis of the system Hamiltonian. The answer may be expressed with only the average energy, $`H_s=\frac{ϵ}{2}\sigma _z+\frac{\mathrm{\Delta }}{2}\sigma _x`$ as a sum of delta functions at the system energies $`\pm \mathrm{}\mathrm{\Omega }/2`$ with weights of the diagonal density matrix elements,

$$\rho _{++,}=(1/2)\left[1\pm \frac{H_s}{\mathrm{}\mathrm{\Omega }/2}\right].$$
(3)
Clearly, if the spin is isolated from the environment, $`H_s=\mathrm{}\mathrm{\Omega }/2`$, the ground state energy, the probability weight to be in an excited state vanishes. This distribution may also be found from knowledge of the isolated eigenenergies, the fact that $`H_s=_jE_j\rho _{jj}`$, and that Tr$`\rho =1`$. This later argument may be extended to $`n`$-state systems given the first $`n1`$ moments of the Hamiltonian and the $`n`$ eigenenergies.

Connection with Real Qubits. The probability weights depend on the energy parameters $`ϵ`$ and $`\mathrm{\Delta }`$, and the expectation values of the Pauli matrices. For real qubits produced in the lab, these will depend on the environment spin-boson . Often, we can link the basic phenomena we have been describing to physical measurements other than energy. For example, in a mesoscopic ring threaded by an Aharonov-Bohm flux $`\mathrm{\Phi }`$, or for the Cooper pair box, the tunneling matrix element depends on flux $`\mathrm{\Delta }(\mathrm{\Phi })`$. The operator one is interested in measuring is given by the normalized persistent current pc , or the expectation of $`\sigma _x`$ (for the symmetric case of $`ϵ=0`$),

$$\sigma _x=\frac{I(\mathrm{\Phi })}{I_0(\mathrm{\Phi })},$$
(4)
where $`I_0(\mathrm{\Phi })`$ is the uncoupled value of the persistent current. A common model for environmental effects is given by coupling the two-state system to a series of harmonic oscillators, the spin-boson model rmp ; pc ; weiss ; spin-boson . In Fig.(1) we have plotted the upper and lower occupation probabilities for the spin-boson model as a function of the coupling constant $`\alpha `$. For the symmetric case ($`ϵ=0`$), we have used the Bethe ansatz solution pc ; bethe , while for $`ϵ`$ finite, we have used the perturbative solution in $`\mathrm{\Delta }/\omega _c`$ which is valid only for larger $`\alpha `$ or $`ϵ`$ pc . Thus the plot is cut off at a small $`\alpha `$. A computational approach calculating the expectation values of the Pauli matrices over the whole parameter range was given in Ref. spin-boson . For $`ϵ=0`$, Eq. (4) determines also the expectation value of $`H_s`$, and thus the weights of the energy distribution function, Eq. (3), are directly related to the persistent current.

Experiments are always carried out at finite temperature, and it is important to demonstrate that there exists a cross-over temperature to the quantum behavior discussed here. In the low temperature limit, the thermal occupation probability is $`p_{up}=e^{(E_2E_1)/kT}`$. In the weak coupling limit for the symmetric spin boson problem, the probability to measure the excited state scales as $`p_{up}=\alpha \mathrm{log}(\mathrm{\Delta }/\omega _c)`$ note2 . Setting these factors equal and solving for $`T^{}`$ yields

$$kT^{}=\frac{E_2E_1}{\mathrm{log}(\alpha \mathrm{log}\frac{\omega _c}{\mathrm{\Delta }})}.$$
(5)
Since $`T^{}`$ scales as the inverse logarithm of the coupling constant, it is experimentally possible to reach a regime where thermal excitation is negligible.

As an order of magnitude estimate, we compare with the Cooper pair box cpbox1 ; cpbox2  which is among the most environmentally isolated solid state qubits qubit . From cpbox2  which found a $`Q10^4`$, we estimate the quantum probability for the box to be measured in the excited state as $`p_{up}10^310^4`$, which is of same order or larger than the thermal probability, $`p_{th}10^4`$. Experimentally, $`p_{up}`$ and $`p_{th}`$ may be confused by fitting data with an effective temperature, $`\rho _{th}=\mathrm{exp}(\beta _{\mathrm{eff}}H_s)`$ note3 . However, one may distinguish true thermal behavior from the effect described here because $`p_{up}`$ and $`p_{th}`$ depend differently on tunable system parameters such as $`\mathrm{\Delta }`$. In fact, $`\beta _{\mathrm{eff}}`$ is an entanglement measure. This chain of reasoning may be inverted to provide an estimate for $`Q`$ given only $`p_{up}`$.

The Harmonic Oscillator. We now consider the entanglement energetics of a harmonic oscillator, $`H_s=p^2/(2m)+(1/2)m\omega ^2q^2`$. Since there are an infinite number of states, the problem is harder. To simplify our task, we assume a linear coupling with a harmonic oscillator bath. This implies that the density matrix is Gaussian so that environmental information is contained in the second moments $`q^2`$ and $`p^2`$ weiss ; ms ,

$$q|\rho |q^{}=\frac{1}{\sqrt{2\pi q^2}}\mathrm{exp}\left\{\frac{(\frac{q+q^{}}{2})^2}{2q^2}\frac{p^2(qq^{})^2}{2\mathrm{}^2}\right\}.$$
(6)
Expectation values of higher powers of $`H_s`$ are non-trivial because $`q`$ and $`p`$ do not commute. The purity of the density matrix Eq. (6) is

$$\mathrm{Tr}\rho ^2=𝑑q𝑑q^{}q|\rho |q^{}q^{}|\rho |q=\frac{\mathrm{}/2}{\sqrt{q^2p^2}}.$$
(7)
The uncertainty relation, $`\sqrt{q^2p^2}\mathrm{}/2`$, guarantees that $`\mathrm{Tr}\rho ^21`$, with the inequality becoming sharp if the oscillator is isolated from the environment. As the environment causes greater deviation from the Planck scale limit, the state loses purity.

The generating function $`Z`$ may be calculated conveniently by tracing in the position basis and inserting a complete set of position states between the operators,

$$Z(\chi )=𝑑q𝑑q^{}q|\rho |q^{}q^{}|e^{\chi H_s}|q.$$
(8)
The first object in Eq. (8) is the density matrix in position representation, given by Eq. (6). The second object may be interpreted as the uncoupled position-space propagator of the harmonic oscillator from position $`q`$ to $`q^{}`$ in time $`i\mathrm{}\chi `$. We find

$$Z=\left\{2E\frac{\mathrm{sinh}\epsilon \chi }{\epsilon }+2A(\mathrm{cosh}\epsilon \chi 1)+\frac{1+\mathrm{cosh}\epsilon \chi }{2}\right\}^{\frac{1}{2}}$$
(9)
where $`\epsilon =\mathrm{}\omega `$, $`2E=m\omega ^2q^2+p^2/m`$ and $`A=q^2p^2/\mathrm{}^2`$. $`E`$ is the average energy of the oscillator, while $`A1`$ is a measure of satisfaction of the uncertainty principle. Eq. (9) has a pleasing limit for the free particle $`\omega 0`$,

$$Z(\chi )_{free}=\left\{1+\chi p^2/m\right\}^{\frac{1}{2}},$$
(10)
which is just the generating function for Wick contractions, $`p^{2n}=(2n1)!!(p^2)^n`$. Thus, in Eq. (9), the inverse square root generates the right combinatorial factors under differentiation, and the nontrivial $`\chi `$ dependence accounts for the commutation relations between $`q`$ and $`p`$. The first few harmonic oscillator energy cumulants may now be straightforwardly found via Eq. (1),

$`H_s^2`$ $`=`$ $`(1/2)[(\epsilon ^2/2)+4E^22\epsilon ^2A],`$ (11)
$`H_s^3`$ $`=`$ $`(E/2)[16E^2+\epsilon ^2(1+12A)],`$ (12)
$`H_s^4`$ $`=`$ $`48E^44\epsilon ^2E^2(1+12A)`$ (13)
$`+`$ $`\epsilon ^4[(1/8)+2A+6A^2].`$
After inserting the mean square values for an ohmic bath (see the discussion above eqs. (15,16)), Eq. (11) is identical to the main result of Ref. nb .

Alternatively, we now consider the diagonal matrix elements $`\rho _{nn}`$. An analytical expression for the density matrix in the energy basis may be found by using the wavefunctions of the harmonic oscillator, $`\psi _n(q)e^{\gamma ^2q^2/2}H_n(\gamma q)`$ where $`\gamma =\sqrt{m\omega /\mathrm{}}`$ and $`H_n(x)`$ is the $`n^{th}`$ Hermite polynomial. In the energy basis, the density matrix is given by $`\rho _{nm}=𝑑q𝑑q^{}\psi _n^{}(q)q|\rho |q^{}\psi _m(q^{})`$. The position space integrals may be done using two different copies of the generating function for the Hermite polynomials. The diagonal elements may be found by equating equal powers of the generating variables. We first define the dimensionless variables $`x=2\gamma ^2q^2`$, $`y=2p^2/(\gamma ^2\mathrm{}^2)`$, and $`D=1+x+y+xy`$. $`x`$ and $`y`$ are related to the major and minor axes of an uncertainty ellipse. The isolated harmonic oscillator (in it’s ground state) obeys two important properties: minimum uncertainty (in position and momentum) and equipartition of energy between average kinetic and potential energies. The influence of the environment causes deviations from these ideal behaviors which may be accounted for by introducing two new parameters, $`a=(yx)/D,b=(xy1)/D`$ with $`1a1`$ and $`0b1`$. The deviation from equipartition of energy is measured by $`a`$, while the deviation from the ideal uncertainty relation is measured by $`b`$. We find

$$\rho _{nn}=\frac{n!}{2^n}\sqrt{\frac{4}{D}}\underset{m=0}{\overset{n}{}}\underset{p=0}{\overset{n}{}}\frac{(a^2)^p(2b)^m}{m!(p!)^2}\delta _{n,m+2p}.$$
(14)
All terms are positive (because $`a`$ appears only squared), less than one, and obey $`\mathrm{Tr}\rho =1`$. The summation over $`m`$ and $`p`$ define a polynomial $`P_n`$ of order $`n`$ in $`a`$ and $`b`$. The probability for the lone oscillator to be measured in an excited state clearly decays rapidly with level number. These probabilities also reveal environmental information. For example, $`P_1=2b`$ and is thus only sensitive to the area of the state, while $`P_2=a^2+2b^2`$ depends on both the uncertainty and energy asymmetry. Additionally, if we expand the first density matrix eigenvalue ms ; weiss  with respect to small deviations of $`x`$ and $`y`$, we recover $`\rho _{11}`$ in agreement with our general argument.

Although $`x`$ and $`y`$ have been treated as independent variables, the kind of environment the system is coupled to replaces these variables with two functions of the coupling constant. For example, with the ohmic bath weiss ; nb  (in the under-damped limit), the variables are

$`x(\alpha )={\displaystyle \frac{1}{\sqrt{1\alpha ^2}}}\left(1{\displaystyle \frac{2}{\pi }}\mathrm{arctan}{\displaystyle \frac{\alpha }{\sqrt{1\alpha ^2}}}\right),`$ (15)
$`y(\alpha )=(12\alpha ^2)x(\alpha )+{\displaystyle \frac{4\alpha }{\pi }}\mathrm{ln}{\displaystyle \frac{\omega _c}{\omega }},`$ (16)
where $`\alpha `$ is the coupling to the environment in units of the oscillator frequency and $`\omega _c`$ is a high frequency cutoff. This bath information is shown in Fig. (2) with $`\omega _c=10\omega `$. The trajectory of the line over the surface shows how the probabilities evolve as the coupling $`\alpha `$ is increased from 0 to 1. Other kinds of environments would trace out different contours on the probability surface.

In conclusion, we have shown that projective measurements of the system Hamiltonian at zero temperature reveals entanglement properties of the many-body quantum mechanical ground state. Consequently, repeated experiments on simple quantum systems give information about the nature of the environment, the strength of the coupling and entanglement. The larger the energy fluctuations, the greater the entanglement. There are several possibilities for experimental implementations. We have mentioned measurement of persistent current pc ; pc2  as well as projecting on the system’s energy eigenstates. Another measurement possibility is a zero temperature activation-like process act  where the dominant mechanism is not tunneling, but the same quantum effects of the environment which we have discussed here.

This work was supported by the Swiss National Science Foundation.


format: text

text:
# Mapping the Cosmic Web with Ly𝛼 Emission

## 1. Introduction

One of the most striking images to emerge from cosmological simulations is the “cosmic web”. Matter collapses into a web of moderately overdense sheets and filaments, with galaxies and galaxy clusters forming through continued collapse at their intersections. This picture explains both the qualitative distribution of galaxies in redshift surveys (e.g., de Lapparent et al. 1986) and many of the characteristics of the Ly$`\alpha `$ forest (see e.g. Rauch 1998). However, mapping the cosmic web remains a difficult venture. Ly$`\alpha `$ forest spectra probe only one-dimensional lines of sight, while galaxy redshift surveys find only the small fraction of baryons embedded inside galaxies, offering an indirect picture of the gaseous filaments.

In this *Letter*, we use a cosmological simulation to show that surveys of hydrogen Ly$`\alpha `$ ($`\lambda `$1216 Å) emission can provide a powerful method to map the cosmic web. Existing analytic estimates suggest that the surface brightness of dense portions of the IGM can be substantial (Hogan & Weymann, 1987; Gould & Weinberg, 1996). Such regions generally lie within filaments but outside of galaxies, so Ly$`\alpha `$ emission offers a more direct map of the gas distribution in filaments than do galaxy surveys. In fact, Ly$`\alpha `$ emission from galaxies has been used to detect a filament at high redshift (Möller & Fynbo, 2001); here we focus on intergalactic emission from low redshifts, where the distances and sky background are relatively small and where Ly$`\alpha `$ maps of the cosmic web can be compared to existing galaxy surveys. We also show that Ly$`\alpha `$-emitting gas constitutes an intrinsically interesting phase: gas that is cooling around collapsed halos (Haiman et al., 2000; Fardal et al., 2001). The distribution of Ly$`\alpha `$ emission therefore provides a window onto the growth of bound objects.

## 2. Simulation and Analysis

We perform our analysis using the G5 cosmological simulation of Springel & Hernquist (2003b). This smoothed-particle hydrodynamics simulation included a multiphase description of star formation which incorporates a prescription for galactic winds (Springel & Hernquist, 2003a). The simulation assumed a $`\mathrm{\Lambda }`$CDM cosmology consistent with the most recent cosmological observations (e.g., Spergel et al. 2003). It has a box size of $`100h^1`$ comoving Mpc (where the Hubble constant $`H_0=100h\text{ km s}\text{-1}\text{ Mpc}\text{-1}`$ and $`h=0.7`$ in the simulation), a spatial resolution of $`8h^1`$ comoving kpc, and a baryonic particle mass of $`3.26\times 10^8h^1\text{ M}\text{}`$.

Furlanetto et al. (2003, hereafter F03), who focus on UV emission from metals in the IGM, describe the main components of our analysis procedure; we briefly summarize them here. We first compute two grids of the Ly$`\alpha `$ emissivity $`\epsilon _\alpha `$ as a function of hydrogen density $`n_\mathrm{H}`$ and gas temperature $`T`$ using CLOUDY 96 (beta 5, Ferland 2003), assuming ionization equilibrium. For the first grid ($`\epsilon _\alpha ^{\mathrm{pi}}`$), we apply the Haardt & Madau (2001) ionizing background, which includes radiation from galaxies and quasars (our main results are insensitive to the background choice; F03). For the second grid ($`\epsilon _\alpha ^{\mathrm{coll}}`$), we neglect all photoionization processes; we use this grid to calculate the emissivity of self-shielded gas.111Note that the electron density (and hence $`\epsilon _\alpha ^{\mathrm{coll}}`$) becomes sensitive to the metallicity for $`T<10^4\text{ K}`$. However, this has negligible effects on our results because at these temperatures the emissivity is vanishingly small. We compare the two grids in Figure 1. The solid curve shows $`\epsilon _\alpha ^{\mathrm{coll}}`$ (note that $`\epsilon _\alpha n_\mathrm{H}^2`$ when collisional processes dominate), and the other curves show $`\epsilon _\alpha ^{\mathrm{pi}}`$ for several different densities. We see that $`\epsilon _\alpha ^{\mathrm{coll}}`$ drops rapidly when the gas recombines at $`T10^{4.2}\text{ K}`$. In contrast, $`\epsilon _\alpha ^{\mathrm{pi}}`$ is still large at low temperatures because the gas remains photoionized. As the density and temperature increase, collisional processes become dominant and $`\epsilon _\alpha ^{\mathrm{pi}}\epsilon _\alpha ^{\mathrm{coll}}`$.

We wish to compute the emission from a slice of the simulation with fixed $`z`$, $`\mathrm{\Delta }z`$, and angular size. We first randomly choose the volume corresponding to this slice (see F03). We classify particles inside the slice as “optically thin” or “self-shielded.” The latter set includes those dense and cool particles with a large enough neutral fraction to shield themselves from the ionizing background. We approximate self-shielding by assuming that all gas with $`n_\mathrm{H}>n_{\mathrm{ss}}`$ and $`T<T_{\mathrm{ss}}`$ is optically thick. We set fiducial thresholds at $`n_{\mathrm{ss}}=10^3\text{ cm}\text{-3}`$ and $`T_{\mathrm{ss}}=10^{4.5}\text{ K}`$; radiative transfer calculations for self-gravitating clouds like those described in Schaye (2001) find that self-shielding becomes important at about this density, assuming that it remains cool (see also Katz et al. 1996). As shown in F03, most particles in this regime lie on a well-defined curve in the $`n_\mathrm{H}`$$`T`$ plane, beginning at moderate overdensity and $`T10^{4.25}\text{ K}`$ and asymptotically approaching $`10^4\text{ K}`$ as the density increases to $`n_\mathrm{H}0.1\text{ cm}\text{-3}`$. This locus represents gas that is cooling and falling into the centers of halos. A small but important fraction of particles (those that have recently been heated by shocks) have slightly higher temperatures. Once we have made this classification, we assign emissivities to optically thin particles using the $`\epsilon _\alpha ^{\mathrm{pi}}`$ grid and (in our fiducial model) to self-shielded particles using the $`\epsilon _\alpha ^{\mathrm{coll}}`$ grid.222Note that this is conservative in the sense that we neglect photoionization by hard UV/X-ray background photons to which the gas is still optically thin. Unfortunately, the emissivity of self-shielded particles is rather uncertain because the simulation does not include radiative transfer, metal line cooling, or local ionizing sources and because we assume ionization equilibrium. For this reason, we examine their emission more closely in §3.2.

Active star formation also produces substantial Ly$`\alpha `$ emission when gas in the host galaxy absorbs ionizing radiation from stars and subsequently recombines. In the simulation, star formation occurs in any gas particle whose density exceeds $`n_\mathrm{H}=0.129\text{ cm}\text{-3}`$ (Springel & Hernquist, 2003a). For each such particle, we convert the star formation rate (SFR) reported by the simulation to a Ly$`\alpha `$ luminosity via $`L_\alpha =10^{42}(\mathrm{SFR}/\mathrm{M}_{}\mathrm{yr}^1)\text{ erg s}\text{-1}`$, which is accurate to within a factor of a few for metallicities between $`0.05\text{ Z}\text{}<Z<2\text{ Z}\text{}`$ (Leitherer et al., 1999), assuming that $`2/3`$ of ionizing photons are converted to Ly$`\alpha `$ photons (Osterbrock, 1989). The assumption of case-B recombination should be valid provided that most ionizing photons are absorbed in the dense interstellar medium of the galaxy. The actual Ly$`\alpha `$ luminosity of a galaxy depends strongly on the distribution of ionizing sources, the escape fraction of ionizing photons, the presence of dust, and the kinematic structure of the gas (e.g., Kunth et al. 2003), so this should be taken as no more than a representative estimate. Note that any star formation unresolved by the simulation could significantly change our results. However, Springel & Hernquist (2003b) show that the global SFR in the simulation matches observational estimates well, so such difficulties are probably not severe.

Finally, we locate each particle on a pixelized map, convert its emissivity to surface brightness, and smooth the result with a Gaussian filter of full width at half maximum $`\mathrm{\Delta }\theta `$ (F03).

## 3. Results

Figure 2*a* shows a thin slice ($`\mathrm{\Delta }z=10^3`$) of our simulated universe at $`z=0.15`$. The panel shows a large volume with coarse resolution; the colorscale includes the full dynamic range of surface brightness $`\mathrm{\Phi }`$. Note how filaments stand out clearly in the map. Figure 2*b* shows part of a filament with a finer resolution ($`13h^1`$ physical kpc). We include only those pixels with $`\mathrm{\Phi }>10\text{ photons cm}\text{-2}\text{ s}\text{-1}\text{ sr}\text{-1}`$; this is a factor $`10`$ smaller than some observational (Brown et al., 2000) and theoretical (Haardt & Madau, 2001) estimates of the diffuse background at $`1200`$$`2000`$ Å. While the mean emission from filaments lies well below the background, we find that numerous discrete regions trace their course: the space density of these Ly$`\alpha `$ emitters is $`0.1h^3\text{ Mpc}^3`$. Wide-field Ly$`\alpha `$ observations are thus a promising technique to map out the cosmic web. Note that decreasing the threshold to $`\mathrm{\Phi }>1\text{ photons cm}\text{-2}\text{ s}\text{-1}\text{ sr}\text{-1}`$ enlarges the coronae but does not significantly increase their number density.

Ly$`\alpha `$ emission occurs over a wide range of physical scales: generally, emitting regions have a central, high surface brightness ($`\mathrm{\Phi }10^3\text{ photons cm}\text{-2}\text{ s}\text{-1}\text{ sr}\text{-1}`$) core of size $`30`$$`50h^1\text{ kpc}`$ surrounded by an irregular, low surface brightness “corona” that can extend to scales $`10^2\text{ kpc}`$. Thus, the Ly$`\alpha `$ emission is much more extended than the galaxies, indicating that this view of the cosmic web traces truly intergalactic gas.

Figure 3 shows histograms of the pixel flux probability distribution function (PDF) $`\mathrm{d}P/\mathrm{d}\mathrm{log}\mathrm{\Phi }`$. The solid curve in each panel shows our fiducial model (using $`\epsilon _\alpha ^{\mathrm{coll}}`$ for self-shielded gas). Figure 3*a* shows the PDF to very small $`\mathrm{\Phi }`$; we see a broad peak centered at $`\mathrm{\Phi }10^3\text{ photons cm}\text{-2}\text{ s}\text{-1}\text{ sr}\text{-1}`$ and a tail extending to much higher surface brightness. As described in Paper I, the location of the peak is determined simply by ionization equilibrium in the mean-density IGM and can be estimated analytically (Gould & Weinberg, 1996). Note that the PDF depends strongly on $`\mathrm{\Delta }\theta `$, because the emitting regions have structure on small scales.

Figure 3*b* shows that the PDF evolves very little from $`z=0.07`$ to $`z=0.45`$. The fraction of bright pixels increases slowly with redshift because both the density and temperature of the IGM increase with redshift.

### 3.1. Ly$`\alpha `$ Photons From Star Formation

Figure 2*d* shows the Ly$`\alpha `$ emission from star formation. While all star-forming regions are surrounded by intergalactic Ly$`\alpha `$ emission, the converse is not true. We do find that Ly$`\alpha `$ emission nearly always surrounds star particles, but often these are dwarf galaxies or lack ongoing star formation. The exceptions are galaxies inside hot clusters, which do not always host Ly$`\alpha `$ emission. Note that Ly$`\alpha `$ emission produced through star formation is much more compact than the coronae (partly because we assume that all the ionizing photons are absorbed locally). Figure 3*b* compares the mechanisms in a more quantitative form. The emission from star-forming gas in the simulation (dot-dashed curve) dominates the brightest pixels but is unimportant for $`\mathrm{\Phi }10^5\text{ photons cm}\text{-2}\text{ s}\text{-1}\text{ sr}\text{-1}`$. If a large fraction of the ionizing photons escape their host galaxy, the number of pixels affected by local star-formation will increase, but their mean brightness will also decrease. This could be particularly important near active galactic nuclei (AGN), which have large ionizing fluxes and large escape fractions. AGN could “light up” large volumes around their host galaxies in Ly$`\alpha `$ emission (Haiman & Rees, 2001). We note that the SFR, and hence Ly$`\alpha `$ emission from star-forming gas, increases with redshift (Hernquist & Springel, 2003).

### 3.2. Self-Shielded Gas

The largest uncertainty in our calculation is probably the treatment of self-shielded gas. One problem is identifying which particles are shielded from the ionizing background; we find, however, that our results do not change significantly in the range $`10^4\text{ cm}\text{-3}<n_{\mathrm{ss}}<10^2\text{ cm}\text{-3}`$, because most gas of these densities has $`T10^{4.1}\text{ K}`$, where $`\epsilon _\alpha ^{\mathrm{pi}}\epsilon _\alpha ^{\mathrm{coll}}`$ (see Fig. 1). A more important concern is that we may have overestimated the temperature of self-shielded particles because the simulation erroneously allows the ionizing background to heat such particles and underestimates the cooling rates in optically thick regions. Radiative transfer calculations with CLOUDY suggest that the simulation may have overestimated the temperature of particles with $`n_\mathrm{H}10^3\text{ cm}\text{-3}`$ by $`25\%`$, which can be significant given the steep dependence seen in Figure 1. The error increases as the gas density and metallicity increase. In a worst case scenario, *all* self-shielded gas would have $`T<10^{4.1}\text{ K}`$, producing essentially no emission. In Figure 2*c*, we approximate this case by setting $`\epsilon _\alpha =0`$ for all self-shielded gas. Comparing to the fiducial model, we find that the core brightness declines dramatically, but the extended coronae are unaffected. We conclude that the large-scale emission is due to optically thin gas while the brightest regions consist of gas at or near to the self-shielding threshold. The resulting PDF is shown as the short-dashed curve in Figure 3*a*; it essentially eliminates the bright tail of the distribution. Figure 3*c* shows that our results *do* depend on $`n_{\mathrm{ss}}`$ if self-shielded gas has $`\epsilon _\alpha =0`$. The PDF changes little for $`10^4\text{ cm}\text{-3}<n_{\mathrm{ss}}<10^3\text{ cm}\text{-3}`$, because such gas has relatively low density and hence small emissivity. However, it approaches the fiducial model if we set $`n_{\mathrm{ss}}=10^2\text{ cm}\text{-3}`$, because gas with $`n_\mathrm{H}10^2\text{ cm}\text{-3}`$ has $`T10^{4.0}\text{ K}`$ in the simulation and hence a small $`\epsilon _\alpha ^{\mathrm{coll}}`$.

On the other hand, the simulation also neglects local ionizing sources. We have already seen that most self-shielded particles reside near galaxies; we would therefore expect the radiation field around these particles to exceed the background. In this case the emissivity is limited by the ionizing flux reaching the gas. We have estimated the number of ionizing photons from young stars in §3.1, although we did not attempt to model the spatial distribution of the resulting emission. In addition, dense gas may no longer be self-shielded, so the diffuse background will also contribute to the Ly$`\alpha `$ emission. We estimate this extra component by assuming that no particles are self-shielded (dotted curve in Figure 3*a*). The number of bright pixels increases because $`\epsilon _\alpha ^{\mathrm{pi}}\epsilon _\alpha ^{\mathrm{coll}}`$ for $`T<10^{4.1}\text{ K}`$. In either case, the net effect of local ionizing sources is to *increase* the Ly$`\alpha `$ emission.

## 4. Discussion

We study Ly$`\alpha `$ emission from gas in the IGM using a high-resolution cosmological simulation. We find that galaxies are typically surrounded by “coronae” of Ly$`\alpha `$ emission with bright, central cores (of diameter $`40h^1\text{ kpc}`$) surrounded by much larger ($`10^2h^1\text{ kpc}`$) regions with surface brightness near the background. While the number of Ly$`\alpha `$ photons due to star formation in the host galaxies can exceed that emitted by the surrounding gas, the spatial extent of the star formation is typically much smaller. High angular-resolution observations (probing physical scales $`15h^1\text{ kpc}`$) should be able to separate the two components. Detection of the bright cores is feasible with existing technology, although the wide-field ultraviolet spectrographs most useful for these studies have not yet been built. The *Galaxy Evolution Explorer* (*GALEX*), with a large field of view but relatively low spectral resolution ($`10`$ Å), may detect the brightest cores. Furthermore, the strong correlation between Ly$`\alpha `$ emission and galaxies may enable a statistical detection of the signal in deep exposures. A project underway to construct a balloon-borne wide-field UV spectrograph will have a limiting sensitivity $`\mathrm{\Phi }500\text{ photons cm}\text{-2}\text{ s}\text{-1}\text{ sr}\text{-1}`$ in a single night’s observation (D. Schiminovich, private communication).

We emphasize that the emission from self-shielded gas is uncertain, because the simulation does not include all of the relevant physics. We have considered a range of models bracketing these uncertainties (see §3.2): in our fiducial model, dense, cool gas produces Ly$`\alpha `$ photons only through collisional processes, but we also consider cases in which it has zero emissivity or in which we allow photoionization. The choice has little effect on the large-scale, low surface brightness emission, which comes primarily from optically thin gas, but it does strongly affect the luminous cores., where most of the gas is shielded from the ionizing background.

We have also neglected dust, which can efficiently destroy Ly$`\alpha `$ photons. Because the emitting region is often considerably larger than the associated galaxy, much of the emitting gas may be relatively pristine. In the simulation, we find that the majority of particles with large Ly$`\alpha `$ emissivity contain no metals (or dust), though note that the simulation does not include mixing between simulation particles. At higher redshifts, several large-scale Ly$`\alpha `$-emitting “blobs” have already been observed (e.g., Steidel et al. 2000), limiting the amount of dust in these environments.

The structure and relatively large spatial size of Ly$`\alpha `$ emitters, together with their locations, suggest that they are “cooling flows” of gas onto halos. Detailed, high resolution Ly$`\alpha `$ observations therefore offer a means to study the accretion and cooling of gas onto galaxies. Galactic winds will also affect the distribution of Ly$`\alpha `$ emission. However, we find that the winds in our simulations have only minor effects on the statistics of Ly$`\alpha `$ emitters.

On larger scales, the Ly$`\alpha `$ emitters are distributed along sheets and filaments. Because these objects are numerous (with a space density $`0.1h^3`$ Mpc-3), Ly$`\alpha `$ line observations offer a new and powerful way to map the cosmic web. Figure 2 shows that such observations do not require high resolution; wide-field observations could efficiently locate (in three dimensions) a large number of Ly$`\alpha `$ emitters. In contrast to galaxy surveys, this approach reveals the location of baryons *outside* of galaxies. Because Ly$`\alpha `$ emission traces the cool gas accreting onto galaxies, it also selects a different population of dark matter halos, including many with few (or even no) associated stars. Follow-up observations at high resolution would allow us to compare Ly$`\alpha `$ emitters to galaxy surveys and learn how galaxies accrete gas at the present day.

We thank D. Schiminovich, C. Martin, and G. Ferland for helpful discussions and B. Robertson for assistance with the computational facilities. This work was supported in part by The W.M. Keck foundation, NSF grants PHY-0070928 and AST 00-71019 and NASA ATP grant NAG5-12140. The simulations were performed at the Center for Parallel Astrophysical Computing at the Harvard-Smithsonian Center for Astrophysics.


text:
# The Optimal Faucet

## Abstract

The production of small fluid droplets relies on an instability of solutions to the Young-Laplace equation. We ask whether smaller droplets can be produced by changing the shape of the nozzle. At a given critical pressure, the circular nozzle actually produces the largest droplet. The droplet volume can be decreased by up to 18% using a triangular nozzle with stretched corners.

A standard protocol for producing small droplets is as follows: a pipette, of circular cross-section, is pressurized at one end, pushing out a small fluid droplet. If the nozzle is sufficiently small, force balance requires that the droplet has constant mean curvature. At a critical pressure, this equilibrium shape becomes unstable, ultimately leading to the droplet detaching from the nozzle.

The volume of fluid entrained during this process is set by the total fluid volume contained in the critical droplet. This volume scales like $`r^3`$, where $`r`$ is the nozzle radius. On the other hand, the critical pressure for ejecting this droplet scales like $`\gamma /r`$, where $`\gamma `$ is the liquid surface tension. Thus, ejecting smaller droplets requires higher pressures. The smallest size droplet that can be ejected is thus determined by the highest pressure that can be reliably applied to the nozzle, without material failure, etc.

One strategy for creating smaller droplets than those dictated by the instability of a static droplet is to use a time varying forcing at the nozzle. This method has achieved an order of magnitude decrease in droplet volume Chen and Basaran (2002).

However, typical nozzles use a circular cross section. It is not unreasonable to imagine that changing the shape of the cross section to be some other shape may decrease the ejected droplet volume, while maintaining the same applied pressure. For example, imagine that we have a circular nozzle with a pendant droplet just below the critical volume: by “squeezing” the shape of the nozzle cross section into an elliptical shape, one might cause the droplet to detach at a lower volume.

In this paper we address the question: what is the shape of a nozzle for which the ejected droplet volume is minimized, for a given applied pressure? We demonstrate that circular nozzles do not eject the smallest droplets; instead, the optimal nozzle more closely resembles an equilateral triangle, albeit with “stretched” corners. The best nozzle shape that we have found has an ejected droplet volume about twenty percent smaller than the circular nozzle with the same critical pressure. Our method is inspired by and extends J. Keller’s classic treatment of the Euler buckling problem with a beam of nonuniform cross section Keller (1960). Recently, the method has been applied to the optimization of a bistable switch Brenner et al. (2003). For a detailed mathematical treatment of capillary surfaces in general, see Finn (1986).

This Letter is organized as follows. We first explain the origin of the pendant droplet instability. Then we describe our method for reducing droplet size. Lastly we provide numerical calculations implementing the method, and present the candidate optimal nozzle.

Pendant Droplet Instability. \- The instability of a droplet protruding from a nozzle is due to a bifurcation, most easily seen in the case of a circular nozzle that is much smaller than the capillary length, which allows us to neglect gravity. The shape of the droplet is then determined by the Young-Laplace equation $`p=\gamma K`$, where $`p`$ is the pressure difference across the liquid/air interface, $`\gamma `$ is the surface tension, and $`K`$ is the mean curvature of the droplet surface. This equation describes a surface of constant mean curvature $`p/\gamma `$ with the nozzle edge as its boundary. If the boundary is a circle, then the solution must be a section of the sphere with mean curvature $`p/\gamma `$. From the familiar relation

$$K_{\text{sphere}}=\frac{2}{\text{sphere radius}}$$
(1)
we deduce that the radius of curvature of the droplet is $`2\gamma /p`$. For small $`p`$, such that the sphere radius is much greater than the nozzle radius, the solution is a shallow spherical cap. But note that its complement, the rest of the sphere, is also a solution. As $`p`$ is increased, these two solutions approach each other until both become a hemisphere with the nozzle at the equator. The pressure at which the two solutions meet is the critical pressure $`p^{}`$, and the corresponding degenerate solution is unstable. Note that the critical pressure is also the maximum pressure, for the nozzle cannot support a sphere smaller than itself.

For a noncircular nozzle, we no longer have such a simple geometric picture, however key features remain. The unstable solution is still characterized by a bifurcation at which two solutions meet, corresponding to the maximum pressure achievable for the given nozzle. The critical pressure for a general nozzle can be computed as follows: let the droplet surface be parameterized as a function $`\stackrel{}{R}(u,v)`$ over a domain $`D`$ in the $`uv`$-plane, which takes value in three dimensional physical space. The boundary of the domain $`D`$ corresponds to a closed curve $`C`$ which represents the nozzle. The curvature is a nonlinear functional of the surface and its derivatives up to second order, hence the equation for the droplet shape has the form

$$\gamma K[\stackrel{}{R},\stackrel{}{}\stackrel{}{R},\stackrel{}{}\stackrel{}{}\stackrel{}{R}]=p,$$
(2)
where $`\stackrel{}{}`$ is the gradient operator in the $`uv`$-plane.

Upon increasing the pressure $`pp+\delta p`$, the surface changes: $`\stackrel{}{R}\stackrel{}{R}+\delta \stackrel{}{R}`$. Equation (2) implies that the variation $`\delta \stackrel{}{R}`$ and $`\delta p`$ are related by

$$\gamma \widehat{L}\delta \stackrel{}{R}=\delta p,$$
(3)
where $`\widehat{L}\delta \stackrel{}{R}`$ is the change in mean curvature induced by the surface change. $`\widehat{L}`$ is a differential operator acting on $`\delta \stackrel{}{R}`$.

At the critical solution, the pressure is at a maximum; therefore, there must be a solution $`w=\delta \stackrel{}{R}`$ to equation (3) with $`\delta p=0`$. The solution $`w`$ satisfies

$$\widehat{L}w=0$$
(4)
with boundary condition $`w=0`$ at $`D`$. Note that the pressure dependence in this formula arises because $`\widehat{L}=\widehat{L}[\stackrel{}{R}]`$ depends implicitly on the pressure $`p`$ through $`\stackrel{}{R}`$. Hence, the existence of a nonzero $`w`$ is a diagnostic for finding the critical solution to (2) and the corresponding critical pressure $`p^{}`$.

Optimization Method. \- Now, to find the optimal nozzle, we need to derive a relation between the change in critical pressure and change in nozzle shape. Since pressure and volume are conjugate variables, increasing critical pressure is tantamount to decreasing critical volume. By iteratively changing the nozzle shape to increase critical pressure, we will thus arrive at a nozzle which produces smaller droplets. We compare the critical volume of the deformed nozzle with that of the circular nozzle that corresponds to the same critical pressure, since pressure is the control variable in practical situations.

Suppose that a given nozzle shape $`C`$ has a critical pressure $`p^{}`$, a critical droplet shape $`\stackrel{}{R}^{}`$, and a corresponding $`w`$. All of these quantities change when the nozzle shape $`CC+\delta C`$. The change in the droplet shape $`\delta \stackrel{}{R}`$ is linearly related to the pressure change $`\delta p`$ by equation (3) with the boundary condition $`\delta \stackrel{}{R}=\delta C`$ at $`D`$. On the other hand, since the critical solution maximizes the critical pressure, $`w`$ does not change to leading order in $`\delta C`$.

The change in critical pressure induced by $`\delta C`$ can therefore be computed by taking the inner product of both sides of (3) with $`w`$:

$`\gamma w,\widehat{L}\delta R`$ $`=`$ $`\gamma \delta R,\widehat{L}w+\gamma {\displaystyle b(\delta R,w)}`$
$`=`$ $`0+\gamma {\displaystyle b(\delta C,w)}`$
$`=`$ $`w,\delta p.`$
Therefore

$$\delta p=\frac{\gamma b(\delta C,w)}{w,1}.$$
(5)
Here $`b(,)`$ denotes the boundary integrand from integrating by parts. The derivation also uses the self adjointness of $`\widehat{L}`$, which is readily demonstrable by explicit computation 111 In computing an inner product $`f,\widehat{L}g𝑑u𝑑vf\widehat{L}g,`$ where $`f`$ and $`g`$ are arbitrary functions of $`u,v`$, and the integration is over $`D`$, we can undo the differentiation on $`g`$, and instead let the adjoint operator $`\widehat{L}^{}`$ act on $`f`$. Self-adjointness ($`\widehat{L}=\widehat{L}^{}`$) can be verified explicitly. In the derivation, this allows us to simply interchange $`f`$ and $`g`$, while introducing the boundary terms from integration by parts.. Equation (5) is an explicit relation between a change in the nozzle shape ($`\delta C`$) and the resulting change in critical pressure.

Explicit Formula for $`\delta p`$. \- We choose the nozzle $`C`$ to lie in the $`xy`$-plane, enclosing the origin. Then the droplet surface may be given in spherical coordinates by the distance from the origin ($`R`$) as a function of the two angles $`\theta [0,\pi /2]`$ and $`\varphi (0,2\pi ]`$. To avoid the coordinate singularity at the pole ($`\theta =0`$) we use $`u,v`$ given by $`u=\mathrm{tan}(\theta /2)\mathrm{cos}(\varphi )`$ and $`v=\mathrm{tan}(\theta /2)\mathrm{sin}(\varphi )`$. Hence the surface is a scalar function $`R(u,v)`$; its domain $`D`$ is the unit disk in the $`uv`$-plane. We retain $`\varphi `$ to denote the polar angle in the $`uv`$-plane.

An appealing feature of this coordinate system is that the line element remains diagonal:

$$ds^2=dR^2+\mathrm{\Gamma }(du^2+dv^2),$$

where $`\mathrm{\Gamma }=4R^2/(1+u^2+v^2)^2`$. It is then straightforward to compute the free energy $`E=(\gamma dApdV)`$ which yields, upon variation, the Young-Laplace equation

$$\stackrel{}{}(C\stackrel{}{}R)+AR=F,$$
(6)
where $`\stackrel{}{}`$ is the usual gradient operator in the $`uv`$-plane. The coefficients are

$`C`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1+(\frac{1+\rho ^2}{2R})^2(\stackrel{}{}R)^2}}},`$
$`A`$ $`=`$ $`C\left({\displaystyle \frac{(\stackrel{}{}R)^2}{R^2}}+{\displaystyle \frac{8}{(1+\rho ^2)^2}}\right),`$
$`F`$ $`=`$ $`p{\displaystyle \frac{4R^2}{(1+\rho ^2)^2}},`$
where $`\rho ^2u^2+v^2`$ is the radial coordinate in the $`uv`$-plane. $`D`$ corresponds to $`\rho =1`$.

In our coordinate system, the pressure change is

$$\delta p=\frac{1}{\delta _wV}𝑑\varphi \delta C\frac{w_\rho R(R^2+R_\varphi ^2)}{(R^2+R_\rho ^2+R_\varphi ^2)^{3/2}},$$
(7)
where $`\delta _wVd^2\rho w\frac{4R^2}{(1+\rho ^2)^2}`$. Here and in the following we use subscripts to denote partial differentiation.

We can recast this expression into a form that is more geometric. First, the contact angle $`\alpha `$ between the drop and the plane of the nozzle is given by $`\mathrm{cot}\alpha (\varphi )=R_\rho /(R^2+R_\varphi ^2)^{1/2}|_{\rho =1}`$ where the right hand side is evaluated at the boundary. Second, we define $`w_{}w_\rho /(R^2+R_\varphi ^2)^{1/2}|_{\rho =1}`$ which can be understood as follows - note that $`w`$ is the difference between the outer and inner solutions as the pressure approaches bifurcation. Using the contact angle given above, this expression is the difference between the slopes (with respect to the vertical) of the outer and inner solutions at the boundary. This is a coordinate independent quantity. Third, we observe that

$$d\varphi \delta cR=\left(d\varphi \sqrt{R^2+R_\varphi ^2}\right)\left(\delta c\frac{R}{\sqrt{R^2+R_\varphi ^2}}\right)=dl\delta N,$$

where $`dl`$ is the line element, and $`\delta N`$ is the change of the nozzle in the direction locally normal to the nozzle. Lastly, the denominator $`\delta _wV`$ in (7) is just the change in volume from changing the surface by $`w`$. Putting these facts together, the pressure change is

$$\delta p=\frac{1}{\delta _wV}𝑑l\delta Nw_{}\mathrm{sin}^3\alpha ,$$

which leads to the prescription for changing the nozzle

$$\delta N\frac{1}{\delta _wV}w_{}\mathrm{sin}^3\alpha .$$
(8)
Clearly, for the circular nozzle, symmetry implies that $`\delta N`$ should be constant. But this amounts to a mere reduction in the size of the nozzle; the shape remains a circle. So the circular nozzle is at an extremum, in fact a minimum of critical pressure for fixed nozzle area.

For a noncircular nozzle, the contact angle isn’t constant, and hence the change according to the above formula cannot be constant. So one may change the critical pressure while fixing the nozzle area. Moreover, since the circular nozzle is the only one (except the infinite strip) with a constant contact angle, the process of deformation does not end.

We apply (8) iteratively to a perturbed circular nozzle to see how the shape evolves away from the circle.

Figure 1 shows the result of iterations starting with a circle deformed by a perturbation with a three-fold symmetry. The perturbation grows with each iteration, and eventually the nozzle shape becomes concave. With each iteration, we have applied a rescaling in order to maintain the nozzle area. Without the area constraint, the nozzle would become arbitrarily small in accordance with (8). We are interested in the shape of the nozzle, not its size. We also apply the Savitzky-Golay filter Press et al. (1992) at each iteration to smooth out the mesh noise. The solutions to the Young-Laplace equations are obtained using the nonlinear PDE solver in the MATLAB$`^\text{®}`$ PDE Toolbox, which implements the finite element method for elliptic equations with variable coefficients, exactly of the form in (6). For each nozzle shape, we start at a pressure below the bifurcation and by choosing different trial solutions obtain both solutions. Then we bring both solutions to just below the critical pressure by stepping up the pressure, using the solution at each step as the trial solution for the next step. We then use the average of the two solutions for our surface, and their difference for $`w`$. The validity of this procedure can be rigorously shown for a circular nozzle, and we expect it to remaind valid for noncircular nozzles as long as the pressure is brought close to critical.

In order to compare and select among nozzle shapes, we need a measure of optimality independent of size. For every nozzle, we rescale its critical volume by the critical volume corresponding to the circular nozzle with the same critical pressure. This dimensionless volume is given by

$$\stackrel{~}{V}=\frac{v^{}}{\frac{2\pi }{3}\left(\frac{2}{p^{}}\right)^3}.$$
(9)
Figure 2 shows a particular sequence of critical properties obtained through our iteration procedure. We see that the critical pressure begins to increase rapidly about the fifth iteration, after which the decrease in $`\stackrel{~}{V}`$ slows down, and the nozzle shape becomes stretched out (see Figure 1). This means that in order to decrease droplet size at a given pressure, one should use a nozzle shape that is roughly triangular, perhaps with somewhat stretched out corners; but further deformation does not lead to significant improvement. Moreover, gravitational instabilities will inevitably become relevant if the “arms” become too long Pitts (1974); Michael and Williams (1976).

It should be emphasized that we have shown a particular example of an improved nozzle, generated by a choice of the initial perturbation. We have tried other perturbations, leading to shapes with, say, four-fold symmetry or without any symmetry, but the three-fold perturbation has yielded the biggest reduction in the normalized critical volume.

So far we have ignored the effects of gravity, but our formalism applies just as well to the problem with gravity. Including gravity means that the pressure would no longer be constant throughout the drop surface, but rather a linear function of height: $`pp\rho _mgh(u,v)`$, where $`\rho _m`$ is the mass density of the liquid, $`g`$ is the gravitational acceleration, $`h`$ is the distance below the nozzle, and $`p`$ now denotes the pressure at the nozzle ($`h=0`$). Although (6) acquires a new term as a result, this term does not contain derivatives and thus does not contribute to the boundary integral. So our formula for the pressure change remains the same in the presence of gravity. To be sure, the nozzle evolution would differ because the contact angle and $`w_{}`$ will be affected by gravity. Moreover, if the nozzle is too large relative to the capillary length, then gravity destabilizes all solutions: it is not possible to suspend a water drop from a meter wide faucet. It would be interesting to examine the case of the intermediate sized nozzle, small enough to have stable solutions, yet large enough to be affected by gravity.

We thank Eric Lauga and Daniel Podolsky for useful discussions at various stages of the work. This research is funded by NSF DMS-0305873 and Harvard MRSEC.


text:
# Possible mechanism for glass-like thermal conductivities in crystals with off-center atoms

(Draft: )

## Abstract

In the filled Ga/Ge clathrate, Eu and Sr are off-center in site 2 but Ba is on-center. All three filler atoms (Ba,Eu,Sr) have low temperature Einstein modes; yet only for the Eu and Sr systems is there a large dip in the thermal conductivity, attributed to the Einstein modes. No dip is observed for Ba. Here we argue that it is the off-center displacement that is crucial for understanding this unexplained difference in behavior. It enhances the coupling between the “rattler” motion and the lattice phonons for the Eu and Sr systems, and turns on/off another scattering mechanism (for 1K $`<`$ T $`<`$ 20K) produced by the presence/absence of off-center sites. The random occupation of different off-center sites produces a high density of symmetry-breaking defects which scatters phonons. It may also be important for improving our understanding of other glassy systems.

An important, universal, characteristic of glass-like systems is a low thermal conductivity $`\kappa `$; $`\kappa `$(T) varies roughly as T2 for T $`<`$ 1K, forms a plateau/dip region somewhere in the range 5-40K, and then increases slowly at higher temperatures. The low T behavior ($`<`$ 1K) is well understood, and is attributed to a broad distribution of tunneling centers (or two level systems), both in glasses and in some disordered crystals that exhibit glass-like behavior such as (KBr)1-x:(KCN)xGrannan88 , while at high temperatures, Rayleigh scattering becomes important. However the intermediate plateau/dip region is less well understood. Several explanations for the plateau have been proposed, but the model of Grannan etalGrannan88  for the mixed KBr:KCN system has an important general mechanism and appears to be applicable for both glassy and disordered crystalline materials. In their model, the plateau/dip is produced by a nearly localized mode \[the libration modes of the CN ion about its center of mass, at THz frequencies\] that resonantly scatters phonons very effectively. Similar models have been applied to glassy seleniumBermejo94  and disordered garnetsGiesting02 . The same approach has also been used to describe the glass-like thermal conductivity in the off-center Eu and Sr filled clathratesNolas00 ; in this case a “rattler” mode with a low Einstein temperature in the range 30-150K (i.e. at THz frequencies) plays the role of the nearly localized mode. However this model fails in explaining the difference in behavior between the Ba and Eu/Sr filled clathrates - no dip is observed for the on-center Ba system. It also cannot explain the $`\kappa `$ $``$ T dependence observed from roughly 1-10KSales01 . Here we will argue that it is the off-center displacement of Eu and Sr that leads to the plateau/dip region in these systems.

In several compounds with large unit cells (skutterudites and type I clathrates), large cages or voids occur in the structure which can be “filled” with several types of atoms. This dramatically alters their physical propertiesBrown80 ; Stetson91 ; Evers94 ; Danebrock96 . When the “filler” ion is considerably smaller than the void, it is loosely bound and can “rattle” around - hence its name. The weak binding leads to a nearly local modeDong01  described by a low Einstein temperature. A low thermal conductivity has been found in a number of such systemsNolas99 ; Nolas00 ; Sales01  and a glass-like model for $`\kappa `$(T), used for the Eu and Sr clathratesNolas00 . Note that a very low value of $`\kappa `$ is crucial for thermoelectric applications because the figure of merit ZT = TS2$`\sigma _e`$/$`\kappa `$ (S is the Seebeck coefficient, $`\sigma _e`$, the electrical conductivity) can be significantly increased. Hence understanding the mechanisms that lead to a small $`\kappa `$, and particularly their ranges of validity, is very important.

Here we consider $`\kappa `$(T) for the (type I) Ga/Ge clathrates whose structure is formed of two cages with internal voids large enough to house another atom (sites 1 and 2); atoms placed in them, such as Eu, are referred to as Eu1 and Eu2. In the larger site 2 cage, Eu2, Sr2, and Ba2 all form rattlers with quite low Einstein temperaturesSales01 ; however Ba2 is on-center while both Eu2 and Sr2 move off-centerNolas98a ; Nolas01 ; Chakoumakos00 ; Chakoumakos01 ; Zhang02 ; Bridges03 .

In Fig. 1a, we re-plot the thermal conductivity data of Nolas etalNolas00 ; $`\kappa `$ is small for Eu and Sr, and there is a well defined dip associated with the Einstein temperatures of the rattlersNolas00 . Nolas etal model the low thermal conductivity below 20K for the Eu and Sr compounds using a broad distribution of off-center tunneling states. At slightly higher temperatures, they used two resonant scattering terms (two Einstein temperatures) to explain the dip (near 25K), plus a Rayleigh scattering term which dominates at high T. Using this sum of terms, they were able to fit $`\kappa `$(T) for both Eu8Ga16Ge30 and Sr8Ga16Ge30 over the range 5- 100K. However a problem arises in using a tunneling model with tunneling splittings up to $``$ 15K;- first the inferred distribution of tunneling states from ultrasonic measurements is much narrower than in a glassKeppens00 , and second (more importantly), very recent Mossbauer measurements indicate that the tunneling frequency is at most 0.44 Ghz - i.e. $``$ 0.02KHermann04 . If one cuts off the tunneling density of states with a step function even at 1K, then $`\kappa `$ rises too rapidly below 15K.

Similar experimental results for $`\kappa `$(T) of these filled clathrates have also been observed by Sales etalSales01  (See Fig. 1b); however, at low T, $`\kappa `$ is considerably lower and the dip is less pronounced. More importantly, for the Eu and Sr data, $`\kappa `$(T) decreases nearly linearly with T below 10K, which can’t be fit with the above model.

The T-dependence for Ba8Ga16Ge30 in which Ba2 is on-centerSales01 , is more perplexing (See Fig. 1b). For this system, $`\kappa `$ is much higher overall and there is no clear dip in $`\kappa `$(T), although the Einstein temperatures are only slightly larger than those of Eu and SrSales01 ; thus the rattler-phonon coupling must be greatly reducedSales01b . Also, the T-dependence below 10-12K is much faster ($`\kappa `$ $``$ T3) for the Ba clathrate. These results raise two important issues - what determines the strength of the coupling between the rattlers and the phonons and how can one understand the T-dependence of $`\kappa `$ from 1-15K? We propose another mechanism to explain this temperature range; it may also contribute to $`\kappa `$ in other glassy systems.

The off-center displacement appears to play two additional, crucial roles in increasing the phonon scattering (in addition to producing tunneling centers). First from XAFS, the off-center atom is partially bonded to the nearest neighbor Ga/Ge atoms on the side of the cageBridges03 , forming a random array of symmetry-breaking defects. We show below that a high density of such defects can lead to a $`\kappa `$ $``$ T dependence at low and medium temperatures. A likely more important consequence of this bonding is that the motions of the off-center rattler atom are coupled more strongly to just a few of the closest atoms in the surrounding Ga/Ge framework and hence to the phonon modes - i.e. the motions of the Eu2 and nearest Ga/Ge atoms are positively correlated.

We first address the use of the tunneling model and the range of reasonable tunneling splittings. Sethna and coworkersGrannan88  suggest that for the light CN molecule the tunneling splittings (for reorientations of the CN axis) might extend to about 1K; measured tunneling splittings on similar, more dilute systems (CN- and OH-) are less than 3KBridges75 ; Narayanamurti70 . For heavy off-center ions such as Ag+ and Cu+ in the alkali halides, the observed tunneling energies are also low, much less than 1KBridges75 ; only for the light off-center 7Li+ ion does the tunneling splitting reach 1.2K at normal pressuresHerendeen69 . The tunneling behavior can be finely tuned by applying hydrostatic pressure to reduce the potential barrier between equivalent off-center sitesMorgan87 ; Wang92 , but even in such cases the tunneling splittings for heavy atoms such as Ag+ are about 1K when the off- to on-center transition begins to occurMorgan87 .

Within the simple two-well tunneling model, having a tunneling splitting ($``$ 10-15K) that is $``$ 10% of the attempt frequency (the rattler frequency) is inconsistent with the assumptions for the model; for a heavy mass such as Eu, it leads to an unreasonably low potential barrier between the wells ($`<`$ 1 meV (12K); - i.e. less than the tunneling splitting). On the experimental side, a broad distribution of tunneling states should also lead to a linear contribution to the heat capacity at low T when the tunneling states are excited; however the heat capacity for the Ba, Sr, and Eu systems at low T is exponential and can be fit to an Einstein modelSales01 ; thus there is no evidence for tunneling states extending up to 10-15K in such data. With the recent observations that the tunneling frequencies are at most 0.02KHermann04  for the Eu system, another mechanism is needed to describe $`\kappa `$ at low T.

Since atomic forces are generally a very strong function of the distance between the atoms/ions, the force between the central atom (e.g. Sr2 or Eu2) and each of the nearest equivalent Ga/Ge atoms in the cage, will increase substantially for an off-center displacement. The resulting vibrations of the off-center rattler/cage system then depend on the rigidity of the cage (the Ga/Ge cage is quite stiffChakoumakos00 ; Sales01 ), the Eu2-Ga/Ge bond strength, and the number of Eu2-Ga/Ge bonds. In the limit of an idealized rigid cage, there would be no motion of the cage atoms (i.e. no phonon-coupling) as the rattler atom vibrates; the reduced mass of the local mode would be that of the rattler atom as assumed previouslyChakoumakos00 ; Chakoumakos01 ; Bridges03 . For the opposite extreme of a very soft cage, the local mode reduced mass for motion along the bond direction would approach the reduced mass of the atom-pair involved, e.g. Eu and Ge in Eu8Ga16Ge30. In this case, there would be large motions of the nearest cage atoms (Ge) as a result of a local mode vibration, which would couple directly to the phonons of the clathrate framework. Because the clathrate cage is quite stiff, the reduced mass is expected to be close (but not equal) to the free rattler mass.

To further understand this coupling, consider a 1-D model of a cage (mass $`M`$) connected to a rattler (mass $`m`$), and let the respective displacements be $`x`$ and $`X`$. Then $`Mx=mX`$; we expect the matrix element for scattering to be proportional to $`x/X`$=$`m/M`$, and the coupling to $`(x/X)^2`$. The effective mass of the cage decreases when the rattler moves off-center (because only a small fraction of the cage atoms are directly bonded to the rattler) and thus the coupling to phonons increases. Also note that for the off-center case, rattler vibrations could be either $``$ radial (along the bond) or perpendicular to the bond. The latter are similar to libration modes and would have a lower $`\mathrm{\Theta }_E`$ than for radial motion which was not considered for the KCl:KCN system.

To investigate the various contributions to $`\kappa `$ we first reproduced the calculation of Nolas et al.Nolas00 . $`\kappa `$ is given by:

$`\kappa `$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle _0^{\omega _D}}c_sC(\omega ,T)l(\omega ,T)𝑑\omega `$ (1)
$`l(\omega )`$ $`=`$ $`(\mathrm{\Lambda }_R^1+\mathrm{\Lambda }_{res}^1+\mathrm{\Lambda }_{TS}^1)^1+l_{min}`$ (2)
$`\mathrm{\Lambda }_R^1`$ $`=`$ $`D(\mathrm{}\omega /k)^4`$
$`\mathrm{\Lambda }_{res}^1`$ $`=`$ $`{\displaystyle \underset{i}{}}C_i\omega ^2T^2/[(\omega _{E_i}^2\omega ^2)^2\gamma _{E_i}(\omega _{E_i}\omega )^2]`$
$`\mathrm{\Lambda }_{TS}^1`$ $`=`$ $`A(\mathrm{}\omega /k_B)tanh(\mathrm{}\omega /2k_BT)`$ (3)
$`+(A/2)(k_B/\mathrm{}\omega +B^1T^3)^1`$
where $`c_s`$ is the speed of sound, C($`\omega `$,T) is the heat capacity of phonons with frequency $`\omega `$ and $`l(\omega )`$ is the total mean free path which has three components - Rayleigh scattering ($`\mathrm{\Lambda }_R`$), resonant scattering from the Einstein modes ($`\mathrm{\Lambda }_{res}`$) and a contribution from a broad distribution of tunneling states ($`\mathrm{\Lambda }_{TS}`$); the lower limit is constrained to $`l_{min}`$. The Einstein frequencies $`\omega _{E_i}`$ were obtained from structural data and the constants A, B, C1, C2, D, $`\gamma _{E_i}`$ and $`l_{min}`$ are given in ref. Nolas00, . We obtained similar fits although the constants varied slightly. We then reduced the coupling to the resonant modes - the Ci constants - by a factor of 10 and left the other terms unchanged. This result is shown by the dashed line (Fig. 1a) and is very similar in shape to the data of SalesSales01  for Ba8Ga16Ge30 (Fig 1b). Thus a major part of the difference in the T dependence of $`\kappa `$ between Ba8Ga16Ge30, and the Eu or Sr compounds can indeed be explained by a decreased coupling between on-center Ba2 and phonons in the Ga/Ge framework as described above.

An important parameter for probing the rattler-phonon coupling is the reduced mass $`\mu `$, of the local mode oscillator: however, $`\mu `$ is not an easy parameter to obtain. One recent study of rattlers in the skutteruditesCao03 , found the reduced mass to be somewhat less than the atomic mass of the rattler (error $``$ 20%), but more measurements are needed.

A second consequence of the rattler being partially bonded to a side of the cage is the formation of a high-density, random network of symmetry-breaking mass defects - i.e. one of the four off-center orientations is occupied, but the off-center site varies randomly throughout the crystal. The disorder introduced by the random occupation of one of the off-center minima provides insight as to why the crystal exhibits glass-like behavior. This model has similar features to the disorder introduced via irradiation of a quartz crystalBerman50 ; Keppens96  but there is a key difference - the atoms are randomly displaced in the irradiated sample which leads to a range of bond lengths and hence to the possibility of a broad distribution of tunneling states. For the Eu/Sr clathrates the off-center displacement is well defined; the disorder arises from a random occupation of the four off-center sites. Locally such disorder changes the speed of sound as the phonon wave passes a rattler atom; the resulting phase changes must then be included. Note that in this case there is an interference in the scattering between different rattler sites, whereas Rayleigh scattering assumes a single scattering process.

Ziman has considered this type of scattering in some detailZiman60  \- in his model, the effective mean free path for this contribution is

$$\mathrm{\Lambda }_{disorder}(q)=\frac{A}{L}\frac{c_s^2}{(\delta c_sq)^2}$$
(4)
where $`q`$ is the wavenumber, $`\delta c_s`$ is the variation of the sound speed, $`L`$ is the smallest distance between scatterers, and $`A`$ is a numerical constant of order 1. For very short wavelengths $`\lambda `$ $`<<`$ $`L`$, $`\mathrm{\Lambda }_{disorder}`$ goes to a constant - but that is already included in the model of Nolas et al.Nolas00  by the $`l_{min}`$ term. Using the Debye model ($`\omega =c_sq`$) we then replaced the tunneling contribution $`\mathrm{\Lambda }_{TS}^1`$ in Equ. Possible mechanism for glass-like thermal conductivities in crystals with off-center atoms by

$$\mathrm{\Lambda }_{disorder}^1=G\omega ^2.$$
(5)
In Fig. 1a we show that with GEu = 3.7 10-19s2/m (GSr = 6.4 10-20s2/m), and a slight change of the Ci parameters this model describes the T behavior of $`\kappa `$ (of Ref. Nolas00, ) just as well as by using the tunneling model. Taking $`L`$ $``$ 5 10-10m, these numbers imply $`\delta c_s`$/$`c_s`$ $``$ 7% (Eu) and 3%(Sr). Note that this disorder mechanism is important at moderately long wavelengths - i.e. for temperatures up to 20-30K, and alone would yield $`\kappa `$(T)$``$ T.

Attempts to fit the Sales et al data for Eu and SrSales01  using the model of Nolas etalNolas01  were not successfulBaumbach04  \- the data from 2-10K vary nearly linearly with T, and cannot be fit by the $``$ T2 dependence for the tunneling model at low T. However, we can fit the data quite well (Fig. 1b) using $`\mathrm{\Lambda }_{disorder}^1`$ but larger values for G and the Ci are needed to fit the lower thermal conductivity.

For Ba, very much smaller values for Ci are needed, which confirms the earlier assertion that the coupling between phonons and the rattler vibrations is reduced. In addition, the $`\mathrm{\Lambda }_{disorder}^1`$ term cannot fit the low T data for Ba below 12K; $`\kappa `$ is higher and the approximately $`T^3`$ dependence suggests that a finite sample size contribution is needed at low T \- details will be given in a separate paperBaumbach04 . The lack of a $`\mathrm{\Lambda }_{disorder}^1`$ term is consistent with Ba being on-center, as in that case the $`\mathrm{\Lambda }_{disorder}^1`$ term should disappear.

In summary, we have pointed out that an off-center rattler enhances the scattering of phonons in two ways. First, the off-center displacement increases the coupling of the rattler motion to the phonons and hence increases a dip in the thermal conductivity near $``$30K. This provides an explanation for the the lack of a dip in $`\kappa `$(T) for the on-center Ba system. Second, the off-center atom on the side of the cage, introduces a quasi-random set of “extra” atoms that will produce local changes in the sound velocity. This in turn leads to a mean free path that varies as $`\omega ^2`$, which provides an alternative explanation for the reduction of the thermal conductivity for 1K $`<`$ T $`<`$ 20K. It results in a linear T dependence in this range which fits the data of Sales etalSales01  well.

Several questions remain to be answered - for the clathrates, what are the effective masses of the local Einstein modes and what are the actual tunneling splittings (there will be 3 levels) in these systems? Can the Sr system be moved on-center via high hydrostatic pressure and thus allow an investigation of variations in the rattler-phonon coupling? Can these systems be optimized for thermoelectric applications over a particular T-range, by varying the off-center rattler’s mass and/or Einstein temperatures? More generally does the off-center displacement play a role in producing the plateau/dip regions in other glassy systems?

###### Acknowledgements.

We thank J. Rudnick and S. Shastry for helpful discussions. The work at UCSC was supported in part by NSF grant DMR0071863. The experiments were performed at SSRL, which is operated by the DOE, Division of Chemical Sciences, and by the NIH, Biomedical Resource Technology Program, Division of Research Resources.


text:
# Statistically Locked-in Transport Through Periodic Potential Landscapes

## Abstract

Classical particles driven through periodically modulated potential energy landscapes are predicted to follow a Devil’s staircase hierarchy of commensurate trajectories depending on the orientation of the driving force. Recent experiments on colloidal spheres flowing through arrays of optical traps do indeed reveal such a hierarchy, but not with the predicted structure. The microscopic trajectories, moreover, appear to be random, with commensurability emerging only in a statistical sense. We introduce an idealized model for periodically modulated transport in the presence of randomness that captures both the structure and statistics of such statistically locked-in states.

Objects driven through periodic potential energy surfaces face a myriad of choices: either they follow the driving force or they become entrained along any of the commensurate directions through the landscape. Variants of this problem appear in areas as diverse as driven charge density waves Brown et al. (1984), electronic energy states in two-dimensional electron gases Wiersig and Ahn (2001), atom migration on crystal surfaces Pierre-Louis and Haftel (2001), chemical kinetics, and flux flow in type-II superconductors Reichhardt and Nori (1999). Quite recently, this problem was investigated Korda et al. (2002) using a monolayer of colloidal spheres in flowing water as a model system and a square array of holographic optical tweezers hot  to provide the periodic potential energy surface. Depending on the array’s orientation with respect to the driving force, the spheres were observed to trace out a Devil’s staircase hierarchy of commensurate directions through the array, with particular directions being preferentially selected over certain ranges of orientations Korda et al. (2002). Trajectories deflected by a preference for commensurability are said to be kinetically locked in to the lattice.

This anticipated result Ott (1993) was accompanied by two surprises. In the first place, not all kinetically locked-in states were centered on simple commensurate directions. Still more surprisingly, particles’ microscopic trajectories in high-order locked-in states did not consist of sequences of commensurate jumps, but rather consisted of seemingly random lower-order hops whose combination, nonetheless, were commensurate. The appearance of statistical rather than deterministic commensurability suggests an unexpected role for thermal randomness in structuring transport through periodic potentials, and has been dubbed statistical lock-in Korda et al. (2002).

This Letter presents a highly simplified model for statistically locked-in transport through mesoscopic potential energy landscapes that nonetheless accounts for the emergence of combination jumps and their statistical commensurability. In particular this model reveals how the potential energy landscape’s structure and extent establish the discrete spectrum of travel directions selected by biased random walkers. While our discussion is directed toward the purely classical behavior of flowing colloids, similar results should emerge for biased quantum mechanical hopping through arrays of potential wells or barriers. Finally, we argue that the kinetically locked-in state selected at a given orientation can depend sensitively on particle size. Statistical lock-in therefore presents new opportunities for continuously sorting heterogeneous materials into multiple fractions simultaneously.

We first consider an array of optical traps, as shown in Fig. 1. The centers of the traps are on a square lattice with lattice constant $`b`$. The array is taken to extend indefinitely in the $`y`$ direction, while its extent in the $`x`$ direction is $`N`$ lattice constants. Rather than attempting to account for the detailed structure of a practical trap array, we instead focus on the simplest model that captures the experimental phenomenology. For this reason, we model the traps as circular regions of radius $`a<b/2`$ centered on grid points. Once a particle enters a trap’s domain, it is translated directly to the center. A particle passing outside the circle is assumed to be unaffected. A “trapped” particle eventually escapes and is carried along to its next encounter by the driving force. This is consistent with the observed dynamics of colloidal particles flowing through holographic optical trap arrays Korda et al. (2002).

In such a situation there are two distinct regimes we may consider. The first is the “ballistic” regime where most particles pass through the lattice without encountering traps. This corresponds to the limit of small trap density, $`\eta =a/b`$. For sufficiently large values of $`\eta `$, traps are packed closely enough together that particles’ trajectories consist mostly of jumps from one trap to another. This we term the “lattice gas” regime. We shall first derive the condition on $`\eta `$ for being in the lattice gas regime and consider only the case that meets this condition for the rest of the discussion. Consider the array as shown in Fig. 1. Consider a line drawn from one of the lattice points on the left edge which we will take to be the origin. Let the line have a slope $`m`$. Now the shortest distance, $`s`$, between the line and a generic point on the lattice, represented by $`(p,q)`$, with $`p\{1,\mathrm{},N\}`$ and $`q`$, is given by

$$s=\frac{b|qpm|}{(1+m^2)^{\frac{1}{2}}}$$
(1)
We now wish to compute the largest possible value, $`s_m`$, of the distance from the line to the lattice point nearest to the line for all possible slopes. If the trap radius $`a`$ is larger than this distance then all straight line trajectories originating at a lattice point on the left side will necessarily intersect a trap’s region of influence. Now the least upper bound on the minimum value of $`|qpm|`$ can be shown to be $`1/(N+1)`$ Jensen (2003). Using the smallest allowed value of $`m`$, which is the slope of the line through the origin and tangent to the circle centered at $`(1,0)`$, we estimate for an upper bound on $`s_m`$

$$s_m\frac{1}{N+1}(b^2a^2)^{\frac{1}{2}}$$
(2)
Thus a sufficient (and necessary) condition that any straight line trajectory starting at the origin must intersect a trap before exiting the lattice is $`s_m=a`$, or

$$\eta >\frac{1}{\sqrt{(N+1)^2+1}}.$$
(3)
For the remainder of the discussion we shall assume that this condition is satisfied and that every trajectory involves a sequence of inter-trap jumps.

Every trajectory starting at a lattice point on the left side of the lattice and exiting at the right side can be decomposed into a discrete sequence of steps from one lattice point to another, followed by exit from the trap array. For ease of modeling we treat the exit also as a step described by the lattice vector to the trap that would have captured the exiting particle had the trap array been extended to infinity in the $`x`$ direction.

In the physical system Korda et al. (2002), colloidal particles are driven through an array of optical traps by a steady flow in the surrounding medium oriented at angle $`\theta `$ to the trap array’s  axis. These particles are also subject to random thermal forces causing them to diffuse. We will first consider the case where there is no diffusion. In the absence of traps, the particles then would travel in straight lines in the $`\theta `$ direction. Within the trap array, however, a particle’s displacement in each step is determined by the position of the center of the next trap it encounters along its straight line trajectory. Thus, even in the absence of diffusion, the effective direction of motion can differ from the direction dictated by the flow. The effective direction of motion will be described by an angle $`\psi `$ given by $`\mathrm{tan}\psi =\mathrm{}_y/\mathrm{}_x`$ where $`\stackrel{}{\mathrm{}}(\mathrm{}_x,\mathrm{}_y)`$ is the lattice vector describing the step. The functional dependence of $`\psi `$ on the flow direction, $`\theta `$, is given by $`\mathrm{tan}\psi =q/p`$, where $`(p,q)`$ are the smallest integers satisfying

$$|qp\mathrm{tan}\theta |\eta (1+\mathrm{tan}\theta ^2)^{\frac{1}{2}}$$
(4)
with $`p\{1,\mathrm{},N\}`$. Figure 2 shows the result.

Corresponding to a given value of $`\theta `$ there exists a unique value of $`\psi `$ and hence a lattice vector, $`\stackrel{}{\mathrm{}}`$, describing that step. In general, a trajectory $`T`$ consists of a sequence of steps, $`T\{\stackrel{}{\mathrm{}}^1,\stackrel{}{\mathrm{}}^2,\mathrm{},\stackrel{}{\mathrm{}}^n\}`$, where $`\stackrel{}{\mathrm{}}^i`$ is a lattice vector describing the $`i`$-th step of an $`n`$-step trajectory. Here $`\stackrel{}{\mathrm{}}^i=(\mathrm{}_x^i,\mathrm{}_y^i)`$, with $`\mathrm{}_x^i\{1,\mathrm{},N\}`$ and $`\mathrm{}_y^i`$, such that $`\mathrm{}_x^i`$ and $`\mathrm{}_y^i`$ are relatively prime, $`_{i=1}^n\mathrm{}_x^iN`$, and $`_{i=1}^{n1}\mathrm{}_x^i<N`$. These conditions ensure that we do not count multiples of lattice vectors as distinct and that each trajectory terminates at the right side of the trap array. We can now define a *mean* effective direction, $`\overline{\psi }`$, for the entire trajectory by

$$\mathrm{tan}\overline{\psi }(T)=\frac{_{i=1}^n\mathrm{}_y^i}{_{i=1}^n\mathrm{}_x^i}$$
(5)
This value characterizes the particles’ overall transport through the array, and would be reflected in measurements of bulk transport properties, such as a Hall voltage in the case of a periodically modulated two-dimensional electron gas Wiersig and Ahn (2001). Without diffusion or some other randomizing process, all the steps in every trajectory would be identical, with $`\overline{\psi }\psi `$ for any orientation $`\theta `$.

We now consider the case in which diffusion or scattering randomizes the particles’ trajectories. The natural assumption is that randomness would smear out the noise-free transport characteristic in Fig. 2, until, on average, the trajectories simply follow the driving direction, $`\overline{\psi }=\theta `$. In fact, far more interesting behavior results.

We model trajectory fluctuations by assuming that a particle leaving a trap is equally likely to travel along any direction in a wedge of opening angle $`2\delta \theta `$ centered on the flow direction. In terms of experimental parameters for colloidal transport, $`\delta \theta `$ depends on the temperature and the particles’ mobility. The path a particle takes through the array is no longer deterministic, and there are many possible trajectories for any given driving direction.

We first compute the probability of a step being a lattice vector $`\stackrel{}{\mathrm{}}=(\mathrm{}_x,\mathrm{}_y)`$ with $`\mathrm{}_x\{1,\mathrm{},N\}`$ and $`\mathrm{}_y`$ such that $`\mathrm{}_x`$ and $`\mathrm{}_y`$ are relatively prime. This probability corresponds to the fraction of angles in the range $`[\theta \delta \theta ,\theta +\delta \theta ]`$ whose effective direction $`\psi `$ corresponds to lattice vector $`\stackrel{}{\mathrm{}}`$. The probability is given by

$$P_{\mathrm{}}(\stackrel{}{\mathrm{}})=\frac{1}{2\delta \theta }_{\theta \delta \theta }^{\theta +\delta \theta }\delta _{\frac{\mathrm{}_y}{\mathrm{}_x},\mathrm{tan}\psi (\theta ^{})}𝑑\theta ^{},$$
(6)
where $`\delta _{j,k}`$ is the Kronecker delta function. If the condition in Eq. (3) is satisfied then there will only be a finite number of lattice vectors with non-zero probability as defined above, and the sum of the probabilities will be unity. Consequently there is only a finite number of possible trajectories consisting of sequences of these lattice vectors. Let there be $`M`$ distinct possible trajectories given by the set $`\{T_i\}`$. Here $`T_i=\{\stackrel{}{\mathrm{}}^1,\stackrel{}{\mathrm{}}^2,\mathrm{},\stackrel{}{\mathrm{}}^{n_i}\}`$, where $`n_i`$ is the number of steps in trajectory $`i`$. The probability for a particle to take a particular trajectory $`T_i`$ is given by

$$P_T(T_i)=\frac{_{j=1}^{n_i}P_{\mathrm{}}(\stackrel{}{\mathrm{}}^j)}{_{i=1}^M_{j=1}^{n_i}P_{\mathrm{}}(\stackrel{}{\mathrm{}}^j)}.$$
(7)
The probability that the particle’s trajectory will carry it in direction $`\overline{\psi _0}`$, is then given by

$$P(\overline{\psi _0})=\underset{i=1}{\overset{M}{}}P_T(T_i)\delta _{\mathrm{tan}\overline{\psi _0},\mathrm{tan}\overline{\psi }(T_i)}.$$
(8)
This distribution is nonzero for only a few discrete directions and thus is better characterized by its most probable value, $`\overline{\psi ^{}}`$, with $`P(\overline{\psi ^{}})P(\overline{\psi })\overline{\psi }`$, rather than its mean.

Figure 3 shows the results for an $`N=9`$ lattice with $`\eta =0.26`$ when the dispersal angle $`\delta \theta `$ that models disorder is increased from the deterministic limit, $`\delta \theta =0`$ to $`\delta \theta =0.09`$. These values were selected to mimic the experimental conditions in Ref. Korda et al. (2002). Data from that study are plotted as squares in Fig. 3 for comparison. The calculated results demonstrate that disorder does not necessarily wipe out the structure of kinetically locked-in transport, but rather reconfigures the pattern of plateaus. In particular, not all statistically locked-in plateaus are centered on commensurate directions. This was one of the principal surprises to emerge from experimental studies of colloid flowing through optical trap arrays Korda et al. (2002) and is a clear feature of the experimental data in Fig. 3. Now this restructuring can be explained as the emergence of high-order locked-in plateaus from a statistical sampling of low-order jumps. Indeed, as the bottom trace in Fig. 3 demonstrates, many of the commensurate directions corresponding to statistically locked-in plateaus are not microscopically accessible at $`\eta =0.26`$. These plateaus are absent, therefore, from the underlying pattern of deterministically locked-in steps.

By contrast to the higher-order states, the principal plateaus at $`q/p=0`$ and $`q/p=1`$ are both microscopically and deterministically locked-in. These correspond to transport along the $`[10]`$ and $`[11]`$ lattice directions, respectively, which cannot be decomposed into lower-index jumps. This distinction also is clear in the experimental data Korda et al. (2002). Under some conditions, other directions such as $`[21]`$ also can become deterministically locked-in, particularly when both $`\eta `$ and $`\delta \theta `$ are small.

The calculated transport characteristics do not agree with the experimental data in all detail. Qualitative differences, such as the jump in the experimental data at $`q/p=1/2`$ can be ascribed to the $`10\times 10`$ structure of the experiment’s optical trap array, which differs from the extended geometry we considered here. Flowing colloidal particles also are not drawn directly to optical traps’ centers as they jump through a trap array, but rather follow more subtle and complicated trajectories. Clearly, though these measurements more closely resemble our model’s predictions than they resemble the simple deterministic spectrum for the same conditions.

The distribution of statistically locked-in states changes dramatically with the size of the array. Accessible states proliferate as $`N`$ increases, first smearing out the spectrum of statistically locked-in states, and then erasing even the deterministic plateaus. In the large-$`N`$ limit and for $`\delta \theta >0`$, the array has no overall effect on transport: $`\overline{\psi ^{}}=\overline{\psi }=\theta `$. The hierarchically structured transport characteristics discovered experimentally and explained here are features of mesoscopic systems.

Being able to tune transport properties by adjusting the number, distribution, and size of discrete potential wells suggests a practical method for continuously separating heterogeneous fluid-borne samples into multiple fractions simultaneously. Consider, for example, objects driven through arrays of discrete optical tweezers Ashkin et al. (1986), each of which is formed by a strongly focused beam of light. Different objects passing through this optical intensity field experience different potential energy landscapes depending on their size, shape, and composition. This effect already has been exploited to sort mixtures into two fractions with a single line of optical traps Ladavac et al. (2003). In terms of our model, the disk-like wells would appear to be larger for larger or more strongly trapped objects. Larger objects also would diffuse less vigorously, and so would sample a simpler spectrum of statistically locked-in states. Consequently, the same array of potential wells can have a substantially different effect on different objects’ transport.

Figure 4 demonstrates that smaller particles can be deflected into high-order statistically locked in states under the same conditions that larger objects are deterministically locked in to the $`[10]`$ direction. The resulting angular separation is large enough that the two resulting fractions could be collected, for example in microfluidic channels. Indeed, this figure makes clear that several fractions might be continuously selected out of a mixed sample by a single static array of optical traps.

We are grateful for extensive enlightening conversations with Tom Witten. This work was supported primarily by the MRSEC program of the National Science Foundation through Grant Number DMR-0213745 and in part by NSF Grants DMR-0304906 and DBI-0233971.