Escape dynamics through a continuously growing leak

TL;DR

This paper models escape dynamics in a leaky chaotic system with a leak size that depends on the number of particles, relevant to astrophysical planetary accretion, and provides both numerical and analytical insights.

Contribution

It introduces a model where the leak size varies with particle number and derives an analytical solution, extending understanding of escape dynamics in such systems.

Findings

Early phase shows deviation from exponential decay

Analytic solution matches classical results in limiting cases

Model applicable to astrophysical planetary accretion

Abstract

We formulate a model that describes the escape dynamics in a leaky chaotic system in which the size of the leak depends on the number of the in-falling particles. The basic motivation of this work is the astrophysical process which describes the planetary accretion. In order to study the dynamics generally, the standard map is investigated in two cases when the dynamics is fully hyperbolic and in the presence of KAM islands. In addition to the numerical calculations, an analytic solution to the temporal behavior of the model is also derived. We show that in the early phase of the leak expansion, as long as there are enough particles in the system, the number of survivors deviates from the well-known exponential decay. Furthermore, the analytic solution returns the classical result in the limiting case when the number of particles does not affect the leak size.

Tables

Table 1. Table 1: The integral on the left-hand side of ( 20 ) expressed by elementary functions.

$\gamma$	$\int \frac{1}{x^{\gamma }}\cdot \frac{1}{1-x}dx$
0	$-\ln \left(1-x\right)$
1/2	$2{\tanh }^{-1}\left(\sqrt{x}\right)$
3/4	$\ln \left(\frac{1+\sqrt[4]{x}}{1-\sqrt[4]{x}}\right)+2{\tanh }^{-1}\left(\sqrt[4]{x}\right)$
1	$\ln \left(\frac{x}{1-x}\right)$
4/3	$\ln \left(\frac{\sqrt{1+\sqrt[3]{x}+\sqrt[3]{x^2}}}{1-\sqrt[3]{x}}\right)-\sqrt{3}{\tanh }^{-1}\left(\frac{1+2\sqrt[3]{x}}{\sqrt{3}}\right)-\frac{3}{\sqrt[3]{x}}$
3/2	$2{\tanh }^{-1}\left(\sqrt{x}\right)-\frac{2}{\sqrt{x}}$
2	$\ln \left(\frac{x}{1-x}\right)-\frac{1}{x}$

Equations

$$M_i=M_0+(N_0-N_i)\cdot m=\mathcal{M}-N_i\cdot m,$$

$$S_{\mathrm{leak}}(M_i)=C_S\cdot M_i^{\gamma },$$

$$p(M_i)=C_p\cdot M_i^{\gamma }.$$

$$\mathcal{N}=\frac{\mathcal{M}}{m},{\kappa }_{\infty }=C_p{\mathcal{M}}^{\gamma },$$

$$x_i=\frac{M_i}{\mathcal{M}},y_i=\frac{N_i}{\mathcal{N}}=N_i\frac{m}{\mathcal{M}},$$

$$x_i+y_i=1$$

$$\begin{array}{rl}{I}_{i+1}& =I_i+K\sin {\Theta }_i,\\ {\Theta }_{i+1}& ={\Theta }_i+I_{i+1}.\end{array}$$

$$\frac{\mathrm{d}N}{\mathrm{d}t}\approx \Delta N=-p\cdot N=-C_p\cdot M^{\gamma }\cdot N$$

$$\frac{\mathrm{d}M}{\mathrm{d}t}\approx C_p\cdot M^{\gamma }\cdot N\cdot m.$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}={\kappa }_{\infty }\cdot x^{\gamma }\cdot \left(1-x\right),$$

$$x_{i+1}=x_i+\Delta x_i$$

$$t(x)=\frac{x^{1-\gamma }}{{\kappa }_{\infty }\cdot (1-\gamma )}{}_2{F}_1\left(1-\gamma ,1;2-\gamma ;x\right)-\tau$$

$$\tau =\frac{x_0^{1-\gamma }}{{\kappa }_{\infty }\cdot (1-\gamma )}{}_2{F}_1\left(1-\gamma ,1;2-\gamma ;x_0\right).$$

$$\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^2}={\kappa }_{\infty }\cdot x^{\gamma -1}\cdot \frac{\mathrm{d}x}{\mathrm{d}t}\cdot \left[\gamma -(1+\gamma )\cdot x\right],$$

$$x_{\mathrm{PoI}}=\frac{\gamma }{1+\gamma }.$$

