2. MARKOV CHAIN MODEL FOR THE NUMBER OF ATOMS IN A MAGNETO-OPTICAL TRAP

Suppose the door to the MOT could be opened the instant it became empty or closed at the instant it became occupied. Then the number of atoms in the trap can be modeled as a Markov chain. Because of physical limitations, the number of atoms in the trap cannot be monitored continuously: The trap can only be checked every T seconds. During that interval T, more than one atom can sneak inside the trap before the door is closed.

Let R_n be a Poisson random variable with parameter λ, representing the load (i.e., the random number of atoms entering the trap at time n) and denote by X_n the number of atoms in the MOT at time n,n = 1,2,…. For a fixed X_n = x, the random loss variable Y_n has a Poisson distribution with parameter μx. A simple recursive model for the number X_n of atoms present in the MOT at step n under a feedback regime is

with R_n and (X_n−1,Y_n−1) independent. As a simple example, suppose that at step n, there are three atoms in the MOT, so X_n = 3; according to the model, the trap is opened. If we lose two atoms, then at step n + 1, there is only one atom left in the trap (X_n+1 = 1), but if we lose all three atoms, X_n+1 is determined by a new load; notice that in the latter case, we only consider one step, although we have two time intervals. Our goal is to get a single atom with high probability; if we know that the trap empties, there is no point in stopping the process.

Clearly, X_n,n = 1,2,…, forms a countable state Markov chain in discrete time. If X_n → X, with the distribution of X being the steady-state (stationary) distribution for the chain, then the following equality of distributions holds:

Here, R is Poisson distributed with parameter λ, and for a fixed X = x, the random variable Y is Poisson distributed with parameter μx.

Our goal is to find the distribution of X and to find the values of λ and μ that maximize the value π₁ = P(X = 1). For this purpose, we determine the probability vector π = (π₀,π₁,…) with π_k = P(X = k), k = 0,1,2,….

Notice that the transition probabilities of the above Markov chain are

In particular,

with d₀ = 1.

Observe that

which is maximized when μ = 0 (the value of λ being immaterial.)

3. STEADY-STATE DISTRIBUTION FOR THE NUMBER OF ATOMS

Theorem 1: For any λ,μ > 0, the Markov chain determined by transition probabilities given in Eq. (2.3) is ergodic. In fact, the chain is geometrically ergodic; that is, for some R(x) and ρ < 1,

n = 1,2,…. The vector π of stationary probabilities has the form given in Eq. (3.13) with the vector u satisfying Eq. (3.8).

Proof: We will prove that the Markov chain determined by Eq. (2.3) has a stationary distribution π by constructing it. For the stationary distribution to exist, we need

so π₀ ≠ 0, and for k = 1,2,…, with v_k = π_k /π₀,

Also,

. Introduce the lower triangular matrix B with the entries

For μ > 0, we have

, and for 0 < μ ≤ 1, the elements of B are in fact the probabilities of the Borel–Tanner distribution, so that, with e^T = (1, 1,…), e^TB = e^T.

According to Eq. (2.3), the reduced transition probabilities matrix P formed by the elements {P(x|y),x,y = 1,2,…} can be written in the form

Here, D is the diagonal matrix with elements 1,2,…,d^T = (d₁,d₂,…), and [ell ]^T = (λ,λ²,…,λ^j/[(j − 1)!],…). Since P(0|y) = e^−λλ^y/y!, one has for v^T = (v₁,v₂,…),

On the other hand,

Combining these formulas, one obtains for u = Dv,

We will prove that the norm of B, as an operator in L₂, is smaller than e^−μ/2 (<1). Hence, I − B is invertible in L₂. To this end, we show that for any i = 1,2,…,

From definition (3.4),

so to prove Eq. (3.9), it suffices to demonstrate that

However, this inequality follows directly from comparison of the coefficients i^k − ki^k−1 ≤ (i − 1)^k, k = 1,…,i − 1 . Since

