text:
# Effective interactions between inclusions in complex fluids driven out of equilibrium

## I INTRODUCTION

A prominent issue in soft condensed matter physics is the understanding of the equilibrium phase behavior of mesoscopic particles immersed in complex fluids: colloidal suspensions colloid , liquid droplets in liquid crystals inclusionsnem , inclusions in lipid membranes prot , charged particles in electrolytes electrolytes ,… A standard and fruitful procedure to describe the large scale properties of inclusions in a soft medium consists in integrating out the numerous ”solvent” degrees of freedom via ensemble averaging colloid ; onsager . The interaction between the embedded particles is then described by “effective potentials”. The latter can modify the genuine interparticles interactions (Coulomb, van der Waals,…) or give rise to entirely new effects. For instance, in some (equilibrium) cases, the external objects do not modify the ground state energy of the medium but only alter its thermal-fluctuations spectrum. The resulting entropic effective interaction has consequently an amplitude proportional to the thermal energy $`k_\mathrm{B}T`$ and a range comparable to that of the correlations of the fluctuations of the medium raminrmp ; aa . Such fluctuation-induced (i.e. entropic) interactions are commonly referred to as thermal Casimir interactions in analogy with their famous quantum equivalent casimir .

More recently, many experimental studies have reported the organization of particles embedded in fluids when the latter are driven out of equilibrium by the application of external fields: vibrated ; spinner ; pompes ; grier . Extensions of the concepts of ”effective potential” and ”entropic forces” to out of equilibrium situations have been scarce grier ; toddbrenner ; hastings . Indeed effective potentials can not be simply derived from a free energy in an out-of-equilibrium context, and only the instantaneous force acting on the external objects for a given configuration of the medium can be defined. The effective interactions between the host objects should then strongly depend on the dynamics ruling the temporal evolution of the medium (in contrast to the equilibrium Casimir case).