$$\displaystyle\lim_{x\rightarrow 0}t=\left\{\begin{array}[]{ll}-\tau&\;\;\;\textrm{for $0\leq\gamma<1$,}\\ -\infty&\;\;\;\textrm{for $1\leq\gamma$.}\end{array}\right.$$

$$t_{\mathrm{cross}}\sim 1/{\kappa }_{\infty }$$

$$t_{\mathrm{cross}}\sim \frac{1}{{\kappa }_{\infty }}\left[1+\frac{x^{1-\gamma }}{1-\gamma }+\frac{(1+\gamma {)}^{1+\gamma }}{{\gamma }^{\gamma }}\frac{\gamma }{1+\gamma }\right],$$

$$\begin{array}{rl}{M}_i& =M_0+(N_0-N_i)m=\mathcal{M}-N_{i}m,\\ S_{\mathrm{leak},i}& =C_S{M}_i^{\gamma },\\ p_i& =C_p{M}_i^{\gamma },\end{array}$$

$$r_{\mathrm{H}}=a{\left(\frac{\mu }{3}\right)}^{1/3}$$

$$v_{\mathrm{esc}}=\sqrt{\frac{4G{m}_{\mathrm{p}}}{r_{\mathrm{H}}}}$$

$$S_{\mathrm{leak}}=A_{\mathrm{r}}{A}_{\mathrm{v}}\sim r_{\mathrm{H}}^2{v}_{\mathrm{esc}}^2\propto m_{\mathrm{p}}^{4/3}.$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}={\kappa }_{\infty }\cdot x^{\gamma }\cdot \left(1-x\right),$$

$$\frac{1}{{\kappa }_{\infty }}\int \frac{1}{x^{\gamma }}\cdot \frac{1}{1-x}\mathrm{d}x=\int 1\mathrm{d}t.$$

$$\int \frac{1}{x}\cdot \frac{1}{1-x}\mathrm{d}x=\ln \frac{x}{1-x},$$

$$x(t)=\frac{1}{1+e^{-{\kappa }_{\infty }(t+\tau )}},$$

$$\tau =\ln \frac{x_0}{1-x_0}$$

$$\frac{1}{1-x}=\sum _{i=0}^{\infty }{x}^i.$$

$$\int \frac{1}{x^{\gamma }}\cdot \frac{1}{1-x}\mathrm{d}x=\int \sum _{i=0}^{\infty }{x}^{i-\gamma }\mathrm{d}x=$$

$$x^{1-\gamma }\sum _{i=0}^{\infty }\frac{x^i}{i-\gamma +1}=\frac{x^{1-\gamma }}{1-\gamma }\sum _{i=0}^{\infty }\frac{(1-\gamma {)}^{(i)}\cdot 1^{(i)}}{(2-\gamma {)}^{(i)}}\frac{x^i}{i!}=$$

$$\frac{x^{1-\gamma }}{1-\gamma }\cdot {}_2{F}_1\left(1-\gamma ,1;2-\gamma ;x\right)$$

Full text

PACS NUMBERS, AND KEYWORDS

leaky systems

††preprint: APS/123-QED

I Introduction

Simple nonlinear dynamical systems in which trajectories may escape through an artificial leak 111Artificial means in this context that if the leak is not present, escape cannot occur at that part of the phase space. placed in the phase space play an important role in recent studies. Various fields of physics deal with either the escape dynamics of the particles or the decay rate of other physical quantities such as sound intensity, light rays, or fractal eigenstates Schneider et al. (2002); Jung et al. (1993); Ernst and Peters (2014); Nagler (2005); Altmann and Tél (2009); Haszpra and Tél (2013); Portela et al. (2008). It has been pointed out that the escape dynamics strongly depends on the leak size, position, and orientation Lai et al. (1999); Życzkowski and Bollt (1999); Afraimovich and Bunimovich (2010); Bunimovich and Dettmann (2007); Dettmann and Leonel (2012); Dettmann (2013); Dettmann and Georgiou (2011) as well as on other pre-defined properties of the leak, for instance, the reflection coefficient Altmann et al. (2013a). Probably the most interesting question is how the escape dynamics changes if the size of the leak varies. Altmann et al. Altmann et al. (2013b) presented numerical results about the relation between the escape rate and the leak size. In their study, however, the measure of the leak was adjusted manually in each case. Recently, Livorati et al. Livorati et al. (2014) studied the escape in case of periodically driven holes. The main results of their work show parameter (amplitude, initial phases, and period of the oscillations) dependent fluctuations superimposed to the classical exponential decay.

Although mathematicians are interested mostly in the limiting case of vanishing small leaks Haydn et al. (2005); Keller and Liverani (2009); Nándori and Szász (2012), in this work we present the decay dynamics through a continuously growing leak, where the size of the leak depends on a given physical property of the escaping particles. The motivation of this study comes from the application of leaky chaotic systems de Assis and Terra (2014); Kovács and Regály (2015); Morrison and Malhotra (2015) and crash tests Zotos (2015, 2016) in dynamical astronomy discussed in details below.

The model of the growing leak introduced here results in a survival probability of non-escaped trajectories that is different from the well-known classical exponential decay Lai and Tél (2011); Tél and Gruiz (2006). Moreover, we found a simple analytical solution describing the escape dynamics until the leak’s expansion stops. A comprehensive numerical investigation is also performed to confirm our analytic results.

The paper is organized as follows. After the Introduction, in Section II, the motivation as an astrophysical application is described. Then, we give a detailed description of the model of a growing leak and its simple numerical implementation to the standard map. The mathematical background is presented in Sections III.1. Section III.2 is devoted to numerical calculations in order to compare analytic results and simulations. Finally, we discuss our results and draw some conclusions in Section IV.

II Model

II.1 Motivation

