text:
# Phenomenological framework for nonequilibrium steady states of molecular motors

## Abstract

An expression for the energetic efficiency of a molecular motor is presented in terms of an effective temperature, which was defined based on the ratio of the correlation function to the susceptibility of its velocity. We also present a numerical result regarding this temperature for a thermal ratchet model. Furthermore, our expression of efficiency using known values for kinesin can adequately explain the experimental results. The framework presented here has a closed form phenomenologically, and is independent of any detailed models.

What is the essential difference between our muscle and the Carnot cycle as an energy transducer? The recent development of single-molecule-detection techniques Block  for molecular motors led to new perspectives for this question, by directly visualizing their noisy way of chemo-mechanical energy transduction at the single-molecule level. Many theoretical studies have been made based on such experiments to present various models Reimann , upon which energetics have also been discussed Parrondo . However, since most of them are bottom-up type theories, they are in a dilemma that they cannot be compared with experiments until a model and parameters for a molecular motor have been determined completely. This is the very difference from a top-down type theory, namely a phenomenology such as thermodynamics, which deals with energetics based on macroscopic observables alone. It thus sounds meaningful to construct a phenomenology for a molecular motor.

In this Letter, we present an expression for the energy efficiency of a molecular motor based on the susceptibility and the correlation function of its velocity, both of which are observables. For this purpose, we introduce an effective temperature as the ratio of the correlation function to the susceptibility. According to Fluctuation Dissipation Theorem (FDT) Kubo , this effective temperature coincides the temperature of the heat bath when the motor is in equilibrium, while we show that it varies in general with respect to the timescale of measurement, if the motor is at a nonequilibrium steady state.

First, we define the effective temperature. The velocity $`v(t)`$ of a colloidal particle in solution responds linearly to a small external force $`f(t)`$ as

$$v(t)v_0=_{\mathrm{}}^t\chi (tt^{})f(t^{})dt^{},$$
(1)
where $`\mathrm{}`$ is the statistical average, $`v_0`$ is the mean velocity of a reference system (a steady system without perturbation), and $`\chi (t)`$ is called a response function. The Fourier-Laplace transformation of $`\chi (t)`$, given as

$$\chi [\omega ]=_0^{\mathrm{}}\chi (t)e^{i\omega t}dt,$$
(2)
is called the susceptibility. The 1st type of FDT states that if the reference system is in equilibrium, the real part of the susceptibility $`\chi ^{}[\omega ]\mathrm{Re}\chi [\omega ]`$ is proportional to the correlation function of the fluctuation of the reference system.

$$\mathrm{\Phi }[\omega ]=T_0\chi ^{}[\omega ],$$
(3)
where

$$\mathrm{\Phi }[\omega ]=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}(v(0)v_0)(v(t)v_0)e^{i\omega t}dt.$$
(4)
$`T_0`$ is the temperature of the reference equilibrium system. We set the Boltzmann constant to unity.

We proceed to the case where the reference system is not in equilibrium but rather in a nonequilibrium steady state. In this case, since FDT may not hold in a normal sense, the coefficient in Eq. (3) differs from usual temperature in general even when the system is in contact with a single heat bath. We can adopt Eq. (3) as a definition of temperature for the nonequilibrium steady state. After calculating the susceptibility (Eq. (2)) and the correlation function (Eq. (4)) for a nonequilibrium steady state as a reference system, the effective temperature is obtained as

$$T[\omega ]\frac{\mathrm{\Phi }[\omega ]}{\chi ^{}[\omega ]}.$$
(5)
We call this a ‘dissipation temperature’. This formulation for a degree of freedom asymmetric with respect to a time reversal, such as velocity, is a counterpart of that introduced in spin glass systems Peliti  for a symmetrical degree of freedom such as magnetization. It may be obvious that the dissipation temperature varies depending on the timescale $`\omega `$ of observation in general.

Before an energetics argument, we compute the dissipation temperature for a typical model of thermal ratchet, known as a flashing ratchet Astumian ; Prost , in a nonequilibrium steady state, to gain insight into this newly defined temperature. A flashing ratchet is one of the well known models for molecular motors Okada , in which Brownian motion of a particle can be rectified in one direction at a nonequilibrium steady state. In this paper, we suppose that a particle has two chemical states, each of which provides a periodic potential profile $`U_1(x)`$ and $`U_2(x)`$, respectively (see the inset in Fig. 1 (a)). The transition between these two states is assumed to be stochastic at the same and constant transition rate $`\alpha `$. The Langevin equation and the transition dynamics for the particle are

$$\begin{array}{c}m\ddot{x}+\gamma \dot{x}=_xU_i(x)+f(t)+\xi (t),\\ \xi (t)=0,\xi (t)\xi (t^{})=2\gamma T_0\delta (tt^{}),\\ \left[1\right]\genfrac{}{}{0pt}{}{{\scriptscriptstyle \genfrac{}{}{0pt}{}{\alpha }{}}}{\genfrac{}{}{0pt}{}{}{\alpha }}\left[2\right].\end{array}$$
(6)
We assume that the mass $`m`$ and viscosity $`\gamma `$ of the particle are the same for each state. $`\xi (t)`$ is white Gaussian thermal noise with a magnitude given by the temperature of the heat bath $`T_0`$.

We numerically solved the equations (6) with some initial conditions. After some time lag, the system reaches a steady state. First, the double-sided power spectrum density is calculated by Fourier transformation for the velocity $`v(t)=\dot{x}(t)`$ of the particle under no perturbation, i.e. $`f=0`$. By the Wiener-Khinchine theorem Kubo , this amount is equal to $`\mathrm{\Phi }[\omega ]`$. Next, the velocity is determined in the presence of a small oscillating force $`f(t)=ϵ\mathrm{cos}(\omega t)`$. The average linear power loss $`f(t)(v(t)v_0)`$ results in the in-phase susceptibility $`\chi ^{}[\omega ]=af(t)(v(t)v_0)/ϵ^2`$, where $`a=1`$ for $`\omega =0`$ and $`a=1/2`$ for $`\omega 0`$ (see Fig. 1 (a)). Finally, we can obtain the dissipation temperature by Eq. (5). We confirmed that the dissipation temperature calculated in this procedure equals $`T_0`$ in the case of $`\alpha =0`$, i.e. in equilibrium. Figure 1 (b) shows the numerically computed dissipation temperature $`T[\omega ]`$ in a nonequilibrium case. As shown, the dissipation temperature has different values with a long time limit ($`\omega 0`$) and a short time limit ($`\omega \mathrm{}`$), and undergoes a sigmoid-like transition at an intermediate frequency $`\omega ^{}`$. The dissipation temperature equals the temperature of the heat bath in the short time limit and increases with a long time limit Harada . The transition timescale $`2\pi /\omega ^{}`$ coincides the mean passage time $`\tau `$ of the particle in this ratchet system, i.e. the mean time to pass along a single period of the potential.

We now address the energetics in terms of this dissipation temperature. We assume that the motion of a ratchet system or a molecular motor can be approximated by the following generalized Langevin equation:

$$_{\mathrm{}}^t\mathrm{\Gamma }(tt^{})v(t^{})dt^{}=f+\xi (t)+\mathrm{\Xi }(t).$$
(7)
$`f`$ is a small external load. $`\xi (t)`$ is a thermal noise which satisfies normal FDT as $`\xi (t)=0`$ and $`\xi (t)\xi (t^{})=T_0\mathrm{\Gamma }(tt^{})`$. Suppose that $`\mathrm{\Xi }(t)`$ is the nonequilibrium portion of noise with the mean value $`\mathrm{\Xi }(t)=\mathrm{\Xi }_0`$ and the correlation $`(\mathrm{\Xi }(0)\mathrm{\Xi }_0)(\mathrm{\Xi }(t)\mathrm{\Xi }_0)=M(tt^{})`$. We assume no correlation between $`\xi (t)`$ and $`\mathrm{\Xi }(t)`$.

By solving Eq. (7) with Fourier transformation and comparing the results with the definition of $`\chi [\omega ]`$, we obtain

$$\chi [\omega ]^1=\mathrm{\Gamma }[\omega ]_0^{\mathrm{}}\mathrm{\Gamma }(t)e^{i\omega t}dt,$$
(8)
where $`\mathrm{\Gamma }[\omega ]`$ is a Fourier-Laplace transformation of the friction function $`\mathrm{\Gamma }(t)`$. Similarly, the following relation is obtained,

$$\mathrm{\Phi }[\omega ]=(T_0\mathrm{\Gamma }^{}[\omega ]+M[\omega ])/|\mathrm{\Gamma }[\omega ]|^2,$$
(9)
where $`M[\omega ]1/2_{\mathrm{}}^{\mathrm{}}M(t)e^{i\omega t}dt`$, and $`\mathrm{\Gamma }^{}[\omega ]\mathrm{Re}\mathrm{\Gamma }[\omega ]`$. From Eq. (5) (note that $`\chi ^{}[\omega ]=\mathrm{\Gamma }^{}[\omega ]/|\mathrm{\Gamma }[\omega ]|^2`$),

$$T[\omega ]=T_0+\frac{M[\omega ]}{\mathrm{\Gamma }^{}[\omega ]}.$$
(10)
This relation corresponds to the 2nd kind of FDT. We define $`\mathrm{\Delta }T[\omega ]T[\omega ]T_0=M[\omega ]/\mathrm{\Gamma }^{}[\omega ]`$.

We now perform long-time averaging on Eq. (7) after multiplying by $`v(t)`$ and obtain the expression

$`fv(t)={\displaystyle _{\mathrm{}}^t}\mathrm{\Gamma }(tt^{})v(t)v(t^{})dt^{}`$
$`+{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{\Gamma }^{}[\omega ]}{\mathrm{\Gamma }[\omega ]}}T_0{\displaystyle \frac{\mathrm{d}\omega }{2\pi }}+\mathrm{\Xi }_0v(t)+{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{\Gamma }^{}[\omega ]}{\mathrm{\Gamma }[\omega ]}}\mathrm{\Delta }T[\omega ]{\displaystyle \frac{\mathrm{d}\omega }{2\pi }}.`$
This expression represents a balance of energy flux. The left-hand side represents the work against the external load. On the right-hand side, the 1st term corresponds to dissipation into the heat bath, the 2nd term is the energy input from the bath, and the 3rd and 4th terms reflect external energy input. Therefore, the efficiency of energy transduction from external energy to work can be defined as

