Nonuniversality and the role of tails in reaction-subdiffusion fronts

Recently there has been a certain controversy about the scaling properties of reaction-subdiffusion fronts. Some works seem to suggest that these fronts should move with constant speed, as do classical reactiondiffusion fronts, while other authors have predicted propagation failure, i.e., that the front speed tends asymptotically to zero. In the present work we confirm by Monte Carlo experiments that the two situations can actually occur depending on the way the reaction process is implemented. Also, we present a general analytical model that includes these two different behaviors as particular cases. From our analysis, we reach two main conclusions. First, the differences found in the scaling properties show the lack of universality of reactionsubdiffusion fronts. Second, we prove that, contrary to the widespread belief, the tail of the waiting time distributions is not always decisive to determine the speed of these fronts, but sometimes it plays just a marginal role in the front dynamics.


I. INTRODUCTION
The field of anomalous diffusion has attracted a great interest in the last years, as many systems out of equilibrium in nature have been proved to exhibit power-law decays in the distribution of their characteristic waiting times.This implies a fundamental deviation from canonical ͑equilibrium͒ statistics ͓1͔.Of particular interest is the case when the particles in these systems experience a reaction process too, so one can try to analyze the influence of anomalous transport on usual reaction-diffusion phenomena as wave propagation ͓2-5͔, Turing patterns ͓6-8͔, segregation ͓9͔, etc. Mancinelli et al. ͓10͔, for example, provided an interesting compilation of formal situations where anomalous transport is encountered and then explored the traveling front solutions obtained by introducing reaction as an additive term.Fractional diffusion equations have also been used to explore the behavior of traveling fronts for Lévy flights ͓11,12͔ finding that these fronts exhibit in some cases exponential acceleration although recent works have shown the limitations of that result ͓13͔.Finally, several attempts have been done to characterize analytically reaction-subdiffusion fronts.These attempts have been directed to determine either the speed of these fronts ͓2,3,5,14͔ or their asymptotic scaling under different initial conditions and different chemical kinetics ͓15,16͔.
Recently, the appropriate way to implement mathematically a reaction process for a system under subdiffusion has been the subject of an important discussion ͑see ͓14,17,18͔ and references there in͒.Sokolov et al. ͓17͔ first detected that if local conservation of particles is imposed in a process A → B and the ͑anomalous͒ transport is assumed to be completely independent of reaction, then the evolution equation for species B contains a term where the Laplacian of species A appears.This means that nontrivial coupled effects arise due to subdiffusion.According to that, reaction-subdiffusion processes cannot be simply described by the independent contribution of transport and reaction terms, but this coupling must be considered.This idea has important consequences, for example, on the form of stationary solutions for reaction-diffusion in finite domains ͓19͔.
In the present paper we focus on the effects that couplings between anomalous transport and reaction have on the dynamics of traveling fronts for the standard reaction A + B → 2A.Some recent works have predicted for this case a propagation failure ͓14͔, i.e., the velocity of the fronts will tend to zero for t → ϱ.This, however, is in contrast with the results found in ͓2,5͔.There, the Hamilton-Jacobi method ͓20͔ has been used to obtain an explicit expression for the speed of reaction-subdiffusion fronts, which seems to indicate the existence of a well-behaved front with constant speed.In order to explain this contradiction, in ͓14͔ it has been argued that the details of the chemical kinetics considered are different in each case.More specifically, these differences are in the way that the waiting time between jumps gets modified as a consequence of the reaction process.So that we can differentiate the following: Case I ͑Refs.͓2,5͔͒: for the random-walk process considered in these references, the internal clock of a particle ͑which measures the time elapsed since the last jump of that particle͒ is set to 0 after it reacts.According to that, one could think that the average time between consecutive jumps gets increased, but in fact the contrary happens ͓18͔.When we set the internal clock to zero, we take the particle back to the body of the waiting time distribution ͑see Fig. 1͒.It means that we are giving the particle a higher probability to jump.In consequence, jumps occur more frequently and the traveling front can be sustained albeit the waiting time distribution exhibits a heavy tail.
Case II ͑Refs.͓14,17͔͒: in these works it is being assumed that the reaction events that a particle can experience do not modify the internal clock of the particle, i.e., the random time between jumps will be exactly the same if the particle has reacted meanwhile or not.Then the transport properties of the system in the regime t → ϱ will be completely governed by the tail of the distribution.
Despite all this discussion and all the efforts made to understand reaction-subdiffusion fronts, we stress that the results reported in ͓2,5͔ for case I have not been verified numerically yet, and the results for case II have been numerically tested only in a very recent work ͓21͔ which was published at the moment that the current paper was in elaboration.So our main objective here is to propose some numerical ͑Monte Carlo͒ experiments able to confirm that ͑i͒ in some situations reaction-subdiffusion fronts do exhibit propagation failure while in other cases constant propagation holds, and ͑ii͒ the Hamilton-Jacobi method can fit correctly the values for the propagation speed.The recent work in Ref.
͓21͔ will be also discussed here, as we have detected some important discrepancies between the numerical results reported there and ours.On the other side, we will propose a general analytical model which includes the two cases I and II as particular cases.Our analysis leads us to some important conclusions regarding the role of the tails in these processes.We shall prove that it is mainly the body ͑not the tail͒ of the distribution what determines the front speed in case I.As a consequence, the expressions obtained in ͓2,5͔ for the speed of reaction-subdiffusion fronts are of limited utility since only the effect of tails was considered there.