The motivation of the present study Goldreich and Ward (1973) is the so-called planetary accretion process which is one of the two competing planet formation scenarios in these days Matsuo et al. (2007). In this process the forming planetary embryo accretes particles from its vicinity until this region – the feeding zone 222The planetary feeding zone is basically the basin of attraction of a given leak where the leak in phase space can be considered as the forming planetesimal. – becomes empty. The increase of the planet depends on the mass of the particles hitting its surface. Obviously the smaller the embryo at the beginning of this process, the more significant the growth by the accretion. As a very simple model of this process one might consider the gravitational planar circular restricted three body problem (RTBP). In RTBP two point masses (star and planet) orbiting their barycenter on a circle and a third mass-less body (test particle) moves in their gravitational potential in the same plane. Although the planet (and also the star) is considered as a point mass, one can define the Hill radius ($r_{\mathrm{H}}$) in which its gravitational influence is dominant. The particles entering the Hill radius with an appropriate velocity, i.e. slower than the escape velocity from this domain, can be removed from the dynamics and marked as escaped. In addition, $r_{\mathrm{H}}$ grows with the mass of the forming planet, see Eq. (16). Therefore, the growth of the planetary embryo can be considered as a growing leak in the phase space. Thus, from dynamical point of view, the accretion stage of the planet formation can be described via leaky chaotic systems. We give an estimate how the leak size depends on the mass in RTBP, see Appendix A.

To illustrate the leaky RTBP, we plot the evolution of a large number of non-interacting test particles initially placed around the planet’s orbit (see Figure 1). Different colors denote different end-states of particles. Trajectories starting from light gray (green online) points remain the part of the system during the whole integration (1000 orbits of the planet). Gray (red online) points represent test particles whose destination is the planet, more precisely, the half of the Hill radius with proper velocity 333Since we consider the planet as a point mass, crash of the particles and the planetary embryo is difficult to calculate numerically. Consequently, half of the Hill radius is chosen as a region wherein the particles are thought to be accreted by the forming planet. This is, obviously, more rigorous criterion than one Hill radius.. Dark gray (blue online) points indicate trajectories scattered out from the system by the planet.

Although the effect of the planet’s mass and size evolution in the RTBP is dominant only in very early stages of the planet formation, the idea of a growing leak, particularly when the size of the leak depends on a physical property of the leaving particles, might shed light on a new kind of escape dynamics generally in leaky chaotic systems.

II.2 Growing leak model

The discrete dynamical system we are to consider here consists a large number of particles and a leak, where under certain conditions, the particles can escape from the system.

The particles are point masses with the same mass $m$, their initial number is $N_{0}$, while after $i$ iterations we denote the number of surviving particles by $N_{i}$. The leak also has an initial and an instantaneous mass, $M_{0}$ and $M_{i}$, respectively. When a particle falls into the leak, its mass is added to that of the leak, thus

[TABLE]

where $\mathcal{M}=M_{0}+N_{0}\cdot m$ is the total mass of the system.

According to the RTBP (Appendix A) a reasonable choice is that the volume of the leak depends on its mass $M_{i}$ in the form of

[TABLE]

where $\gamma$ is a positive constant. The coefficient $C_{S}$ can be written as $C_{S}=C_{P}\cdot S_{\mathrm{total}}$. Here $C_{P}>0$ denotes a normalization constant while $S_{\mathrm{total}}$ is the volume of the ergodic part of the phase space. The factor $C_{P}$ allows us to control the final size of the leak, (a leak of moderate size avoids excessive restructuring of the phase space). Let $p$ be the escape probability that a particle leaves the system (through the leak) in the next iteration. We suppose that the escape probability is proportional to the actual size of the leak compared to the whole phase space, that is, $p=S_{\mathrm{leak}}/S_{\mathrm{total}}.$ That is, the escape probability (see Eq. (2)) is given by

[TABLE]

Generally, the escape probability is changing as the mass (and size) of the leak is increasing.

At this point, it is useful to introduce some new constants and variables:

[TABLE]

[TABLE]

where $\mathcal{N}$ is the number of particles corresponding to the total mass $\mathcal{M}$, $\kappa_{\infty}$ is the asymptotic escape rate when all the mass of the system is in the leak, $x$ is the ratio of the mass of the leak and the total mass (mass ratio), $y$ is the ratio of the number of the particles which are outside the leak to the total number of the particles $\mathcal{N}.$ It is obvious that

[TABLE]

for all time instant.

We will use these dimensionless quantities through the rest of the paper.

The assumption of a small leak in our model corresponds to the pure exponential survival probability, i.e. when the system shows strong chaotic properties. That is, if a static leak with size equal to the final size of the evolving leak (set by $C_{\mathrm{P}}$) produces exponential decay, we consider that this measure of the leak is small enough to our purposes and fits to the zero order approximation $p=S_{\mathrm{leak}}/S_{\mathrm{total}},$ widely used in the literature, see for example Dettmann and Georgiou (2009). In addition, the exponential decay can also be observed in weakly chaotic systems for short times until the hyperbolic dynamics dominate.