$$e\frac{fv}{Q_{\mathrm{ir}}+\mathrm{\Xi }_0v},$$
(12)
where $`Q_{\mathrm{ir}}_{\mathrm{}}^{\mathrm{}}\left(\mathrm{\Gamma }^{}[\omega ]/\mathrm{\Gamma }[\omega ]\right)\mathrm{\Delta }T[\omega ]d\omega /2\pi `$. It is clear from this expression that a dissipation temperature greater than the heat bath temperature leads to irreversible dissipation and decreases the efficiency. On the other hand, in the ideal case where $`\mathrm{\Delta }T[\omega ]=0`$ and $`f\mathrm{\Xi }_0`$, which corresponds to a quasistatic process, we get $`e1`$. Thus, Eq. (12) represents how much work can be extracted from thermodynamical free energy Derenyi . This definition of efficiency must be identified from what usually appears in the theory for the Carnot cycle in thermodynamics, $`\eta `$, which reflects how much work can be extracted from an absolute amount of input energy. In the present context, $`\eta `$ is expressed as $`\eta fv/\left(Q+\mathrm{\Xi }_0v\right)`$, where $`Q_{\mathrm{}}^{\mathrm{}}\left(\mathrm{\Gamma }^{}[\omega ]/\mathrm{\Gamma }[\omega ]\right)T[\omega ]d\omega /2\pi `$.

It is also clear that the condition $`Q_{\mathrm{ir}}0`$ is necessary because the efficiency $`e`$ must be less than unity. Thus, the fact that the dissipation temperature is equal to or greater than the heat bath temperature is associated with the stability of the nonequilibrium steady state.

Let us compare this result with the experimental results for the conventional kinesin. Although neither the correlation function nor susceptibility have been measured for kinesin, information is available regarding the diffusion coefficient, mobility and so on. Thus, we adopt the following approximation based on the assumption that the dissipation temperature of kinesin has the same characteristics as that of a flashing ratchet (Fig. 2 (b)). First, we suppose that the dissipation temperature of kinesin is also a sigmoidal function with a transition around the mean passage time $`\tau `$. We also suppose that the imaginary part of $`\mathrm{\Gamma }[\omega ]`$, which is related to the inertia of the motor, is negligible above the timescale of inertia $`m/\gamma `$, which is much smaller than $`\tau `$. Under these assumptions, $`Q_{\mathrm{ir}}`$ can be approximated as $`Q_{\mathrm{ir}}2\mathrm{\Delta }T[0]\omega ^{}/2\pi =2(D/\mu T_0)/\tau `$, where $`D`$ and $`\mu `$ are the diffusion coefficient and mobility, respectively. Thus, efficiency can be expressed as

$$e\frac{\mu f(\mathrm{\Xi }_0f)}{2(D/\mu T_0)/\tau +\mu \mathrm{\Xi }_0(\mathrm{\Xi }_0f)}.$$
(13)
We used the values shown in Table 1 from Refs. Svoboda ; Nishiyama .

The value of $`\mathrm{\Xi }_0`$ was determined from a stalling force for kinesin, since the force-velocity curve of kinesin is almost linear. Based on these values, the long-time dissipation temperature of kinesin is $`T[0]=D/k_\mathrm{B}\mu 900`$ K, while $`T_0300`$ K (We only get back Boltzmann constant $`k_\mathrm{B}`$ in this paragraph). Equation (13) can be evaluated as $`ef(8f)/\left(16.8+8(8f)\right)`$. Figure 2 shows this function along with the experimental data for kinesin taken from Ref. Nishiyama . These data agree very well without any fitting parameters.

A remarkable feature of the present theory is its closed formalism only with observables. We assumed no specific model to calculate the efficiency. However we need to discuss the assumption used here: the generalized phenomenological Langevin equation (7). First, we must carefully consider whether the motion of a molecular motor is a Markovian process that can be described with such a linear Langevin equation. This assumption requires that the internal force of the motor $`\mathrm{\Xi }(t)`$ must be small enough. For the same reason, the present theory is limited on the case of a small external load $`f`$ nonlinearity . This assumption seems to be satisfied for kinesin, since its response to a external force is linear. Moreover, although we assumed an equation only for velocity, it is necessary that there is little or no dissipation for the other degrees of freedom, so that the above assumption does not affect the energetics. If this point is violated, an argument for another degree of freedom is conducted in parallel with the argument presented here, by measuring the dissipation temperature for such a degree of freedom and constructing a suitable generalized Langevin equation. In this case, the dissipation temperature will be a tensor (we do not present the details here). However, the fact that the present theory which considers only velocity can reproduce the experimental data of kinesin might mean that dissipation for the other degrees of freedom is negligible compared to that for velocity.

In conclusion, we have presented an expression for the energetic efficiency of a molecular motor using observable quantities alone. A concept of the dissipation temperature is also introduced based on the susceptibility and the correlation function of the motor. The dissipation temperature for a flashing ratchet exhibits two regimes with respect to the observation timescale. In the short-time regime, the dissipation temperature equals the temperature of the heat bath, while it increases in the long-time regime. The two regimes meet around the mean passage time of the ratchet system. The expression of efficiency derived in the present paper showed a good agreement with the experimental data for kinesin.

The dissipation temperature is an observable quantity in principle using recently developed single-molecule-detection techniques for molecular motors, although there has been no previous observation. The measurement of the dissipation temperature of processive motors, such as kinesin, or cooperative systems, such as muscle fibers or flagella, should be an interesting problem in itself, to clarify their energetic properties over a wide range of timescale. The data will be also useful for testing previously described models of molecular motors.

The issue of effective temperature is also related to fundamental problems of thermodynamics in nonequilibrium systems. It has been recently reported that the effective temperature based on FDT actually has the nature of temperature in several nonequilibrium systems Barrat ; Sasa . We need a further investigation to judge whether the dissipation temperature discussed here can also be a natural temperature for molecular motors, which hopefully leads us to an understandng of nonequilibrium thermodynamics for molecular motors.

The author acknowledge useful discussions with Prof. K. Yoshikawa, Dr. S. Sasa, and Ms. Hayashi. This work is supported in part by Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists (No. 05494).


format: text

text:
# 1 Introduction

## 1 Introduction

An important question in the theory of quasirandomness, uniform distribution of points, and diophantine approximation is the following: For which $`n`$ is it true that there exists an integer $`k`$, $`1k<n`$ and $`(k,n)=1`$, so that $`k/n`$ has a continued fraction whose partial quotients are bounded in average by a constant $`B`$? That is, if we write $`k/n=[0;a_1,a_2,\mathrm{},a_m]`$, we wish to find $`k`$ so that

$$t^1\underset{i=1}{\overset{t}{}}a_iB$$

for all $`t`$ with $`1tm`$. Denote by $`(B)`$ the set of all $`n`$ for which such a $`k`$ exists. These sets are discussed at length in  and the related matter of partial quotients bounded uniformly by a constant appears as an integral part of . This latter question is closely connected with Zaremba’s Conjecture (), which states that such a $`k`$ exists for all $`n>1`$ if we take $`B=5`$.

Define the continuant $`K(a_1,a_2,\mathrm{},a_m)`$ to be the denominator of the continued fraction $`k/n=[0;a_1,a_2,\mathrm{},a_m]`$. In , it is proven that, if $`S_n(B)`$ is the number of sequences $`𝐚=(a_1,\mathrm{},a_m)`$ bounded uniformly by $`B`$ with $`K(𝐚)n`$ and $`H(B)`$ is the Hausdorff dimension of the set of continued fractions with partial quotients bounded uniformly by $`B`$, then

$$\underset{n\mathrm{}}{lim}\frac{\mathrm{log}(S_n(B))}{\mathrm{log}n}=2H(B).$$

Then, in , $`H(2)`$ is calculated with a great deal of accuracy: $`H(2)0.53128`$. Therefore, $`S_n(2)`$, and thus the number of $`p/q`$ with $`qn`$ whose partial quotients are bounded by $`2`$, is $`n^{1.0625\mathrm{}+o(1)}`$. (This improves the previous best known lower bound, $`n^{1.017}`$ computed in , slightly.)

Define $`\overline{S}_n(B)`$ to be the number of sequences $`𝐚=(a_1,\mathrm{},a_m)`$ with partial quotients bounded in average by $`B`$ so that $`K(𝐚)n`$. Clearly, $`\overline{S}_n(B)S_n(B)`$, so $`\overline{S}_n(2)n^{1.0625}`$. In the next section, we prove something much stronger, however – an exponent of $`1.5728394`$ – thus providing a lower bound in the first nontrivial case. Section 3 discusses the implications for the density of $`(2)`$ and a few open problems.

## 2 The Proof

###### Theorem 1.

For any $`ϵ>0`$, $`\overline{S}_n(2)n^{2\mathrm{log}2/\mathrm{log}(1+\sqrt{2})ϵ}`$.

###### Proof.

The proof consists of two parts: computing the number of positive sequences of length $`m`$ bounded in average by $`2`$, and then computing the smallest possible $`m`$ so that $`K(a_1,\mathrm{},a_m)>n`$ and the $`a_i`$ are bounded in average by $`2`$.

First, we wish to know how many sequences $`(a_1,\mathrm{},a_m)`$ there are with $`a_j1`$ for each $`j[m]`$ and $`_{j=1}^ra_j2r`$ for each $`r[m]`$. Call this number $`T(m)`$. By writing $`b_j=a_j1`$, we could equivalently ask for sequences $`(b_1,\mathrm{},b_m)`$ with $`b_j0`$ for each $`j[m]`$ and $`_{j=1}^rb_jr`$ for each $`r[m]`$. This is precisely the number of lattice paths from $`(0,0)`$ to $`(m,m)`$ which do not cross the line $`y=x`$, and so $`T(m)`$ is the $`m^{\text{th}}`$ Catalan number, or $`(m+1)^1\left(\genfrac{}{}{0pt}{}{2m}{m}\right)=4^{m(1o(1))}`$.

In the following lemmas, we show that $`K(a_1,\mathrm{},a_m)n`$ if $`m\mathrm{log}n(1o(1))/\mathrm{log}(1+\sqrt{2})`$. Therefore, setting $`m`$ as large as possible, we have at least

$$4^{\mathrm{log}n(1o(1))/\mathrm{log}(1+\sqrt{2})}=n^{2\mathrm{log}2/\mathrm{log}(1+\sqrt{2})o(1)}$$

sequences with partial quotients bounded in average by $`2`$ and continuant $`n`$ . ∎

We must show that the size of a continuant with partial quotients bounded in average by $`B`$ is at most the largest size of a continuant with partial quotients bounded by $`B`$.

###### Lemma 2.

If the sequence $`(a_1,\mathrm{},a_m)`$ of positive integers is bounded in average by $`B>1`$, then $`K(a_1,\mathrm{},a_m)K(\underset{m}{\underset{}{B,\mathrm{},B}})`$.

###### Proof.

We prove the Lemma by a “shifting” argument. That is, we perform induction on the size of the entry $`a_j`$ such that $`a_j>B`$ and $`j`$ is as small as possible. If $`𝐚=(a_1,\mathrm{},a_m)`$ contains no $`a_t>B`$, we are done, because increasing the partial quotients can only increase the continuant. If there is some $`a_t>B`$, let $`t2`$ be the smallest such index. We consider two cases: (i) $`a_tB+2`$ or $`a_{t1}<B`$, and (ii) $`a_t=B+1`$, $`a_k=B`$ for $`skt1`$ for some $`2st1`$, and $`a_{s1}<B`$. (Clearly, $`𝐚(B,B,\mathrm{},B,B+1,a_{t+1},\mathrm{},a_m)`$, since this sequence is not bounded in average by $`B`$. Therefore we may assume $`s2`$.)