II. RESULTS FROM THE NUMERICAL EXPERIMENTS
Our Monte Carlo experiments have been performed by following the prescriptions given in ͓15͔.We set a fixed number N of particles ͓N A ͑x͒ particles of species A and N B ͑x͒ of species B͔ at every xth node of a one-dimensional ͑1D͒ lattice of size L. In order to generate symmetric traveling fronts we choose all the particles in the lattice to be initially of species B ͓it is, N A ͑x͒ = 0 for any x͔ and replace at t = 0 the B particles by A particles at the central node of the lattice, denoted by x = 0. Hence, the front position is defined as the furthest node ͑from the center͒ where the number of particles of species A is above a certain threshold N A ‫ء‬ .Every particle is updated at periodic intervals separated by a time step ⌬t.At each time step every particle of species A is given the opportunity to react ͑become a particle of species B͒ with probability kN A ͑x͒N B ͑x͒⌬t, where k is the reaction constant.Also, if the particle has been waiting at the same node for a time larger or equal than its random time to the next jump, then it jumps with probability 1/2 to the right and 1/2 to the left ͑we only allow jumps to the nearest neighbors for simplicity͒.The random waiting time to the next jump for each particle is chosen from a series of values distributed according to ͓15,22͔ with 0 Ͻ ␣ Ͻ 1, which is a treatable and well-studied expression exhibiting an asymptotic power-law decay ͑͒ ϳ ␣␤ −1 −1−␣ .
In order to compare with continuous models, which is one of our main purposes, we must choose N sufficiently large and ⌬t small.We have taken N =10 4 and ⌬t =2ϫ 10 −3 ͑with arbitrary units͒ in most of the results presented here and have checked that more extreme values of N and ⌬t do not lead to a significant improvement of our results.For small values of N discretization effects would arise; these effects deserve an exhaustive analysis and so they will be addressed by separate in a forthcoming paper.Note also that according to our transport algorithm, the number of particles at each node of the lattice is not constant in time, but it will fluctuate around its mean value N.This could be avoided by introducing a more complex rule for transport.However, for N large the effect of these fluctuations can be neglected, and actually we have verified that they do not affect the results presented in this paper.Also, in order to prevent undesirable boundary effects we always use a lattice with periodic boundary conditions.
Cases I and II can be implemented in the experiments by renewing the internal clock of particles in the same fashion as described above.That is, a new waiting time is chosen for a particle only in case it jumps for case II ͑so its internal clock is set to 0͒, while in case I the waiting time of the particle is renewed also after it reacts.In Fig. 2 we present the results for both cases and for the same values of the parameters k, ␣, and ␤.For the case I there is always a linear dependence between front position and time ͓Fig.2͑a͔͒ so a constant-speed stationary solution is found.On the contrary the case II exhibits, after a short transient, a convex growth which confirms the tendency v → 0 found in ͓14,21͔.To stress the differences between both cases, in Fig. 2͑b͒ we plot the number of jumps performed in the whole lattice as a function of time.For case II a scaling t ␣−1 is found which can be easily justified ͑see Sec.V below͒, while for case I this scaling stops holding for large times and tends to a constant value.This suggests that in case II, as time goes by more particles get "trapped" in the tail of the waiting time distribution and so they stay waiting for very long times, which results in a continuously decreasing number of jumps in the system and a failure in the front propagation.On the contrary, in case I the reset of the waiting times after a reaction event ensures that there are always some particles which "escape" from the tail of the distribution and so the number of jumps do not decrease to 0 for t → ϱ.