Furthermore, in case of weak chaos the growing leak in the model presented should avoid the quasiperiodic domain in the phase space. On the other hand, if the leak intersects the KAM tori during its growth, the survival probability will decay with lower different rate. In other words, since the regular domain behaves as a forbidden region for trajectories originating outside, the leak biting into it will have an unreachable part for those trajectories resulting in a different escape probability. However, this is no longer true when the leak originally contains islands or more precisely when the ratio of the regular islands inside and outside the leak remains constant.

II.3 Simplified numerical experiment

In order to analyze the escape dynamics through a continuously growing leak defined by Eq. (3), we introduce a simple test system. Our numerical experiments are based on the standard map ($\text{mod}\;2\pi$) which describes the Poincaré map of the kicked rotator.

This choice makes it possible to check the leak’s expansion in both co-ordinate and velocity directions, respectively. The standard map (SM) reads as follows

[TABLE]

In Eq. (4) $K$ denotes the strength of the perturbation and allows to study either fully hyperbolic dynamics ($K$=5.19) or mixed phase space structure, e.g. $K=2.7$.

An other reason we consider the SM is that it allows us to mimic the conservative dynamics in the RTBP where regular islands are also embedded in the chaotic sea producing the well-known structure of the phase space similar to that in Fig. 2.

For simplicity, we presume that the leak grows equally in $I$ and $\Theta$ directions, i.e. it conserves its original shape. In order to avoid the early irregular effects in escape rate due to the location and density of the initial conditions, a threshold time is obtained before the leak is opened. Thus, we have a uniform distribution of the trajectories in the ergodic region of the phase space. The threshold time is set to be $i=250$ in all simulations.

Figure 2 shows the phase space portrait of the SM for $K$=2.7. We place a square-shaped leak centered at point $(I,\Theta)=$(5,5) with initial size $(\Delta I,\Delta\Theta)=(0.01,0.01)$ ($S_{leak}^{(0)}=10^{-4}$) 444 Note that the mass and the size of the leak are identical parameters of the problem. The instantaneous size can be obtained from the current mass and vice versa.

and store the number of escaped trajectories at every iteration step. The semi-diagonals indicate the expansion until the leak reaches its final size $(\Delta I,\Delta\Theta)\approx(0.316,0.316)$ ($S_{leak}^{(\infty)}=0.1$). Initial conditions are placed uniformly in the black square ($3.2\leq\Theta\leq 3.7$, $3.2\leq I\leq 3.7$) far from KAM islands as well as the final leak.

The result of a test run is displayed in Figure 3. It is clearly visible that the well-known exponential decay of the non-escaped trajectories starts after $\mathtt{\sim}$3000 iterations (blue squares). Red triangles denote the instantaneous leak size, $S_{\mathrm{leak}}$, which is growing rapidly until it reaches its final ($\mathtt{\sim}$90%) size. One can also observe that the exponential decay starts roughly when the expansion of the leak ceases. We can, thus, presume that the exponential behavior is a consequence of the stationary leak size with escape rate $\kappa_{\infty}$.

The semi-logarithmic plot of the non-escaped trajectories allows one to find the asymptotic escape rate, $\kappa_{\infty}$ as $t\to\infty$ (for strong chaotic regime). This simulation yields $\kappa_{\infty}=0.00254$.

Furthermore, the numerical investigation confirms the naive idea that until the leak’s expansion is present, the instantaneous escape rate, $\kappa(t),$ and also the escape probability is changing in time according to $\mathrm{d}(\ln y_{n})/\mathrm{d}t=-\kappa(t)$. However, when the growth slows down significantly $\kappa(t)$ reaches the asymptotic escape rate $\kappa_{\infty}$ (green asterisks), see Figure 3. This behavior can be explained as follows. At the beginning of the simulation ($t<2000$) a very large number of escaping trajectories feed the small leak in one iteration step and, therefore, its mass (size) growth is accelerating. Beyond a certain limit the mass (or equivalently the number) of escaping particles in one iteration compared to the mass of the leak becomes small, i.e. escape is present with moderate increase of the leak size. In this case ($2000\leq t\leq 5000$), however, there are enough particles in the system to observe the exponential decay.

The reason for the larger dispersion in $\kappa(t)$ and its deviation from $\kappa_{\infty}$ beyond $t\approx 5000$ is twofold. On the one hand, the number of non-escaped trajectories, after 5000 iterations, becomes so small ($\sim$100) that the statistic is unreliable. On the other hand, Figure 3 shows the simulation for $K$=2.7, in which case KAM tori are responsible for stickiness and consequently a power-law decay of trajectories for longer escape times (not shown). In other words, $\kappa(t)$ would follow the horizontal dashed line in case of the fully hyperbolic dynamics, for instance, $K\geq$ 5.19, with an arbitrarily large $N_{0}.$

III Results

III.1 Analytic solution

After having some impression about the escape dynamics from numerical simulations, in this section, we show that a continuous approximation of the temporal behavior of the model can be described by analytic formulae.

We consider the particle number $N_{i}$ and all the other related discrete functions $M_{i}$, $x_{i}$, and $y_{i}$ as being continuous functions $N(t)$, $M(t)$, $x(t)$, and $y(t)$.

Practically, we can do that because the particle number and the typical timescale (number of iterations) of the process is also much higher than unity ($N_{0}\gg 1$).