, a theorem by Schur (Halmos [4, Problem 37], with p_i ≡ 1, implies that the operator defined by the matrix B has norm in L₂ less than e^−μ/2 < 1. Note that b₁₁ = e^−μ, so that ∥B∥₂ ≥ e^−μ.

Since [ell ] belongs to the space L₂, the solution u of Eq. (3.8) is also in L₂ and can be found from the formula

Now we can prove the existence and uniqueness of the stationary distribution π for μ > 0 and λ > 0. Start with u given by Eq. (3.12) and put v = D⁻¹u. By the Cauchy–Schwartz inequality, v ∈ L₁. The stationary distribution π is given by

The uniqueness of the stationary distribution π follows from the uniqueness of solution for Eq. (3.8). By Theorem 6.9 in Kemeny, Snell, and Knapp [7, p.135], the Markov chain with the transition probabilities given by Eq. (2.3) is ergodic.

Our argument shows that the operator ζI − B is invertible if ζ > e^−μ/2. It follows that for such values of ζ, one can find in L₂ the solution u of the equation

for any right-hand side [ell ] in L₂. This fact implies that one can solve the “drift” inequalities

with positive V and C = {0}. It follows from Nummelin [9, p.90] that our chain is geometrically ergodic, which concludes the proof of the theorem. █

The referee suggested a proof of geometric ergodicity by checking Eq. (3.15) with V(x) = β^x for some β > 1:

Notice that as

, and for some 0 < ζ < 1, e^{−μ(1−1/β)} < ζ. It follows that Eq. (3.15) holds for some finite x_ζ and positive β where C = {k : k ≤ x_ζ}.

To get a useful approximation to π₁, we use the form of the matrix B and denote elements of

by α_ij = α_ij(μ), and put

.

Theorem 2: The probability π₁ is given by

with an upper bound (3.24). An approximation (3.26) is valid, and a lower bound (3.31) for π₁ holds, leading to approximation (3.36) for 0 < μ < 1. Small values of μ and λ are optimal.

Proof: Using the identity (I − B)⁻¹ = I + B(I − B)⁻¹, the α_ij can be found from the recurrent formula

for i ≥ j. Note that

One has

so that using Eq. (3.13), we obtain Eq. (3.17). Note that

It is easy to see by induction that

. If for k = 1,…,j − 1, α_k1 ≤ α₁₁ /k, then for j ≥ 2,

Also, as iα_ij ≤ jα_jj,

From Eq. (3.22), we obtain

which shows that for this probability to be close to unity, one must have μ ≈ 0 or λ must be large. For

it follows from Eqs. (3.22) and (3.23) that lim_n→∞ π₁⁽ⁿ⁾ = π₁ and

where

is the incomplete gamma function.

An argument similar to that in Eq. (3.23) can be used to demonstrate that for any fixed j, iα_ij is a nonincreasing sequence in i = j,j + 1,…, so that the limit, δ_j = lim_i→∞ iα_ij, exists. Since ∥(I − B^T)⁻¹∥₂ ≤ [1 − ∥B∥₂]⁻¹ and, for a fixed l, b_lp is a decreasing sequence in p ≤ l, we get, for i ≥ j + 2 from Eq. (3.18),

so that

For 0 < μ < 1, lim_i→∞ i²b_ij = 0, and, therefore, δ_j = δ_j+1 for any j, demonstrating that δ_j ≡ δ. It follows that

Using Eq. (3.23), a lower bound for π₁ is obtained,

To get a formula for δ when 0 < μ < 1, recall that for any positive λ,

As neither α_kj nor d_k depend on λ, it follows that for any

, which implies that

Indeed,

since, as in the derivation of the expected value of the Borel–Tanner distribution (e.g., formula (16) in Haight and Breuer [3]),

These results lead to an approximation for π₁; if one replaces u₁ by λα₁₁ + λ²α₂₁ + δ(e^λ − 1 − λ − λ²/2) and sets u₂ ≈ λ²α₂₂ + δ(e^λ − 1 − λ − λ²/2), u_j ≈ δ(e^λ − 1 − λ − λ²/2),j ≥ 3. With

denoting the integral exponential function and C denoting Euler's constant, the approximate formula for π₁ when 0 < μ < 1 takes the form

For small values of μ,

This ratio is a decreasing function of λ; for small λ,π₁ ∼ λ/(λ + μ), so that π₁ → 1 as λ ↓ 0,μ/λ → 0. Therefore small values of λ and μ = o(λ) are optimal. █

Theorem 2 gives a method for approximating the probability of interest π₁. Given ε, from Eq. (3.26) we determine the level of truncation n. From Eq. (3.25), using the recurrence formulas (3.18) and (3.19), we obtain the approximate value π₁ⁿ of π₁ such that |π₁ − π₁ⁿ| ≤ ε. Figure 1 depicts the probability π₁ for λ,μ ∈ (0,2.5) using ε = 10⁻⁶.

The plot of the probability π1 as a function of λ and μ.

The referee has suggested using the geometric distribution as the loss distribution. Thus, assume that Y_n = min(X_n,G), where G is a random variable with geometric distribution over nonnegative integers with parameter p. Then the inverse of I − B can be obtained explicitly. For any load distribution R with the finite first moment, the same technique as in Theorem 1 shows that the stationary distribution is obtained from Eq. (3.13):

In particular,

which is an increasing function in p for any load distribution such that 1 − P(R = 0) ≥ P(R = 1)E(R). The latter condition holds for the Poisson load distribution. Therefore, values of p close to 1 and λ ≈ 0 are optimal for this distribution.