Case (i):

Let $`𝐛=(b_1,\mathrm{},b_m)`$ = $`(a_1,\mathrm{},a_{t1}+1,a_t1,\mathrm{},a_m)`$. We show that $`K(𝐛)>K(𝐚)`$. First, note that

$$\underset{j=1}{\overset{r}{}}b_j=\{\begin{array}{cc}_{j=1}^ra_j\hfill & \text{ if }rt1\hfill \\ 1+_{j=1}^{t1}a_j\hfill & \text{ if }r=t1.\hfill \end{array}$$

Since $`a_tB+1`$, $`_{j=1}^{t1}a_jtBB1`$, so $`1+_{j=1}^{t1}a_j(t1)B`$, and $`𝐛`$ is bounded in average by $`B`$. Second, note that it suffices to consider the case of $`t=m`$, since, if $`K(b_1,\mathrm{},b_j)>K(a_1,\mathrm{},a_j)`$ for $`1jt`$, then $`K(𝐛)>K(𝐚)`$. (That is, $`K()`$ is monotone increasing.)

Let $`q_j=K(a_1,\mathrm{},a_j)`$ and $`q_j^{}=K(b_1,\mathrm{},b_j)`$. (We use the convention that $`q_j=0`$ when $`j<0`$ and $`q_0=1`$.) Clearly, $`q_j=q_j^{}`$ if $`j<t1`$. When $`j=t1`$, we have $`q_{t1}^{}>q_{t1}`$ by monotonicity. When $`j=t`$,

$$q_t=a_tq_{t1}+q_{t2}=a_t(a_{t1}q_{t2}+q_{t3})+q_{t2}=(a_ta_{t1}+1)q_{t2}+a_tq_{t3},$$

and

$`q_t^{}`$ $`=(b_tb_{t1}+1)q_{t2}^{}+b_tq_{t3}^{}`$
$`=((a_t1)(a_{t1}+1)+1)q_{t2}+(a_t1)q_{t3}`$
$`=q_t+q_{t2}(a_ta_{t1}1)q_{t3}.`$
Since $`a_ta_{t1}+2`$ and $`q_{t2}>q_{t3}`$, we have

$$q_t^{}q_t+q_{t2}q_{t3}>q_t.$$

Case (ii).

Now, assume that $`a_t=B+1`$, $`a_k=B`$ for $`skt1`$ for some $`2st1`$, and $`a_{s1}<B`$. Then define $`𝐛=(b_1,\mathrm{},b_m)`$ by letting $`b_j=a_j`$ if $`js1`$ and $`jt`$; $`b_{s1}=a_{s1}+1`$; and $`b_t=a_t1`$. Again, we may assume that $`t=m`$. Then

$$\underset{j=1}{\overset{r}{}}b_j=\{\begin{array}{cc}_{j=1}^ra_j\hfill & \text{ if }r=t\text{ or }r<s1\hfill \\ 1+_{j=1}^ra_j\hfill & \text{ if }s1rt1.\hfill \end{array}$$

For any $`r`$ such that $`s1rt1`$,

$$\underset{j=1}{\overset{r}{}}a_j=\underset{j=1}{\overset{t}{}}a_j\underset{j=r+1}{\overset{t}{}}a_jBt(B(tr1)+(B+1))Br1.$$

Therefore, $`_{j=1}^rb_jBr`$ for all $`r[t]`$, and we may conclude that $`𝐛`$ is bounded in average by $`B`$.

Define $`F_k`$ as follows: $`F_0=0`$, $`F_1=1`$, and, for $`k>1`$, $`F_k=BF_{k1}+F_{k2}`$. Then it is easy to see by induction that

$$K(\underset{k}{\underset{}{B,\mathrm{},B}},x)=F_{k+1}x+F_k.$$

Also,

$$K(y,c_1,\mathrm{},c_r)=yK(c_1,\mathrm{},c_r)+K(c_2,\mathrm{},c_r).$$
(1)
Taking $`k=ts`$, we deduce

$$K(a_{s1},\mathrm{},a_t)=a_{s1}((B+1)F_{k+1}+F_k)+(B+1)F_k+F_{k1},$$

and

$`K(b_{s1},\mathrm{},b_t)`$ $`=(a_{s1}+1)(BF_{k+1}+F_k)+BF_k+F_{k1}`$
$`=K(a_{s1},\mathrm{},a_t)+(Ba_{s1})F_{k+1}`$
$`K(a_{s1},\mathrm{},a_t)+F_{k+1}.`$
If $`s=2`$, we are done. Otherwise, we use that

$`K(b_{s2},\mathrm{},b_t)`$ $`=a_{s2}K(b_{s1},\mathrm{},b_t)+K(b_s,\mathrm{},b_t)`$
$`a_{s2}K(a_{s1},\mathrm{},a_t)+F_{k+1}+b_tK(b_s,\mathrm{},b_{t1})+K(b_s,\mathrm{},b_{t2})`$
$`=a_{s2}K(a_{s1},\mathrm{},a_t)+F_{k+1}+K(a_s,\mathrm{},a_t)K(a_s,\mathrm{},a_{t1})`$
$`=K(a_{s2},\mathrm{},a_t).`$
Now, inductive application of (1) to the continuants $`K(b_{sj},\mathrm{},b_t)`$, $`3js1`$, yields $`K(𝐛)K(𝐚)`$, since $`a_{sj}=b_{sj}`$ in this range.

By repeating cases (i) and (ii) as appropriate, we will eventually reach a sequence of partial quotients bounded by $`B`$, and at each stage we never decrease the corresponding continuant. The result therefore follows. ∎

It remains to find a bound on $`K(B,\mathrm{},B)`$.

###### Lemma 3.

If $`B1`$, $`K(\underset{m}{\underset{}{B,\mathrm{},B}})\left(\frac{1}{2}(B+\sqrt{B^2+4})\right)^{m+1}`$.

###### Proof.

We proceed by induction. The case $`m=0`$ is trivial. Suppose it is true for all $`m<M`$. Then, by (1),

$`K(\underset{M}{\underset{}{B,\mathrm{},B}})`$ $`=BK(\underset{M1}{\underset{}{B,\mathrm{},B}})+K(\underset{M2}{\underset{}{B,\mathrm{},B}})`$
$`B\left({\displaystyle \frac{1}{2}}(B+\sqrt{B^2+4})\right)^M+\left({\displaystyle \frac{1}{2}}(B+\sqrt{B^2+4})\right)^{M1}`$
$`\left({\displaystyle \frac{1}{2}}(B+\sqrt{B^2+4})\right)^{M1}\left({\displaystyle \frac{1}{2}}B^2+{\displaystyle \frac{1}{2}}B\sqrt{B^2+4}+1\right)`$
$`=\left({\displaystyle \frac{1}{2}}(B+\sqrt{B^2+4})\right)^{M+1}.`$
## 3 The Density of $`(2)`$

###### Corollary 4.

There is a constant $`C`$ and a subset $`S`$ of the positive integers such that $`\mathrm{log}|S[n]|/\mathrm{log}n\mathrm{log}2/\mathrm{log}(1+\sqrt{2})o(1)0.786`$ so that, for each $`nS`$, there exists a $`k[n]`$, $`(k,n)=1`$ so that $`k/n`$ has partial quotients bounded in average by $`2`$.

###### Proof.

Let $`U`$ be the set of all reduced fractions $`p/q`$, $`1<p<q`$, whose partial quotients $`𝐚=(a_1,a_2,\mathrm{},a_m)`$ are bounded in average by $`2`$ and such that $`𝐚^{}=(a_2,\mathrm{},a_m)`$ is bounded in average by $`2`$. The number of such $`𝐚`$ with $`K(𝐚)n`$ is at least twice the number of sequences $`𝐚^{}=(a_2,\mathrm{},a_m)`$ bounded in average by $`2`$ with $`K(𝐚^{})n/3`$, because, if $`[𝐚^{}]=p/q`$, then $`K(𝐚)=a_1q+p3K(𝐚^{})n`$. (The fact that $`𝐚^{}`$ is bounded in average by $`2`$ implies that $`[1,𝐚^{}]`$ and $`[2,𝐚^{}]`$ are also.) Then, since every rational has at most two representations as a continued fraction, the number of elements of $`U`$ whose denominator is $`n`$ is at least $`\overline{S}_{n/3}(2)`$, which is at least $`n^{2\mathrm{log}2/\mathrm{log}(1+\sqrt{2})o(1)}`$. Let $`S`$ be the set of denominators of fractions appearing in $`U`$. If $`p/q=[𝐚]`$ is in $`U`$, then $`[𝐚^{}]=(qa_1p)/p`$, so $`p`$ is the continuant of a sequence whose partial quotients are bounded in average by $`2`$. Therefore, $`\overline{S}_{n/3}(2)|S[n]|^2`$, and we may conclude that $`\mathrm{log}|S[n]|/\mathrm{log}n\mathrm{log}2/\mathrm{log}(1+\sqrt{2})o(1)`$. ∎

Attempts by the author to find a generalization of the above result to $`(B)`$ by applying much more careful counting arguments when $`B>2`$ have failed thus far. It would also be interesting to (i), calculate the Hausdorff dimension of the set of reals in $`[0,1)`$ whose partial quotients are bounded in average by $`B`$, and (ii), draw a connection, similar to that of the “uniform” case, between this quantity and the asymptotic density of $`(B)`$.

AMS Classification Numbers: 11K50; 11K38


text:
# Migration of extrasolar planets to large orbital radii

## 1 Introduction

With one exception (Konacki et al. 2003), all confirmed extrasolar planets have been discovered by the Doppler velocity technique. The selection effects inherent to radial velocity surveys (Cumming, Marcy & Butler 1999) favor the detection of planets at small orbital radii. To date, about half of the known planets have semi-major axis $`a<1`$ AU, while the most distant - 55 Cnc d - lies at $`5.9`$ AU from its parent star. 111From the online Extrasolar Planets Encyclopedia, at http://cfa-www.harvard.edu/planets/catalog.html (Schneider 2003) as of March 26th, 2003.. Indirect evidence, however, suggests that there could be a sizable population of massive planets at much greater radii. Recent observations of dusty debris around Vega have been interpreted as suggesting the presence of a planet of a few Jupiter masses with $`a>30`$ (Wilner et al. 2002). Further, simulations modeling circumstellar dust discs suggest a planet lies at a distance of $`5565`$ AU from Epsilon Eridani (Ozernoy et al. 2000).

Forming planets in these outer locations is difficult. Gas giants must form before the disc dissipates, at timescales no greater than $`510`$ Myr (Haisch, Lada, & Lada 2001). In standard core accretion models (Safranov 1969), the timescale for building the core of a giant planet increases rapidly with radius, with a $`t_{form}`$ scaling approximately as $`a^2`$ (Pollack et al. 1996). Although such models are undoubtedly oversimplified (Pollack et al. 1996; Bryden, Lin & Ida 2000), it is hard to avoid the conclusion that forming massive planets at radii of several tens of AU within $`10`$ Myr is difficult. Indeed, this has led to the suggestion that Uranus and Neptune may have formed at smaller radii in our own Solar System (Thommes, Duncan & Levison 1999, 2002). Motivated by these issues, we investigate the possibility of forming massive planets at small $`a`$, followed by outward migration. In Sections 2 and 3, we consider sequentially the two mechanisms that have been extensively studied in the context of inward migration: planet-disc interactions (Goldreich & Tremaine 1980; Lin & Papaloizou 1986; Lin, Bodenheimer & Richardson 1996; Trilling et al. 1998) and gravitational scattering after disc dissipation (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997; Ford, Havlickova & Rasio 2001; Terquem & Papaloizou 2002). Our conclusions are briefly summarized in Section 4.

## 2 Migration via gas disc interactions

### 2.1 Methods

We calculate the orbital evolution of massive planets embedded within an evolving protoplanetary disc using a variant of the approach described by Armitage et al. (2002). We use a simple, one-dimensional (i.e. vertically averaged) treatment to model the evolution of a protoplanetary disc evolving under the action of both internal viscous torques and external torques from one or more embedded planets (Goldreich & Tremaine 1980; Lin & Papaloizou 1986; Trilling et al. 1998; Trilling, Lunine & Benz 2002). For a disc with surface density $`\mathrm{\Sigma }(R,t)`$, the governing equation is,

$$\frac{\mathrm{\Sigma }}{t}=\frac{1}{R}\frac{}{R}\left[3R^{1/2}\frac{}{R}\left(\nu \mathrm{\Sigma }R^{1/2}\right)\frac{2\mathrm{\Lambda }\mathrm{\Sigma }R^{3/2}}{(GM_{})^{1/2}}\right]+\dot{\mathrm{\Sigma }}_w.$$

Here, $`\nu `$ is the kinematic viscosity which models angular momentum transport within the disc gas, and $`\dot{\mathrm{\Sigma }}_w`$ is a term which allows for mass to be lost from the disc – for example as a consequence of photoevaporation. The second term within the brackets describes how the disc responds to the planetary torque, $`\mathrm{\Lambda }(R,a)`$, where this function is the rate of angular momentum transfer per unit mass from the planet to the disc. For a planet in a circular orbit at radius $`a`$, we take,

$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{q^2GM_{}}{2R}}\left({\displaystyle \frac{R}{\mathrm{\Delta }_p}}\right)^4R<a`$
$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{q^2GM_{}}{2R}}\left({\displaystyle \frac{a}{\mathrm{\Delta }_p}}\right)^4R>a`$ (1)
where $`q=M_p/M_{}`$, the mass ratio between the planet and the star,

$$\mathrm{\Delta }_p=\mathrm{max}(H,|Ra|),$$
(2)
and $`H`$ is the scale height of the disc. Guided by detailed protoplanetary disc models (Bell et al. 1997), we adopt $`H=0.05R`$.

The transfer of angular momentum leads to orbital migration of the planet at a rate,

$$\frac{\mathrm{d}a}{\mathrm{d}t}=\left(\frac{a}{GM_{}}\right)^{1/2}\left(\frac{4\pi }{M_p}\right)_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}R\mathrm{\Lambda }\mathrm{\Sigma }dR,$$
(3)
if the only torque on the planet comes from the gravitational interaction with the disc.