The time derivative of $N(t)$ is approximately the negative of the average number of escaping particles $\Delta N$ during one iteration which is $p\cdot N$, so we can write

[TABLE]

where we used Eq. (3). As $\Delta M=-\Delta N\cdot m$, the time derivative of $M(t)$ is

[TABLE]

Combining Eq. (6), $M(t)=x(t)\cdot\mathcal{M}$, and $N(t)=y(t)\cdot\mathcal{N}=(1-x(t))\cdot\mathcal{M}/m$, we get a first-order separable ordinary differential equation for $x(t)$:

[TABLE]

Derivation of the solution can be found in Appendix B.

Equation (7) is a continuous approximation of the recursive difference equation

[TABLE]

where $\Delta x_{i}=\kappa_{\infty}\cdot x_{i}^{\gamma}\cdot\left(1-x_{i}\right)$, which gives the exact description of the discrete-time problem.

The implicit solution of (7) can be given by

[TABLE]

where the constant of integration $\tau$ follows from the initial value $x_{0}$ as

[TABLE]

The solution of Eq. (7), $x(t),$ has a point of inflection (PoI) for all $\gamma>0.$ The second derivative of $x$ from (7)

[TABLE]

from which the $x$ coordinate of the inflection point ($x_{\mathrm{PoI}}$) can be obtained

[TABLE]

We further elaborate on the error properties of the above solution in Appendix C.

We can distinguish two parts of the leak-growing process. The separatrix is the point of inflection of the $x(t)$ function. Figure 4 shows the functions $x(t)$ for different $\gamma$s. For the sake of comparison the graphs are shifted leftward, thus, the inflection points are placed exactly above a row at $t=0.$

The mass growth $x(t)$ beyond the point $x_{\mathrm{PoI}}$ (or $t=0$) has the same characteristic for different $\gamma$s. The reason is that in the limit $t\rightarrow\infty$, $x\rightarrow 1,$ Eq. (7) can be written as $\mathrm{d}x/\mathrm{d}t\approx-\kappa_{\infty}x$ which means that function $x(t)$ approximates 1 exponentially with exponent $-\kappa_{\infty}$ and the process does not depend on $\gamma.$

This is, however, not the case to the left of the point of inflection. In the limit of $x\rightarrow 0$, Eq. (7) can be written as $\mathrm{d}x/\mathrm{d}t\approx\kappa_{\infty}x^{\gamma}$ which means that the solution $x(t)\approx\left[\kappa_{\infty}(1-\gamma)(t+\tau)\right]^{\frac{1}{1-\gamma}}$ follows a power-law and contains both $\kappa_{\infty}$ and $\gamma.$

Furthermore, in this regime $\gamma$ defines two different behaviors. Considering the case of $\gamma<1$ we have a point where $x(-\tau)=0.$ That is, the integration constant $\tau$ is suitable to determine a time instant in the past when the mass of the leak was zero, i.e. when the whole growing process began. While in the case of $\gamma\geq 1$ the function $x(t)$ approaches zero only in the limit of $t\to-\infty.$ In summary

[TABLE]

Nevertheless, it is obvious from Eq. (7) that $\kappa_{\infty}$ is inversely proportional to the timescale of the process. The condition that the timescale have to be much higher than unity is equivalent to $1/\kappa_{\infty}\gg 1$. This fact is important to ensure that the continuous time approximation, Eqs. (5) and (6), is valid in our model.

The adopted model of growing leak defines a stochastic process, whose complete description is possible only by using the probability theory. The question arises naturally, how the probability mass function of the particle number can be calculated after the $i$th iteration if the initial one is known? The question is important because if the standard deviations are considerable, then we need the probability mass functions in order to have a complete description. Otherwise, the averaged behavior, studied previously, describes the process well. In Appendix D we derive the probability mass functions, and study its properties this problem.

We should mention that during the calculation we assumed that $\gamma>0.$ However, it is obvious that solutions of Eq. (7) can also be found for negative exponents in a similar way. The discussion of the case $\gamma<0$ is beyond the scope of the present study.

III.2 Numerical tests

After discussing the analytic description of the survival probability, we confirm the validity of our calculations by running several numerical simulations. In order to demonstrate the general phenomenon of escape dynamics, we use different $\gamma$ values in our calculations.

First, the results of the hyperbolic and mixed dynamics are compared. In this calculation we show that for different system parameters $K=2.7\text{ and }5.19$ the analytical solution works very well. Figure 5(a) shows the ratio of non-escaping trajectories $y(t)$ for the $\gamma=1$ case, i.e. the leak size depends linearly on mass. One can easily see that the analytical solution (dashed and dotted dashed lines) fit the numerical data fairly accurately, especially for small iteration numbers, $t<2000.$ In order to be able to compare the accuracy of the results quantitatively, we calculate the relative difference between the simulated data (S) and the analytic solution (C). The difference $S-C$ in percentages is plotted in Figure 5(b). It shows the same tendency what we can observe by naked eye in panel (a). The $S-C$ diagram remains under 4% level until $t\approx 2500.$ In addition, $S-C$ shows that in the case of $\gamma=1$ the analytic solution is more accurate for fully hyperbolic dynamics ($K=5.19$) than for mixed phase space ($K=2.7$) for $t>2500$. The reason of that comes form Eq. (21), since it turns to be purely exponential for $t\gg 1$, that is $y(t)\sim exp(-\kappa_{\infty}t)$. In addition, the decay of $y(t)$ in the latter case starts to deviate from the exponential due to the sticky effect of the KAM tori.