III. MODEL
In order to understand in deeper detail the situation, let us propose a model which is able to capture the essential dynamics of these processes.First, for the reaction kinetics we consider the catalytic conversion A + B → 2A of B particles into A, which leads to the form where R A ͑x , t͒ represents the number of new particles of species A appeared from the reaction.
In order to implement transport, we will follow the general ideas from the continuous-time random walk ͓23,24͔.So the spatiotemporal evolution of the number of particles of species A can be written as where J A ͑x , t͒ represents the number of particles of species A jumping to the xth node at time t.The first term on the right-hand side ͑rhs͒ of Eq. ͑3͒ stands for the contribution from those A particles initially at t = 0 that have not performed their first jump yet.We define ͑t͒ and R ͑t͒ as the probabilities to remain at the current node at least for a given time for those particles appeared as a consequence of a jump ͑͒ or as a consequence of a reaction event ͑ R ͒.For the former case, we have so there is a direct relation with the waiting time distribution ͓Eq.͑1͔͒.On the other hand, the specific form of R ͑͒ can be independent of ͑͒ but depends on the way we assign waiting times to the particles after they have reacted.This will allow us to differentiate between cases I and II.Similarly to Eq. ͑4͒, the relation must hold, where R ͑͒ is defined as the waiting time distribution for particles A appeared from reaction kinetics.
According to the model, the quantity J A ͑x , t͒ + R A ͑x , t͒ determine the total income of new particles A to the xth node at time t.The contribution from jumps will be expressed by The first term in the rhs of this equation represents the firstjump contribution of those particles that have been waiting at their initial position up to time t.The second term follows from the individuals jumping to the xth node from its neighbors in the lattice.Finally, the third term represents those particles of species A, having appeared at x −1 or x + 1 as a consequence of a reaction event, that jump to x after a random time distributed according to R ͑t͒.Now, we can try to solve system ͑2͒-͑6͒ by transforming from the real space ͑x , t͒ to the Fourier-Laplace space with coordinates ͑q , s͒.By doing this, one obtains the explicit expression for the number of particles, where we use the variables s and q to specify that the functions have been transformed to the Fourier-Laplace space.The corresponding generalized master equation ͑GME͒ from Eq. ͑7͒ is found by inverting back to the real space

‫ץ‬N
where the functions M i ͑t͒ ͑i =1,2,3͒ are defined through their Laplace transforms Note that M 1 corresponds to the so-called memory kernel, and so the GME of the standard CTRW is recovered from Eq. ͑8͒ if the reaction process is obviated ͑it is, for R A =0͒. Also, the explicit dependence of M 2 and M 3 on the transport functions and R is a clear trace of the couplings between reaction and diffusion mentioned above.So, in the present model we find that the coupling is made evident within the reaction term, while in the aforementioned approaches ͓17͔ it is found that the master equation for species A depends explicitly on the transport properties of species B. At the end, the idea is that memory effects due to anomalous transport always make transport and reaction terms nonseparable within the master equation.