In the core accretion model for giant planet formation (e.g. Pollack et al. 1996), the accretion of the gaseous envelope is predicted to take longer than any other phase of the formation process (several Myr in the baseline model for Jupiter presented by Pollack et al. 1996). In particular, the time scale for accretion is much longer than the time scale on which a sufficiently massive planet can establish a gap in the protoplanetary disc, since numerical simulations show that an approximately steady-state gap can be set up by a massive planet within $`10^2`$ orbital periods (e.g. Lubow, Seibert & Artymowicz 1999). A consequence of the inequality of these time scales is that massive planets – those of several Jupiter masses – probably accrete most of their envelopes subsequent to the development of a gap in the protoplanetary disc. Numerical simulations show how this accretion may occur. Gas from the outer disc penetrates the leaky tidal barrier created by the planet, and flows inward to form a small circumplanetary disc around the growing planet (Lubow, Seibert & Artymowicz 1999; D’Angelo, Henning & Kley 2002; Bate et al. 2003).

The existence of mass flow across gaps onto planets is intrinsically a two (or three) dimensional phenomenon (e.g. the discussion in Artymowicz & Lubow 1996). The torque function (Eq. 1) used in our one dimensional code establishes a clean gap for all planet masses above about $`0.1M_J`$, and this gap precludes any mass flow across the gap, or onto the planet. To allow for the mass growth of planets, we have therefore modified the one dimensional treatment to explicitly include mass flow from the outer disc on to the planet. We begin by making an approximate fit to the results of two-dimensional numerical simulations (Lubow et al. 1999; D’Angelo et al. 2002). We define the efficiency of mass accretion across the gap via a parameter $`ϵ`$, which is the planetary accretion rate as a fraction of the disc accretion rate at large radii (away from the location of the planet). The results of the aforementioned numerical simulations can then be approximated by the formula,

$$\frac{ϵ}{ϵ_{\mathrm{max}}}1.668\left(\frac{M_p}{M_J}\right)^{1/3}e^{\frac{M_p}{1.5M_J}}+0.04,$$
(4)
where $`M_J`$ is the mass of Jupiter and $`ϵ_{\mathrm{max}}`$ is an adjustable parameter which can be used to test how the results depend upon the overall efficiency of planetary accretion. We use the above equation to calculate at each timestep the appropriate planetary accretion rate. We then remove the required amount of mass from the first zone on the outer edge of the gap, and add it to the mass of the planet. Note that we assume that all the mass flow onto the planet originates from the outer disc, and do not permit any material to ‘bypass’ the planet and flow directly from the outer disc to the inner disc.

Mass accretion across the gap onto a planet may also be expected to lead to accretion of angular momentum, given that the specific angular momentum of gas at the outer gap edge exceeds that of the planet. This effect – which it is easy to show can have a significant influence on the migration rate – cannot be straightforwardly measured from existing numerical simulations222Bate et al. (2003), for example, explicitly exclude this advective torque from their estimates of the migration rate, due to difficulties in measuring the net torque from the small scale circumplanetary disc formed in their simulations.. For this paper, we adopt the simplest approach, and assume that the accreted gas has the same specific angular momentum as gas at a radius $`R_{gap}=1.6a`$. This fixed radius approximates the location of the outer edge of the gap throughout most of the calculation.

Eq. (2.1) is solved on a fixed, non-uniform mesh using standard explicit numerical methods (e.g. Pringle, Verbunt & Wade 1986). The mesh is uniform in a scaled variable $`X\sqrt{R}`$. Typically, 300 grid points are used, with an inner boundary at 0.1 AU and an outer boundary at 200 AU. A zero-torque ($`\mathrm{\Sigma }=0`$) boundary condition is applied at $`R_{in}`$. For the protoplanetary disc model adopted, the outer boundary is at sufficiently large radius that the choice of boundary condition there has no influence on the results.

### 2.2 Protoplanetary disc model

We model the protoplanetary disc as a viscous accretion disc (Lynden-Bell & Pringle 1974) suffering mass loss at large radius as a consequence of photoevaporation (e.g. Shu, Johnstone & Hollenbach 1993). Motivation for considering models of this form is provided first and foremost by observations of photoevaporative flows in Orion (Johnstone, Hollenbach & Bally 1998), and is discussed further by Clarke, Gendrin & Sotomayor (2001), Matsuyama, Johnstone & Hartmann (2003), and Armitage, Clarke & Palla (2003). Briefly, we write the viscosity as a fixed (in time) power-law in radius. For most of the calculations we have adopted the form,

$$\nu =3\times 10^{13}\left(\frac{R}{1\mathrm{AU}}\right)\mathrm{cm}^2\mathrm{s}^1.$$
(5)
This yields a steady-state surface density profile $`\mathrm{\Sigma }R^1`$, which is similar to that derived from more detailed protoplanetary disc models over the radii of interest (Bell et al. 1997). To test how sensitive the conclusions are to this assumption, we have also run one model with a different form for the viscosity,

$$\nu =1.3\times 10^{13}\left(\frac{R}{1\mathrm{AU}}\right)^{3/2}\mathrm{cm}^2\mathrm{s}^1.$$
(6)
Mass loss from the disc scales with radius as,

$`\dot{\mathrm{\Sigma }}_w`$ $`=`$ $`0,R<R_g`$
$`\dot{\mathrm{\Sigma }}_w`$ $``$ $`R^1,R>R_g,`$ (7)
with $`R_g=5\mathrm{AU}`$. We express the normalization of the mass loss via a parameter $`\dot{M}_{\mathrm{wind}}`$, which is defined as the total mass loss from the disc for a disc with an outer edge at 25 AU. The instantaneous rate of mass loss will therefore differ from this value depending upon the extent of the disc.

The initial surface density profile for the $`\nu R`$ disc model is,

$$\mathrm{\Sigma }=\mathrm{\Sigma }_0\left(1\sqrt{\frac{R_{in}}{R}}\right)\frac{1}{R}e^{R/R_0},$$
(8)
while the $`\nu R^{3/2}`$ model is identical except for the replacement of $`1/R`$ by $`1/R^{3/2}`$. Here, $`\mathrm{\Sigma }_0`$ is a constant used to define the initial accretion rate, while $`R_0`$ is a truncation radius which sets a smooth exponential cut off to the surface density at large radius. For $`RR_0`$, a disc described by this initial condition has a constant accretion rate, so we specify the initial surface density of our models via this inner accretion rate $`\dot{M}_{\mathrm{init}}`$.

### 2.3 Results

Fig. 1 shows how the evolution of the accretion rate and disc mass varies with the strength of photoevaporative mass loss. We have computed models with our standard viscosity ($`\nu R`$) that have $`\dot{M}_{\mathrm{wind}}=10^9M_{}\mathrm{yr}^1`$, $`\dot{M}_{\mathrm{wind}}=2.5\times 10^9M_{}\mathrm{yr}^1`$, and $`\dot{M}_{\mathrm{wind}}=5\times 10^9M_{}\mathrm{yr}^1`$. For consistency with observational determinations of protoplanetary disc parameters in nearby star-forming regions (e.g. Gullbring et al. 1998), we adopt for all of these models an initial accretion rate of $`\dot{M}_{\mathrm{init}}=5\times 10^8M_{}\mathrm{yr}^1`$, and a truncation radius of $`R_0=10\mathrm{AU}`$. This yields an initial disc mass of $`0.066M_{}`$. An additional model (shown as the dashed curve in Fig. 1) was calculated with the $`\nu R^{3/2}`$ viscosity law and a mass loss of $`\dot{M}_{\mathrm{wind}}=5\times 10^9M_{}\mathrm{yr}^1`$. With identical choices of $`\dot{M}_{\mathrm{init}}`$ and $`R_0`$, the initial disc mass for this model was $`0.067M_{}`$.