Physically more interesting cases are when $\gamma\neq 1$ but rational. Let us recall our motivation, the planet formation analogy in the planar RTBP. The size of the leak in the phase space in this particular case is proportional to $m_{p}^{4/3},$ see Eq. (18) in the Appendix.

Figure 6(a) shows the number of surviving particles, $x(t),$ and the mass growth of the leak, $y(t),$ for $\gamma=4/3$ (squares and triangles, respectively). The analytic solution goes together with the numerical simulation also for this value of $\gamma.$ As is well seen in panel (b) the $S-C$ diagram remains under the 5% level until the leak reaches its final mass, $t\approx 8000.$ This is not true, however, at the very beginning of the iteration, $t<10$ after opening the leak. In this regime sudden changes in the number of escaping trajectories appear. Trajectories situated exactly ’above’ the leak and its pre-images disappear immediately from the system. This rapid change in the number of particles is, however, not covered by the analytic solution and, consequently, large differences may show up in the first phase of the $S-C$ diagram.

In the previous two examples we considered particles with equal masses, $m=1.$ A more realistic scenario is when the particles in various physical problems have different masses corresponding to a certain distribution. The log-normal distribution is a good choice to describe the particle size (and/or mass). We present a simulation for $\gamma=2/3$ with different kind of mass distributions, see Figure 7. The numerical results in panel (a) show what can also be derived directly from Equations (5) and (6): the mass growth of the leak does not depend on the mass of the individual particles but only on the mean value of the distribution. Consequently, the leak’s mass changes in time with the same rate for both equal mass particles (pink squares) and log-normal distribution (red triangles), and also for other distributions such as uniform and normal (stars and circles in Figure 7(a), respectively). The statistical fluctuations in leak’s mass, smaller than 15%, disappear after  200 iterations, panel (b).

IV Summary and Discussion

The model Equations (5) and (6) describe the escape dynamics in a leaky chaotic system when the size of the leak is growing in time and the expansion depends on the particles’ mass. Consequently, the escape probability is time-dependent. The analytic solution to the problem provides a power-law behavior at the very early stage ($x\approx 0$) of the dynamical evolution. This phase depends on the exponent $\gamma$ in Eq. (2). However, for larger $t,$ when the feeding of the leak diminishes, the survival decay turns to be exponential. Between these two limits the escape rate is time dependent.

The qualitative picture is the following. After the leak reaches roughly the 90% of its final measure, or more precisely, beyond the point of inflection of $x(t)$, the speed of the growth slows down. After this point the growth of the leak is so slow that it can be thought of as a static leak, and the decay rate turns to be exponential, see Figure 8(a). Numerical simulations verify that the escape rate $\kappa_{\infty}$ (short thick solid line) for a static leak (red triangles) of size 0.1 is the same as in the case of a growing leak (blue squares) when it reaches 90% of its final size (also 0.1), panel (a).

In addition, this behavior is in a very good agreement with the analytical solution describing the early stage escape dynamics. The effect is considerable for relatively short times only as long as enough number of particles are in the system, therefore, the presence of the well-known power-law decay of stickiness (tail of the distribution) in mixed phase space is not affected by the size variation of the leak. However, the crossover time, when the nonhyperbolic part of the chaotic saddle starts to dominate, can be updated.

The crossover time $t_{\mathrm{cross}}$ in weakly chaotic regime is written as follows (Eq. (89) in Altmann et al. (2013a))

[TABLE]

with the assumption that the leak size is small. The growing leak model provides a simple generaliztion of this naive approximation in $\gamma\leq 1$ case

[TABLE]

where the second and the third terms in the bracket define the shift ($t_{\mathrm{shift}}$) the crossover experiences, see the shematic view in Fig. 8b. The second term is the time of the growth until the leak mass is moderate, see the approximation of Eq. (7) when $x\ll 1$, while the third term can be derived from the slope of the function $x(t)$ at point $x_{\mathrm{PoI}}$, Fig. 4. It can be easily obtained that $t_{\mathrm{shift}}\to 0$ when $\gamma\to 0$ and $x\ll 1$.

Equation (25) properly describes also the limit case $m_{i}\to 0.$ Namely, if the mass of the particles tends to zero, i.e. the growth of the leak is fairly slow, one recovers the classical exponential decay for the surviving trajectories. We note that the same effect can be seen when the initial mass of the leak $x_{0}$ is set so large that even the massive particles ($m_{i}>0$) falling into it do not have any effect on the leak’s mass and, therefore, it can be considered as a static leak.