IV. APPLICATION TO CASE I
Now we are in position to analyze cases I and II as particular cases of the model above.Case II has been defined by the condition that the internal clock of particles is set to 0 after the particles react.It leads to the identity R ͑t͒ = ͑t͒, i.e., the random waiting times are distributed identically to every newcome, no matter if it has appeared from a jump or a reaction event.It is straightforward to see that for this case the GME ͑8͒ turns into which corresponds to the reaction-random-walk case already studied in many previous works ͑see, for example, ͓2,20͔͒.For this specific case, it has been shown that wave-front solutions can be analytically characterized by means of the Hamilton-Jacobi method ͓2,5͔.So that Monte Carlo experiments are expected to yield a stationary front solution traveling with constant speed as confirmed by our results in Fig. 2. If one applies the Hamilton-Jacobi method ͑see ͓2,20͔͒ to determine the front speed from Eq. ͑10͒ together with Eq. ͑2͒, the expression found is After inserting Eq. ͑1͒ into Eq.͑11͒ we obtain v = min where ⌫͓•,•͔ is the incomplete gamma function.So that we can compare now the value of v from this expression ͑the minimum must be computed numerically͒ with that found from Monte Carlo experiments.In Fig. 3 we show this comparison-the circles represent the numerical results and the solid lines are obtained from Eq. ͑12͒-which yields an excellent agreement.Likewise, note that the results from Eq. ͑12͒ are in disagreement with those found in Refs.͓2,5͔ ͑represented by dashed lines in Fig. 3͒.There the front speed was calculated by using the approximation ͑s͒ϳ1−Cs ␣ ͑with C constant͒ corresponding to an asymptotic decay ͑͒ ϳ −1−␣ for large .This disagreement seems to suggest that the tail of the waiting time distribution is not necessarily the responsible for the dynamics of the front.Actually, the contrary is closer to truth in this case.The dynamics of the front is governed by those particles that jump sooner, it is, those which are in the body of the distribution ͑Fig.1͒.As the new A particles appeared from the reaction kinetics are assigned a new random waiting time, these are the ones that are expected to determine the front speed.These ideas have been confirmed by our Monte Carlo experiments.We have found that very similar front speeds are found for waiting time distributions ͑͒ with the same body but very different tails, while different speeds arise if the body is different albeit the tails are exactly the same.This can also be seen in Fig. 3, where we have chosen the parameter C from ͑s͒ϳ1−Cs ␣ in such a way that the tail of this distribution coincides with that of Eq. ͑1͒.Despite the tail of the distribution is so exactly the same, different values for v arise.Obviously, the differences become notorious for ␤ ӷ 1, as for that regime the distribution ͓Eq.͑1͔͒ has a larger body.On the contrary, for ␤ Ӷ 1 the distribution tends to its asymptotic form ͑͒ ϳ ␣␤ −1 −1−␣ very fast and so the front speeds become very similar.
These ideas serve to solve the paradox reported in ͓2͔ and later discussed in ͓5͔.The authors there found that reactionsubdiffusion fronts could be faster than classical reactiondiffusion fronts.They suggested that this result was not possible since subdiffusive particles move asymptotically slower than diffusive ones and so should do the corresponding front solutions.Our present analysis show that in fact the paradox does not exist.The speed of reaction-subdiffusion fronts arising from Eq. ͑10͒ must not necessarily be lower than that of classical reaction-diffusion fronts, since the tail of the distribution does not govern the front dynamics.

V. APPLICATION TO CASE II
Case II is a more complex situation, as we have to differentiate the waiting time distributions for newcomes appeared from a jump or from a reaction event.If the particles of species A and B have exactly the same transport properties ͑as assumed in ͓14,17͔͒, then R ͑͒ can be interpreted in terms of the distribution of waiting times of arbitrary age ͑DWTAA͒ defined in ͓25͔.This refers to the stationary versus non stationary state condition for the CTRW reported in ͓24,26,27͔.Consider that the CTRW process started at time t =−t a but we start our observation ͑it is, we start counting the waiting times͒ at t = 0; the corresponding distribution of random times to the first jump is called DWTAA and denoted by t a ͑͒.This formal definition coincides with the situation considered in our case II, except that here the CTRW started at t = 0 and we want to determine the distribution of waiting times to the first jump after a particle has reacted at time t.So that we can write the distribution R ͑͒ as ͓see Eq. ͑5͒ in ͓25͔͔

