*3.1. Distribution over the Two Reservoirs in Case of Lévy Noise*

If the noise in the setup of Figure 3 is Gaussian, then the system will relax to an equilibrium with equal concentration in the two reservoirs. Each particle then has a probability *P*<sup>1</sup> = *R*<sup>2</sup> 1/(*R*<sup>2</sup> <sup>1</sup> + *<sup>R</sup>*<sup>2</sup> <sup>2</sup>) to be in the larger reservoir and a probability *<sup>P</sup>*<sup>2</sup> = *<sup>R</sup>*<sup>2</sup> 2/(*R*<sup>2</sup> <sup>1</sup> + *<sup>R</sup>*<sup>2</sup> <sup>2</sup>) to be in the smaller reservoir. The probability to be in a certain reservoir is in that case, simply proportional to the volume of that reservoir. In 2D, the "volume" is the area *Vi* = *πR*<sup>2</sup> *<sup>i</sup>* /2.

Next consider Lévy particles. The distribution will now be different. As was shown in the previous section and in Appendix A, Lévy particles tend to accumulate near the walls and in the smaller "nooks and corners". With Lévy particles, the probability to be in the smaller reservoir will be larger than the reservoir's fraction of the total volume.

For a stochastic simulation, we let the semicircular walls be "sticky" again, i.e., the particle comes to a standstill upon hitting the wall and only displaces again if a subsequent computed step leads to motion inside the system. If the linear vertical wall in the middle is hit, an elastic collision occurs. Thus, that wall is "bouncy".

We will use the 2D solution for a circle, *pst*(*r*)=(*α*/(2*π*)) <sup>1</sup> − *<sup>r</sup>*<sup>2</sup> *<sup>α</sup>*/2−<sup>1</sup> , to come to an estimate of the steady-state distribution for the setup in Figure 3. We move to a description where *ρi*(*ri*), with *i* = 1, 2, denotes the normalized particle density in reservoir *i* at a distance *ri* from the opening. With:

$$\rho\_i(r\_i) = \frac{a}{\pi R\_i^2} \left( 1 - \left(\frac{r\_i}{R\_i}\right)^2 \right)^{a/2 - 1} \tag{5}$$

it is readily verified that:

$$\int\_{r\_i=0}^{R\_i} \int\_{\phi=-\pi/2}^{\pi/2} \rho\_i(r\_i) r\_i \, dr\_i \, d\phi = 1. \tag{6}$$

With a large number of particles in the setup, there will be a relaxation to a distribution with a fraction *ϕ*<sup>1</sup> in reservoir 1 and a fraction *ϕ*<sup>2</sup> in reservoir 2. Obviously we have *ϕ*<sup>1</sup> + *ϕ*<sup>2</sup> = 1. For any distribution over the two reservoirs we have:

$$
\rho(r\_1, r\_2) = \varphi\_1 \rho\_1(r\_1) + \varphi\_2 \rho\_2(r\_2). \tag{7}
$$

It is easy to see that *ρ*(*r*1,*r*2) = 1, where the integration is over the entire 2-semicircle system in the figure.