Due to the leak expansion we can consider an instantaneous chaotic saddle in our model at every time step. This object is reducing as the leak is growing and converges to that invariant set which corresponds to the final leak size. This process results in a temporally changing chaotic saddle and a non-stationary exponent of the survival probability (escape rate $\kappa(t)$). A similar phenomenon can be found in Motter et al. (2013) where the exponent is also time dependent (see Eq. (1) in Motter et al. (2013)). In contrast of the similarity, the temporarily changing chaotic saddle should not be confused with the transient chaotic saddle introduced in Motter et al. (2013).

In summary, we have presented an analytic description of the escape of the trajectories through a continuously growing leak both in fully hyperbolic and in mixed phase space.

We stress, however, that during the whole calculation we did not utilize explicitely the fact that $m$ is the mass of the particle, though the basic motivation is related to the mass growth of a planetary embryo. Therefore, one can reformulate the model in a more general way. Let us write Eqs. (1), (2), and (3) together as follows

[TABLE]

where now $m$ is a physical property of the particles, $M_{n}\text{ and }M_{0}$ are the evolved and initial additive property of the leak, and $\mathcal{M}=M_{0}+N_{0}m.$ Other quantities are the same as given in the introduction of the model in Section II.2. This means that the analytical method presented in this paper might be suitable to predict the characteristics of the escape dynamics in different kinds of systems where the leak size depends on some specific physical property of the particles (charge, spin, energy level, chemical composition, etc.).

We also would like to draw the attention to the limitation of present model. In fact, the dynamics in the standard map does not depend on the size of the leak. In other words, the leak affects only the escape rate but not the individual survival trajectories themselves. This is not the case, for instance, in the restricted three body problem, where the growing planetary mass governs the dynamics of the surviving particles and, therefore, should also modify the escape dynamics. Considering such an extension in the SM, a natural choice could be the introduction of a variable nonlinearity parameter $K(M)$ whose value could also depend on the leak mass/size. Studying this effect is postponed to future studies.

Appendix A The exponent $\gamma$ in the planar RTBP

In this section we show a short derivation for how the size of the leak in the RTBP depends on the mass of the planetary embryo. First, we can introduce a four-dimensional leak in the phase space of the RTBP. Two dimensions out of four cover the physical extent of the planet ($0.5r_{\mathrm{H}}$) in the configuration space, i.e. the small gray circle at the position (-1,0) in Figure 1. The remaining two components whose absolute value is the escape velocity at half of $r_{\mathrm{H}}$ describe the size of the leak in the velocity space. In fact, the Hill radius and the escape velocity, as described above, can be written as a function of the planet’s mass. Hence, the size of the 4 dimensional leak in phase space depends only on the mass of the planet $(m_{p}).$

The Hill radius $r_{\mathrm{H}}$ is defined

[TABLE]

where $\mu=m_{\mathrm{p}}/M_{\mathrm{s}}$ is the planet-to-star mass ratio and $a$ is the planet’s semi-major axis. In addition, a particle must have a smaller velocity than the escape velocity in order to be trapped in a pre-defined region, e.g. in one half of the Hill radius. The escape velocity from $0.5r_{\mathrm{H}}$ reads

[TABLE]

where $G$ denotes the gravitational constant and $m_{\mathrm{p}}$ is the planetary embryo’s mass.

Thus, the size of the leak ($S_{\mathrm{leak}}$) in the phase space of the RTBP is obtained as the product of the spatial ($A_{\mathrm{r}}=\pi r_{\mathrm{H}}^{2}$) and velocity extensions ($A_{\mathrm{v}}=\pi v_{\mathrm{esc}}^{2}).$ That is, we have a leak size with $\gamma=4/3$

[TABLE]

Appendix B Solution of Eq. (7)

Let us recall Eq. (7)

[TABLE]

After arrangement and integration we have

[TABLE]

In the special case of $\gamma=1$

[TABLE]

and

[TABLE]

where

[TABLE]

is the constant of integration. In the case of $\gamma\neq 1$, first, we consider the fact that

[TABLE]

Now, the integral on the LHS of Eq. (20) can be written as

[TABLE]

where $q^{(i)}$ is the rising Pochhammer symbol

[TABLE]

and

[TABLE]

is the Gaussian hypergeometric function Abramowitz and Stegun (1970); Zwillinger and Moll (2014). Taking $1/\kappa_{\infty}$ on the LHS and performing the integration on the RHS, the solution as given by Eq. (9) is obtained.

For certain rational $\gamma$ values the implicit solutions of Eq. (7) (corresponding to the integral on the left-hand side of (20)) are summarized in Table 1.

Interestingly, in addition to $\gamma=1,$ the solutions for $\gamma=0$ and $\gamma=1/2$ can also be given in explicit forms as follows

[TABLE]

and

[TABLE]

respectively. Equation (25) provides the classical exponential decay for $\gamma=0,$ when the leak is stationary.

Appendix C Error analysis

During the simulation, the sequence of the averaged mass ratios $(x_{i})^{\infty}_{i=0}$ is governed by the recursive formula (8). The differential equation (7) and its implicit solution (9) give only a continuous approximate solution of the original discrete problem. The question arises naturally, how good the approximation (9) is?