͑13͒
Here, ͑n ͉ ͒ is the probability that n jumps occur during the interval ͑0,͒ provided that the last jump is performed exactly at time ; so we can identify ͑1 ͉ ͒ = ͑͒.This n-jump function can be determined from the convolution ͑n͉͒ = ͵ 0 ͑n − 1͉tЈ͒͑ − tЈ͒dtЈ.͑14͒ Equation ͑13͒ can then be understood as follows.The first term in the rhs represents the contribution from particles that have not performed their first jump when observation starts ͑i.e., at time t͒, while the second term gives us the contribution from particles that have already performed 1,2,3… jumps before.For a more detailed description of the DWTAA and its connection to the CTRW the reader is addressed to Refs.͓25,28͔.
Despite we are able to propose an exact expression for R ͑͒, note that the resolution of the problem comes out to be very difficult if not impossible.First of all, it is not possible to find an exact analytical expression of R ͑͒ from Eq. ͑13͒ with Eq. ͑1͒.What is more, we have that R ͑͒ depends explicitly not just on the waiting time but also on the absolute time t.To stress this dependence we use in the following the notation R ͑ , t͒ instead of R ͑͒.According to that, the Fourier-Laplace transform of model ͑2͒-͑6͒ cannot be analytically computed, and so a GME cannot be explicitly found for this case.Anyway, we can still try to find some approximated expressions in the asymptotic regime t → ϱ.For this limit, Eq. ͑13͒ together with Eq. ͑1͒ leads to the expression It means that the jump probability of a particle after a reaction event tends to 0 for t → ϱ.So that we can assume that in the asymptotic regime the contribution of newcomes arrived from jumps is much larger than the contribution from particles arrived from reaction events.Accordingly, we assume R Ϸ 0 in expression ͑8͒.From that we can try to apply again the Hamilton-Jacobi method to find the asymptotic value for the front speed; this leads to the expression Since ͑s͒ must be always a monotonically decreasing function of s, it follows that the function s / cosh −1 ͑1 / ͑s͒͒ grows monotonically with s and so the minimum in expression ͑16͒ will be always at s = 0.In consequence, the asymptotic front speed predicted is v = 0, in accordance with the result found in ͓14͔ by alternative methods.The interpretation of this propagation failure is simple.In the asymptotic regime all the particles are trapped in the tail of the distribution and they just wait without performing any jumps, so the front speed eventually tends to 0. The scaling behavior t ␣−1 found in Fig. 2͑b͒ for the number of jumps versus time shows this tendency too and can be justified as follows.According to the notation from our model in Sec.III, the number of jumps performed is equivalent to ͚ x ͓J A ͑x , t͒ + J B ͑x , t͔͒, where the sum is performed over all the nodes of the lattice.The evolution of the total number of particles ͑A plus B͒ is not affected by the reaction so we have from the standard CTRW that which, after inserting Eq. ͑1͒, shows an asymptotic scaling J A ͑x , t͒ + J B ͑x , t͒ϳt ␣−1 , in agreement with the numerical results in Fig. 2͑b͒.Finally, from the results shown in Fig. 2͑a͒ one may conclude that a scaling v ϳ t −1 holds too, at least for a certain region of times.However, according to our Monte Carlo experiments, this scaling is quite poor for a wide range of parameter values, specially for ␣ small.Alternatively, we have tried to determine from our model in Sec.III and by standard scaling techniques ͓15,16͔ whether any simple scaling could be analytically deduced, without success.

Comparison with Ref. [21]
As mentioned above, while the present manuscript was in elaboration a new article ͓21͔ has been published where the results predicted in ͓14͔ for the case II have been tested numerically.Except for some minor details, which are not expected to affect the results, the numerical experiments proposed in that work are equivalent to ours.However, the range of parameters explored by the authors there is quite different to that used here.Specifically, we have chosen N large and ⌬t small in order to make sure that a comparison with models based on the CTRW is feasible, while in ͓21͔ they used ⌬t = 1 and quite small values of N. Two main conclusions reached in that work were: ͑i͒ propagation failure for t → ϱ occurs as a consequence of a decelerating frontwhich is the same result we find here-and ͑ii͒ a scaling v ϳ t −1 holds, with = ͑␣ +1͒ / 2. Since this last idea contradicts our results above, we have decided to use our Monte Carlo experiments to reproduce the results presented in Fig. 3 of Ref. ͓21͔.Note that the scaling v ϳ t −1 is found there for less than two decades, typically for times from 10 3 to 2 ϫ 10 4 .This is not sufficient in general to ensure that a power law holds, so we have extended these experiments up to t =10 6 .The corresponding results are shown in Fig. 4, where we plot N AT ͑t͒ as a function of time.Here, N AT ͑t͒ represents the total number of A particles in the lattice N AT ͑t͒ = ͚ x=−L/2 L/2 N A ͑x , t͒, which will be proportional to the front position, and so it should fulfill N AT ϳ t .However, from our results it is clear that these data do not fit asymptotically a single power law, as reported in Table I.For example, for the case ␣ = 0.6 we obtain an exponent = 0.785 in the time interval 10 2 Ͻ t Ͻ 2 ϫ 10 4 , very similar to the value found in ͓21͔, but = 0.727 is found for 10 3 Ͻ t Ͻ 10 5 and = 0.706 for 10 4 Ͻ t Ͻ 10 6 .Only for larger values of ␣ the fit seems to be better, but it could be due to the fact that t =10 6 is still too small to observe the real asymptotics for that case.So, our results suggest that the scaling v ϳ t ͑1−␣͒/2 reported numerically in Ref. ͓21͔ is not robust.