As expected from previous calculations (Clarke et al. 2001; Matsuyama et al. 2003), all four models show qualitatively similar evolution. There is an initial phase in which the disc mass and accretion rate decline slowly, due primarily to mass accretion onto the star. Subsequently, the mass and accretion rate drop more rapidly as the evolution becomes dominated by mass loss via the wind (Clarke, Gendrin & Sotomayor 2001). Higher rates of mass loss reduce the disc lifetime, but all models have observationally acceptable lifetimes in the range between 4 Myr and 7 Myr. Most importantly for our purposes, in all four models mass loss is at least reasonably important (relative to accretion) for the overall evolution of the disc. The fraction of the initial disc mass that is lost in the wind varies between 40 percent and 52 percent for the $`\nu R`$ models, and is 39 percent for the $`\nu R^{3/2}`$ model. We note that the fraction lost in the wind is only a weak function of $`\dot{M}_{\mathrm{wind}}`$, because the interval of time over which the wind acts is reduced for higher instantaneous mass loss rates.

To investigate how planets migrate within the evolving disc, we run the protoplanetary disc models repeatedly. In each run, we add an initially low mass planet ($`0.1M_J`$) into the disc at a specified time and radius. We then allow the planet to grow and migrate within the evolving disc, and record the final planet mass and orbital radius once the disc has been either accreted or lost in the wind. By varying the formation time $`t_{\mathrm{form}}`$, and the formation radius $`a_{\mathrm{form}}`$, we study how the final outcome depends upon when and where in the disc massive planets form.

Fig. 2 shows the results for the disc with the standard viscosity law and varying rates of mass loss. We considered planet formation radii of 5 AU and 10 AU, safely outside any estimate of the snow line (Sasselov & Lecar 2000), and took $`ϵ_{\mathrm{max}}=1`$. The final planet mass and orbital radius are plotted as a function of the formation time. For the low and intermediate rates of mass loss, the sense of orbital migration is predominantly inward. Planets formed near the end of the disc lifetime end up in orbits close to where they formed, accrete relatively little disc gas, and remain as low mass objects. Planets formed earlier migrate inwards under the action of gravitational torques, and have time to grow to several Jupiter masses. These results are consistent with previous studies of migration (Trilling et al. 1998; Armitage et al. 2002; Trilling, Lunine & Benz 2002).

For the highest rate of mass loss from the disc, however, qualitatively different evolutionary tracks, shown in the right-hand panels of Fig. 2, are obtained. The enhanced mass loss means that there is a larger range of disc radii across which the radial velocity of the gas is outward. This can drive significant outward migration. For our choices of parameters, we find that inward migration persists for all planets formed at 5 AU, while outward migration is the rule for planets formed at 10 AU. For formation times between about 1 Myr and 2 Myr, migration approximately stalls (similar to the behaviour reported by Matsuyama, Johnstone & Murray 2003, though for slightly different reasons), while for earlier formation times the gas can drive these outer planets to radii of around 20 AU or greater. Significant accretion onto the planet occurs throughout this time, so the planets stranded at larger radii are all predicted to be massive objects.

Any attempt to distil the inherently multi-dimensional physics of planetary migration into a fast one dimensional scheme is bound to be approximate, and there are particularly obvious uncertainties in our models for planet growth and disc evolution. We have already demonstrated, as our main result, that for planets forming at radii of around 10 AU a switch between inward and outward migration occurs when $`\dot{M}_{\mathrm{wind}}`$ is varied by a factor of a few. The mass loss rate via photoevaporation is clearly a vital control parameter. To check how important some of the other parameters are, we have recalculated the migration of planets formed at 10 AU in two different models. In one, we altered the assumed disc viscosity (to $`\nu R^{3/2}`$ rather than $`\nu R`$), with the new viscosity chosen to produce a steady-state surface density profile $`\mathrm{\Sigma }R^{3/2}`$. This scaling is one used often in studies of the Solar nebula (Weidenschilling 1977). A second model was computed with the standard viscosity law, but with an accretion efficiency parameter $`ϵ_{\mathrm{max}}=0.5`$. This change halves the rate of growth of planets via accretion. Both models used the high rate of mass loss previously found to be conducive to outward migration.

Fig. 3 shows the final radii attained by planets in the three models. Outward migration occured in all three models, reflecting the primary importance of the assumed mass loss rate in determining the fates of the model planets. Changing the efficiency of accretion on to the planet made negligible difference to the final planetary radii (though it reduced the final masses of the model planets by approximately a factor of two). Substantially greater migration, however, was obtained in the calculation with the different viscosity law, despite the fact that a smaller fraction of the initial disc mass was actually lost to the wind in this case. We interpret this as being a side effect of the different surface density profile. Outward migration occurs when the torque from the inner disc exceeds that from the disc at radii beyond the planet’s orbit. The steeper surface density profile of the $`\mathrm{\Sigma }R^{3/2}`$ means that there is less mass initially exterior to the planet. As this mass is lost in the wind, the now unbalanced torque from the inner disc is more effective in driving the planet outward.

### 2.4 Observational implications

What are the implications of our gas disc migration calculations for the origin of planets at large orbital radii? We believe that three general conclusions can be drawn. First, mass loss from the outer disc can drive substantial outward migration, even when the mass loss is modest enough that the disc can survive for several Myr. Higher rates of mass loss would lead to more dramatic migration, but the resulting short disc lifetimes might preclude planet formation. In the central regions of the Orion Nebula, for example, Johnstone et al. (1998) infer mass loss rates between $`2\times 10^8M_{}\mathrm{yr}^1`$ and $`6\times 10^7M_{}\mathrm{yr}^1`$, with correspondingly small estimated disc lifetimes. Our results suggest that photoevaporation can have a major impact on planetary migration even in substantially more benign environments. Second, outward migration driven by the gas disc is favoured in systems where photoevaporative mass loss is stronger. The predicted rate of mass loss due to photoevaporation for a Solar mass star is rather small (Shu, Johnstone & Hollenbach 1993), so for relatively isolated stars we would expect outward migration only for stars significantly more massive than the Sun. An alternative possibility is that the mass loss is driven by external irradiation, as in the case of Orion (Johnstone, Hollenbach & Bally 1998). Finally, outward migration via this mechanism is a relatively slow process, which occurs on the viscous time of the protoplanetary disc. There is ample time for initially low mass planets to accrete substantial gaseous envelopes, so we would expect planets at large radius to be massive objects.

## 3 Migration in multiple planet systems

### 3.1 Introduction

If the initial outcome of the planet formation process is a system of several massive planets, subsequent gravitational interactions can lead to many possible outcomes. Planets may collide, be ejected from the system, or settle into quasi-periodic or periodic orbits. Theoretical investigations of the orbital evolution of two-planet systems (Ford, Havlickova, & Rasio 2001) and few-planet systems (Chambers, Wetherill, & Boss 1996) have been conducted with detailed descriptions of special cases, such as equal mass planets on coplanar orbits.

Gladman (1993) established that systems with two close planets exhibit chaotic but quasi-periodic behavior given the appropriate initial conditions. Ford, Havlickova, and Rasio (2001) explored the evolution of such systems when the mass of each of two planets revolving around a solar-type star in nearly circular, coplanar orbits is equal to $`10^3M_{}`$. We will expand on this study by considering different mass ratios for the two planets revolving around a central star. In doing so, we will show that outward migration is possible, for either planet.

### 3.2 Methods

#### 3.2.1 Motivation

The setup for our set of simulations is motivated primarily by the gravitational scattering experiments performed by Ford, Havlickova, & Rasio (2001). That study considered the interaction between two planets of mass $`10^3M_{}`$ (on the order of a Jupiter mass) and in an initially close configuration. The resulting branching ratios of system outcomes was explored. Systems became unstable by ejecting a planet or by a collision, one between planets or one between a planet and a star. Stable systems remained or settled into quasi-periodic orbits over $`2`$ Myr.

The initial configuration of both planets in Ford, Havlickova, & Rasio’s (2001) simulations is motivated by the “Hill Stability Criterion”, an analytic result first applied to planetary systems by Gladman (1993). A system is said to be Hill stable if the planets cannot approach each other closely for all time. We adopt Gladman’s (1993) notation: $`\mathrm{\Delta }`$ represents the least separation for which both planets, for sure, will be Hill stable. Nothing can be said about the Hill stability of a system for a separation less than $`\mathrm{\Delta }`$. $`\mathrm{\Delta }_{ch}`$ represents the greatest separation at which “global chaos” will occur. A system that is “globally chaotic” might produce collisions or ejections. For a separation greater than $`\mathrm{\Delta }_{ch}`$, the planets might exhibit stable quasi-periodic orbits. The geometry is illustrated in Fig. 4. Gladman’s (1993) analytic second order expansion for small values of $`\mu _1`$ and $`\mu _2`$ yields an approximate expression for $`\mathrm{\Delta }`$,

$$\mathrm{\Delta }23^{\frac{1}{6}}(\mu _1+\mu _2)^{\frac{1}{3}}+23^{\frac{1}{6}}(\mu _1+\mu _2)^{\frac{2}{3}}\frac{11\mu _1+7\mu _2}{3^{\frac{11}{6}}(\mu _1+\mu _2)^{\frac{1}{3}}}.$$

Unlike $`\mathrm{\Delta },`$ an approximation for $`\mathrm{\Delta }_{ch}`$ can only be found empirically. By developing overlap resonance criteria for two-planet systems, Wisdom (1980) derived the following approximation when $`\mu _1=\mu _2=\mu `$:

$$\mathrm{\Delta }_{ch}2\mu ^{\frac{2}{7}}$$

Ford, Havlickova, & Rasio (2001) ran simulations where the initial semimajor axis separation of the planets lay between $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }_{ch}`$, so that the planets would neither be in a Hill Stable configuration nor become unstable immediately. In this work, we sample the entire initial separation range spanned by $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }_{ch}`$ in order to best detect planets that may migrate outwards and remain on quasi-periodic orbits.

#### 3.2.2 Simulation Setup

We denote the semimajor axis of the initially inner planet $`a_1`$, the initially outer planet $`a_2`$, and the respective planet mass/central mass ratios as $`\mu _1`$ and $`\mu _2`$. Other orbital elements will distinguish the planets with a subscript of “1” or “2”. All runs were performed with $`a_1=1`$ AU, so that the results may be scaled easily to any ratio of semimajor axes, such as for $`a_1=510`$ AU.