In this paper, we attempt to extend the paradigm of Casimir effective interactions to objects immersed in soft systems ”shaken” by external energy sources that can not be a priori modeled by heat baths. We use a model, detailed in section II, where both the medium and the objects are very simple. (i) The hosting medium is described by a scalar free field $`\varphi `$ living in a $`d`$ dimensional space ($`d=1,2,3`$): for instance, $`\varphi `$ could represent the height profile of a fluid interface ($`d=2`$) or a contact line ($`d=1`$interface , or the angular deviation of the director of a nematic liquid crystal ($`d=3`$prostdegennes . (ii) The external objects are two identical rigid parallel plates that enforce a zero field on their surfaces (Dirichlet boundary conditions). Among the numerous possible choices for the shaking sources, we explore three different cases of experimental relevance (c.f. Fig. 1): (A) a localized monochromatic shaking, (B) an uniform monochromatic shaking, (C) a spatially uncorrelated colored noise. All these shaking induce effective interactions between the plates, which we characterize in section III by the corresponding average forces, discussing their sign, amplitude and range. We compare these features with those of the usual thermal Casimir effect. The detail of the calculations and a more thorough analysis of the shaking-induced interactions are gathered in section IV. This more technical section also includes some comments on the time dependence and fluctuations of these shaking-induced forces. Section V ends the paper with a synthetic summary of our main results and an outlook on possible applications.

## II Model and notations.

### II.1 Energetics and Force on the plates :

We consider two plates in a soft medium separated by a distance $`L`$ much shorter than their lateral extension $`L_{}`$ as depicted in Fig. 2. Their thickness plays no role in all that follows and will be set to $`0`$. The soft medium is modeled by a scalar field $`\varphi `$ associated to the elastic Hamiltonian :

$$=\frac{\kappa }{2}d^d\stackrel{}{R}\left[\varphi (\stackrel{}{R})\right]^2,$$
(1)
with $`\kappa `$ the elastic modulus and $`\stackrel{}{R}=(x,𝐫)`$ where $`x`$ is the coordinate normal to the plate, see Fig. 2. $`\varphi `$ is taken dimensionless.

We restrict ourselves to strong plates-field interactions, modeled by Dirichlet boundary conditions (DBC) on the plates: $`\varphi (x=0,𝐫)=\varphi (x=L,𝐫)=0`$.

We first focus our attention on the force $`F^>`$ that the soft medium exerts on the right hand side of the plate located at $`x=0`$. For a given configuration of the elastic field, $`F^>`$ is given by the integral over $`𝐫`$ of the stress tensor component $`T_{xx}(𝐫)`$ landau . The components $`T_{ij}`$ of this stress tensor can be obtained by the virtual displacement method taking into account the DBC aa . In the present geometry, $`T_{xx}(𝐫)=\frac{1}{2}\kappa \left[_x\varphi (0^+,𝐫)\right]^2`$, so that

$$F^>=\frac{\kappa }{2}𝑑𝐫[_x\varphi (0^+,𝐫)]^2.$$
(2)
Note that this force pushes the plate toward the negative $`x`$ direction whatever the field conformation.

The net force $`F`$ on the plate is the algebraic sum of the contributions of the two sides: $`FF^>+F^<`$. The force $`F^<`$ is given by a formula similar to Eq. (2), so that:

$$F=\frac{\kappa }{2}𝑑𝐫[_x\varphi (0^+,𝐫)]^2+\frac{\kappa }{2}𝑑𝐫[_x\varphi (0^{},𝐫)]^2.$$
(3)
### II.2 Dynamics :

For sake of clarity all calculations are performed within the simplest local and strongly dissipative dynamical scheme: the $`\varphi `$ field evolves under the application of an external ”shaking” source $`\eta `$ as :

$`\gamma _t\varphi `$ $`=`$ $`{\displaystyle \frac{\delta }{\delta \varphi }}+\eta ,`$ (4)
$`=`$ $`\kappa ^2\varphi +\eta ,`$
$`\varphi (0,𝐫,t)`$ $`=`$ $`\varphi (L,𝐫,t)=0.`$ (5)
Here $`\gamma `$ is a generalized friction coefficient.

A description of the dynamics of the plates is beyond the scope of this paper, and we assume that they are fixed. $`L`$ is thus a constant. After partial Fourier transform with respect to $`t`$ and $`𝐫`$, equations (4,5) can be recast into:

$$\varphi _{\omega ,𝐪}(x)=_0^L𝑑x^{}_{\omega ,𝐪}^>(x,x^{})\eta _{\omega ,𝐪}(x^{}),$$
(6)
where the response function $`^>`$ corresponds to the diffusion kernel in the slab geometry with DBC and with a diffusion constant $`\kappa /\gamma `$. The Fourier transforms used in this paper are defined by $`f_𝐪d𝐫f(𝐫)\mathrm{exp}(i𝐪.𝐱)`$ and $`g_\omega 𝑑tg(t)\mathrm{exp}(i\omega t)`$.

A generalization of our results obtained with (4,5) to other slow dynamics that obey an equation of the form (6) with a different kernel will be commented upon throughout the text. This generalization is simple provided the dynamic linear response of $`\varphi `$ relates spatial and temporal scales through an algebraic relation $`l_\omega \omega ^{1/z}`$, where $`l_\omega `$ is the spatial extension of the elastic distortion resulting from a localized periodic shaking of pulsation $`\omega `$, and $`z`$ is the so-called dynamical exponent. The diffusive model described above by Eq. (4,5) indeed fits in this picture with $`z=2`$ and:

$$l_\omega =\left[(\gamma /\kappa )\omega \right]^{1/2}$$
(7)
## III Average Forces induced by three kinds of shaking

To set a reference for further comparison, we first recall the expression of the average thermal force $`F_{\mathrm{Casimir}}`$ on the plate at $`x=0`$ if the whole medium is thermally excited by a heat bath imposing a temperature $`T`$ kardarlee . The simplest derivation consists in computing the total free energy of the system for a given interplate distance $`L`$, and taking its derivative with respect to $`L`$. An alternative approach is to consider the model Langevin dynamics for the field given by Eq.(4,5), with $`\eta `$ a thermal Gaussian white noise of zero mean and of variance $`\eta (\stackrel{}{R},t)\eta (\stackrel{}{R}^{},t^{})=2\gamma k_\mathrm{B}T\delta (\stackrel{}{R}\stackrel{}{R}^{})\delta (tt^{})`$. Averaging over realizations of the noise leads to the same universal expression for the net thermal Casimir force aa :

$$F_{\mathrm{Casimir}}=A_dk_\mathrm{B}T\frac{L_{}^{d1}}{L^d},$$
(8)
with $`A_d=(d1)\mathrm{\Gamma }(d/2)\zeta (d)/(4\pi )^{d/2}`$. This universal mean force describes a long range attraction that depends on the temperature of the heat bath but not on specific material properties such as $`\kappa `$ and $`\gamma `$.

We now calculate the average forces generated by the three non-thermal shaking processes mentioned in the introduction, with a field $`\varphi `$ that evolves according to Eq. (4,5).

### III.1 Average force induced by a localized monochromatic shaking

We start with the simplest choice for the shaking in Eq. (4), namely a point-like monochromatic shaker oscillating at a pulsation $`\omega _\mathrm{s}`$, located at an equal distance $`L/2`$ from the two plates:

$$\eta (\stackrel{}{R},t)=\eta _\mathrm{s}a^d\mathrm{cos}(\omega _\mathrm{s}t)\delta (xL/2)\delta (𝐫),$$
(9)
The noise amplitude has been arbitrarily split into the product of a microscopic volume $`a^d`$ where the medium is shaken, times a shaking amplitude $`\eta _\mathrm{s}`$ that has the dimension of an energy density.

Such a process could be, for example, the result of the forcing of a microscopic ”active” particle located between the two plates (e.g. a magnetic particle under the influence of an external oscillating magnetic field), in which case $`a`$ is typically fixed whereas $`\eta _\mathrm{s}`$ and $`\omega _\mathrm{s}`$ can be externally tuned.

The net force $`F`$ felt by the plate at $`x=0`$, reduces to its right hand side contribution $`F^>`$ since the outer medium remains undistorted ($`\varphi (x<0)=0`$). Moreover the geometry is here totally symmetrical so that the force felt by the plate at $`x=L`$ is: $`F`$. Denoting $`F_\mathrm{A}^>`$ the time averaged force for this process, and anticipating on the detailed discussion of section IV-A, we use Eq. (23) to write:

$$F_\mathrm{A}^>=\frac{(\eta _\mathrm{s}a^d)^2}{\kappa }\frac{1}{L^{d1}}\mathrm{{\rm Y}}_A(\frac{1}{2},\frac{L}{l_{\omega _\mathrm{s}}}).$$
(10)
Again $`l_{\omega _\mathrm{s}}=[(\gamma /\kappa )\omega _\mathrm{s}]^{\frac{1}{2}}`$ is the dynamical length associated with the shaking pulsation $`\omega _\mathrm{s}`$. The adimensional function $`\mathrm{{\rm Y}}_\mathrm{A}(1/2,L/l_{\omega _\mathrm{s}})`$ is plotted in Fig. 3, and its general expression (including the shaker position dependence) is given below in Eq. (23).

Two different limit regimes can be identified: If $`L/l_{\omega _\mathrm{s}}1`$, the shaking period is much larger than the relaxation of any interplate field excitation, so that the elastic deformations generated at $`x=L/2`$ propagate along the $`x`$ axis up to the plates surfaces. In this quasistatic limit, the time averaged force $`F_\mathrm{A}^>`$ felt by the plate at $`x=0`$ is comparable to the elastic force on the plates under the application of a constant static “effective” perturbation $`\frac{1}{\sqrt{2}}\eta _\mathrm{s}a^d\delta (xL/2)\delta (𝐫)`$. The latter induces static deformations of lateral extension $`L`$ on the plate. Thus, $`F_\mathrm{A}^>`$ scales as $`\eta _\mathrm{s}^2a^{2d}/(\kappa L^{d1})`$ independently of the shaking period. In contrast, the dynamic length is much smaller than the plate-shaker distance if $`L/l_{\omega _\mathrm{s}}1`$. The elastic deformations are exponentially screened before reaching the plates surfaces. Therefore, the force amplitude decays exponentially as shown in Fig. 3.

In summary, a localized periodic shaking induces an effective interplate interaction. Whereas $`F_{\mathrm{Casimir}}`$ is universal and leads to an attraction between the two plates, the interaction obtained here is repulsive and does depend both on the elastic modulus of the surrounding“fluid” and on its dynamics (through $`l_{\omega _\mathrm{s}}`$). Finally, both the range and the amplitude of the present shaking-induced force can be externally tuned by varying the pulsation $`\omega _\mathrm{s}`$ and the amplitude $`\eta _\mathrm{s}`$ of the shaking. The dependence of the force on the shaker position as well as its fluctuations will be discussed in section IV.

### III.2 Average force induced by a uniform monochromatic shaking

We now suppose that the elastic medium is everywhere periodically and homogeneously shaken:

$$\eta (\stackrel{}{R},t)=\eta _\mathrm{s}\mathrm{cos}(\omega _\mathrm{s}t).$$
(11)
Posible experimental realizations include: the homogeneous and periodic shaking of a metallic elastic line using a magnetic field couder  or the driving of an interface between two dielectic liquids using AC electric fields xx . One may also use flexoelectric effects in liquid crystal to generate such a shaking process prostdegennes2 .

By symmetry the two plates feel opposite instantaneous forces. The exact expression for the time averaged force $`F_\mathrm{B}`$ on the $`x=0`$ plate is computed in section IV:

$$F_\mathrm{B}=\frac{\eta _\mathrm{s}^2L_{}^{d1}}{\kappa }l_{\omega _\mathrm{s}}^2\mathrm{{\rm Y}}_\mathrm{B}\left(\frac{L}{l_{\omega _\mathrm{s}}}\right).$$
(12)
The dimensionless quantity $`\mathrm{{\rm Y}}_\mathrm{B}`$, see Eq. (30), is plotted in Figure 4: it is non-monotonic and changes sign periodically. Given the plates geometry and the shaking translational invariance, the net normal stress $`F_\mathrm{B}/L_{}^{d1}`$ depends neither on the plates lateral extension $`L_{}`$ nor on the dimension $`d`$. The fast and slow shaking asymptotics can thus be qualitatively inferred from dimensional analysis. If the medium is slowly shaken ($`L/l_\omega 1`$), the net force amplitude diverges as $`l_{\omega _\mathrm{s}}^21/\omega _\mathrm{s}`$. First, notice that a constant uniform source $`\eta _\mathrm{s}`$ between the plates induces an homogenous normal stress on the inner side of the plate ($`x=0^+`$) that scales necessarily as $`F_\mathrm{B}^>/L_{}^{d1}\eta _\mathrm{s}^2\kappa ^1L^2`$. Second, on the left hand side of the plate ($`x<0`$), $`l_{\omega _\mathrm{s}}`$ is the only length scale available to construct the normal stress $`F_\mathrm{B}^</L_{}^{d1}\eta _\mathrm{s}^2\kappa ^1l_{\omega _\mathrm{s}}^2`$. In this slow shaking regime this latter contribution dominates, so that $`\mathrm{{\rm Y}}_\mathrm{B}`$ plateaus to a finite value. Turning now to the $`L/l_\omega 1`$ limit, if the interplate distance goes to infinity, the two sides of the plate face identical semi infinite media. Consequently, they are pushed in opposite directions with the same amplitude and thus feel on average no net force. So, in this limit $`L/l_\omega 1`$ that also corresponds to fast shaking at fixed $`L`$, $`F_\mathrm{B}`$ decays to zero, as observed in Fig. (4).

Returning to the most remarkable features, Fig. 4 and Eq. (12) show that the fluctuations of a uniformally shaken medium induce effective interplate forces that strongly differ from $`F_{\mathrm{Casimir}}`$. They are indeed completely tunable unlike their equilibrium analog. Their sign, amplitude and range can all be controlled externally by varying the shaking period and amplitude.

### III.3 Average force induced by a stochastic shaking: colored noise

Finally the case of stochastic shaking sources is considered. As mentioned above, if $`\eta `$ is the white Langevin noise modeling the coupling to a heat bath, the resulting field fluctuations induce thermal Casimir interactions between the plates. Any other noisy shaking drives the system out of equilibrium. Numerous examples of elastic materials shaken by “active” stochastic processes are provided by biological systems, e.g. active membranes pompes , actin-myosin gels gels , etc.

Without attempting to describe accurately the specifics of the noise in a given system, we focus here on the simple case of a spatially uncorrelated colored noise of zero mean, seeking for the influence of its color on the interplate force generated. We describe the fluctuations of this noise by:

$`\eta (\stackrel{}{R},t)\eta (\stackrel{}{R}^{},t^{})=\eta _\mathrm{s}^2a^dS(tt^{})\delta (\stackrel{}{R}\stackrel{}{R}^{}),`$ (13)
where the correlation $`S(t)`$ is defined by its Fourier transform $`S_\omega `$, also referred to as the noise power spectrum:

$$S_\omega =\frac{\mathrm{\Omega }_\mathrm{s}}{\mathrm{\Omega }_\mathrm{s}^2+\omega ^2}.$$
(14)
In Eq. (13) $`a`$ is a microscopic distance chosen so that the mean squared noise satisfies $`\eta ^2(\stackrel{}{R},t)=\eta _\mathrm{s}^2`$. Again $`\eta _\mathrm{s}`$ has the dimension of an energy density. The white noise limit corresponds to $`\mathrm{\Omega }_\mathrm{s}\mathrm{}`$, with $`\eta _\mathrm{s}^2a^d\mathrm{\Omega }_\mathrm{s}^1`$ kept constant.

The mean net force $`F_\mathrm{C}`$ on the $`x=0`$ plate is exactly computed in section IV-C:

$$F_\mathrm{C}=\frac{\eta _\mathrm{s}^2a^d}{\kappa }\frac{L_{}^{d1}}{L^{d2}}\mathrm{{\rm Y}}_\mathrm{C}\left(\frac{L}{l_{\mathrm{\Omega }_\mathrm{s}}}\right).$$
(15)
The dimensionless ratio $`\mathrm{{\rm Y}}_\mathrm{C}/A_d`$ is plotted in Fig. 5, using Eq. (IV.3). We recall that $`A_d`$ measures the amplitude of the thermal Casimir force ($`A_d=(d1)\mathrm{\Gamma }(d/2)\zeta (d)/(4\pi )^{d/2}`$).

Postponing a more involved discussion to next section, we emphasize here three points. First, the color of the noise has no influence on the sign of the force: the two plates attract each other whatever $`\mathrm{\Omega }_\mathrm{s}`$ and $`\eta _\mathrm{s}`$. Second, $`\mathrm{{\rm Y}}_\mathrm{C}`$ decreases monotonically with $`L/l_{\mathrm{\Omega }_\mathrm{s}}`$, in other words the larger the noise correlation time $`\mathrm{\Omega }_\mathrm{s}^1`$, the stronger the force felt by the two plates. Third, Fig. (5) indicates that $`\mathrm{{\rm Y}}_\mathrm{C}(u)/A_d`$ decays as $`1/(2u^2)`$ when $`u1`$. This implies that the attraction is long-ranged and decays only algebraically with the interplate distance. A connection with the Casimir result appears in this regime: if one sets $`\eta _\mathrm{s}^2a^d\mathrm{\Omega }_\mathrm{s}^1=2\gamma k_\mathrm{B}T`$, then in the limit white noise limit $`\mathrm{\Omega }_\mathrm{s}\mathrm{}`$, $`F_\mathrm{C}`$ converges as expected to $`F_{\mathrm{Casimir}}=A_dk_\mathrm{B}TL_{}^{d1}/L^d`$.

To conclude, we have shown that Casimir-like attractions can be induced by colored noisy shaking. These attractions decay algebraically at long distance as in the thermal case and contrary to the two former monochromatic cases. Their amplitude depend explicitly on the material properties and their scaling on the noise color.

## IV Shaking-Induced Forces : Detailed analysis.

In this section we provide a more detailed analysis of the interactions induced by the three kind of shakings: in addition to the explicit derivation of the average of the corresponding forces, we also comment on their temporal fluctuations.

We start with the derivation of the tool that will allow us to compute the forces on the $`x=0`$ plate in all situations, namely the linear response function $`_x`$ of the field’s gradient at the surface of the plate to a point-like shaker in the medium (by definition, the kernel $``$ is the solution of the dynamic equations (4,5), with $`\eta (\stackrel{}{R},t)=\delta (\stackrel{}{R}\stackrel{}{R}^{})\delta (tt^{})`$). We compute this quantity using an image method which conveys an intuitive picture MF : as in electrostatics, these equations can be solved replacing the constraint at the boundaries (5) by the addition of an ad hoc distribution of image sources outside the integration domain.

We consider first a shaking source located in the interplate region $`0<x^{}<L`$ and compute the response $`_x^>`$ on the $`x=0^+`$ side of the plate. In the present slab geometry, the images are the mirror images of the original source through virtual reflecting planes located at $`x=nL`$ as depicted in Fig. 6. Consequently, the response $`_{\omega ,𝐪}^>(x,x^{})`$ can be expressed as an infinite sum over all the images contributions:

$`_{\omega ,𝐪}^>(x,x^{})`$ $`=`$ $`{\displaystyle \underset{n}{}}_{\omega ,𝐪}\left(2nL+x^{}x\right)`$
$``$ $`{\displaystyle \underset{n}{}}_{\omega ,𝐪}\left(2nLx^{}x\right),`$
where $``$ the diffusion kernel in an infinite domain :

$`_{\omega ,𝐪}(x,x^{})={\displaystyle \frac{\mathrm{exp}\left(\sqrt{q^2+il_\omega ^2}\left|xx^{}\right|\right)}{2\kappa \sqrt{q^2+il_\omega ^2}}},`$ (17)
with $`\sqrt{i}1`$. Performing the sum (IV) and taking the derivative with respect to $`x`$ yields:

$$_x_{\omega ,𝐪}^>(0,x^{})=\frac{\mathrm{sinh}\left[\sqrt{(qL)^2+i(L/l_\omega )^2}(1x^{}/L)\right]}{\kappa \mathrm{sinh}\left[\sqrt{(qL)^2+i(L/l_\omega )^2}\right]}.$$
(18)
A simpler calculation can be carried out if the source is in the negative $`x`$ region ($`x^{}<0`$). In this case the only image source is the symmetric in the $`x=0`$ plane of the original source, so that the response of the field gradient on the $`x=0^{}`$ side of the plate is described by:

$$_x_{\omega ,𝐪}^<(0,x^{})=\frac{e^{\sqrt{(qx^{})^2+i(x^{}/l_\omega )^2}}}{\kappa }.$$
(19)
### IV.1 Localized Periodic Shaking.

If the only driving is a single monochromatic shaker located between the plates at $`x=x_\mathrm{s}`$,

$$\eta =\eta _\mathrm{s}a^d\delta (xx_\mathrm{s})\delta (𝐫)\mathrm{cos}(\omega _\mathrm{s}t),$$
(20)
then the general expression of the net force acting on the $`x=0`$ plate is simply obtained from the definition of the response function $`_x^>`$:

$`F(t)={\displaystyle \frac{\kappa (\eta _\mathrm{s}a^d)^2}{4}}\left[{\displaystyle _𝐪}\right|_x_{\omega _\mathrm{s},𝐪}^>(0,x_\mathrm{s})|^2]\times `$ (21)
$`\left[1+\mathrm{cos}(2\omega _\mathrm{s}t)\right],`$
with the short hand notation $`_𝐪𝑑𝐪/(2\pi )^{d1}`$ phase .

There is obviously no force on the left side of the plate, so that using (18) we obtain the exact expression for the instantaneous force on this plate. Its time average $`F_\mathrm{A}^>`$ describes a repulsion between the plates:

$$F_\mathrm{A}^>(x_\mathrm{s},\omega _\mathrm{s})=\frac{(\eta _\mathrm{s}a^d)^2}{\kappa L^{d1}}\mathrm{{\rm Y}}_\mathrm{A}(\frac{x_\mathrm{s}}{L},\frac{L}{l_{\omega _\mathrm{s}}})$$
(22)
with,

$$\mathrm{{\rm Y}}_\mathrm{A}(u,v)=\frac{1}{4}_𝐪\left|\frac{\mathrm{sinh}\left[\sqrt{q^2+iv^2}(1u)\right]}{\mathrm{sinh}(\sqrt{q^2+iv^2})}\right|^2.$$
(23)
We now give a closer look at the asymptotics corresponding to the fast and slow shaking limits.

Quasi-static shaking: We first consider the regime where the driving period is much larger than the relaxation time of any interplate field deformation, i.e. the dynamical length is much larger than the interplate distance ($`L/l_{\omega _\mathrm{s}}1`$).

The limiting value $`\mathrm{{\rm Y}}_\mathrm{A}(u,0)`$ can be exactly computed by integrating over the in plane $`𝐪`$ modes before summing up all the contributions from all images in Eqs. (IV,21). If $`d=3`$ $`\mathrm{{\rm Y}}_\mathrm{A}(u,0)=_u\left[\pi \left(1u\right)\mathrm{cot}\left(\pi u\right)2\mathrm{\Gamma }u+2\left(1u\right)\psi \left(u\right)\right]/(32\pi )`$, with $`\psi `$ the Euler $`\psi `$-function, and $`\mathrm{\Gamma }`$ the Euler constant gradstein , if $`d=2`$, $`\mathrm{{\rm Y}}_\mathrm{A}(u,0)=\left[1+\pi \left(1u\right)\mathrm{cot}(\pi u)\right]/(8\pi )`$, and if $`d=1`$, $`\mathrm{{\rm Y}}_\mathrm{A}(u,0)=(1u)^2/4`$. The scaling function $`\mathrm{{\rm Y}}_\mathrm{A}(x_\mathrm{s}/L,0)`$ is plotted in Figure 7 (log-log plot). The fast decay for values of $`x_\mathrm{s}/L`$ larger than $`\frac{1}{2}`$ witnesses that the elastic distortions are essentially occurring in the $`L/2<x<L`$ region. Conversely a shaker close to $`x=0`$ induces strong variations of the elastic field in the vicinity of the plate, leading to an algebraic divergence of the force when $`x_\mathrm{s}/L0`$ (if $`d>1`$).
The force dependence on the shaker position, can also be qualitatively understood thanks to an electrostatic analogy. In the present quasi-static limit, the kernel $`_{\omega _\mathrm{s}}^>`$ corresponds to the inverse of the Laplace operator (with DBC on the plates). Consequently $`\varphi `$ plays the role of an electrostatic potential, the plates mimic a grounded planar capacitance, and the shaking source a point like particle bearing a charge $`\eta _\mathrm{s}a^d`$. The force on the plate is then analogous to the electrostatic force produced by a point like charge jackson . From the image expansion (IV) we know that this force is identical to the one produced on a virtual plane at $`x=0`$ by an infinite number of $`(+\eta _\mathrm{s}a^d)`$ charges located at $`x=2nL+x_\mathrm{s}`$ and of $`(\eta _\mathrm{s}a^d)`$ charges at $`x=2nLx_\mathrm{s}`$. Hence:
–If $`x_\mathrm{s}<L`$ (c.f Fig. 6): the electric field on the plate surface is dominated by the image charges at $`|x|=x_\mathrm{s}`$, so $`F(\eta _\mathrm{s}a^d)^2/(\kappa x_\mathrm{s}^{d1})`$. Higher order reflections produce image sources that form dipoles viewed from $`x=0`$. Their subdominant contributions to the electrostatic force actually reduce the above estimation. This implies that the further the second plate, the bigger the force produced by a single shaker.
–if $`x_\mathrm{s}L`$ (c.f Fig. (6)): the image charges form dipoles located at $`x=(2n+1)L`$. The electric field on the plate is dominated by the influence of the dipole located at $`|x|=L`$ (the closest from the $`x=0`$ plate). The force therefore scales $`F_\mathrm{A}^>(\eta _\mathrm{s}a^d)^2(1x_\mathrm{s}/L)^2/(\kappa L^{d1})`$ .

Fast shaking:

If $`L/l_{\omega _\mathrm{s}}1`$, the electrostatic analogy only provides a partial description of the force generation. In this case, only the first reflexion contribute substantially to the image expansion (IV) as if the second plated (at $`x=L`$) was pushed to infinity. The only remaining length scales are $`x_\mathrm{s}`$ and $`l_{\omega _\mathrm{s}}`$ and the average force can be approximated by

$`F_\mathrm{A}^>{\displaystyle \frac{(\eta _\mathrm{s}a^d)^2}{\kappa x_\mathrm{s}^{d1}}}\stackrel{~}{\mathrm{{\rm Y}}}_\mathrm{A}(x_\mathrm{s}/l_{\omega _\mathrm{s}}),`$ (24)
with $`\stackrel{~}{\mathrm{{\rm Y}}}_\mathrm{A}(u)=\frac{1}{4}_𝐪\left|\mathrm{exp}(\sqrt{q^2+iu^2})\right|^2`$. This dimensionless screening factor is plotted in figure 8. If the dynamic length $`l_{\omega _\mathrm{s}}`$ is larger than the distance separating the shaking source from the plate, $`x_\mathrm{s}/l_{\omega _\mathrm{s}}1`$, the amplitude of the elastic field excitations are exponentially reduced before reaching the plate’s surface. More precisely, in this regime :

$`\stackrel{~}{\mathrm{{\rm Y}}}_\mathrm{A}(u)(u)^{\frac{d1}{2}}\mathrm{exp}(\sqrt{2}u).`$ (25)
Oppositely, $`\stackrel{~}{\mathrm{{\rm Y}}}_A(x_\mathrm{s}/l_{\omega _\mathrm{s}})`$ has a non vanishing finite limit if $`x_\mathrm{s}/l_{\omega _\mathrm{s}}`$ goes to zero. $`F_\mathrm{A}^>`$ then corresponds to the electrostatic force acting on a conducting plate facing a single “effective” charge $`\eta _\mathrm{s}a^d`$ located at $`(x=x_\mathrm{s},𝐫=0)`$, so that, as in the quasi-static case above, $`F_\mathrm{A}^>(\eta _\mathrm{s}a^d)^2/(\kappa x_\mathrm{s}^{d1})`$.

Fluctuations of the shaking-induced force: Beyond the definition and the characterization of the average interplate interaction, it is crucial to notice that the present shaking induced force is a time varying quantity. Rewriting (21) as $`F(t)=F_\mathrm{A}^>\left[1+\mathrm{cos}(2\omega _\mathrm{s}t)\right]`$, it is obvious that the force experienced by the immersed plates can not be a priori reduced to its sole average value. As a matter of fact the fluctuations of the instantaneous force are comparable to its average amplitude.

### IV.2 Homogenous Monochromatic Shaking.

We now consider that the soft medium is uniformally shaken by $`\eta (\stackrel{}{R},t)=\eta _\mathrm{s}\mathrm{cos}(\omega _\mathrm{s}t)`$. $`\eta `$ is the superposition of in phase shakers of the –A– type uniformally distributed in space. All the elementary shakers produce periodic elastic deformations that “interfere” and give rise to the instantaneous forces $`F^>(t)`$ (resp. $`F^<(t)`$) at $`x=0^+`$ (resp. $`x=0^{}`$) obtained from Eqs. (2,6):

$`F^>(t)={\displaystyle \frac{\kappa }{4}}L_{}^{d1}|{\displaystyle _0^L}_x_{\omega _\mathrm{s},\mathrm{𝟎}}^>(0,x^{})dx^{}|^2\times `$ (26)
$`\left[1+\mathrm{cos}(2\omega _\mathrm{s}t+\theta ^>)\right],`$
with $`\theta ^>=2\mathrm{arg}\left[_0^L_x_{\omega _\mathrm{s},\mathrm{𝟎}}^>(0,x^{})dx^{}\right]`$, and

$`F^<(t)={\displaystyle \frac{\kappa }{4}}L_{}^{d1}|{\displaystyle _0^{\mathrm{}}}_x_{\omega _\mathrm{s},\mathrm{𝟎}}^>(0,x^{})dx^{}|^2\times `$ (27)
$`\left[1+\mathrm{cos}(2\omega _\mathrm{s}t+\theta ^<)\right],`$
with $`\theta ^<=2\mathrm{arg}\left[_0^{\mathrm{}}_x_{\omega _\mathrm{s},\mathrm{𝟎}}^<(0,x^{})dx^{}\right]`$. At every instant these two forces compete. They push the plate in opposite directions to give rise to the net force $`F(t)`$. The latter is exactly computed using Eqs. (3,18,19,26,27) , and after some elementary algebra phase :

$`F(t)`$ $`=`$ $`F_\mathrm{B}\left[1+{\displaystyle \frac{\mathrm{cos}(2\omega _\mathrm{s}t)}{\mathrm{cos}\left[L/(\sqrt{2}l_{\omega _\mathrm{s}})\right]}}\right].`$ (28)
The average $`F_\mathrm{B}`$ is given by:

$`F_\mathrm{B}`$ $`=`$ $`{\displaystyle \frac{\eta _\mathrm{s}^2L_{}^{d1}}{\kappa }}l_{\omega _\mathrm{s}}^2\mathrm{{\rm Y}}_\mathrm{B}\left({\displaystyle \frac{L}{l_{\omega _\mathrm{s}}}}\right),`$ (29)
$`\mathrm{{\rm Y}}_\mathrm{B}(u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{cos}(u/\sqrt{2})}{\mathrm{cos}(u/\sqrt{2})+\mathrm{cosh}(u/\sqrt{2})}}.`$ (30)
The function $`\mathrm{{\rm Y}}_\mathrm{B}`$ is plotted in section III-B, see Fig. 4, where its fast and slow asymptotics are qualitatively studied.

Mean force tunability: The tunability of the present shaking-induced force is clearly revealed by Eqs. (29) and (30). Precisely for fixed $`L`$: (i) The shaking pulsation $`\omega _\mathrm{s}`$ sets the sign of $`F_\mathrm{B}`$ to $`\mathrm{sign}(\mathrm{cos}[L/(\sqrt{2}l_{\omega _\mathrm{s}})])`$, and its range to $`l_{\omega _\mathrm{s}}`$ (via the $`\mathrm{cosh}`$ term in (29)),(ii) then $`\eta _\mathrm{s}`$ allows to adjust the mean force amplitude arbitrarily. This second assertion is actually true only if $`L\frac{\pi }{\sqrt{2}}(2n+1)l_{\omega _\mathrm{s}}`$, with $`n`$ integer. In such cases the plates experience no mean net force.

Fluctuations of the shaking-induced force: Eq. (28) indicates actually that $`F(t)`$ “fluctuates” periodically around $`F_\mathrm{B}`$ with a pulsation $`2\omega _\mathrm{s}`$. Moreover the amplitude of the $`2\omega _\mathrm{s}`$ components of the force is $`1/\mathrm{cos}[L/(l_{\omega _\mathrm{s}}\sqrt{2})]`$ times larger that its mean value $`F_\mathrm{B}`$ whatever the shaking parameters $`\omega _\mathrm{s}`$ and $`\eta _\mathrm{s}`$.

Shaking-induced organization: The two preceeding comments strongly suggest that the homogenous shaking may be used to dynamically induce spacial “localization” of the plates. As a matter of fact, when mesoscopic objects are immersed in a soft medium, their intrinsic dynamics (ignored in the present work) necessarily filters out high frequency effects of external forcing. If their response time is much larger than $`1/(2\omega _\mathrm{s})`$, these “slow” moving objects only respond to the stationary component $`F_\mathrm{B}`$ of the shaking-induced force. The corresponding steady state interplate distance in the present simple geometry is then controllable by tuning $`\omega _\mathrm{s}`$.

### IV.3 Colored Noisy Shaking.

The medium surrounding the plates is now supposed to be driven by stochastic shaking sources with zero mean and mean squared fluctuations given by Eq. (13). The noise averaged net force $`F_\mathrm{C}`$ experienced by the plate at $`x=0`$ is the sum of the mean values of the two independent forces $`F^>`$ and $`F^<`$. Using Eqs. (12,18,20) this force is written:

$`F_\mathrm{C}={\displaystyle \frac{2L_{}^{d1}}{a^d}}{\displaystyle _\omega }S_\omega [{\displaystyle _0^L}(F_\mathrm{A}^>(x,\omega )+F_\mathrm{A}^<(x,\omega ))dx`$ (31)
$`+{\displaystyle _L^{\mathrm{}}}F_\mathrm{A}^<(x,\omega )dx].`$
Exact results: For a Lorentzian power spectrum, see Eq. (14), the integration over frequency and position of all the elementary contributions $`F_\mathrm{A}(x,\omega )`$ in (31) yields:

$`F_\mathrm{C}`$ $`=`$ $`{\displaystyle \frac{\eta _\mathrm{s}^2a^dL_{}^{d1}}{\kappa L^{d2}}}\mathrm{{\rm Y}}_\mathrm{C}\left({\displaystyle \frac{L}{l_{\mathrm{\Omega }_\mathrm{s}}}}\right),`$ (32)
$`\mathrm{{\rm Y}}_\mathrm{C}\left(u\right)`$ $`=`$ $`{\displaystyle \frac{1}{4u^2}}{\displaystyle _𝐪}q\left(\mathrm{coth}q1\right)`$
$``$ $`{\displaystyle \frac{1}{4u^2}}{\displaystyle _𝐪}\sqrt{q^2+u^2}\left[\mathrm{coth}\left(\sqrt{q^2+u^2}\right)1\right].`$
The scaling function $`\mathrm{{\rm Y}}_\mathrm{C}`$ is plotted in Fig. 5. To achieve the exact computation of $`F_\mathrm{C}`$, the integral over the in-plates $`𝐪`$ modes in (IV.3) is performed. After some algebra we obtain: If $`d=1`$, $`\mathrm{{\rm Y}}_\mathrm{C}=(1u[\mathrm{coth}(u)1])/u^2`$, if $`d=3`$, $`\mathrm{{\rm Y}}_\mathrm{C}=[3\zeta (3)\pi ^26(u^2u\mathrm{log}(1e^{2u})+3\mathrm{L}\mathrm{i}_2(e^{2u})]/(48\pi u^2)`$, where $`\zeta `$ is the Rieman Zeta function and $`\mathrm{Li}_2`$ the dilogarithm function : $`\mathrm{Li}_2(x)=_{k=1}^{\mathrm{}}x^k/k^2`$. Finally for the $`d=2`$ case $`\mathrm{{\rm Y}}_\mathrm{C}=2_n[K_0(2nu)+K_2(2nu)]`$, where $`K_n`$ is the $`n^{th}`$ K-Bessel function gradstein .

Notice that this last calculation can actually be bypassed, if one already knows the thermal Casimir force the left plate would feel if the soft medium had a finite correlation length $`\xi `$, precisely if the elastic Hamiltonian were: $`=(\kappa /2)[\stackrel{}{}\varphi (\stackrel{}{R})]^2+[\varphi (\stackrel{}{R})/\xi ]^2d^d\stackrel{}{R}`$. Denoting this force $`\stackrel{~}{F}_{\mathrm{Casimir}}(L/\xi )`$, we can show that Eqs. (32,IV.3) can be recast into bafg :

$$F_\mathrm{C}=\frac{\eta _\mathrm{s}^2a^d\mathrm{\Omega }_\mathrm{s}^1}{2\gamma k_\mathrm{B}T}\left[\stackrel{~}{F}_{\mathrm{Casimir}}\left(0\right)\stackrel{~}{F}_{\mathrm{Casimir}}\left(\frac{L}{l_{\mathrm{\Omega }_\mathrm{s}}}\right)\right],$$
(34)
where $`\stackrel{~}{F}_{\mathrm{Casimir}}\left(0\right)`$, of course corresponds to the Casimir force the plate experiences in a scale-free fluctuating medium, c.f. Eq. (8).

Qualitative approach: Eq. (31) shows clearly that $`F_\mathrm{C}`$ is the result of the incoherent sum over frequencies and positions of the forces produced by localized shakers of the –A– type. Note that the sum over the shaking frequencies is weighted by the noise power spectrum $`S_\omega `$. This decomposition is now used to justify first the sign of the mean force, then the algebraic decays of the mean force as observed in Eq.(32) and Fig. (5).

Only the elementary monochromatic shakers localized at distances $`|x|<L`$ on the left hand side of the plate have a challenger in the $`0<x<L`$ inner region. Referring to the discussion in Section IV-A, we know that $`F_\mathrm{A}^<>0`$ and that $`\left|F_\mathrm{A}^>(x,\omega )\right|<\left|F_\mathrm{A}^<(x,\omega )\right|`$. Moreover, $`S_\omega `$ is necessarily a positive quantity for all frequencies. It then turns out that the net mean force is positive, for any choice of $`S_\omega `$: the two objects attract each other. Besides, $`F_\mathrm{A}^>(x,\omega )F_\mathrm{A}^<(x,\omega )`$ for shaker located at $`x<L/2`$. This important property leads to the cancellation of the short distance divergences in the first term of equation (31). Finally, $`\left|F_A^>(x,\omega )\right|\left|F_A^<(x,\omega )\right|`$ if $`|x|>L/2`$, so that Eq. (31) can be approximated by $`F_\mathrm{C}2L_{}^{d1}a^d_\omega S_\omega _{L/2}^{\mathrm{}}F_A^<(x,\omega )𝑑x`$.

To sum up, the elastic field deformations produced by the –A– type shakers located at $`x<L/2`$ dominate the mean net force $`F_\mathrm{C}`$.

Without refering to any specific choice for $`S_\omega `$, the noise power spectrum is here characterized only by its width $`\mathrm{\Omega }_\mathrm{s}`$, see e.g. Eq. (13). Thus, going to the $`(x,l_\omega )`$ variables and defining $`L^{}\mathrm{max}(L,l_{\mathrm{\Omega }_\mathrm{s}})`$ we easily deduce:

$$F_\mathrm{C}\frac{\eta _\mathrm{s}^2a^d\mathrm{\Omega }_\mathrm{s}^1L_{}^{d1}}{\gamma }_L^{}^{\mathrm{}}\frac{dl_\omega }{l_\omega ^3}_L^{l_\omega }\frac{dx}{x^{d1}},$$
(35)
since $`F_\mathrm{A}^<(x,\omega )`$ is exponentially weak if $`x>l_\omega `$ and scales as $`1/x^{d1}`$ if $`x<l_\omega `$, c.f. Fig. (9).

We can now integrate (35) and distinguish the two cases:

(i) $`L/l_{\mathrm{\Omega }_\mathrm{s}}1`$. In this wide noise spectrum limit $`F_\mathrm{C}(\eta _\mathrm{s}^2a^d\mathrm{\Omega }_\mathrm{s}^1/\gamma )L_{}^{d1}/L^d`$, the thermal Casimir force scaling is obtained as anticipated in section III.

(ii) $`L/l_{\mathrm{\Omega }_\mathrm{s}}1`$. In this case the noise correlation time $`\mathrm{\Omega }_\mathrm{s}^1`$ is larger than the field relaxation over distances smaller than L. The number of elementary shakers contributing to the force is substantially increased compared to (i). Depending on the space dimension $`d`$, we obtained the three scaling forms: if $`d=3`$, $`F_\mathrm{C}(\eta _\mathrm{s}^2a^3\kappa ^1)L_{}^2/L`$, if $`d=2`$, $`F_\mathrm{C}(\eta _\mathrm{s}^2a^2\kappa ^1)L_{}\mathrm{log}(l_{\mathrm{\Omega }_\mathrm{s}}/L)`$, an finally if $`d=1`$, $`F_\mathrm{C}(\eta _\mathrm{s}^2a\kappa ^1)l_{\mathrm{\Omega }_\mathrm{s}}`$. In agreement with the above exact calculations and with Fig. 5.

Fluctuations of the shaking-induced force: To conclude this section, we address briefly the question of the fluctuations of the force. In order to asses their relative importance with respect to the mean force value, we introduce the dimensionless ratio: $`F_\mathrm{C}/\mathrm{\Delta }F`$, where the mean squared deviation $`\mathrm{\Delta }F^2`$ is defined by: $`\mathrm{\Delta }F^2F(t)^2F_\mathrm{C}^2`$. Whereas $`F_\mathrm{C}`$ is the result of the competition between the mean forces $`F_\mathrm{C}^>`$ and $`F_\mathrm{C}^<`$ of opposite signs, the forces’ fluctuations on the two sides add up and $`\mathrm{\Delta }F^2=(\mathrm{\Delta }F^>)^2+(\mathrm{\Delta }F^<)^2`$ with obvious notations. Note that the noise’s four points correlation function is a priori required to compute $`\mathrm{\Delta }F`$. We make here one more assumption and consider that the $`\eta (\stackrel{}{R},t)`$ are Gaussian random variables with correlations given by Eq. (13,14). This allows to fully characterize the shaking process and consequently the variance of the shaking-induced force. For sake of clarity we present without additional details the scaling form of the mean over variance ratio (ignoring logarithmic corrections):

$$\frac{F_\mathrm{C}}{\mathrm{\Delta }F}\left(\frac{L_{}}{L}\right)^{\frac{d1}{2}}\{\begin{array}{cc}\left(\frac{L}{l_{\mathrm{\Omega }_\mathrm{s}}}\right)^{\frac{3d}{2}}\hfill & \mathrm{if}L/l_{\mathrm{\Omega }_\mathrm{s}}1\hfill \\ \left(\frac{l_{\mathrm{\Omega }_\mathrm{s}}}{L}\right)^{\frac{d+1}{2}}\hfill & \mathrm{if}L/l_{\mathrm{\Omega }_\mathrm{s}}1\hfill \end{array}$$
(36)
we have also assumed $`d>1`$, since the 1D case deserve a more careful analysis as in the equilibrium context more carefully studied in bafg . The above expressions shows that, the present shaking induced force can not be a priori reduced to its sole $`F_\mathrm{C}`$ value $`F_\mathrm{C}`$. Indeed, depending on the geometrical aspect ratio $`L_{}/L`$ and on the noise correlation time $`\mathrm{\Omega }_\mathrm{s}^1`$ the variance of the force can dominate its mean value by orders of magnitude.

### IV.4 Toward more complex dynamics

To what extent do the results presented above apply to systems evolving according to more complex dynamics? It is noticeable that only the explicit form of the scaling functions $`\mathrm{{\rm Y}}_\mathrm{X}`$ requires a precise description of the field dynamics. Conversely, all the qualitative analysis, as well as the asymptotic expressions for the various forces, involve only the definition of the dynamical length $`l_\omega `$. We thus expect that they should be extendable to any other dynamic scheme relating algebraically the spatial and the temporal relaxation scale.

## V Summary and Outlook

To summarize, a simple model has allowed us to show that the concept of fluctuation-induced interactions can be extended to complex fluids driven to out of equilibrium. Our main result concern identical plate-like inclusions immersed in an homogenous medium externally driven by :–A– a monochromatic localized shaking source located between the plates, –B– a monochromatic and homogenous shaking and –C– a noisy colored shaking. These external processes generate effective forces on the plates without acting directly on them. The main features of the mean forces induced by these three simple shaking are summarized in the table below and compared to their equilibrium thermal Casimir analog.

| Shaking | Amplitude | Range | Sign |
| --- | --- | --- | --- |
| A | $`\eta _\mathrm{s}^2`$ | $`l_{\omega _\mathrm{s}}`$ | repulsion |
| B | $`\eta _\mathrm{s}^2`$ | $`l_{\omega _\mathrm{s}}`$ | tunable |
| C | $`\eta _\mathrm{s}^2`$ | power law | attraction |
| Thermal | $`k_\mathrm{B}T`$ | power law | attraction |

Five points are worth a highlight:

(i) Whereas a shaking at a single frequency leads necessarily to short range interactions, noisy shakings generate effective forces that decays algebraically with the interplate distance.

(ii) If the external shaking occurs only between the two plates the resulting interaction is always repulsive. Conversely, when distributed on the whole sample the external shakers generate forces on both sides of the plates. The resulting net force is attractive for spatially uncorrelated shaking. On the contrary, with coherent shaking, interference phenomena between the elastic medium distortions can lead to attractive or repulsive interactions.

(iii) The soft medium-plates coupling has been here modeled by Dirichlet boundary conditions. As in the equilibrium case our results are expected to hold for other static boundary conditions representing strong interactions. In contrast the description of situations of weaker coupling and/or cases where dynamic boundary conditions are required could yield to new phenomenology.

(iv) The generalization of our results to an elastic medium characterized by an intrinsic relaxation scale $`\xi `$ (see the discussion above Eq. (34)) can be done without much effort following blindly the procedures described in this paper. The only qualitative difference in the results summarized in the above table is the forces’ range: For shakings of the –A– and –B– type the range of the force becomes: $`\mathrm{min}(\xi ,l_\omega )`$. In the noisy shaking case –C– the force decay is not algebraic anymore but exponential for interplate distances larger than $`\xi `$.

(v) Beyond the derivation of their average values, we have also shown that the shaking-induced forces are in general expected to be strongly fluctuating quantities, whatever their precise origin.

To conclude, let us recall that our description is clearly extremely simplified and of course is not supposed to model any real system. Even within this simple framework the extension of our result to many moving objects remains complex. But, the analysis of simplified situations such as the one described here could certainly help understanding ordering and dynamic behavior of inclusions in out-of equilibrium complex fluids.


format: text

text:
# Microlensing limits on numbers and orbits of extra-solar planets from the 1998-2000 OGLE events.

## 1 INTRODUCTION

In 1995, very precise radial velocity (Doppler) measurements resulted in the first detection of an extra-solar planet orbiting a main-sequence star (?). This discovery spawned a series of campaigns aimed at finding new exo-planets, placing limits on their abundance and eventually understanding their formation processes.

A variety of techniques is currently being used to search for planetary signatures. The list includes pulsar timing (?), transit searches (?) and gravitational microlensing (??). This paper considers the latter.

Photons emitted by a background stellar source are deflected as they come near the influence of the gravitational field of a massive foreground object which acts as a gravitational lens. This results in two images of the source, one on either side of the lens. If the source, lens and observer are perfectly aligned, we observe a ring image around the lens whose radius is defined as the Einstein ring radius and is characteristic of each lens.

The Einstein ring radius is given by:

$$\theta _E=R_E/D_d=\sqrt{\frac{D_{ds}}{D_sD_d}\frac{4GM(\theta _E)}{c^2}},$$
(1)
where $`D_{ds}`$, $`D_d`$ and $`D_s`$ are the lens-source, observer-lens and observer-source distances respectively. The Einstein ring provides a natural angular scale to describe the lensing geometry for several reasons. If multiple images of the source are produced, the typical angular separation of the images is of the order of 2$`\theta _E`$. Furthermore, sources that are closer than about $`\theta _E`$ to the optical axis experience strong lensing and are therefore significantly amplified. On the other hand, sources which are located well outside this ring receive very little amplification.

In microlensing, the separation of the images created by the lensing effect ($`\theta 10^3`$arcsec) is too small to be resolved by current telescopes and one can observe only the combined flux. The resulting lightcurve is symmetric in time with its maximum amplification $`A_0`$ at the time of closest approach $`t_0`$ between the projected position of the source on the lens plane and the lens itself. Its shape is described by:

$$A(u)=\frac{u^2+2}{u(u^2+4)^{1/2}},$$
(2)
where $`A(u)`$ is the total amplification and $`u`$ is the unlensed angular separation of source and lens in units of the angular Einstein ring radius $`\theta _E`$ (?):

$$u=\frac{\theta }{\theta _E}=\left[u_0^2+\left(\frac{2(tt_0)}{t_E}\right)\right]^{1/2},$$
(3)
where $`t_E=2\theta _ED_d/v_{}`$ and $`v_{}`$ is the transverse velocity between the source and lens.

This simple point-source point-lens (PSPL) model can describe most of the cases well enough. However, there are cases when the shape of the lightcurve is not symmetric and exhibits significant deviations. These so-called anomalies of the lightcurve can be due to several factors and have been extensively examined in recent literature (?????). The most interesting of these are the anomalies which can be attributed to the binary nature of the lens. The possibility that the binary lens system may be a star-planet system of extreme mass ratio, has spawned dedicated observing campaigns to reveal such planets (???).

In section 2 we discuss our sample of events and fitting method. Section 3 deals with the notion of ‘detection zones’. These are zones on a $`\mathrm{\Delta }\chi ^2`$ map as a function of planet position where the presence of a planet is excluded by our analysis. In section 4 we explain our treatment of the fit residuals and discuss what threshold value we should adopt for the analysis. The best planet candidates are presented in section 5, followed by a discussion on planet detection probabilities in section 6. In section 7 we provide the theoretical background to getting limits on the number of planets per star and we conclude with a summary of the analysis in section 8.

## 2 ANALYSIS OF OGLE LIGHTCURVES

### 2.1 THE SAMPLE OF EVENTS

The datasets presented here were obtained from the DoPhot analysis of the events available at the OGLE website. For the purposes of this analysis we only considered the events for which the PSPL model provided adequate fits.

Out of a total of 162 events observed by the Optical Gravitational Lensing Experiment (?) collaboration during a three year period (1998-2000) we found that for 145 of these the reduced $`\chi ^2`$ of the PSPL model fit, adjusting 6 parameters, approached a value of unity indicating a good fit. We discuss our method of fitting the lightcurves in the next section.

### 2.2 POINT-SOURCE POINT-LENS FITS

In comparison with our previous work (?), for this paper we have improved our numerical treatment of the PSPL model by replacing the approximate (?) formalism valid only for extreme mass ratios, $`q<<1`$, with a binary lens model valid for an arbitrary mass ratio. Our neglect of finite-source effects (?) renders our PSPL lightcurves inaccurate for extreme mass ratios ($`qu_s^2`$), and high amplification ($`Au_s^1`$), where $`u_s`$ is the source size in units of the Einstein ring. In addition, for the discussion in Section 2.4 we allow the lensed source to be blended with a constant source, $`A(t)A(t)+b`$, with $`b>0`$, so that $`f_bb/(1+b)`$, with $`0<f_b<1`$, is the fraction of the apparent unlensed source flux contributed by nearby constant sources.

To illustrate our analysis we present our fits for events 1998bul14 and 2000bul31 in Figure 1. These are high magnification events and the fitted event parameters are indicated on the top left corners of the plots. Our estimates of the parameters for the event 1998bul14 are in agreement with those published by the PLANET collaboration (?).

Our initial fits were accomplished by adjusting four parameters describing the shape of the PSPL lightcurve, the time of maximum amplification $`t_0`$, event timescale $`t_E`$ (time to cross Einstein ring diameter), maximum amplification $`A_0`$ and the baseline magnitude $`I_0`$. The residuals of the 4-parameter fits were generally consistent with the estimated error bars during the bright phases, but often larger than expected near the baseline magnitude $`I_0`$. This probably reflects the difficulty in obtaining accurate photometric measurements in crowded fields (see section 2.3).

To assess evidence of lightcurve deviations at a first glance, we impose a criterion of $`\chi ^2/(N4)<3.5`$, for the 4-parameter fit, which helps to identify the obvious non-PSPL lightcurves ($`N`$ is the number of data points fitted). This criterion was not met by 46 events. Upon further examination of these lightcurves, we concluded that for 17 of these, the effects of stellar binary lensing, lensing of binary sources, parallax motion and stellar variability were responsible for anomalies in the lightcurve (see appendix A). We removed these 17 events from our list since our model does not account for these effects. The remaining 29 events were kept for further analysis despite the inferior quality of the lightcurves. These were events that were either sampled only during one part of the amplification phase, were imaged very sparsely for the duration of the amplification and/or showed significant scatter in the baseline.

### 2.3 ADJUSTMENT OF THE ERROR BARS

Since our assessment of the evidence for lensing by planets is based on the significance of the residuals to the fitted PSPL lightcurve, we cannot ignore the systematic discrepancy between the fit residuals and the published error bars. When we simply scale the reported error-bars by a factor $`f`$, we find that the residuals are often larger than expected in the flat part of the lightcurve, and smaller than expected in the amplified parts of the lightcurve. This suggests that an additive flux error, in addition to an overall scaling of the error bars, may be needed to account for the rms residuals in all parts of the lightcurve. To deal with this we decided to include a crowded field error in our model, $`\sigma _0`$. This parameter accounts for increased scatter in the measurements of the unlensed flux. Thus, if blending and variable seeing effects cause the baseline flux to vary by more than the theoretical limit from photon-counting noise due to the background photons, the increase in $`\sigma _0`$ allows for that extra variance. We also scale the reported error bars by a factor $`f`$. Thus the error bar $`s_i`$ on the flux assumes the form:

$$s_i=\left(\sigma _{0}^{}{}_{}{}^{2}+f^2\sigma _{i}^{}{}_{}{}^{2}\right)^{1/2},$$
(4)
where $`\sigma _0`$ is the additive flux error intended to account for the crowded field effects and $`\sigma _i`$ is the original error bar. The initial estimates for $`\sigma _0`$ and $`f`$ are 0 and 1 respectively. This 2-parameter noise model gives a more satisfactory description of the residuals for all events.

With these two extra parameters to adjust the size of the error bars, we cannot use $`\chi ^2`$ minimization to optimize the fit since the fit can achieve $`\chi ^2=0`$ by making the error bars infinitely large. Therefore we use a maximum likelihood criterion. Assuming Gaussian error distributions, maximizing the likelihood is equivalent to minimizing

$$\chi ^2+2\underset{i=1}{\overset{N}{}}\text{ln}(s_i).$$
(5)
This introduces an appropriate penalty for making the error bars large since that increases the value of the second term in the equation.

Tables 1, 2, 3, 4 summarise the best-fit parameters we obtained by maximizing the likelihood of the data over 6 parameters. There are three columns in the Tables that give the reduced $`\chi ^2`$ values for four, five and six parameters respectively, the first being evaluated with the published error bars, the latter using the method outlined above to adjust them. The last method represents a substantial improvement with $`\chi ^2/(N6)1`$.

### 2.4 BLENDING

An effect that we expect to be affecting all lightcurves to some extent is blending. Blending is common in the photometry of crowded fields such as the Galactic Bulge. Its influence might lead to misinterpretation of the baseline magnitude of the source and thus inaccurate estimation of the true maximum amplification and timescale of the lensing event. As blended events may be chromatic, multi-band photometry can help estimate this effect (??).

We have performed fits to the events using 6 (no blending) and 7 (with blending) parameters. We chose to use the results of the 6 parameter fits in this paper. For some individual cases we also consider the 7 parameter results (also see ?). In a poorly sampled lightcurve, a few points high or low can make the best fit indicate substantial blending even when there is no real evidence for it. Including a blending flux in the model results in larger values for the $`A_0`$ and $`t_E`$ parameters. This effect boosts the detection probabilities. Since we do not want to artificially enhance these, we choose the 6 parameter model thereby underestimating the probability of detection.

For all 145 events, we calculate the ratio of amplifications ($`A_0(7)/A_0(6)`$) and event timescales ($`t_E(7)/t_E(6)`$) obtained by fitting the data with (7 parameter fit) and without (6 parameter fit) allowing for a blend fraction. The results are plotted in Figures 3(a) and 3(b). For high blend fractions $`f_b>10`$% the effects are dramatic on the values of $`A_0`$ and $`t_E`$. We estimate that $`23`$% of our events are high amplification events ($`A_0>5`$) that suffer from serious blending ($`f_b>10`$%). By not modeling this effect we are underestimating the true values of $`A_0`$ and $`t_E`$ and thus the planet detection probability for the affected events. It follows that the limits we are quoting in section 6 for the total detection probability are lower limits. The underestimate of the detection probability by not including blending in the analysis is of the order of $`13`$%. We present a selection of events fitted with and without a blend fraction in Table 5.

Some of the events from the 1998-1999 seasons have been re-identified and re-analyzed using the method of difference image analysis (DIA) by the OGLE team (?). Figure 4(a) shows how our derived values for the maximum amplification (using the 6 parameter fit) compare with the difference image analysis results. Figure 4(b) compares the values for $`A_0`$ derived using the 6 and 7 parameter fits. Both the DIA analysis and the 7 parameter fits (including blending) return a higher value for the maximum amplification.

## 3 PLANET DETECTION ZONES

Microlens event monitoring in the OGLE dataset is not consistently intensive enough to detect deviations due to small planets, having sampling gaps of many hours, occasionally days. Therefore any evidence of a planetary perturbation will be confined to very few data points. In the case that a candidate deviation is noted, the corresponding frame(s) must be thoroughly examined since cosmic ray hits affecting the CCD data might be mistaken for a planetary signature. In this situation, only a very significant departure from the PSPL lightcurve could be accepted as evidence of a planet.

Figure 2 presents a $`\mathrm{\Delta }\chi ^2`$ map as a function of planet position for events 1998bul14 and 2000bul31, and for a planet/star mass ratio $`q=10^3`$. These show how the $`\chi ^2`$ of the fit changes from the PSPL fit when we add a planet at different places on the lens plane. The gray zones on the plots indicate regions where the $`\mathrm{\Delta }\chi ^2`$ is zero, i.e. placing the planet here does not affect the lightcurve in any way. White zones are negative detection zones, or exclusion zones, which show regions of $`\mathrm{\Delta }\chi ^2`$ $`+`$25. These are the regions where the planet can be excluded at the time of observation since at these places the planet is close enough to one of the images to cause a large perturbation of the lightcurve near one of the data points. Note that the detection zones close to the Einstein ring of the lens are larger because the planet is perturbing a more highly amplified image. The detection zones also depend on the size of the error bars, with smaller error bars resulting in larger detection zones.

The appearance of black zones on the plots would signify that a better fit has been achieved with a planet at that position interacting with one of the two images of the source star ($`\mathrm{\Delta }\chi ^2`$ $``$25). The black zones are positive detection zones. If a lightcurve provides evidence for a planet, the $`\mathrm{\Delta }\chi ^2`$ map will have a clear black spot corresponding to the planet’s position. With a poorly sampled lightcurve we expect two black spots appearing at both possible image positions since we cannot tell which image the planet is interacting with. With good sampling, a black spot would show only at the correct position of the planet at the time of observation.

We have examined such $`\mathrm{\Delta }\chi ^2`$ maps for all 145 events and find that for the majority of the events there is no compelling evidence for a planet. However, each data point in each lightcurve rules out the presence of a planet near the major and minor image positions. Based on this method, we quantify in section 6 our non-detection results by calculating the detection probability for each event for two different mass ratios.

## 4 ANALYSIS OF RESIDUALS

The 3-year OGLE dataset has 145 lightcurves with $`40`$ data points each (within the range $`3<a/R_E<3`$) for which the PSPL model provides an adequate fit. This gives a total of $`N=5840`$ data points. If we consider $`N`$ data points with Gaussian errors, the largest residual (largest false alarm) is approximately

$$S_F(N)2.14\sqrt{\mathrm{log}_{10}N},$$
(6)
in units of $`\sigma `$. We therefore expect the largest single point residual in the 3-year OGLE microlensing dataset to be approximately $`S_F4.15\sigma `$. If the data point with the largest residual is well separated in time from its neighbors, in comparison with the duration of the planet lens anomaly, then it should be possible for the planet lens anomaly to fit the largest residual almost exactly without being held back by other nearby data points. The largest reduction in $`\chi ^2`$ when a planet is included in the lens model should therefore be approximately $`S_{F}^{}{}_{}{}^{2}(4.15)^217`$.

The $`\mathrm{\Delta }\chi ^2`$ distribution we obtain from our fits to the 145 lightcurves with star+planet models for $`q=10^3`$ is shown in Figure 5. The median $`\mathrm{\Delta }\chi ^2`$ value is 6.27 and the largest is 75.67. The latter is much larger than the expected value found above. We make a Monte-Carlo estimate for the theoretical distribution by generating $`N_{dp}`$ standard Gaussian errors, where $`N_{dp}`$ is the number of data points in a light curve in the range $`3<x/R_E<3`$, and taking the largest outlier $`S_F(N_{dp})`$ for each of the lightcurves. We then plot the values of $`[S_F(N_{dp})]^2`$ and compare them with the observed distribution in Figure 5. Although we have already corrected the error bars as described in section 2.3, the $`\mathrm{\Delta }\chi ^2`$ values (median 6.27) are larger than expected (median 4.18). We have stretched the theoretical distribution by a factor $`f^2=1.5`$ to match the median of the observed $`\mathrm{\Delta }\chi ^2`$ values. In Figure 5 we show the expected distribution before (dot-dashed) and after (dashed) scaling the $`\mathrm{\Delta }\chi ^2`$ values up. Even after this scaling, a significant tail of events with high $`\mathrm{\Delta }\chi ^2`$ values remains, extending up to $`\mathrm{\Delta }\chi ^2`$=75.

To identify events that may have anomalies due to planets, we impose a cut on events with $`\mathrm{\Delta }\chi ^2>`$ 10 (last column in Tables 1, 2, 3, 4) in the non-Gaussian tail of Figure 5. This cut leaves us with 38 events to examine. We re-fit these events including a blending parameter (7 parameter fit) and notice that for 5 of these the $`\mathrm{\Delta }\chi ^2`$ is now below 10 and we therefore reject them as candidates. Upon examination of this list of events we selected 5 as the most likely planetary candidates. The remaining ones were not examined any further because the deviations seemed to be due to noise or too far away from the event peak to be considered good candidates. Plots of the event lightcurves, residuals and $`\mathrm{\Delta }\chi ^2`$ maps for the selected events can be seen in Figures 6, 7 and 8. We briefly discuss these next.

## 5 DISCUSSION OF BEST PLANET CANDIDATES

These events were re-analyzed including a blend fraction parameter to assess the significance of the residuals.

1998bul38: This event (Fig 6(a)) has a baseline $`I`$ magnitude of 14.703 and was imaged only during the rising amplification phase. This results in a wider space for the parameters to vary, thus the constraints on such events are less stringent than on events that have been imaged on both the rise and decline of amplification. There is a single point on this lightcurve that deviates by $`+9\sigma `$ at $`x/R_E2`$. OGLE has re-analyzed this event with image subtraction and the deviation is not observed anymore so we must disqualify this event as a candidate (priv. communication with A. Udalski).

1999bul37: Sparsely monitored about the peak, this event (Fig 6(b)) shows a $`>+4\sigma `$ residual near the peak with $`2\sigma `$ residuals on either side. The $`\mathrm{\Delta }\chi ^2`$ plot shows a clear dark spot. This event is disqualified as a candidate after image subtraction analysis by OGLE.

2000bul03: There are several data points deviating by $`3\sigma `$ about the peak and the $`\mathrm{\Delta }\chi ^2`$ plot shows a dark spot corresponding to a deviation caused by a perturbation of the minor image. This event (Fig 7(a)) is highly blended. It has not been imaged during the rise and there may be doubts whether this is due to microlensing.

2000bul12: This is a high amplification event (Fig 7(b)) where the best PSPL model without blending fails to fit the wings and the peak data points of the lightcurve. We refit the data using our seven parameter fit which allows for a blend fraction. This fit suggests that 2000bul12 is highly blended having a blend fraction of $`f_b=0.68`$. Figure 7(b) shows several black zones, where the 7-parameter blend model fit is significantly improved by including a planet. We must therefore retain this event as a possible planet anomaly.

Having reanalyzed OGLE data for this event by image subtraction, the MOA team recently claimed (?) that they derive a particularly high peak magnification of $`160`$ while their analysis raises interesting possibilities for a planetary presence. Our best-fit values for the maximum amplification are 50.7 and 159.0 for 6 and 7 parameters respectively. Note that the reported OGLE value for the amplification of this event using the standard Paczynski fit is 50.678. There is a strong correlation between the blend fraction and the derived maximum amplification which leads to uncertainty in $`A_0`$ for events that are significantly blended and where the blend fraction is not very well determined. HST images may be able to resolve individual sources and reduce the uncertainty.

2000bul41: The most interesting deviations from the PSPL model are seen in Fig 8. Several data points obtained by OGLE in the period HJD 2451745.683 to HJD 2451746.738 deviate from the fitted curve by several sigma and seem to suggest the presence of an anomaly in the region. There is also a single point close to the peak that deviates by $`+7.7\sigma `$. Figure 8 shows a clear black zone so we retain this event as a potential planet anomaly. However, the PLANET collaboration observed this event and disregard the point close to the peak as spurious but they have obtained no data over the 2 day period before the peak. OGLE confirms that the deviation at 2451773.59746 is fully artificial since the image is bad due to a tracking problem. Larger than usual deviations at HJD 245174 are also rejected since they seem to have been caused by bad weather. This event is therefore not retained as a candidate.

To summarise the results of this section, of the 145 events which passed our original rejection criteria, we identify 38 events with $`\mathrm{\Delta }\chi ^2>`$ 10 in the non-Gaussian tail of the residuals. Of these, only 5 events show anomalous structure in the lightcurve that could plausibly be attributed to planets.

For the 5 candidates presented, none of the anomalies is well enough sampled to be securely identified as a planet. For three of these events, the causes of the deviations have been confirmed to be not of planetary nature. This leaves us with only 2 possible candidates, 2000bul03 and 2000bul12.

## 6 PLANET DETECTION CAPABILITIES

The analysis above has identified a few events with anomalies that could be due to planets. In this section we focus our attention on the vast majority of events that are consistent with the PSPL model. The data points obtained on these events, exclude planets in certain regions near each lens. We combine the results from all the observed events to place upper limits on the number of planets per star as a function of orbit size and planet mass.

### 6.1 PLANET EXCLUSION ZONES VS $`a/R_E`$ AND $`q=m_p/m_{}`$

Following the lightcurve reduction and fitting, we calculate the net detection probability for each of the sampled events using the method outlined in ?.

We calculated the probability for two different planet/star mass ratios $`q=10^3`$ and $`q=10^4`$ for all 145 OGLE events. Briefly, the detection probability of finding a planet at any position on the lens plane for an orbital radius $`a`$ is calculated as:

$$P(\text{det}|a,q)=P(\text{det}|x,y,q)P(x,y|a)𝑑x𝑑y.$$
(7)
We take the first term on the right-hand side to be

$$P(\text{det}|x,y,q)=\{\begin{array}{cc}1\hfill & \text{if }\mathrm{\Delta }\chi ^2>\text{ 25}\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$$
(8)
Here, $`\mathrm{\Delta }\chi ^2`$ is the change in $`\chi ^2`$ for a planet at $`x,y`$ relative to the no-planet model. This term becomes significant when the planet at $`x,y`$ lies close to one of the images of the source at the time of one of the data points in the lightcurve. The second term $`P(x,y|a)`$ is obtained by randomly orienting the planet’s assumed circular orbit of radius $`a`$ and then projecting it onto the $`x,y`$ plane of the sky. This term may be written as:

$$P(x,y|a)=\{\begin{array}{cc}\frac{1}{2\pi a\sqrt{a^2r^2}}\hfill & \text{for }r=\sqrt{x^2+y^2}<a\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$$
(9)
For the majority of the events the detection limits are below 5%, while well sampled, high amplification events are dominant in imposing the most useful constraints. The total probability is given at the end of this section for detection of planets with mass ratios $`q=10^3`$ and $`q=10^4`$. Planets can be most easily detected when their orbit size matches the Einstein ring radius. For $`a<<R_E`$ the planet is too close to the lens star, well inside the main detection zones that straddle the Einstein ring. For large orbits, $`a>>R_E`$, the detection probability drops off as $`P(a/R_E)^2`$, equation (9), because such planets spend only this small fraction of their time projected on the sky in the region of the star’s Einstein ring.

Figure 9 examines the detection probability for $`a=R_E`$ as a function of peak amplification $`A_0`$ and event timescale $`t_E`$. We expect the detection probability to scale with the size and number of the detection zones arising from individual measurements. The detection zone sizes scale roughly as the area of the planet’s Einstein ring times the source amplification, while the number of detection zones scales as the event duration divided by the typical sampling time. Thus we expect the planet detection probability from the OGLE lightcurves to scale roughly as

$$P0.03\left(\frac{A_0}{3}\right)\left(\frac{t_E}{40\mathrm{d}}\right)\left(\frac{q}{10^3}\right).$$
(10)
The results in Fig. 9 are consistent with this expectation. The strong increase in detection probability with $`A_0`$ is clear, while the dependence on $`t_E`$ is less clear. If we plot residuals relative to the above model, the remaining scatter arises from details of the actual timing of observations with respect to the lightcurve peak.

Fig 10 shows the result of summing the detection probability $`P_s(\text{det}|a,q)`$ for star $`s`$ over all stars:

$$P(a,q)=\underset{s}{}P_s(\text{det}|a,q).$$
(11)
This quantifies the planet detection capability of our analysis of the OGLE lightcurves. If all stars have $`\eta `$ planets with orbit radius $`a`$ and mass ratio $`q`$, then the number of planets we expect to detect is

$$<n>=\eta P(a,q).$$
(12)
The highest detection probability occurs as expected for $`aR_E`$. The expected $`(a/R_E)^2`$ scaling is evident for large orbits and the probability is 10 times lower for $`q=10^4`$ than with $`q=10^3`$, confirming the expected linear scaling with $`q`$. For comparison we also plot the curves for $`q=10^3`$ and $`q=10^4`$ but using a different detection threshold value of $`\mathrm{\Delta }\chi ^260`$ and 100. We see that $`<n>`$ scales roughly as $`(\mathrm{\Delta }\chi ^2)^{1/2}`$, so that raising this $`\mathrm{\Delta }\chi ^2`$ threshold from 25 to 100 lowers the expected number of detections from $`7\eta `$ to $`4\eta `$.

The detection probability curve for $`q=10^3`$ peaks at $`a/R_E=1.1`$ where $`<n>7.4\eta `$. Similarly, the curve for $`q=10^4`$ also peaks at $`aR_E`$ with $`<n>1.3\eta `$. For $`\eta =1`$ and $`q=10^3`$ we expect a maximum of 7.42 detections. If the planet anomaly candidates discussed in section 5 are in fact not due to planets, $`n<1`$, then the number of planets per lens star is:

$$\eta <\frac{1}{P(a,q)}0.14\left(\frac{q}{10^3}\right)^1\left(\frac{\mathrm{\Delta }\chi ^2}{25}\right)^{1/2}$$
(13)
for planets with $`aR_E`$.

Following ? we can define the ‘lensing zone’ as $`a_1<a<a_2`$ where $`a_1=0.6R_E`$ and $`a_2=1.6R_E`$. This is the region where detection of planets by microlensing effects is reasonably efficient. If we let each lens star have $`\eta _{\mathrm{LZ}}`$ planets in the lensing zone (LZ), and distribute the orbits uniformly in $`\mathrm{log}a`$, then the expected number of planet detections is:

$$n=\eta _{\mathrm{LZ}}\frac{_{a_1}^{a_2}P(a,q)d\mathrm{log}a}{\mathrm{log}(a_2/a_1)}.$$
(14)
Evaluating this using the results in Fig 10, we find that a detection of zero planets, $`n<1`$, corresponds to

$$\eta _{\mathrm{LZ}}<0.18\left(\frac{q}{10^3}\right)^1\left(\frac{\mathrm{\Delta }\chi ^2}{25}\right)^{1/2}.$$
(15)
### 6.2 PLANET EXCLUSION ZONES VS $`a`$(AU) AND $`m(m_{JUP})`$

In Figure 10 we have presented a general description of the total detection probability in terms of the mass ratios and the projected separation of the companion from the lens in units of the Einstein ring radius. Here we convert these to physical units using a similar method to that followed by ? to give a clearer indication of planetary constraints placed by the OGLE dataset. The results are shown in Fig 11.

To convert to physical units we assume

$$R_E=1.9\text{AU}\left(\frac{t_E}{45\mathrm{d}}\right),$$
(16)
and

$$M=0.3M_{}\left(\frac{t_E}{45\mathrm{d}}\right)^{1/2},$$
(17)
where 45d is the median event timescale for the OGLE events, 1.9 AU is the median Einstein ring radius, and 0.3 $`M_{}`$ is the median lens mass. These scalings are consistent with a model in which most of the OGLE events arise from sources and lenses that are both galactic Bulge stars, with $`D_d6`$kpc, $`D_s8`$kpc.

This results in the combined detection probability for all OGLE events as a function of planetary orbital radius, as shown in Figure 11. Returning to Figure 10, we can see that both mass ratios peak at $`aR_E`$. In particular, for $`q=10^3`$, the average total detection probability for the lensing zone is $`P(\text{det}|a_1<a<a_2)=5.6`$. From Figure 11, the average detection probability for the region $`1\text{AU}<a<4\text{AU}`$ is $`P(\text{det}|1\text{AU}<a<4\text{AU})=3.2`$ for $`q=10^3`$ and the probability peaks at around 3 AU.

The curves in Fig 11 are broader and more symmetric than those in Fig 10. This arises because the detection curves all peak at $`aR_E`$, but each event has a different $`R_E`$ with the slower events corresponding on average to larger Einstein rings (eqn16). This horizontal blurring in Fig 11 makes the curve wider than in Fig 10 and also reduces the expected number of detections to

$$n4\eta \left(\frac{q}{10^3}\right)\left(\frac{\mathrm{\Delta }\chi ^2}{25}\right)^{1/2},$$
(18)
for planets with $`a3`$AU, where the detection efficiency is highest.

For zero detections, $`n<1`$, the upper limit on the number of planets per lens star is

$$\eta 0.23\left(\frac{q}{10^3}\right)^1\left(\frac{\mathrm{\Delta }\chi ^2}{25}\right)^{1/2},$$
(19)
for planets with $`a3`$ AU, and

$$\eta _{14\mathrm{A}\mathrm{U}}0.31\left(\frac{q}{10^3}\right)^1\left(\frac{\mathrm{\Delta }\chi ^2}{25}\right)^{1/2},$$
(20)
for planets uniformly distributed in $`\mathrm{log}a`$ over $`a=14`$ AU.

The solid curves in Fig 10 show the expected number detections for Jupiter-mass planets. With the star mass given by eqn 17, the mass ratio for a planet of mass $`m`$ is

$$q=3\times 10^3\left(\frac{t_E}{45\mathrm{d}}\right)^{1/2}\left(\frac{m}{m_J}\right).$$
(21)
The results for Jupiter-mass planets, $`m=m_J`$, are therefore found by a small linear extrapolation of those for $`q=10^3`$ (dotted) and $`q=10^4`$ (dashed). For zero detections, $`n<1`$, the upper limit on the number of planets per lens star is

$$\eta 0.16\left(\frac{m}{m_J}\right)^1\left(\frac{\mathrm{\Delta }\chi ^2}{25}\right)^{1/2},$$
(22)
for planets with $`a3`$ AU, and

$$\eta _{14\mathrm{A}\mathrm{U}}<0.21\left(\frac{m}{m_J}\right)^1\left(\frac{\mathrm{\Delta }\chi ^2}{25}\right)^{1/2},$$
(23)
for planets uniformly distributed in $`\mathrm{log}a`$ over $`a=14`$ AU.

## 7 UPPER LIMITS ON NUMBER OF PLANETS PER STAR

Our analysis of the OGLE events provides evidence for $`n`$ planets, where $`n2`$, since only 2 candidate planetary lens anomalies were identified, and some or both of these may be due to noise glitches rather than real planets. If we assume that all stars have $`\eta `$ planets with orbit radius $`a`$ and mass ratio $`q`$ uniformly distributed in log($`\eta `$), then the expected number of events is

$$<n>=\eta P(q,a),$$
(24)
where

$$P(a,q)=\underset{s}{}P_s(\text{det}|q,a).$$
(25)
From Bayes theorem, the probability distribution on $`\eta `$ given a measurement of $`n`$ is

$$P(\eta |n)=P(n|\eta )P(\eta )/P(n).$$
(26)
Since we have no knowledge on the distribution of $`P(\eta )`$ we shall assume $`P(\eta )\eta ^1`$. The last term in the equation,

$$P(n)=P(n|\eta )P(\eta )\text{d}\eta ,$$
(27)
is a normalization constant to ensure that $`P(\eta |n)\text{d}\eta =1`$. $`P(n|\eta )`$ is a likelihood function, i.e. the probability that we will detect $`n`$ planets given that there are $`\eta `$ per star. The number of detected planets $`n`$ is a Poisson random variable with expected value $`<n>=\eta P(q,a)=E`$. Then we have that

$$P(n|\eta )=\frac{E^n\text{exp}(E)}{n!}.$$
(28)
If the data hold no conclusive evidence of planets, i.e. $`n=0`$, then the constraint on $`\eta `$ is an exponential distribution:

$$P(\eta |n=0)=\frac{P(n=0|\eta )P(\eta )}{_0^{\mathrm{}}P(n=0|\eta )P(\eta )\text{d}\eta }.$$
(29)
The general form of $`P(\eta |n)`$ can easily be calculated to be:

$$P(\eta |n)=\frac{f^nP(q,a)^{n+1}\text{exp}(\eta P(q,a))}{n!}.$$
(30)
This is shown on Figure 12. If there are $`n`$ deviations which we assume are due to planets, then the expected number of planets per star of type $`q,a`$ is:

$$<\eta >=\eta P(\eta |n)\text{d}\eta =\frac{n}{P(q,a)}.$$
(31)
Under our assumptions and for a detection threshold value of $`\mathrm{\Delta }\chi ^225`$, we can state from the results presented in Figure 11 that the OGLE dataset indicates that less than 21 ($`n`$) %, where $`n2`$, of the lens stars have Jupiters orbiting them at an orbital radius of $`1<a<4`$ AU.

Our result uses a smaller $`\mathrm{\Delta }\chi ^2`$ threshold than the results presented by the PLANET (??) collaboration. After searching for planetary signatures in 43 intensively monitored microlensing events, they concluded that less than 33% of the $`0.3M_{}`$ stars that typically comprise the lens population have Jupiter-mass companions with semi-major axes in the range of $`1<a<4`$ AU. We must note however that their detection threshold value was very conservatively set at $`\mathrm{\Delta }\chi ^260`$. For comparison, we calculate the detection probability for the same region ($`1.5<a<4`$ AU) for $`\mathrm{\Delta }\chi ^2=60`$ to be $`28`$%. This is lower than the limit imposed by PLANET (also note that our model assumptions lead to an underestimate of the true detection probability of the order of $``$13%).

## 8 SUMMARY

We have analyzed three years (1998-2000) of OGLE observations of microlensing events to place limits on the abundance of ‘cool Jupiters’. We fitted a total of 145 events using a maximum likelihood fit that adjusts 6 parameters. We computed detection probability maps for each event using $`\mathrm{\Delta }\chi ^2`$ threshold values of 25, 60, 100. If we assume that all lenses are of the same star type and all have one Jupiter analogue at $`a/R_E1`$, then from the combined detection probability results we infer that a maximum of 7 planets with $`q=10^3`$ should have been detected. Our selection criteria returned 5 candidate events for an assumed mass ratio of $`q=10^3`$. Three of these were then confirmed to be due to other effects so only two remain as plausible candidates. Our results suggest that less than 14% of the lens stars have ‘cool Jupiters’ orbiting them at $`a/R_E1`$. We translate that in AU, recalculate the probability for Jupiter-mass planets, and from Figure 11 we can state that the OGLE dataset indicates that less than 21 ($`n`$)% of the lens star population has Jupiters orbiting them at an orbital radius of $`1<a<4`$ AU. $`n2`$ is the number of planet anomaly candidates that are actually due to planets. We conclude that observing time is more efficiently allocated by observing many events with sampling intervals that produce non-overlapping detection zones than using intensive sampling on a few events.

## 9 ACKNOWLEDGMENTS

We thank the OGLE team for placing their data on these events in the public domain and in particular A. Udalski and P. Wozniak for useful feedback. Keith Horne was supported by a PPARC Sr. fellowship and by Beatrice Tinsley visiting Professorship at the University of Texas, Austin.

## Appendix A Excluded events

We give a brief description of the events excluded from the PSPL fits and briefly discuss the possible nature of the deviation.

1998bul12: Having a baseline $`I`$ magnitude of 12.932, this event shows magnitude variations not consistent with the PSPL model. The magnitude starts increasing at JD$`2450000920`$ and reaches a maximum value of $`12.67`$mag at JD$`2450000`$ 950 . Then it drops off but not following the standard Paczynski bell-shaped curve. The lightcurve observed could be that of a variable star.

1998bul28: This event displays a double peaked structure. It could be interpreted as a binary lens or a binary source. The PLANET collaboration has obtained data on this event and they find the best binary-lens model in the range $`q=10^410^2`$ has $`\mathrm{\Delta }\chi ^219`$ (?).

1998bul32: This event is sparsely monitored about the peak which slackens our interpretation of its likely nature. It exhibits a magnitude change of $`\mathrm{\Delta }I>6`$mag and clearly deviates from the PSPL model. We can make no definite claims about the cause of this variability.

1999bul11: This is a binary lens event with clear signs of caustic crossings. The amplification increases rapidly as the source crosses the first caustic, then falls off temporarily only to increase again as the source crosses the second caustic. The magnitude then returns to its baseline value.

1999bul17: (also MACHO 99-BLG-28) The data on this event are consistent with binary microlensing.

1999bul19: A superposition of two lightcurves resulting from binary source lensing can explain the effects seen in this lightcurve.

1999bul23: Anomalies expected from binary microlensing events are obvious in this lightcurve as well. The PLANET collaboration have published a paper presenting their fits to this event. The source star is a G/K sub-giant in the Galactic Bulge. Their best-fit parameters are a mass ratio $`q0.39`$ and projected separation $`d2.42`$ (?).

1999bul25: (also MACHO 99-BLG-35) Based only on the OGLE public available data, it remains unclear whether this event is a binary lens or a binary source. The deviations around the time of peak amplification are short lived and can be accounted for by both models. However, extra data obtained by the PLANET collaboration, but unavailable to the community, strongly suggest that this event is caused by a binary lens.

1999bul28: A very short timescale event with $`t_E`$ (radius) of 5.972 days. Sparsely monitored about the peak, it is not obvious what the slight deviations from the PSPL model are due to.

1999bul29: The PLANET collaboration excludes event 1998BUL29 from their analysis since they find that a point-lens finite-source model fits their data better. We exclude it too.

1999bul32: There seem to be a number of effects present in this lightcurve. Since this is is a long timescale event with $`t_E155`$ (radius) we expect the parallax effect due to the Earths motion to be detectable. The event also seems to be severely blended.

1999bul40: A faint event with baseline magnitude $`I`$=19.780. OGLE data do not help distinguish what causes the observed deviation. PLANET have unpublished data on this event and claim this is another binary lensing event.

1999bul42: For 1999BUL42 the binary nature of the lens is seen in the lightcurve. The deviation lasts for $`21`$ days. No fits are available for the data.

2000bul28: Recently reduced PLANET data from two observing sites appear to suggest that this event is due to multiple lens microlensing. Two caustic crossings may have occurred at JD-2450000=1677 and JD-2450000=1682.

2000bul38: This is a caustic crossing binary lens event.

2000bul43: A long timescale event with $`t_E174`$ (radius) days. The event is observed only during the increasing amplification phase and the parallax effect due to the Earths orbital motion is obvious in the lightcurve.

2000bul46: A binary lens event.