The steady state occurs if there are as many 1 → 2 transitions as there are 2 → 1 transitions. We will next derive what values of *ϕ*<sup>1</sup> and *ϕ*<sup>2</sup> lead to steady state. In the above figure, imagine a semicircular strip of width *dri* at a distance *ri* from the opening. The number of particles in the strip is *ρi*(*ri*)*πridri* (*i* = 1, 2). We assume that for *r* > *r*0, we are in the region where the power-law-description of the tail of the Lévy distribution (*pα*(*r*) <sup>∝</sup> *<sup>r</sup>*−(*α*+1) as *<sup>r</sup>* <sup>→</sup> <sup>∞</sup>) applies. The probability that the displacement during one timestep is *larger than r* is then proportional to *r*−*α*. For small *d* and sufficiently large *r*0, the angle *θ*, cf. Figure 3, will be small and we have *d* = *θr*. For a Lévy jump to lead to a particle transiting to the other reservoir, the jump must also be in the right direction. This brings in a factor (*d*/*r*) cos *φ*, where *φ* is the indicated angle of the position on the semicircle with the horizontal. Integrating over *φ* from −*π*/2 to *π*/2, the full direction factor is found to be 2*d*/*r*. All in all, during one timestep we have for the number of cross-reservoir transitions from a distance between *r* and *dr*:

$$dn\_i^{tr}(r\_i, r\_i + dr\_i) \propto q p\_i \frac{d}{r\_i} r\_i^{-a} \rho\_i(r\_i) r\_i \, dr\_i. \tag{8}$$

Integrating from *r*<sup>0</sup> to the boundary *Ri*, we obtain for the number of Lévy-jump-associated transitions from reservoir *i*:

$$N\_i^{tr} \propto \frac{q\_i}{R\_i^2} \int\_{r\_i = r\_0}^{R\_i} r\_i^{-a} \left(1 - \left(\frac{r\_i}{R\_i}\right)^2\right)^{a/2 - 1} dr\_i. \tag{9}$$

The proportionality constant (associated with the ∝) and the *r*<sup>0</sup> (the radius from which the power law is taken to describe the Lévy-stable distribution) are the same for both reservoirs. At this point, it is also important to realize that for the Lévy jumps to dominate the number of 1 → 2 and 2 → 1 transitions, *R*<sup>1</sup> and *R*<sup>2</sup> must both be much larger than *r*0.

*Mathematica* will readily give an analytical solution for the integral Equation (9). The solution involves the hypergeometric function [32]. Before working out Equation (9) in its full generality, we make a further simplification that will not affect the solution too much: As *Ri r*<sup>0</sup> for both *i* = 1 and *i* = 2, we take *r*<sup>0</sup> = 0 to be the lower limit of the integral. With 0 < *α* < 1, the integral will not diverge with *ri* → 0. Next, the all-important reservoir radius *Ri* can be scaled out of the actual integral and incorporated in the prefactor:

$$N\_i^{tr} \propto \frac{q\_i}{R\_i^2} R\_i^{-a} R\_i \int\_{r\_i = 0}^{R\_i} \left(\frac{r\_i}{R\_i}\right)^{-a} \left(1 - \left(\frac{r\_i}{R\_i}\right)^2\right)^{a/2 - 1} d\left(\frac{r\_i}{R\_i}\right). \tag{10}$$

Upon taking *u* = *ri*/*Ri* and *v* = *u*<sup>2</sup> (so *dv* = *du*<sup>2</sup> = 2*u du* and thus *du* = 1/(2 <sup>√</sup>*v*) *dv*), further simplification is achieved:

$$N\_i^{tr} \propto \varphi\_i R\_i^{-1-a} \int\_{v=0}^1 v^{-a/2-1/2} (1-v)^{a/2-1} \, dv. \tag{11}$$

The integral on the right-hand side is the well-known Euler integral, which is also known as the beta function [32]. Ultimately, this integral depends only on *α*. It is finite for 0 < *α* < 1 and as it is the same for both reservoirs, we find:

$$N\_i^{tr} \propto \varphi\_i R\_i^{-1-\alpha}.\tag{12}$$

The steady state condition is *<sup>ϕ</sup>*1*R*−1−*<sup>α</sup>* <sup>1</sup> <sup>≈</sup> *<sup>ϕ</sup>*2*R*−1−*<sup>α</sup>* <sup>2</sup> . With *<sup>ϕ</sup>*<sup>1</sup> <sup>+</sup> *<sup>ϕ</sup>*<sup>2</sup> <sup>=</sup> 1 we then get:

$$\varphi\_1 \approx \frac{R\_1^{1+a}}{R\_1^{1+a} + R\_2^{1+a'}}, \quad \varphi\_2 \approx \frac{R\_2^{1+a}}{R\_1^{1+a} + R\_2^{1+a'}} \quad \text{and} \quad \frac{\varphi\_1}{\varphi\_2} \approx \left(\frac{R\_1}{R\_2}\right)^{1+a}.\tag{13}$$

The better approximation is obtained by not fully carrying through the *r*<sup>0</sup> = 0 simplification of the last paragraph. That the simple approximation according to Equation (13) fails for larger values of *α* is partly due to scaling issues. For the analytic approximation to be consistent with the numerics, we need Δ*t* 1/*<sup>α</sup>* to be significantly smaller than *r*<sup>0</sup> (cf. Equation (3) with *σ* = 1). Setting *r*<sup>0</sup> = 0 leads to a range where this is no longer true. As *α* becomes larger, this range becomes larger. Keeping *r*<sup>0</sup> > 0 in Equation (9), we find after some algebra and use of *Mathematica* for the equivalent expression of Equation (12):

$$\boldsymbol{N}\_{i}^{tr} \propto \boldsymbol{\varrho}\_{i} \left[ \boldsymbol{R}\_{i}^{-1-a} - \frac{\sqrt{\pi}}{\Gamma(\frac{a}{2})\Gamma(\frac{3-a}{2})} \frac{\boldsymbol{r}\_{0}^{1-a}}{\boldsymbol{R}\_{i}^{2}} \,\_{2}F\_{1} \left( \frac{1-a}{2}, 1 - \frac{a}{2}; \frac{3-a}{2}; \left( \frac{\boldsymbol{r}\_{0}}{\boldsymbol{R}\_{i}} \right)^{2} \right) \right],\tag{14}$$

where <sup>2</sup>*F*1(*a*, *b*; *c*; *z*) is the aforementioned hypergeometric function. It is readily verified that the second term in the square brackets dominates for *α* → 2 and small *r*0. This is due to the *<sup>r</sup>*1−*<sup>α</sup>* <sup>0</sup> term. The hypergeometric function is defined as a power series [32] and under the *<sup>r</sup>*<sup>0</sup>  1 condition we can still take (*r*0/*Ri*)<sup>2</sup> ∼ 0 and hence <sup>2</sup>*F*1(.) ≈ 1. The ratio of particles in the two reservoirs is then:

$$\frac{\varrho p\_1}{\varrho\_2} \approx \left(\frac{R\_1}{R\_2}\right)^{1+a} \frac{\sqrt{\pi}}{\sqrt{\pi}} \frac{R\_2^{a-1} - r\_0^{a-1} \Gamma(\frac{3-a}{2}) \Gamma(\frac{a}{2})}{\sqrt{\pi} \, R\_1^{a-1} - r\_0^{a-1} \Gamma(\frac{3-a}{2}) \Gamma(\frac{a}{2})},\tag{15}$$

which reduces to (*R*1/*R*2)1+*<sup>α</sup>* (cf. Equation (13)) if we take *<sup>α</sup>* <sup>&</sup>lt; 1 and *<sup>r</sup>*<sup>0</sup> <sup>→</sup> 0 concurrently. Note, furthermore, that the equilibrium distribution, i.e., *ϕ*1/*ϕ*<sup>2</sup> = (*R*1/*R*2)2, is properly approached if we concurrently take *α* → 2 and *r*<sup>0</sup> → 0. Both the approximations according to Equations (13) and (15) are depicted in Figure 4 and compared with the results of a stochastic simulation. Finally, it is worth pointing out that Equation (15) is still an approximation. The power law, Equation (1), that characterizes the Lévy-stable distribution is not valid for small values of *ξ*. For values of *ξ* near zero, the distribution is Gaussian-like and this is what is relevant for the behavior of particles close to the opening, i.e., *r* → 0. Gaussian diffusion near the opening will lead to a continuous and differentiable steady-state concentration profile near the opening. This is also what Figure 5 shows.

**Figure 4.** For the setup of Figure 3 with *R*<sup>1</sup> = 10 and *R*<sup>2</sup> = 1, we let *ϕ*<sup>1</sup> and *ϕ*<sup>2</sup> represent the fraction of particles in reservoir 1 and 2, respectively, at steady state. The curves depict the analytic approximations, Equation (13) (dashed) and Equation (15) (solid), of *ϕ*1/*ϕ*2. Each dot is the result of a stochastic simulation of 40,000 particles for 4 <sup>×</sup> <sup>10</sup><sup>5</sup> timesteps (with <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> 0.001) following a 2 <sup>×</sup> 105 timestep relaxation period. For the approximation according to Equation (15), we let *r*<sup>0</sup> = 0.05 and find good agreement with the result of the stochastic simulation.

Figure 4 shows the ratio *ϕ*1/*ϕ*<sup>2</sup> as a function of *α* and compares the result of a stochastic simulation with Equations (13) and (15). We took *R*<sup>1</sup> = 10 and *R*<sup>2</sup> = 100. For *α* → 0, the simple approximation according to Equation (13) leads to *ϕ*1/*ϕ*<sup>2</sup> = 10. For the more sophisticated approximation according to Equation (15), the *ϕ*1/*ϕ*<sup>2</sup> value at *α* → 0 can be brought arbitrarily close to 10 by taking *R*<sup>1</sup> and *R*<sup>2</sup> much larger than *r*0. There is 10 times as much "sticky wall" in the large reservoir and this result tells us that for *α* → 0, all particles are concentrated at the sticky walls as would intuitively be expected.

The result that is derived in Appendix B hints at the reason that *α* = 1 is "almost like" *α* = 2. As we move away from the opening, the probability to hit the opening *decreases* as *r*−*α*. However, with a homogeneous distribution of particles, the number of particles at a distance between *r* and *r* + *dr increases* proportional to *r*. For an *n*-dimensional setup, the increase is proportional to *rn*−<sup>1</sup> (for *n* = 2 we have circular strips and for *n* = 3 we have spherical shells). All in all, we find that the number of "hits" from a distance *r* is proportional to *rn*−*α*−1. Note that for *n* = 3, the entire range of *α* leads to an increase of "hits" with *r*. We have not done any further investigation of the 3D case. We see that for *n* = 2, an increase of "hits" with *r* only occurs if *α* < 1. For 1 < *α* < 2, the number of "hits" decreases with *r* and in that case transitions mostly happen from the region around the

opening. This decrease with 1 < *α* < 2 also means that the particle exchange through the opening does not "feel" the different radii of the different reservoirs anymore.

**Figure 5.** The figure on the left depicts a steady-state distribution for 50,000 Lévy particles in a two-reservoir confinement as depicted in Figure 3 after 10<sup>5</sup> timesteps. We have *R*<sup>1</sup> = 2, *R*<sup>2</sup> = 1, and the opening has a width *d* = 0.1. For the figure on the right we started with a steady-state distribution and ran the simulation for another 10<sup>5</sup> iterations. We took a horizontal strip through the center with a width of 0.02 and partitioned it into 300 bins. Particles in each bin were counted and the results of the subsequent 10<sup>5</sup> iterations were added. The solid line represents the resulting normalized 1D histogram. The dashed reference curve is the solution Equation (4). For the left reservoir, the domain was scaled to a length 2. Normalization of the combination of analytic solutions was done such that the probability to be in the left reservoir is 2/3. It is readily verified that this leads to continuity at the location of the opening.

Equation (4) describes and Figure 5 shows a nonhomogeneous distribution: As we move away from the opening, the concentration actually increases. This should add to the exponent *n* − *α* − 1 that we derived in a previous paragraph. Some of this effect is incorporated in the approximation that led to Equation (15). Both that approximation and the simulations show an asymptotic approach to (*R*1/*R*2)<sup>2</sup> as *<sup>α</sup>* → 2 and *<sup>r</sup>*<sup>0</sup> → 0.