VI. CONCLUSIONS
We have presented a reaction-random-walk model in order to compare and confirm some results obtained by different authors previously about the dynamics of reaction-subdiffusion fronts.Our work provides numerical evidence that the dynamics of these fronts is very different depending on the way that reaction is implemented ͑case I or II͒.This shows that the scaling properties of reaction-subdiffusion fronts are not universal, contrary to classical ͑Markovian͒ reaction-diffusion fronts.In the classical case, no memory effects are present, which always results in the same scaling of wave-front solutions.
At the same time, our numerical experiments yield some interesting results about the role of the tails of the waiting time distributions considered in reaction-subdiffusion systems.In our case II, studied previously by Sokolov et al. ͓17,14,21͔, the long tail of the distribution ͑1͒ is the responsible for the propagation failure, as more particles get gradually trapped in the tail.This results in a frozen state in which all particles remain waiting and so front propagation is not possible.On the other side, for the case I we have found that the Hamilton-Jacobi method is able to predict the constant speed found numerically, provided one takes into account the whole distribution of waiting times and not just the tail of the distribution.This is because, contrary to the widespread belief, the tail of the distribution is not necessarily the responsible for the front dynamics.
Let us conclude by mentioning that in order to get a realistic implementation of reaction-subdiffusion processes for the case of chemical reactions we do not think neither case I nor case II are completely satisfactory.The model presented here in Sec.III represents just an attempt to unify and  generalize these approaches.At practice, we think that the reaction kinetics may probably affect somehow the waiting times to the next jump of particles, but it is not expected in general to make the internal clock of the particle be set to 0 ͑as claimed in case II͒, except maybe for some catalytic isomerization reactions and other simple cases.At the end, all these ideas are based on simplified analytical models, while greater efforts at the level of single-molecule experiments are needed before we can elucidate the real statistical properties of individual reaction events.

FIG. 1 .
FIG. 1. Schematic representation of a general probability distribution function.

FIG. 2 .
FIG. 2. Results obtained from the Monte Carlo experiments for the cases I and II in the text.Two different values of ␣ are shown, 0.3 ͑circles and triangles for cases I and II, respectively͒ and 0.7 ͑squares and inverted triangles for cases I and II, respectively͒.The rest of parameter are set as fixed, with ␤ = 0.01 and k = 0.2.͑a͒ The front position is computed as the furthest node from the center where N A Ͼ N ‫ء‬ = 0.1N, with N =10 4 .The lines represent a scaling proportional to t, which confirms the constant front speed found only for case I. ͑b͒ The total number of jumps includes jumps from all particles ͑A and B͒ in the whole lattice.The lines represent the scaling t ␣−1 .

FIG. 3 .
FIG. 3. Comparison between the numerical values of the front speed ͑Monte Carlo experiments, points͒ and the Hamilton-Jacobi results derived here ͓Eq.͑12͒, solid lines͔ and in Ref. ͓2͔ ͑dashed lines͒.The value of the parameter C is chosen so that the distribution ͓Eq.͑1͔͒ coincides in the Laplace space with ˆϳ 1−Cs ␣ in the asymptotic regime.Two different values of k are shown, while ␣ = 0.5 is set as fixed.

FIG. 4 .
FIG. 4. Plot of the total number of A particles N AT versus time for three different values of ␣ ͑see legend͒.To reproduce the cases studied in ͓21͔, the other parameters have values ␤ =1, k = 0.006, N = 10, and ⌬t =1.

TABLE I .
Values of the exponent obtained from a power-law fit N AT ϳ t from the data in Fig.4for different time intervals.␣10 3Ͻ t Ͻ 2 ϫ 10 4 10 4 Ͻ t Ͻ 2 ϫ 10 5 5 ϫ 10 4 Ͻ t Ͻ 10 6