All numerical runs of the three-body system were performed with a Bulirsch-Stoer routine from the HNbody integration package (Rauch & Hamilton 2002). The routine ran for 2 Myr with an accuracy parameter of $`10^{12}`$, and was rerun with uniformly smaller initial timesteps when necessary until this accuracy was achieved. For each run, the Bulirsch-Stoer routine began with an initial timestep of $`0.05`$ yr, and the orbital elements of each body were output every $`0.01`$ Myr. In most cases, energy and angular momentum errors, expressed by $`(EE_0)/E_0`$ and $`(L_ZL_{Z_0})/L_{Z_0}`$, where $`E_0`$ is the initial total system energy and $`L_{Z_0}`$ is the initial total angular momentum in the direction perpendicular to the orbital plane, did not exceed $`10^7`$. In the most pathological cases, energy and z-angular momentum errors were conserved to within $`10^4`$. Angular momenta in the other two directions were typically conserved to two order of magnitudes better than the z-angular momentum.

We are most interested in the evolution of stable systems in which both planets remain bound. We define stable systems as systems where both planets have semimajor axes which never exceed $`10^3`$ AU and eccentricities which never exceed $`0.99`$. Further, we deemed a system unstable if at any time the HNbody code outputted a negative value for the semimajor axis or eccentricity of either planet.

A preliminary exploration of phase space revealed that of the systems that become unstable, most did so within 0.05 Myr. Thus, each system underwent several checks at 0.05 Myr; unstable systems were terminated, and stable systems were evolved for a total of 2 Myr. In order to reduce computer time, we imposed additional conditions on systems after 0.05 Myr. Systems not satisfying these conditions were terminated. Given that $`a_1=1`$ AU for each run, the conditions are

1) 1.2 AU $`<a_2<`$ 3.0 AU,

2) $`a_1<1.05`$ AU,

3) If $`a_2/a_1<1.5`$, then $`e_2<0.1`$,

4) If $`a_2/a_1>1.5`$, then $`e_2<0.3`$.

These conditions were chosen based on a preliminary exploration of the properties of unstable systems. For example, we found that after 0.05 Myr, if the semimajor axis of the outer planet is within 50% of $`a_1`$, then a high ($`>0.1`$) eccentricity of the outer planet implies that the system will become unstable.

By definition, long-term chaotic behavior may differ drastically due to an infinitesimal change in initial conditions. Therefore, because computers use finite-precision arithmetic, individual runs are largely irreproducible from machine to machine. In this context, one can only make statements about the probability of a system behaving in a certain manner (Quinlan & Tremaine 1992). Thus, for each particular dimension of phase space explored, we performed 300 runs, each with randomly chosen initial orbital parameters that lay within specific ranges.

We expanded on the results of Ford, Havlickova, and Rasio (2001) by investigating the effects of altering the planet/star mass ratio for each planet. For each pair of mass ratios, 300 runs were performed, each run with both planets having a randomly chosen initial eccentricity between 0 and 0.01 and initial inclination between $`0^{}`$ and $`5^{}`$, along with an argument of perihelion, longitude of ascending node, and mean anomaly each randomly chosen between $`0^{}`$ and $`360^{}`$.

We performed 9 sets of 300 runs for each of the following four outer planet mass ratios $`\mu _2`$: $`1\times 10^4,`$ $`5\times 10^4,`$ $`1\times 10^3,`$ and $`5\times 10^3.`$ Each set of runs corresponds to a different inner planet mass ratio. We varied the mass of the inner planet within one order of magnitude of the mass of the outer planet. As mentioned, the initial semimajor axis of the inner planet for every run was set at 1.0 AU. The initial semimajor axis of the outer planet was chosen in a regime that exhibited unpredictable behavior but allowed for quasi-periodic orbits. Explanation of this regime follows.