Let us consider two successive terms of the original sequence $x_{i}$ and $x_{i+1}$ (see the inset of Fig. 9). According to the approximation $t(x)$, the time interval between the two states is $t(x_{i+1})-t(x_{i})$ instead of 1. The difference $\Delta t(x_{i})=1-\left[t(x_{i+1})-t(x_{i})\right]$ is the (relative) error of the approximation caused by one iteration. As $x_{i+1}=x_{i}+1/t^{\prime}(x_{i})$, function $\Delta t$ can be expressed as

[TABLE]

Figure 9 shows the functions $\Delta t(x)$ for different $\gamma$s. In the cases of $\gamma=4/3$ and $\gamma=1$ the relative errors remain unter $1.2\%$ (in general under $\kappa_{\infty}/2$).

Unfortunately, in the third case ($\gamma=2/3$), $\lim_{x\to 0}t(x)=1$ (100% relative error). In small $x$ approximation, more precisely if $x_{i}\ll\kappa_{\infty}^{\frac{1}{1-\gamma}}$ , the recursive formula (8) can be approximated by $x_{i+1}\approx\kappa_{\infty}\cdot x_{i}^{\gamma}$. This recursive sequence can be written in explicit form as

[TABLE]

This sequence is increasing really fast from any astronomically small value $x_{0}$ to $x_{i}\approx\kappa_{\infty}^{\frac{1}{1-\gamma}}$. If $x_{0}=10^{-a}\cdot\kappa^{\frac{1}{1-\gamma}}_{\infty}$ and $x_{i}=10^{-b}\cdot\kappa^{\frac{1}{1-\gamma}}_{\infty}$ then the time period of the growing is

[TABLE]

For example, in the case of $x_{0}=10^{-100}\cdot\kappa^{\frac{1}{1-\gamma}}_{\infty}$ ($a=100$) and $x_{i}=0.98\cdot\kappa^{\frac{1}{1-\gamma}}_{\infty}$ ($b\approx 0.01$), $i\approx 23$. The continuous approximation does not describe this fast growing process.

Although the relative error decreases under $3.7\%$ at $x\approx 0.01$ (see Fig. 9), according to our numerical results the global error is acceptable if the initial mass ratio $x\gtrapprox\kappa_{\infty}^{\frac{1}{1-\gamma}}$. For example, in the corresponding case of Fig. 11 in spite of the initial mass ratio ($x_{0}=10^{-6}$) is slightly smaller than $\kappa_{\infty}^{\frac{1}{1-\gamma}}$, the global error remains under $4\%$.

Appendix D Distributions

Let $\xi_{i}$ be discrete random variables associated with the number of particles after the $i$th iteration. Here we derive the probability mass function $P(\xi_{i+1}=k),\;k=0\dots N_{0}$ by assuming that it is known from the earlier iterations $P(\xi_{i}=j),\;j=0\dots N_{0}.$

The number of escaping particles during one iteration follows binomial distribution. Let us suppose that there are $j$ particles in the system ($\xi_{i}=j$) and after one iteration the number of particles is $k\leq j$ ($\xi_{i+1}=k$), then the number of escaping particles is $j-k.$ Using the formula of the binomial distribution, we can write the following conditional probability

[TABLE]

where $p_{j}$ is the escape probability which corresponds to the particle number $N=j,$ namely

[TABLE]

According to the law of total probability, we can write

[TABLE]

thus we get a recursive formula for $P(\xi_{i+1}=k).$ If the initial number of particles is set to be $N_{0}$ then the initial distribution reads

[TABLE]

and any $P(\xi_{i+1}=k)$ probability can be calculated recursively by using (30)–(35).

In order to check whether the analytic model is valid, several calculations of distribution series were carried out. Figure 10 shows the 1st, 5th, and 9th deciles (10-quantiles) of the series of distributions for $\gamma=4/3.$ This calculation is suitable to test the accuracy of the particle number ratio $y(t)=1-x(t)$ calculated in the section III.1. The analytic solution is also plotted (black curve) together with the statistical results. One can see that the function $y(t)$ is close the decile curves which means that the analytic solution is suitable to approximate the discrete process.

We also verified these result by analyzing distributions for different $\gamma$s. Figure 11 illustrates the results for $\gamma=2/3,\;1$, and $4/3$. The other parameters were $N_{0}=10^{6}$, $m=1$, and $C_{p}=(2\pi)^{-2}\cdot N_{0}^{-\gamma}$ in all three cases (the initial escape probabilities were the same, $p_{0}=10^{-4}\cdot(2\pi)^{-2}$ and $\kappa_{\infty}\approx 1/(2\pi)^{2}=2.53\cdot 10^{-2}$). The distributions were calculated until their averages decreased under 0.1 percent of the initial particle number ($E(\xi_{i})<10^{-3}\cdot N_{0}$).

Figure 11(b) shows the standard deviations in all three cases. In general, the standard deviations are not negligible but remain relatively small.

Acknowledgements.

We are indebted to G. Kovács and T. Tél for useful discussions. The authors also thank the anonymous referees their valuable comments and suggestions that helped to improve the text significantly. This work was partially supported by the OTKA Grant No. NK100296, K119993, and PD121223. TK also thanks for the support for the Fulbright Commition and the Hungary Initiatives Foundation.