As mentioned previously, the region bounded by $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }_{ch}`$ describe systems that $`\mathrm{𝑚𝑎𝑦}`$ exhibit Hill instability and $`\mathrm{𝑚𝑎𝑦}`$ exhibit quasi-periodic motion. When both are exhibited, $`\mathrm{𝑠𝑜𝑚𝑒}`$ planets exhibit outward migration, which is the topic of this section. So as to include all regions from which migrating behavior might arise, we took $`\mathrm{\Delta }_{ch}=2[min(\mu _1,\mu _2)]^{\frac{2}{7}}.`$ Because $`a_1=1`$ AU for each run, values of $`a_2`$ were randomly chosen between $`1+\mathrm{\Delta }_{ch}`$ and $`1+\mathrm{\Delta }`$ for the simulations.

Our initial conditions for the scattering experiments are similar to those adopted by most previous authors (e.g. Ford, Havlickova & Rasio 2001; Marzari & Weidenschilling 2002; Chambers, Wetherill & Boss 1996), in that we start with fully-formed planets in relatively close proximity to each other, and consider only their mutual gravitational interactions. Such initial conditions are chosen primarily for simplicity and compatibility with previous work, and can only partially be justified on physical grounds. In particular, since many of the trial systems prove to be unstable over time scales comparable to the lifetime of protoplanetary discs (Haisch, Lada & Lada 2001), there is likely to be an earlier epoch during which both planet-planet gravitational interactions and planet-disc interactions are important. This earlier stage of evolution could lead to an unstable multiple planet system in at least two ways. First, two planets might form in well-separated orbits, but then subsequently migrate into an unstable configuration as a consequence of disc-driven migration (Goldreich & Tremaine 1980; Lin & Papaloizou 1986; Lin, Bodenheimer & Richardson 1996). For this to happen, the planets would have to avoid becoming trapped into resonance during the migration process (Snellgrove, Papaloizou & Nelson 2001; Lee & Peale 2002; Murray, Paskowitz & Holman 2002). Alternatively, the planets might form close together, and be stabilized against immediate violent instability by the presence of a surrounding gas disc (Lin & Ida 1997; Nagasawa, Lin & Ida 2003)333In related work, Agnor & Ward (2002) considered the damping of terrestrial planet eccentricity by a remnant gas disc.. Further, unstable planetary systems may form after the disc has dissipated. Planet-planet interactions alone in crowded systems of Jovian-mass planets most often leaves only two survivors in close, quasiperiodic but ultimately unstable orbits (Adams & Laughlin 2003). Although any of these pathways could plausibly lead to initial conditions that are similar to those which we (and other authors) have assumed, a full treatment will obviously need to model the more complex interactions that are possible during the phase when the gas disc is being dispersed.

### 3.3 Results

Table 1 and Figs. 6-9 displays the data we took. Table 1 summarizes the data in terms of mass ratio and extent of migration, while the figures provide a more detailed look at migration behavior. Fig. 5 provides a representative sample of the final semimajor axes and eccentricities of both planets when significant outward migration occurs, and provides a perspective on how data is plotted in Figs. 6-9. For example, the 43 stable systems displayed in Fig. 5 are plotted as the solid lines in Fig. 8. Fig. 5 illustrates that planets that migrate both inward and outward tend to increase their eccentricity. Although the extent of migration in both directions is similar, those planets that migrate inward tend to increase their eccentricities at a higher rate.

Figs. 6-9 each represent a set of six cumulative probability distributions that display the fraction of stable systems vs. semimajor axis ratios, eccentricities, and inclinations of the outer and inner planet after exactly 2 Myr. The upper panel in each figure represents the initially outer planet, and the lower panel represent the initially inner planet. Each figure keeps the mass of the outer planet fixed, but varies the mass of the inner planet. Although we obtained nine curves for each graph on each figure, we show only the curves that exhibit significant radially outward migration, plus some that do not, and we do not show any curves where less than 10% of the 300 runs were stable.

Detecting migration for the outer planets is difficult because their semimajor axes were randomly chosen for each run. To aid in determining the net radial movement of these planets, we define $`\alpha _1a_1(\mathrm{final})/a_1(\mathrm{initial})`$ and $`\alpha _2a_2(\mathrm{final})/a_2(\mathrm{initial})`$. Each panel in Figs. 6-9 contain distributions of semimajor axis ratios, rather than absolute semimajor axes, because those are scalable in this study.

Table 1 and Figs. 6-9 illustrate that outward migration of both the inner and outer planets occurs more frequently with smaller values of $`\mu _2`$. Fig. 6 displays extensive migration. In the lowest mass case, almost 60% of the stable inner planets migrate at outward to at least 150% of their original semimajor axis. As $`\mu _1`$ is increased, migration occurs less frequently; for $`\mu _1=2.5\times 10^4`$ no inner planets migrated out to 150% of their original semimajor axis. The correlation of mass ratio to final state behavior is less apparent but still present in the final eccentricity and inclination curves. In both graphs, the least-mass case prompts the highest final eccentricity and inclination values for the orbit of the inner planet.

The final state of the outer planet, represented by the lower panel of Fig. 6, reflects the exchange of angular momentum in this system. One sees that when an inner planet migrates outward, it leaves the initially outward planet back so that $`\alpha _2<1`$. The one instance where the outer planet experiences significant migration is the highest mass case - the same case where $`\alpha _1<1.5`$. The eccentricity and inclination distributions for the outer planet in this massive case differ drastically from the behavior seen for other inner planet masses.

In Fig. 7, $`\mu _2=5\times 10^4`$ and the result is significantly less migration than seen in Fig. 6. Similarly to Fig. 6, the highest mass case exhibits the greatest migration of the outer planet. In this case, $``$20% of the outer planets and $``$15% of the inner planets migrate beyond 150% of their original semimajor axis. Halving the inner mass value reduces the probability for $`\alpha _1<1.5`$ or $`\alpha _2<1.5`$ to $``$5%. In contrast to Fig. 6, the least mass case fails to show significant outward migration of the inner planet, despite the inner planet exhibiting the highest final inclinations.

In Fig. 8, where the outer mass is a full order of magnitude more massive than in Fig. 6, we see even less migration. Inner bodies for the smallest mass case exhibit high eccentricities and inclinations, but only $``$5% satisfy $`\alpha _1>1.5`$. The case where $`\mu _1=\mu _2=1\times 10^3`$ is the same case studied by Ford, Havlickova, & Rasio (2001). Their Figs. 11-12 can be compared to the semimajor axis ratio distributions, eccentricity distributions, and inclination distributions in Fig. 8.

No migration of the inner planet occurs for any case in Fig. 9. However, the outer planet shows significant migration. Further, as shown by the sequence of four curves in the lower panel, the higher the mass of the outer planet, the greater the extent of the migration. For the higher mass cases, about 25% of the initally outer planets satisfy $`\alpha _2>2.5`$. Remarkably, these planets all retain final eccentricities $`<0.4`$, in stark contrast to the large final eccentricities of the migrating planets of Figs. 6-8.

Because of the chaotic nature of the three body problem, one encounters difficulty when predicting the appropriate timescale over which to integrate. We adopted a similar timescale to the one used by Ford, Havlickova, & Rasio (2001), but recognize that systems that appear to be stable at 2 Myr might become unstable at some future time. We extended the running time of one system to sample the consequences. After 2 Myr, 55 out of the 300 systems for the case $`\mu _1=5\times 10^3,\mu _2=1\times 10^3`$ remain stable. By running this system for 10 Myr, we found only 3 out of the 55 systems became unstable, and did so quickly. Further, these three systems became unstable before 3 Myr.

For any set of 300 runs of a chaotic system, the slightest change in any input parameter of the Bulirsch-Stoer algorithm might drastically alter the results of any individual runs, but keep the same global behavior. The extent of the invariance of this global behavior is a function of the code used and the number of systems sampled. To explore this measure of invariance for the runs presented here, we reran two systems with four different initial tolerances ($`5\times 10^{11}`$, $`1\times 10^{11}`$, $`5\times 10^{12}`$, $`1\times 10^{12}`$). One of these systems ($`\mu _1=\mu _2=5\times 10^3`$) exhibited negligible migration, and the other ($`\mu _1=1\times 10^5,\mu _2=1\times 10^4`$) exhibited significant migration. In order to achieve the largest feasible sample size for this error analysis, we tripled the number of runs performed to 900 for the system exhibiting negligible migration. As Fig. 10 illustrates, runs with different tolerances are practically indistinguishable for the case of no migration, but vary up to 20% for the case of significant migration.

### 3.4 Summary

The gravitational interaction of a pair of unequal mass planets that lie in close initial configurations allow either planet to migrate outward to at least twice it’s initial semimajor axis. Planets that migrate outward lie in quasi-stable, slightly inclined eccentric orbits with an eccentricity that spans the entire permissible range and an inclination up to about $`10^{}`$. Although either planet in any given system may drift outward, the smaller the mass of the planets, the higher the tendency for the initially inner planet to migrate outward. Further, less massive giant planets such as Uranus or Neptune are more likely to migrate outward than Jupiter-mass planets. This result, derived from our gravitational scattering simulations alone, is consistent with Thommes, Duncan & Levison’s (1999, 2002) conclusion that Uranus and Neptune’s current location is a result of outward migration amongst Jupiter and Saturn.

## 4 Conclusions

As a first step in explaining the presence of planets at large orbital radii, we have shown that outward migration of protoplanets is possible both by planet-disc interactions and by planet-planet gravitational scattering without the presence of a disc. Strong mass loss in discs coupled with planetary cores that are formed at about $`10`$ AU allow planets to migrate outward in discs to radii that are as much as a factor of several in excess of their initial semimajor axes. Planets that migrate in such a manner are likely to be massive. We predict that gas-driven outward migration should be most likely to occur around more massive stars, whose strong UV flux can drive a powerful photoevaporative outflow. Subsequently, when in the appropriate chaotic regime, planets within a multiple planet system may migrate outward due to gravitational scattering alone. Planets that migrate in this manner may be massive or not, however low mass objects tend to exhibit the most extensive outward migration. Orbital migration is typically accompanied by an increase in eccentricity that spans the allowable range for elliptic orbits.

## Acknowledgments

We thank an anonymous referree for helpful comments, Kevin Rauch and Doug Hamilton for use of their integrator and assistance in it’s operation, as well as the Colorado Rings Group for useful discussions.

This paper is based upon work supported by NASA under Grant NAG5-13207 issued through the Office of Space Science.


text:
# Molecular Discreteness in Reaction-Diffusion Systems Yields Steady States Not Seen in the Continuum Limit

(August 5, 2004)

## Abstract

We investigate the effects of spatial discreteness of molecules in reaction-diffusion systems. It is found that discreteness within the so called Kuramoto length can lead to a localization of molecules, resulting in novel steady states that do not exist in the continuous case. These novel states are analyzed theoretically as the fixed points of accelerated localized reactions, an approach that was verified to be in good agreement with stochastic particle simulations. The relevance of this discreteness-induced state to biological intracellular processes is discussed.

Many systems in nature that involve chemical reactions can be studied with the help of reaction-diffusion equations. For certain processes, a relatively small number of suitably chosen continuous macroscopic variables yields excellent descriptive results. In biological systems, however, not only is the variety of chemicals enormous, the number of molecules of each of the chemical species can range from the relatively very large to the relatively very small. Now, if the species with small numbers of molecules were irrelevant, obviously, their existence could be ignored and one could focus on the species with large numbers of molecules that can effectively be described by a continuous variable. However, it should not really come as a surprise that it was found that, in general, species with small numbers of molecules cannot be neglected and that certain functions in cells can critically depend on very small fluctuations Mikhailov ; Blumenfeld . Indeed, in prior studies on reaction-diffusion systems some effects of fluctuations on pattern formation were found (see e.g., fluctuation1 ; fluctuation2 ). Stochastic differential equations are often used to study effects of fluctuations.

Of course, on a microscopic level chemicals are composed of molecules, and the actual reactions occur between these molecules. Therefore, in principle, reaction events must be integer and change only discretely. In analysis with stochastic differential equations, though, the fluctuations are regarded as continuous changes. Clearly, this approximation can only be valid if applied to fluctuations that involve sufficiently large numbers of molecules and should not be applied when relevant chemical species are very rare.

In order to address this issue, we previously studied the effects of discreteness in simple autocatalytic reaction network systems and reported discreteness-induced transitions as well as drastic effects on concentrations YTKK2001 ; YTKK2003 . A key feature of these systems was, however, that the medium was assumed to be well-stirred.

In contrast, in a system with diffusion in space, the total number of molecules may vary from point to point. By assuming that the reaction is fast and the diffusion is slow, locally, the discreteness of the molecules can become important. In fact, this can even be the case if the total number of molecules is large but spread out over a large area as well.

Therefore, a length scale should be considered such that it can serve as a benchmark for judging whether or not a continuum approximation is applicable. To consider this problem, the ratio between the reaction and diffusion rates is important and a candidate for the length scale is the typical distance over which a molecule diffuses during its lifetime, i.e., before it undergoes reaction as defined by Kuramoto Kuramoto1 ; Kuramoto2 . For reference, let us briefly review the work.

Consider the reaction 111In Kuramoto1 ; Kuramoto2 , this is expressed as $`A+MX+M`$, $`2XE+D`$, where, the concentrations of $`A`$, $`M`$, $`E`$ and $`D`$ are taken to be constant in order to make an analysis of the equations possible.

$$A\stackrel{k}{}X,2X\stackrel{k^{}}{}B.$$

If the concentration of $`A`$ is set to be constant, $`X`$ is produced at a constant rate $`k`$ while decaying by the reaction $`2XB`$ at a rate $`k^{}`$. The average concentration of $`X`$ at the steady state is $`X=\sqrt{kA/2k^{}}`$, where, for simplicity, $`A`$ is the concentration of the chemical $`A`$. Thus the average lifetime of $`X`$ at the steady state is estimated to be $`\tau =1/(2k^{}X)=1/\sqrt{2kk^{}A}`$. Suppose that $`X`$ molecules diffuse with the diffusion constant $`D`$. The typical length over which an $`X`$ molecule diffuses in its lifetime is then estimated to be

$$l=\sqrt{2D\tau },$$
(1)
which is called the Kuramoto length Kampen .

The Kuramoto length $`l`$ represents the relation between the reaction rate and the diffusion rate. When the system size is smaller than $`l`$, its behavior is dominated by diffusion and local fluctuations rapidly spread throughout the system. Contrastingly, if the system size is much larger than $`l`$, fluctuations are localized only in a small part of the system, and distant regions fluctuate independently.

In this reasoning, it is assumed that the average distance between molecules is much smaller than $`l`$. Thus the actual discreteness of the molecules can be ignored, and the concentration of the chemical $`X`$ can be regarded as a continuous variable. However, if the average distance between molecules is comparable to or larger than $`l`$, local discreteness of molecules may not be negligible. Suppose a chemical $`A`$, with very low concentration, produces another chemical $`B`$. The average lifetime of $`B`$ is short, such that the Kuramoto length of $`B`$ is shorter than the average distance between adjacent $`A`$ molecules. With this setting, chemical $`B`$ molecules may be considered as localized around $`A`$ molecules. This is especially so if the reactions involve 2nd or higher orders of $`B`$. Then the localization of chemical $`B`$ may drastically alter the total rate of the reactions, and the effect of the local discreteness of the molecules may thus be rather significant.

In order to systematically investigate the effects of the local discreteness of the molecules, we consider a simple one-dimensional reaction-diffusion system with 3 chemicals ($`X_1`$, $`X_2`$, and $`X_3`$) and the following 4 reactions

$`X_2+X_3\stackrel{k_1}{}X_2+X_1;X_3+X_1\stackrel{k_2}{}2X_3`$
$`2X_2\stackrel{k_3}{}X_2+X_1;2X_1\stackrel{k_4}{}X_1+X_2.`$
Here, we assume that the first two reactions are much faster than the others, i.e., the reaction constants satisfy $`k_1,k_2k_3>k_4`$. To be specific, we take $`k_1=k_2=100r`$, $`k_3=ar`$, and $`k_4=r`$ ($`r>0`$, $`1<a100`$).

In the continuum limit, $`c_i(t,x)`$, the concentration of chemical $`X_i`$ at time $`t`$ and position $`x`$, is governed by the reaction-diffusion equation for the system given by

$`{\displaystyle \frac{c_1}{t}}`$ $`=`$ $`100r(c_1c_2)c_3r(c_1^2ac_2^2)+D_1{\displaystyle \frac{^2c_1}{x^2}}`$ (2)
$`{\displaystyle \frac{c_2}{t}}`$ $`=`$ $`r(c_1^2ac_2^2)+D_2{\displaystyle \frac{^2c_2}{x^2}}`$ (3)
$`{\displaystyle \frac{c_3}{t}}`$ $`=`$ $`100r(c_1c_2)c_3+D_3{\displaystyle \frac{^2c_3}{x^2}}`$ (4)
where $`D_i`$ is the diffusion constant of $`X_i`$. The system is closed and thus the total concentration $`c`$ is conserved. For simplicity, we assume $`D_i=D`$ for all $`i`$.

The reaction-diffusion equation has fixed points at $`(c_1,c_2,c_3)=(0,0,c),(\sqrt{a}c/(\sqrt{a}+1),c/(\sqrt{a}+1),0)`$ for all $`x`$. By performing a straightforward linear stability analysis, it is shown that only the former is stable. Indeed, by starting from an initial condition with $`c_i>0`$, this reaction-diffusion equation always converges to the fixed point $`(0,0,c)`$.

Next, in order to obtain insights into the case when the continuum limit cannot be taken we carry out direct particle simulations. Each molecule diffuses randomly (showing Brownian motion) in a one-dimensional space with periodic boundary conditions (length $`L_x`$). When two molecules are within a distance $`d_r`$ they react with a certain probability and the total number of molecules ($`N`$) is conserved.

First, we investigate the case with $`a=4`$ and show time series of the number of molecules $`N_i`$ of chemical species $`X_i`$ in Fig. 1. As can be seen, $`N_1`$ and $`N_2`$ do not converge to $`0`$ but to relatively large numbers. As can be expected the final concentrations depend on $`r`$ and $`D`$ and for $`X_2`$ it is depicted in Fig. 2. Approximately, the concentration turns out to be proportional to $`\sqrt{r/D}`$ when $`N_1,N_2N`$.

To elucidate the origin of this proportionality, we take a closer look at the Kuramoto length, which, of course, depends on the molecule species. In the case of the $`X_1`$ molecules it is given by $`l_1=\sqrt{D/50rc_3}`$, as the average lifetime of $`X_1`$ is $`1/100rc_3`$. Here we consider the situation $`N_1,N_2N`$, so that $`c_3c`$. In the discussion below, we assume that $`l_1=\sqrt{D/50rc}=\sqrt{DL_x/50rN}`$.

Using this length $`l_1`$, the density of the remaining $`X_2`$ molecules is found to be about $`0.1`$ molecule per $`l_1`$, independent of the parameters, as shown in Fig. 2. After relaxation, this density does not depend on the initial conditions, as long as $`N_i1`$ is satisfied initially. Furthermore, the density is independent of the system size $`L_x`$, if $`L_xl_1`$, so that the number of remaining molecules $`N_2`$ is simply proportional to $`L_x`$. Consequently, in this analysis one obtains a finite $`c_2`$ regardless of the system size or initial conditions which is clearly different from the continuum limit where $`c_2`$ goes to $`0`$.

In this system, $`X_1`$ molecules are produced by $`X_2`$ molecules. If $`\lambda _2`$, the average distance between $`X_2`$ molecules, is smaller than $`l_1`$, the distributions of $`X_1`$ around neighboring $`X_2`$ molecules overlap each other significantly and one can regard $`X_1`$ to be uniformly distributed. In contrast, if $`\lambda _2`$ is much larger than $`l_1`$, molecules $`X_1`$ will localize around the $`X_2`$ molecules (The size $`L_x\lambda _2`$). Then, the reaction $`2X_1X_1+X_2`$ is accelerated when compared to the case that the same total number of $`X_1`$ molecules is uniformly distributed.

We define the acceleration factor $`\alpha (\lambda _2,l_1)`$ as the ratio between the reaction rate with localized $`X_1`$ and the reaction rate with uniformly distributed $`X_1`$. If $`\lambda _2l_1`$, it is expected that $`\alpha 1`$. Assuming that the distribution of $`X_1`$ is continuous and represented by the concentration $`c_1(x)`$ 222Here, only the $`X_1`$ species is relevant to this reaction, so that it is not necessary to consider detailed structures smaller than the typical distance between $`X_1`$ molecules and the total rate of the reaction can therefore be described by a smoothened distribution., the acceleration factor can be expressed as

$$\alpha =\frac{c_1^2}{c_1^2}=\frac{L_x^1c_1^2𝑑x}{\left(L_x^1c_1𝑑x\right)^2}.$$
(5)
For simplicity, we assume that the distribution of the localized $`X_1`$ molecules is Gaussian with a standard deviation $`l_1`$ centered around the $`X_2`$ molecules (which may overlap each other). Suppose that the $`X_2`$ molecules are randomly distributed over the system with an average distance $`\lambda _2`$, we then obtain 333The acceleration factor $`\alpha `$ is estimated as follows. We assume that the distribution of localized $`X_1`$ molecules is Gaussian with a standard deviation $`l_1`$ around the $`X_2`$ molecules. I.e. $`\rho _i(x)=(\sqrt{2\pi }l_1)^1\mathrm{exp}((xx_i)^2/2l_1^2)`$, where $`x_i`$ is the position of each $`X_2`$ molecule. The total distribution (concentration) of $`X_1`$ is $`c_1(x)=_i\rho _i(x)`$, and $`c_1=\rho _i(x)𝑑x/\lambda _2=1/\lambda _2`$. Since the molecules $`X_2`$ are randomly distributed, $`c_1^2=(\rho _i)^2=\left(\rho _i\right)^2+\rho _i^2=c_1^2+(2\sqrt{\pi }l_1)^1c_1`$
($`L_xl_1,\lambda _2`$). Thus, $`\alpha =c_1^2/c_1^2=1+(2\sqrt{\pi }l_1)^1c_1^1=1+\lambda _2/(2\sqrt{\pi }l_1)`$.
Consequently, we obtain $`\alpha =1+{\displaystyle \frac{1}{2\sqrt{\pi }}}{\displaystyle \frac{\lambda _2}{l_1}}`$.

$$\alpha =1+\frac{1}{2\sqrt{\pi }}\frac{\lambda _2}{l_1}=1+\frac{1}{2\sqrt{\pi }l_1c_2}.$$
(6)
On the other hand, the average lifetime of $`X_2`$ molecules is much longer, so that the Kuramoto length for $`X_2`$ molecules is longer than $`\lambda _2`$. Consequently, the reaction $`2X_2X_2+X_1`$ is not accelerated by localization.

Provided that $`N_1,N_2N_3`$, $`N_1N_2`$ due to the fast reactions $`X_2+X_3X_2+X_1`$ and $`X_3+X_12X_3`$. As a result, the ratio between the two reaction rates is given by

$$\frac{\text{The rate of }(X_1X_2)}{\text{The rate of }(X_2X_1)}\frac{\alpha k_4N_1^2}{k_3N_2^2}\frac{\alpha }{a}.$$
(7)
Following eq. (7), the two reaction rates are balanced if $`N_2`$ takes a value such that $`\alpha =a`$ is satisfied. Corresponding to $`\alpha =a`$, a novel fixed point appears at

$$c_1=c_2=\left(2(a1)\sqrt{\pi }l_1\right)^1(=c_s),$$
(8)
provided $`c_1,c_2c_3`$ and $`c_3=c`$. The stability of this fixed point is analyzed, by linearizing eqs. (6) and (8) around the fixed point. Noting that

$$\alpha =1+\frac{(a1)c_s}{c_2}=a\frac{a1}{c_s}\delta c_2+o(\delta c_2),$$
(9)
with $`c_1=c_s+\delta c_1`$ and $`c_2=c_s+\delta c_2`$, and rewriting eqs. (2) and (3) with $`\alpha `$ in eq. (9), we obtain

$`\left(\begin{array}{c}\dot{c_1}\\ \dot{c_2}\end{array}\right)`$ $`=`$ $`r(\begin{array}{cc}2ac_s100c& (3a1)c_s+100c\\ 2ac_s& (3a1)c_s\end{array})\left(\begin{array}{c}\delta c_1\\ \delta c_2\end{array}\right)`$ (17)
$`+o(\delta c_1,\delta c_2).`$
The Jacobi matrix has two negative eigenvalues, and the fixed point is stable (This is natural, since if $`\alpha <a`$, $`N_2`$ decreases, leading to the increase of $`\alpha `$, and vice versa). This fixed point (steady state) is distinct from that of the original reaction-diffusion equation, $`(0,0,c)`$.

From eq. (6), $`\alpha `$ becomes $`4`$ when $`\lambda _2/l_1=6\sqrt{\pi }10.6`$. In our simulation with $`a=4`$, about $`0.1`$ $`X_2`$ molecule per $`l_1`$ remains, as shown in Fig. 2. In other words, $`\lambda _2/l_110`$, in good agreement with the estimate.

By changing $`a`$, we numerically obtain the relation between the $`\lambda _2/l_1`$ and the actual acceleration factor $`\alpha `$, again agreeing well with the above theoretical estimate $`\alpha =1+\frac{1}{2\sqrt{\pi }}\frac{\lambda _2}{l_1}`$, as shown in Fig. 3.

In the estimate above, we consider the case that $`N_1,N_2N`$. On the other hand, if $`N`$ is set to be smaller than the estimated value of $`N_2`$ at the steady state, $`N_2`$ increases to satisfy the balance, and finally reaches the state $`N_1+N_2=N`$, $`N_3=0`$, which corresponds to the unstable fixed point of the reaction-diffusion equation, $`(\sqrt{a}c/(\sqrt{a}+1),c/(\sqrt{a}+1),0)`$, as shown in Fig. 1 (a).

The localization of $`X_1`$ cannot be maintained without the spatial discreteness of $`X_2`$ molecules. In reaction-diffusion equations, any pattern will disappear eventually given a sufficiently long evolution time unless it is somehow sustained. This is even the case when the initial distribution of $`X_2`$ is discrete. But again, it is essential to recall that reaction-diffusion equations are an approximation and in that sense an idealization. In reality, a single molecule itself can of course not be broadened by diffusion and the spatial discreteness of $`X_2`$ molecules is always maintained. By itself, a molecule is a diffusion-resistant pattern.

The alteration of the steady state due to localization is not limited to the present type of reaction network. Provided that the conditions

1. Chemical $`A`$ generates another chemical species $`B`$.
2. The lifetime of $`B`$ is short or the diffusion of $`B`$ is slow so that the Kuramoto length of $`B`$ is much smaller than the average distance between $`A`$ molecules.
3. The localization of molecule $`B`$ accelerates some reactions.

are satisfied, discreteness may alter the dynamics. The last condition is easily satisfied if species $`B`$ is involved in second or higher order reactions. Finally, if

1. The acceleration alters the density of $`A`$ molecules,

the above acceleration mechanism may control the density of $`A`$ to produce a novel steady state.

As for the localization effect by the discreteness of catalytic molecules, Shnerb et al. recently showed that it can amplify autocatalytic reaction-diffusion processes Solomon2000 ; Louzoun . In their model, however, the density of the catalyst is fixed as an externally given value, and the concentration of the product, localized around the catalyst, diverges in time. In our mechanism, the density of the catalyst ($`A`$, or $`X_2`$) changes autonomously and reaches a suitable value to produce the discreteness effect. Hence the effect of discreteness is controlled by the discreteness itself, leading to a novel steady state. Indeed, theoretical estimates for the novel concentrations based on the self-consistent fixed point of acceleration due to the localization agree well with numerical results.

In so far as the conditions (i)–(iv) are met, our result does not depend on the details of the reactions, and should generally be valid for reaction-diffusion systems. We have carried out simulations of similar reaction-diffusion systems, and again the discreteness effect led to novel pattern formation that cannot be accounted for by Turing type mechanisms (with or without noise).

Experimental verification of our results should be possible by suitably designing a reaction system, with the use of, say, microreactors or vesicles. Also, in biological cells, many chemicals work at low concentrations on the order of 1 nM or less. Furthermore, diffusion is sometimes restricted, e.g. due to surrounding macro-molecules, and may be slow. In such an environment, it is probable that the average distance between the molecules of a given chemical species is much larger than the Kuramoto lengths of some of the other chemical species. Indeed, biochemical systems contain various higher order reactions and positive feedback mechanisms that might naturally support the conditions (iii)–(iv) above.

###### Acknowledgements.

This research is supported by grants-in-aid for scientific research from the Ministry of Education, Culture, Sports, Science and Technology of Japan (11CE2006, 15-11161). One of the authors (Y.T.) is supported by a research fellowship from Japan Society for the Promotion of Science.