Dependence of the transportation time on the sequence in which particles with different hopping probabilities enter a lattice
aa r X i v : . [ n li n . C G ] S e p Dependence of the transportation time on the sequence in which particles withdifferent hopping probabilities enter a lattice
Hiroki Yamamoto , ∗ Daichi Yanagisawa , , and Katsuhiro Nishinari , School of Medicine, Hirosaki University, 5 Zaifu-cho Hirosaki city, Aomori, 036-8562, Japan Research Center for Advanced Science and Technology, The University of Tokyo,4-6-1 Komaba, Meguro-ku, Tokyo 153-8904, Japan Department of Aeronautics and Astronautics, School of Engineering, The University of Tokyo,7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan (Dated: September 19, 2019)Smooth transportation has drawn the attention of many researchers and practitioners in severalfields. In the present study, we propose a modified model of a totally asymmetric simple exclusionprocess (TASEP), which includes multiple species of particles and takes into account the sequence inwhich the particles enter a lattice. We investigate the dependence of the transportation time on this‘entering sequence’ and show that for a given collection of particles group sequence in some casesminimizes the transportation time better than a random sequence. We also introduce the ‘sortingcost’ necessary to transform a random sequence into a group sequence and show that when this isincluded a random sequence can become advantageous in some conditions. We obtain these resultsnot only from numerical simulations but also by theoretical analyses that generalize the simulationresults for some special cases.
I. INTRODUCTION
Transportation systems are key topics in social or bi-ological systems [1]. In social systems, researchers havesought to obtain smooth transportation in various situa-tions, such as production flow [2, 3], vehicular traffic [4–7], and pedestrian evacuation [8–11]. On the other hand,for biological systems, intracellular transportation alongmicrotubules has been vigorously investigated [12–14].Among various transportation models, the asymmetricsimple exclusion process (ASEP), pioneered by MacDon-ald and Gibbs [15, 16], has attracted much attention. It isa stochastic process on a one-dimensional lattice in whichparticles move asymmetrically. A derivative of ASEP, inwhich particles are allowed to hop unidirectionally (leftto right in the present study) is called a totally asym-metric simple exclusion process (TASEP). In the fieldof nonequilibrium statistical mechanics, researchers haveapplied TASEP to various transportation problems, suchas molecular-motor traffic [17–20], vehicular traffic [5–7, 21–23], and the exclusive-queuing process [24–26], es-pecially since the TASEP with open boundary conditionshas been solved exactly [27–29].In practice, researchers struggle to achieve smooth op-eration for various tasks, smooth logistics for variousproducts, and an effective evacuation method for pedes-trians in various situations, such as exit plans from sportsstadiums and concert venues. To attain smooth flow insuch situations, we often consider the sequence in whichwe perform tasks and pedestrians move because this se-quence may affect the total performance of the systems.For example, slow pedestrians may block fast ones at theback of a narrow street, which worsens pedestrian flow. ∗ [email protected]; To investigate how the abovementioned sequences affectpedestrian flow in various systems, herein, we proposea modified TASEP comprising a finite number of multi-species particles, in which the entering sequences of theparticles are considered.In the proposed model, we consider the number of par-ticles to be finite and study the transportation times ofthose particles. Note that we do not consider the steadystate of the system itself. Minemura et.al [30] investi-gated the transportation time for a hopping probabilitythat depends upon the lot size, using the single-speciesTASEP with a finite number of particles. Other relatedworks [31–34] also adopted a finite number of particles.In those models, however, particles circulated through asystem comprising a lattice and a particle pool while theinput or output rate was varied.Additionally, the concept of multiple particles has al-ready been extensively studied [35–58]. For example,second class particles were introduced in Refs. [35–45]and more than two-species particles were introduced inRefs. [46–55]. However, most of these studies focusedon mathematically exact solutions to the systems underconsideration by using such as Matrix Product Ansatzand did not much consider the application of the studiedmodel to real-world situations. Furthermore, owing totheir simplicity, periodic-boundary conditions have beenadopted in many studies [35, 39, 47, 49–51, 53, 54, 56–58]. Studies on multi-species ASEP with open bound-aries and random updating were undertaken only re-cently [37, 40, 41, 44, 45, 52, 55]; for example, Ref. [55]obtained the exact phase diagram for a multi-species(more than two-species) ASEP. The present investiga-tion primarily focuses on the problem of minimizing thetransportation time, adopting open-boundary conditionsand parallel updating. With the same boundary condi-tions and updating rules as the present study, Ref. [59]adopted particles with disorder, whereas jumping par-ticles were introduced in Ref. [60]. Note that majorityrelated works considering multi-species particles assumethat swapping between different types of particles, i.e.,bidirectional particle hopping, can occur, whereas ourmodel prohibits swapping [61].In contrast, to the best of our knowledge, no TASEPinvestigations that focus on the entering sequence of theparticles (the key highlight of our model) have been re-ported thus far. Herein, we have considered two specialtypes of sequences in particular: ‘random sequences’ and‘group sequences,’ and we have compared the transporta-tion times for these two types of sequences. In associationwith the entering sequence, we have introduced the sort-ing cost in our model. Without sorting, particles are typ-ically transported at random, i.e., in a random sequence.Therefore, considering the cost of sorting particles from arandom sequence into a group sequence is useful. In thepresent study, we define this sorting cost and comparethe results obtained with and without sorting.We have determined the dependence of the transporta-tion time on the entering sequence of the particles fromnumerical simulations based on our model. Moreover,we find that the optimal sequence can vary, dependingupon choice of parameter set, when the sorting cost isconsidered. In addition, we have succeeded in obtainingmathematical proofs of the simulation results for somespecial cases.The remainder of the present study is organized as fol-lows. Section II describes the details of our proposedmodel and some important parameters, modifying theoriginal TASEP. In Sec. III, we present and discussthe results of numerical simulations using the modifiedTASEP. Section IV presents theoretical analyses of thesimulation results for some special cases. The paper con-cludes in Sec. V.
II. MODEL DESCRIPTIONA. Original (single-species) TASEP withopen-boundary conditions
The original TASEP with open-boundary conditionsis defined as a one-dimensional lattice of L sites, labeledfrom left to right i = 0 , , ......, L − α , andleave the lattice from the right boundary with probabil-ity β . In the bulk of the lattice, if the right-neighboringsite is empty, a particle hops to that site with probability p ; otherwise it remains at its present site. Our modified TASEP differs from this original one in the following fourways. (cid:1) p p (cid:1) (cid:2) L sites (cid:1) (cid:1) i = 0 i = 1 i = L- FIG. 1. (Color online) Schematic illustration of the originalTASEP with open-boundary conditions.
B. Difference 1: Finite number of particles
First, the number of particles N is finite, as illustratedin Fig. 2. The system evolves until the N th particleleaves the lattice. We define the transportation time T as the time gap between the start of the simulation andthe time when the N th particle leaves the lattice. (cid:1) (cid:1) (cid:2) L sites N particles FIG. 2. (Color online) Schematic illustration of TASEP witha finite number of particles. This figure shows the case N = 12. C. Difference 2: Multi-species particles
Second, our model adopts multi-species particles, i.e.,particles with different hopping probabilities, as illus-trated in Fig. 3. Specifically, each of the N particlesis allocated to one of S species, where 1 ≤ S ≤ N . Par-ticles that belong to each species s ( s = 1 , , ......, S ) allhave the same hopping probability p = p s (0 < p s ≤ S = 1 our model reduces to the single-species TASEP, whereas with S = N all particles havedifferent hopping probabilities. The fraction of all the N particles allocated to each species s is defined as r s ,obviously satisfying P Ss =1 r s = 1. (cid:1) p (cid:1) L sites (cid:1) p (cid:1) p (cid:2) p (cid:3) Species Species Species Species FIG. 3. (Color online) Schematic illustration of TASEP withmulti-species particles. In this figure, we show a case with S = 3, where the red particles belong to species 1, greenones to species 2, and the yellow ones to species 3. D. Difference 3: Consideration of enteringsequence of particles
Third, we consider the sequence in which the particlesenter the lattice (i.e., the ‘entering sequence’), which isthe most important feature in our model. Specifically,particles form a queue before the left boundary and enterthe lattice according to the sequence, as illustrated inFig. 4. In the present study, we investigate two types ofsequences: ‘random sequences’ and ‘group sequences,’ asillustrated in Fig. 5.In a random sequence, particles line up randomly re-gardless of their hopping probabilities. A random se-quence thus has N ! / Q Ss =1 ( r s N )! patterns. Note that inreal situations without any controls, random sequencescan be assumed to occur spontaneously.On the other hand, in a group sequence, particles formgroups of the same species and line up group by group.There are S ! possible patterns of group sequences, whichare clearly among the random sequences.For the case S = N , where all hopping probabilitiesare different, we bunch the particles with similar hop-ping probabilities close together with each other as muchas possible, imaginarily considering them as ‘continuousgroups.’ Consistent with this idea, we define a group se-quence with S = N as either an ascending or a descend-ing sequence. Note that we define such a sequence byconsidering the rightmost particle to be the first particlein the sequence.We define the transportation times for the random andgroup sequences to be T R and T G , respectively. E. Difference 4: Introduction of the sorting cost
Finally, we introduce the cost of sorting the particlesand investigate the effect of the sorting cost on the trans-portation time. Here, we define the sorting cost as theminimal number of exchanges K ( τ a , τ b ) necessary to sortthe particles form sequence τ b to sequence τ a , where τ a and τ b represent the sequence after sorting and be-fore sorting, respectively. Note that the arguments of K ( τ a , τ b ) will be abbreviated in obvious cases.In the present study, τ a ( τ b ) correspond to τ G ( τ R ), (cid:1) (cid:1) (cid:2) L sites (cid:1)(cid:1) Sequence Sequence 2Sequence FIG. 4. (Color online) Schematic illustration of the enteringsequences of particles. In this figure, we show threeexamples among all 90 { = 6!/(2!2!2!) } possible sequences forthe case N = 6, S = 3, and r = r = r = 1 /
3. Note thatSequence 1 is one example of a group sequence, whereas theothers are examples of random sequences. (cid:1) (cid:1) (cid:1)(cid:1) (cid:1) (cid:1)(cid:1)(cid:1)
FIG. 6. (Color online) Two examples with K ( τ a , τ b ) = 2 forthe case N = 6 and S = 3. For each sequence τ R , we chosethe one of all 6 (= 3!) possible sequences τ G so that K ( τ a , τ b ) is minimized. III. SIMULATION RESULTS
In this section, we use numerical simulations to inves-tigate the dependence of the transportation time on theentering sequence of the particles.In all the simulations below, we set L = 200 and N =10 , L an N in AppendixA. We determine the value of T for each parameter, andaverage T over 100 trials for Fig. 7 and over 10 trials forFigs. 8, 9, and 10). A. Without sorting cost ( λ = 0 ) In this subsection, we set λ = 0, i.e., we do not includethe sorting cost.In Fig. 7 we plot the simulation values of the numberof particles that have not yet exited the lattice at time t for S ∈ { , , N } . We fix ( α, β ) = (1 ,
1) for (a)–(c) and( α, β ) = (0 . , .
2) for (d)–(f). In the figures, we referto the number of particles that have not yet exited thelattice at time t simply as the ‘remaining particles.’ Thesimulation starts at t = 0, and the number of particlesbecomes 0, i.e., the N th particle exits the lattice, at t = T .We note two important phenomena in Figs. 7 (a)–(c).First, surprisingly, T G is smaller than T R for all threevalues of S when α = β = 1. This result implies thatthe group sequences yield smoother transportation thanthe random ones for the cases ( α, β ) = (1 , T G seems not to depend upon the order of each group inthe group sequence, which can take S ! possible patterns.On the other hand, in Figs. 7 (d)–(f), unlike the casesin Figs. 7 (a)–(c), the difference between T R and T G seems to vanish.In order to compare the difference between T R and T G for various ( α, β ), we define ∆ T as the ratio of the changefrom T R to T G ; that is,∆ T = T G − T R T R . (1) From this definition of ∆ T , ∆ T < T >
0) indicatesthat group (random) sequences are preferable for smoothtransportation. Note that in the following, to calculate∆ T we assume that each group in a group sequence isarranged in ascending order in terms of species number s . The simulation values of ∆ T for various ( α, β ) with(a) S = 2, (b) S = 3, and (c) S = N are plotted inFig. 8. Note that the black lines represent the bound-aries between the low-density/high-density (LD/HD) andthe maximal current (MC) phases of the single-speciesTASEP with hopping probability p (boundary A) in Fig.8 (a), and p (boundary B1) and p (boundary B3) in Fig.8 (b), respectively.Figure 8 shows that for all three values of S , ∆ T issmall in the region where min( α, β ) is relatively large.[In Fig. 8 (b), ∆ T finally yields to a constant value inthe upper-right region beyond boundary B3.] On theother hand, ∆ T is small in the region where min( α, β ) isrelatively small. [In Fig. 8 (a) and (b), ∆ T is almost 0,especially in the lower-left region beyond the boundaryA or B1.] Here, we term the region with ∆ T < T R > T G ), whereas wedesignate the region with ∆ T ≈ T R ≈ T G ), if it exists.These results indicate that group sequences can maketransportation smoother than random sequences whenthe system is mainly governed by the bulk region of thelattice, but the dependence on the type of sequences van-ishes (or decreases) when the system is mainly governedby the boundaries. B. With sorting cost ( λ > ) In this subsection, we consider the sorting cost by vary-ing λ for the same parameter sets in the previous subsec-tion. Appendix B presents specific schemes for obtainingthe minimal number of exchanges necessary to sort theparticles in the simulations.Figure 9 plots ∆ T for (a) S = 2, (b) S = 3, and(c) S = N as functions of λ for various ( α, β ) ∈{ (0 . , . , (0 . , . , (1 , } , which are plotted as blackcrosses in Fig. 8. We emphasize again that in the regionwith ∆ T < T > λ = 0 correspond to those obtained withoutconsidering the sorting cost.As discussed in the previous subsection, we note that∆ T ≤ α, β ) when λ = 0, indicatingthat sorting is almost always beneficial for smooth trans-portation. However, once the sorting cost is considered,the sign of ∆ T can become positive, especially in the re-gion where min( α, β ) is relatively small, indicating thatsorting is not always beneficial. Note that the curves of( α, β ) =(0.6, 0.6) and (1, 1) are observed to overlap eachother in Fig. 9 (b) and (c), unlike Fig. 9 (a). This hap- FIG. 7. (Color online) Simulation values of the number of particles remaining at time t with λ = 0 for (a) S = 2, (b) S = 3,and (c) S = N with ( α, β ) = (1 , S = 2, (e) S = 3, and (f) S = N with ( α, β ) = (0 . , . S = 2, S = 3,and S = N , respectively, we set ( p , p ; r ) = (0 . ,
1; 0 . p , p , p ; r , r , r ) = (0 . , . , .
8; 0 . , . , . p s = 1 − . N − s ) / ( N −
1) ( s = 1 , , ......, N ), respectively, fixing λ = 0. The notation ’ p s → p t ’ means that a group ofspecies s is followed by a group of species t . pens because there is no difference in ∆ T at these twopoints when λ = 0, as we can see in Figs. 8 (b) and (c).Figure 10 plots ∆ T for (a) S = 2, (b) S = 3, and(c) S = N for various ( α, β ) with λ = 1. In this fig-ure, we note the existence of a new region in which∆ T >
0, which we term a ‘random-advantageous region’( T R < T G ). This new region widens as λ increases, finallyresulting in the complete disappearance of the group-advantageous region for large enough λ . Note that Fig.10 (c) exhibits only a random-advantageous region. IV. THEORETICAL ANALYSES
In this section, we show that the simulation resultscan be theoretically reproduced in some special cases.Specifically, we have succeeded in obtaining a mathemat- ical proof of the appearance of the group-advantageousregion for any group number S ( >
1) when λ = 0. A. Approximate flow of a multi-species TASEP
In this subsection, before calculating T , we brieflydiscuss the steady-state flow Q S of the multi-speciesTASEP that corresponds to a random sequence. Wewrite Q S = Q S ( p , ......, p S ; r , ......, r S ), with the argu-ments abbreviated in obvious cases. When the flow Q is simulated for each parameter set, we first evolve thesystem for 10 time steps and then average over the next10 time steps. FIG. 8. (Color online) Simulation values of ∆ T for various ( α, β ) with (a) S = 2, (b) S = 3, and (c) S = N . The parametersother than ( α, β ) are the same as in Fig. 7. Note that three black crosses in each panel represent ( α, β ) =(0.1, 0.2), (0.6, 0.6),and (1, 1), respectively. The color scale at the right of each panel represents the value of ∆ T . -0.2-0.100.10.2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-0.2-0.100.10.2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 ( (cid:1) , (cid:2) ) = ( , ) (cid:1) (cid:1) (cid:1) ( (cid:1) , (cid:2) ) = ( , )( (cid:1) , (cid:2) ) = ( , ) (a) -0.200.20.40.60.8 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 ( (cid:1) , (cid:2) ) = ( , )( (cid:1) , (cid:2) ) = ( , )( (cid:1) , (cid:2) ) = ( , ) (cid:1) (c) ( (cid:1) , (cid:2) ) = ( , ) (cid:1) ( (cid:1) , (cid:2) ) = ( , )( (cid:1) , (cid:2) ) = ( , ) (b) FIG. 9. (Color online) Simulation values of ∆ T with (a) S = 2, (b) S = 3, and (c) S = N as functions of λ for various( α, β ) ∈ { (0 . , . , (0 . , . , (1 , } . The parameters other than ( α, β ) are the same as in Fig. 7.FIG. 10. (Color online) Simulation values of ∆ T for various ( α, β ) with (a) S = 2, (b) S = 3, and (c) S = N for λ = 1. Theother parameters are the same as in Fig. 7. P P ∗ P P ∗ P P ∗ = − α (1 − α ) β − α )(1 − β ) p p rα rαβ − p β rα (1 − β ) 0 1 − β − r ) α (1 − r ) αβ − p β − r ) α (1 − β ) 0 0 0 1 − β P P ∗ P P ∗ P P ∗ (2) FIG. 11. (Color online) Simulation (circles) and theoretical (curves) values of (a) Q , (b) Q , and (c) Q N for L = 2 asfunctions of α for various β ∈ { . , . , } . The other parameters are fixed at (a) ( p , p ; r ) = (0 . ,
1; 0 . p , p , p ; r , r , r ) = (0 . , . , .
8; 0 . , . , . p s = 1 − . S − s ) / ( S −
1) ( s = 1 , , ......, N ) and r s = 1 /N (= 1 / , L = 2 This subsection presents the derivation of an approx-imate Q S based on a Markov chain model. Due to thedifficulty of considering general values of L (the lengthof the lattice) and S (the number of particle species),we consider the simplest case—with L = 2 and S = 2.As two species of particles exist—that is, particles withhopping probability p and particles with p —each sitemay have three states: ‘unoccupied (state 0),’ ‘occupiedby a particle 1 (state 1),’ and ’occupied by a particle2 (state 2).’ This results in 9 possible states; however,noting that it is not necessary to distinguish the particleat site 1 because it always leaves the lattice with prob-ability β , the number of possible states can be reducedto 6. Here, we define the probability distribution P ij ( i = 0 , , , j = 0 , ∗ ), where i and j represent the statenumber of site 0 and 1, respectively. Note that state ∗ indicates either of state 1 or 2.The master equations for the steady state are summa-rized in Eq. (2), using the relation r (= r )+ r = 1. Notethat r and r are replaced with r and 1 − r , respectively,for the case S = 2. In addition, P ij must satisfy thenormalization condition X i =0 P i + X i =0 P i ∗ = 1 . (3)From Eqs. (2) and (3), the flow of the system canbe written as a function of p , p , and r ; that is, Q ( p , p ; r ), is given by the following expression: Q ( p , p ; r ) = p P + p P = p p A { (1 − r ) p + rp } A + p p B , (4)where A = αβ ( α + β − αβ ) (5)and B = α + β − α β − αβ + αβ. (6) The specific forms of the probability distributions aresummarized in Appendix C.For r = 1 and p = p , the system reduces to the single-species TASEP with the flow Q ( p ), where Q ( p ) = pApB + A . (7)Note that the flow of the single-species TASEP for gen-eral L is exactly solved in Ref. [63]. Therefore, as-suming that the value p = p h satisfies the condition Q ( p , p ; r ) = Q ( p ), we can derive p h = p p (1 − r ) p + rp . (8)The quantity p h is termed the harmonic mean of p and p . Accordingly, for L = 2, Q is equivalent to Q ( p = p h ). This relation holds for any species number S ( > β ∈ { . , . , } with (a) S = 2, (b) S = 3 and (c) S = N (= 10 , (cid:1) (cid:1) (cid:1) p (cid:1) p (cid:2) p (cid:2) platoon FIG. 12. (Color online) Schematic illustration of a platoon.In this figure, we set S = 2, with the red particles belongingto species 1 (faster) and the green ones to species 2 (slower).A green particle blocks the red particles behind it, so thatthe trailing red particles cannot hop with probability p butonly with probability p , which is less than p h .
2. General L ( > For general L and S , it is complicated to solve themaster equations. Therefore, in this subsection, we in- FIG. 13. (Color online) Phase diagrams. The color bars indicate the simulation values of (a) Q ( p = p min ), (b) Q ( p , p ; r ),and (c) Q ( p = p h ), respectively, obtained by fixing ( p , p ; r ) = (0 . ,
1; 0 . p min = 0 .
5, and p h = 1 / (0 . / . . /
1) = 2 / α, β ) =(0.1, 0.2), (0.2, 0.1), and (1,1), respectively.FIG. 14. (Color online) Calculated values of ∆ Q for various( α, β ), fixing ( p , p ; r ) = (0 . ,
1; 0 . stead assume an inequality, based on the results in theprevious subsection and the qualitative discussions, andconfirm the validity of the inequality by the simulations.First, for general L and S , Q S is clearly larger than Q ( p = p min ), where p min = min { p , p , ......, p S } .In addition, for L >
2, a platoon can be observed in thebulk of the lattice, in which a slower particle behaves asa bottleneck, and faster particles behind it cannot hopwith a probability larger than that of the smaller one,i.e., less than p h , as shown in Fig. 12. This phenomenonsuppresses the flow, implying that Q S is smaller than Q ( p = p h ).Consequently, Q S satisfies the following inequality; Q ( p = p min ) < Q S < Q ( p = p h ) . (9)In this subsection, we hereafter consider the case S = 2.Figure 13 shows the phase diagrams obtained byplotting the simulation values for (a) the single-species TASEP with p = p min , (b) the two-species TASEP, and(c) the single-species TASEP with p = p h , respectively.Note that Q ( Q ) are the simulation (theoretical) values(and similarly hereafter).Comparing these three figures shows that Eq. (9) ob-viously holds. In addition, as in Figs. 13 (a) and (c), wefind that three different phases—HD, LD, and MC—alsoexist in Fig. 13 (b). Due to Eq. (9), the boundariesbetween the LD (HD) and MC phases of Fig. 13 (b) liebetween those of Figs. 13 (a) and 13 (c). Note that theblack lines in Figs. 13 (a) and (c) are theoretical bound-aries, based on the fact that the boundary between theLD and MC phases of the single-species TASEP withhopping probability p [27] can be written as (cid:26) α = 1 − √ − p ∧ − √ − p < β < β = 1 − √ − p ∧ − √ − p < α < . (10)Here, as for ∆ T , we define ∆ Q as the ratio of the changefrom Q ( p = p h ) to Q ; that is,∆ Q = Q − Q ( p = p h ) Q ( p = p h ) , (11)and we note that ∆ Q = 0 when Q ( p = p h ) = 0.Figure 14 shows ∆ Q for various ( α, β ), for the fixedparameter set ( p , p ; r ) = (0 . ,
1; 0 . p (boundary 1) and the single-species TASEPwith p = p h (boundary 2). Therefore, the lower-left(upper-right) region beyond boundary 1 (boundary 2)corresponds to the LD/HD (MC) phases both for the two-species TASEP and for the single-species TASEP with p = p h . This figure confirms that ∆ Q starts from 0 inthe LD/HD phase, decreases, and finally yields to a con-stant value in the MC phase as ( α, β ) approaches theupper right.Figure 15 plots ∆ Q as a function ( p , p ) for vari-ous ( α, β ) ∈ { (0 . , . , (0 . , . , (1 , } , fixing r = 0 . FIG. 15. (Color online) Calculated values of ∆ Q for various ( p , p ), S = 2, and r = 0 .
5, (a) ( α, β ) = (0 . , .
2) (LD), (b)( α, β ) = (0 . , .
1) (HD), and (c) ( α, β ) = (1 ,
1) (MC). Q (cid:1) / Q (cid:2) (cid:2) p (cid:3) p (cid:3) (cid:4) (cid:1)(cid:1)(cid:1) L = L = L = (b) S = S = S = N (= 10,000) L = L = L = L = L = L = L = L = L = L = L = L = L = L = L = FIG. 16. (Color online) Simulation values of the ratio (a) Q /Q ( p = p h ), (b) Q /Q ( p = p h ), and (c) Q N /Q ( p = p h ) asfunctions of α for various L ∈ { , , , , , } , fixing β = 0 .
6. The other parameters are fixed at (a)( p , p ; r ) = (0 . ,
1; 0 . p , p , p ; r , r , r ) = (0 . , . , .
8; 0 . , . , . p s = 1 − . S − s ) / ( S −
1) ( s = 1 , , ......, N )and r s = 1 /N (= 1 / , Note that both the single-species TASEP with hoppingprobability p h and the two-species TASEP exhibit theLD, HD, and MC phases with ( α, β ) =(0.1, 0.2), (0.2,0.1), and (1, 1), respectively. This is because ∀ ( p , p )(0 . ≤ p , p ≤
1) 0 . < − √ − p min ≤ − √ − p h .For example, Fig. 13 confirms that those three points ex-ist within each corresponding phase for ( p , p ) = (0 . , Q is approx-imately 0, whereas ∆ Q deviates from 0 in Fig. 15 (c),as is also observed in Fig. 14. These phenomena can beexplained as follows.First, in the LD/HD phase, Q is mainly governed bythe input/output probability, leading to ∆ Q →
0, i.e., Q approaches Q ( p = p h ). This is because Q deviatesfrom Q ( p = p h ) mainly due to the existence of platoons,which do not influence the flow much in this phase. Notethat ∆ Q decreases as α or β approaches 0.2, becausethe influence of platoons increases, approaching the MCphase of the two-species TASEP.On the other hand, in the MC phase, Q is mainlygoverned by the bulk region of the lattice. Therefore,the existence of platoons has a more critical influence on Q , causing ∆ Q to deviate from 0; i.e., Q < Q ( p = p h ). Especially as | p − p | increases, the extent of the deviation also increases. This is because the effect ofplatoons increases when there is a large gap between p and p .Figure 16 plots (a) Q /Q ( p = p h ), (b) Q /Q ( p = p h ), and (c) Q N /Q ( p = p h ) for various L ∈{ , , , , , } . Both of Q S ( S = 2 , , , Q ( p = p h ) are obtained by the simulations. In all thefigures, we observe that for L > Q S /Q ( p = p h ) gen-erally becomes less than 1, i.e., Q S < Q ( p = p h ), espe-cially when α increases, i.e., the system approaches andexhibits the MC phase. Note that for S = 10 , L = 2 , , L ≥
10. Those results imply thatEq. (9) and its qualitative discussions can be applicablefor general L and S . B. Relation between T R and T G without the sorting cost Hereafter, we assume p < p < ...... < p S , S > ∀ s r s >
0, and α > λ = 0, i.e., we do not considerthe sorting cost. If for any number of particle species s ,0 r s N is large enough for T G and T R to be determined bythe steady-state flow (see Appendix E), we obtain T R ≈ NQ S (12)and T G ≈ S X s =1 r s NQ ( p = p s ) . (13)Note that this approximation immediately implies the in-dependence of T G from the order of the group sequence.Strictly speaking, T G can differ depending on the orderof each group in the group sequence. However, that dif-ference can be ignored for large N (see Appendix E).In addition, we define the transportation times of theparticles with the same hopping probabilities p h and p min as T H ≈ NQ ( p = p h ) (14)and T M ≈ NQ ( p = p min ) , (15)respectively. From Eqs. (9), (12), (14) and (15), weimmediately obtain the inequality T H < T R < T M . (16)In the following, we show that a general relation be-tween T R and T G can be obtained mathematically forgeneral S ( > T H and T G , and not by comparing T R and T G directly. Here, we introduce the new func-tion f ( α, β ; p , ......, p S ; r , ......, r S ), which is defined asfollows: f ( α, β ; p , ......, p S ; r , ......, r S ) = T H − T G . (17)Because we can assume α < β without lossof generality, we adopt this assumption in thefollowing discussion, writing in abbreviated form f ( α, β ; p , ......, p S ; r , ......, r S ) = f ( α ). Note that forcases with α ≥ β , the theoretical results can be obtainedsimply by replacing α (LD) with β (HD).A contour map of f ( α ) in the ( α, β ) plane exhibits fourlarge regions, which are summarized in Tab. I.In the following subsections, we examine the behaviorof f ( α ) according to this classification.
1. Region 1: α < − √ − p In this region, all the steady-state phases of the single-species TASEP for any p s exhibit the LD phase. Here, TABLE I. Classification of RegionsRegion No. Range1 α < − √ − p − √ − p ≤ α < − √ − p h − √ − p h ≤ α < − √ − p S − √ − p S ≤ α the steady-state flow for the single-species TASEP withparallel updating [27] is given by Q ( p ) = α p − αp − α for LD phase .
12 (1 − p − p ) for MC phase . (18)Therefore, we obtain T G ≈ S X s =1 r s N ( p s − α ) α ( p s − α ) (19)and T H ≈ N ( p h − α ) α ( p h − α ) . (20)From Eqs.(19) and (20), we obtain f ( α ) as f ( α ) = N ( p h − α ) α ( p h − α ) − S X s =1 r s N ( p s − α ) α ( p s − α ) . (21)After some calculations, we obtain f ( α ) <
0; (22)the detailed derivation is given in Appendix F.
2. Region 2: − √ − p ≤ α < − √ − p h This region is further divided into ( u −
1) subregions,as summarized in Tab. II.
TABLE II. Classification of subregions in Region 2Subregion No. Range2–1 1 − √ − p ≤ α < − √ − p − √ − p u ≤ α < − √ − p u +1 ... ...2– v − √ − p v ≤ α < − √ − p v +1 ... ...2–( u −
1) 1 − √ − p u − ≤ α < − √ − p h In Subregion 2– v ( v = 1 , , ......, u − p = p , ......, p v , p h exhibits the MC1phase, whereas that with p = p v +1 , ......, p S displays theLD phase. Therefore, using Eq. (18), we obtain Q ( p ) =
12 (1 − p − p ) for p = p , ......, p v − , p h .α p − αp − α for p = p v +1 , ......, p S . (23)From Eqs. (13), (14), and (23), T G and T H can be writtenas follows; T G ≈ v X s =1 r s N − √ − p s + S X s = v +1 r s N ( p s − α ) α ( p s − α ) (24)and T H ≈ N ( p h − α ) α ( p h − α ) . (25)From Eqs. (24) and (25), we thus obtain f ( α ) in theform f ( α ) ≈ N ( p h − α ) α ( p h − α ) − v X s =1 r s N − √ − p s − S X s = v +1 r s N ( p s − α ) α ( p s − α ) . (26)For 1 − √ − p v ≤ α < − √ − p v +1 , due to α
s = v + 1 , ......, S and p h − α >
0. The function f ( α ) is continuous and differentiable with respect to α ,including at each boundary (see Appendix G). However,the signs of f ( α ) and df ( α ) /dα are not specified. Notethat the following conditionlim α → q h − df ( α ) dα > , (27)where q h = 1 − √ − p h indicates that f ( α ) increasesmonotonically at least near the boundary between Sub-region 2–( u −
1) and Region 3. This is discussed in Ap-pendix H.
3. Region 3: − √ − p h ≤ α < − √ − p S Similarly to Region 2, this region is further dividedinto ( S − u + 1) subregions, as summarized in Tab. III.Note that Subregion 3–1 vanishes in the case p h = p u ,resulting in ( S − u ) subregions.In Region 3– v ( v = u, u + 1 , ......, S ), the single-speciesTASEP with p = p , ......, p v − , p h exhibits the MC phase,whereas that with p = p v , ......, p S displays the LD phase.Therefore, using Eq. (18), we obtain Q ( p ) Q ( p ) =
12 (1 − p − p ) for p = p , ......, p v − , p h .α p − αp − α for p = p v , ......, p S . (28) TABLE III. Classification of subregions in Region 3Subregion No. Range3– u − √ − p h ≤ α < − √ − p u u + 1) 1 − √ − p u ≤ α < − √ − p u +1 ... ...3– v − √ − p v − ≤ α < − √ − p v ... ...3– S − √ − p S − ≤ α < − √ − p S From Eqs. (13), (14), and (28), T G and T H can be writtenas follows; T G ≈ v − X s =1 r s N − √ − p s + S X s = v r s N ( p s − α ) α ( p s − α ) (29)and T H ≈ N − √ − p h . (30)From Eqs. (29) and (30), f ( α ) becomes f ( α ) ≈ N − √ − p h − v − X s =1 r s N − √ − p s − S X s = v r s N ( p s − α ) α ( p s − α ) . (31)For 1 − √ − p v − ≤ α < − √ − p v , due to α < p v <...... < p S , we obtain p s − α > s = v, ......, S . Sim-ilarly to Region 2, f ( α ) is continuous and differentiablewith respect to α including at each boundary.After some calculations,we find that df ( α ) /dα satisfies df ( α ) dα > f ( α ) is a monotonically increasingfunction of α throughout Region 3.
4. Region 4: − √ − p S ≤ α In this region, all the steady-state phases of the single-species TASEP for any hopping probabilities p s and p h exhibit the MC region. Note that this region vanishes inthe case p S = 1 because 1 − √ − p S = 1.Therefore, we obtain T G ≈ S X s =1 r s N − √ − p s (33)and T H ≈ N − √ − p h , (34)2both of which are independent of α . From Eqs. (33) and(34),we thus obtain f ( α ) as f ( α ) = 2 N − √ − p h − S X s =1 r s N − √ − p s . (35)After some calculations, we obtain f ( α ) > , (36)the detailed derivation of which is given in Appendix J.
5. Relation between T R and T G With the results of Subsec. IV B 1–IV B 4, we can ob-tain a general relation between T R and T G for some spe-cial cases. Table IV summarizes the signs of f ( α ) and df ( α ) /dα in each region. Note that ‘U’ indicates thatthe sign is unclear. TABLE IV. Sign of f ( α ) and df ( α ) /dα Region No. f ( α ) df ( α ) /dα − U2 U U3 U +4 + 0
Considering Tab. IV and the continuity of f ( α ) includ-ing at each boundary (see Appendix G), we find from theintermediate value theorem that ∃ α cr such that f satisfies f ( α = α cr ) = 0 ⇔ T H = T G (37)in Region 2 or 3. The specific conditions that α cr mustsatisfy are given in Appendix K.Defining α cr , max as the largest value among the quan-tities α cr , we obtain f ( α ) > T H > T G —in theregion where α > α cr , m . This is because f ( α ) is continu-ous and increases monotonically from a point in Region2 (and through Region 3), to yield f ( α ) > T G < T H < T R (38)in the region α > α cr , m . Eq. (38) means that ∆ T < α, β ) is relatively large. This result indicates thatthe group-advantageous region must appear even in acase with p S = 1, for which Region 4 vanishes.In analogy with the discussion above, we can also pre-dict that a region with ∆ T < S = N . C. Relation between T R and T G with sorting cost In this subsection, we discuss the change in the relationbetween T R and T G when λ >
0, i.e., when the sortingcost is included. In the following, we first obtain a generalformula for the sorting cost and then evaluate upper andlower limits to λ .
1. General formula for the sorting cost
First, we calculate mathematically the averaged mini-mal number of exchanges necessary to sorting the parti-cles from random to group sequences.We here define K as the averaged value of K , usingthe fact that τ R can take N ! / Q Ss =1 ( r s N )! patterns withequal probability. We thus have K = Q Ss =1 ( r s N )! N ! X ∀ τ R K ( τ G , τ R ) . (39)If K ′ ( τ G , τ R ) is the minimal number of exchanges nec-essary to sort the particles from a random sequence τ R toa given fixed group sequence τ G , then K ′ ( τ G , τ R ) satisfies K ′ ( τ G , τ R ) = min { K ( τ G , τ R ) } . (40)Note that the number of elements of { K ( τ G , τ R ) } is equalto that of { τ G } from the definition. Eqs. (39) and (40)indicate that the best group sequence τ G can vary de-pending on the particular random sequence τ R .Due to the difficulty of a general calculation of K , weinstead calculate K ′ , which is defined as follows: K ′ = Q Ss =1 ( r s N )! N ! X ∀ τ R K ′ ( τ G , τ R ) , (41)where τ G is a fixed sequence out of the set { τ G } for allpossible τ R .For S = 2 and S = N , K ′ can be generally calculatedas K ′ = r (1 − r ) N for S = 2 ,N − N X k =1 k for S = N, (42)the detailed derivations of which are discussed in Ap-pendix L.Figure 17 shows the ratio K/K ′ for (a) S = 2 and (b) S = N . Both figures show that K/K ′ ≈
1, i.e., K ≈ K ′ ,indicating that there is no problem in substituting K ′ for K for large enough N .In the following calculations, we therefore use K ′ in-stead of K because K ′ can be represented by a generalformula, whereas K cannot.3 r (cid:1) (cid:1) / (cid:1) (cid:1) (cid:2) (a) N = N = N = N (b) (cid:1) (cid:1) / (cid:1) (cid:1) (cid:2) FIG. 17. (Color online) (a) Simulation values of the ratio
K/K ′ for various N ∈ { , , , , , } for S = 2. (b)Simulation values of the ratio K/K ′ for S = N . Note that K and K ′ are each simulation values obtained byrespectively averaging over 100 trials.
2. Upper and lower limits to λ cr We first define λ = λ cr ≥
0, as the value for which T R = T G . Note that λ cr is defined to be equal to 0if T R ≤ T G when λ = 0. From the definition of λ cr ,the random-advantageous region appears when λ > λ cr .Based on the discussions in Subsec. IV B and IV C 1, wehere evaluate λ cr for S = 2 and α < β .To take into account the sorting cost, we add the term λK ( λ >
0) to T G ; that is, T G ≈ λK ′ + S X s =1 r s NQ ( p = p s ) . (43)Conversely, we do not add that term to T R because arandom sequence means a sequence without sorting.We also define λ H and λ M as the values of λ for which T H = T G and T M = T G , respectively. Note that λ H canhave negative values, because T H can be less than T G .From Eq. (16) and Subsec. IV B, when λ = 0 therelations among T R , T G , T H , and T M must satisfy one ofthe following three inequalities: T H < T R ≤ T G < T M , (44) T H < T G < T R < T M , (45) or T G < T H < T R < T M . (46)Therefore, the relations of λ cr , λ H , and λ M can be writtenas follows: max(0 , λ H ) ≤ λ cr < λ M , (47)where we note that by definition λ cr ≥
0, whereas λ H canbe either negative or positive, while λ M must be positive.In Region 1, i.e., α < −√ − p , where T H < T G with λ = 0, λ H must be negative, while λ M must be positive,and satisfy N ( p − α ) α ( p − α ) ≈ λ M r (1 − r ) N + rN ( p − α ) α ( p − α ) + (1 − r ) N ( p − α ) α ( p − α ) . (48)In Region 2, i.e., 1 − √ − p ≤ α < − √ − p h , λ H can be either negative or positive, whereas λ M must bepositive. The quantities λ H and λ M satisfy rN ( p h − α ) α ( p h − α ) ≈ λ H r (1 − r ) N + 2 rN − √ − p + (1 − r ) N ( p − α ) α ( p − α ) (49)and 2 N − √ − p ≈ λ M r (1 − r ) N + 2 rN − √ − p + (1 − r ) N ( p − α ) α ( p − α ) , (50)respectively.In Region 3, i.e., 1 − √ − p h ≤ α < − √ − p , λ H can be either negative or positive, and λ M must bepositive. Thus, λ H and λ M satisfy2 N − √ − p h ≈ λ H r (1 − r ) N + 2 rN − √ − p + (1 − r ) N ( p − α ) α ( p − α ) (51)and 2 N − √ − p ≈ λ M r (1 − r ) N + 2 rN − √ − p + (1 − r ) N ( p − α ) α ( p − α ) , (52)respectively.4 TABLE VI. Upper and lower limits to λ cr Region No. Upper and lower limits to λ cr ≤ λ < rα p − α p − α − p − α p − α ! ( , r (1 − r ) p h − α α ( p h − α ) − r − √ − p − (1 − r )( p − α ) α ( p − α ) !) ≤ λ cr < r − √ − p − p − α α ( p − α ) ! ( , r (1 − r ) − √ − p h − r − √ − p − (1 − r )( p − α ) α ( p − α ) !) ≤ λ cr < r − √ − p − p − α α ( p − α ) ! r (1 − r ) − √ − p h − r − √ − p − − r − √ − p ! < λ cr < r − √ − p − − √ − p ! FIG. 17. (Color online) Simulation values (black circles) and the theoretical existence range of λ cr (yellow region) as functionsof α . The other parameters are fixed at ( β ; p , p ; r ) = (a) (1;0.5,1;0.5) and (b) (1;0.5,0.6;0.5). In Region 4, i.e., 1 − √ − p ≤ α , due to T G < T H In the present study, we have used a modified TASEPto analyze the dependence of the transportation time onthe entering sequences of particles, using both the nu-merical simulations and theoretical analyses.Here, we summarize a number of important results.In Sec. III, we discovered that there exists an impor-tant ‘group-advantageous region’ where T R > T G whenmin( α, β ) is relatively large and the sorting costs are ne-glected. When sorting costs are introduced, a new re-gion called a ‘random-advantageous region’ appears with T R < T G . In addition, the group-advantageous regionshrinks and finally disappears as λ increases. We ex-plored these phenomena for various S ∈ { , , N } .In Sec. IV, we analyzed the simulation results by5employing mathematical approaches for certain specialcases. Using some approximations, we have showntheoretically that without the sorting cost the group-advantageous region must appear for any parameter sets( S, p s , r s ). Moreover, we have succeeded in deriving theupper and lower limits to the value of λ cr where T R = T G by obtaining a general formula for the sorting cost.Our findings can be applied to real-world situations,such as providing efficient operation for various tasksand smooth logistics for various products and yieldingan effective evacuation method for pedestrians. Specif-ically, for smooth operation, we can determine whetherwe should begin tasks without considering the operationsequence or otherwise. Similarly, for smooth logistics,we can select whether the products should be bunchedwith nearly equal sizes. Furthermore, for ensuring effec-tive evacuation of pedestrians, we can determine whetherthe bunching of pedestrians having nearly equal veloci-ties should be conducted before transportation. The cri-teria for these judgments depend on the magnitude ofthe consideration or bunching cost ( λ ). Note that thesemagnitudes significantly differ from each other, i.e., con-sidering only the sequence of tasks is typically deemedcheaper (have a smaller λ ) than sorting various pedestri-ans and products. ACKNOWLEDGMENTS This work was partially supported by JST-Mirai Pro-gram Grant Number JPMJMI17D4, Japan, JSPS KAK-ENHI Grant Number JP15K17583, and MEXT as ’Post-K Computer Exploratory Challenges’ (Exploratory Chal-lenge 2: Construction of Models for Interaction AmongMultiple Socioeconomic Phenomena, Model Develop-ment and its Applications for Enabling Robust andOptimized Social Transportation Systems) (Project ID:hp180188). Appendix A: Validity of our selection of L and N In this Appendix, we briefly discuss the validity of se-lecting L = 200 and N = 10 , L , we comparethe simulation values of Q for L = 200 and L = 1 , Q /Q ′ , where Q and Q ′ rep-resent the flow of the multi-species TASEP with L = 200and L = 1 , α for various β ∈ { . , . , } . The result that Q /Q ′ ≈ L = 200. Thus, wechoose L = 200 to decrease the simulation time.On the other hand, the assumption that T is deter-mined by a steady-state flow may be inappropriate forsmall N . Therefore, we have compared the results for N = 10 , 000 and N = 20 , L = 200.Figure 19 shows the ratio T /T ′ , where T and T ′ rep-resent the transportation times for N = 10 , 000 and (cid:1) (cid:1) = (cid:1) = (cid:1) = (cid:1) (cid:1) / (cid:1) (cid:1) (cid:2) FIG. 18. (Color online) Simulation values of the ratio Q /Q ′ as a function of α for various β ∈ { . , . , } . Theother parameters are fixed at ( p , p ; r ) = (0 . , 1; 0 . N = 20 , α for various β ∈ { . , . , } . The result that T /T ′ ≈ . 5, i.e., that T is proportional to N , indicates that the assumption canbe regarded as valid for N = 10 , N = 10 , 000 similarly to decrease the simulation time. (cid:1) T / (cid:1) (cid:1) (cid:1) = (cid:1) = (cid:1) = FIG. 19. (Color online) Simulation values of the ratio T /T ′ as a function of α for various β ∈ { . , . , } . The otherparameters are fixed at L = 200 and ( p , p ; r ) = (0 . , 1; 0 . Appendix B: Simulation schemes for obtaining theminimal number of necessary exchanges In this Appendix, we briefly describe the specific simu-lation schemes we used to obtain K . We emphasize thatthe cost of counting or comparing particles and the dis-tances between exchanged particles are both ignored inthe following.First, for S = 2, τ a = τ G can have only one of twopatterns. Once τ G is fixed to be either of these two se-quences, we can immediately obtain the number of par-ticles placed at the wrong areas in sequence τ b = τ R ,which is twice as large as the number of necessary ex-changes (see also Appendix L). Consequently, comparingthe results for the two τ G gives the smaller number as K .6Second, for S = 3, τ a = τ G can have six patterns.Once τ G is fixed at one of these six sequences, we canimmediately obtain the number of particles placed at thewrong areas in any sequence τ b = τ R . After selectingone species, which we first replace at the correct loca-tion, we exchange all particles of that species that areplaced in the wrong areas in sequence τ b = τ R . The sub-sequent procedure is similar to the case for S = 2. Conse-quently, comparing the six results for each τ G again givesthe smallest number as K . Note that we can similarlycalculate the numbers for general S > S = N , τ a = τ G can have one of two pat-terns: either an ascending or a descending sequence. Oneexchange is needed for each particle in τ b = τ R for whichthere exists a particle with a smaller (larger) hoppingprobability than the noted particle. This is termed a‘selection sort.’ This procedure starts from the leadingparticle. Consequently, by comparing the results for thetwo τ G , the smaller number is again selected as the min-imal number of necessary exchanges. Appendix C: Probability distribution with L = 2 and S = 2 Here, we summarize the probability distributions with L = 2 and S = 2, which can be obtained from Eqs. (2)and (3). The specific forms are described as follows: P = p p (1 − α ) β { (1 − r ) p + rp } A + p p B = p h (1 − α ) β A + p h B ,P ∗ = p p αβ { (1 − r ) p + rp } A + p p B = p h αβA + p h B ,P = rp A { (1 − r ) p + rp } A + p p B = rp h Ap ( A + p h B ) ,P ∗ = rp p α (1 − β ) { (1 − r ) p + rp } A + p p B = rp h α (1 − β ) A + p h B ,P = (1 − r ) p A { (1 − r ) p + rp } A + p p B = (1 − r ) p h Ap ( A + p h B ) ,P ∗ = (1 − r ) p p α (1 − β ) { (1 − r ) p + rp } A + p p B = (1 − r ) p h α (1 − β ) A + p h B , (C1)where A = αβ ( α + β − αβ ) , (C2) B = α + β − α β − αβ + αβ, (C3) and p h = p p (1 − r ) p + rp . (C4) Appendix D: Q S for general S with L = 2 In this appendix, we prove that for general S with L = 2, Q S is equal to Q ( p = p h ).From the results with L = 2 and S = 2 (see AppendixC), we can conjecture the probability distributions forgeneral S with L = 2 as P = p h (1 − α ) β A + p h B ,P ∗ = p h αβA + p h B ,P s = r s p h Ap s ( A + p h B ) ,P s ∗ = r s p h α (1 − β ) A + p h B , (D1)where s = 1 , , ......, S .On the other hand, the master equations of the steadystate are summarized as 2( S + 1) equations: P = (1 − α ) P + (1 − α ) βP ∗ ,P ∗ = (1 − α )(1 − β ) P ∗ + S X k =1 p k P k ,P s = r s P + r s αβP ∗ + (1 − p s ) P s + βP s ∗ ,P s ∗ = (1 − β ) P s ∗ + r s α (1 − β ) P ∗ , (D2)where s = 1 , , ......, S . In addition, P ij must satisfy thenormalization condition S X i =0 P i + S X i =0 P i ∗ = 1 . (D3)We can confirm that Eqs. (D1) satisfy Eqs. (D2)and (D3). With Penron-Frobenius theorem regardingstochastic matrix, this indicates that Eqs. (D1) areunique solutions for Eqs. (D2) and (D3).From Eqs. (D1), the flow of the system is given by thefollowing expression: Q S = S X s =1 p s P s = S X s =1 p s r s p h Ap s ( A + p h B )= Q ( p = p h ) . (D4)7 Appendix E: Validity of the approximation for T In this Appendix, we briefly demonstrate the validityof Eq. (13).Figure 20 (a) shows the ratio T G , sim /T G , theo as a func-tion of α for various β ∈ { . , . , } with S = 2 and Fig.20 (b) shows the same ratio for S = 3. Note that T G , sim and T G , theo represent the values of T G from the simu-lations and that given by Eq. (13), respectively. Bothfigures show that T G , sim /T G , theo ≈ 1, indicating that Eq.(13) provides a good approximation for T G . (cid:1) (cid:1) = (cid:1) = (cid:1) = (cid:1) (cid:1) , (cid:3) (cid:4) (cid:5) / (cid:1) (cid:1) , (cid:6) (cid:7) (cid:8) (cid:9) (cid:1) (cid:1) = (cid:1) = (cid:1) = (cid:1) (cid:1) , (cid:3) (cid:4) (cid:5) / (cid:1) (cid:1) , (cid:6) (cid:7) (cid:8) (cid:9) FIG. 20. (Color online) Simulations values of the ratio T G , sim /T G , theo as a function of α for various β ∈ { . , . , } with (a) S = 2 and (b) S = 3. The other parameters arefixed at (a) ( p , p ; r ) = (0 . , 1; 0 . 5) and (b)( p , p , p ; r , r , r ) = (0 . , . , . 8; 0 . , . , . Strictly speaking, T G , sim /T G , theo must be larger than1 on average. This is mainly due to the fact that T G , sim includes T , which is the time required for the first parti-cle to reach the right-hand boundary, whereas T G , theo ig-nores that time. This also indicates that T G can differ de-pending on the order of each group in the group sequence(i.e., the hopping probability of the leading group). How-ever, this difference has little influences on the theoreticalresults, as explained below.First, T can be estimated as T ≈ Lp s , (E1)where s = 1 , , ......, S and the time steps before the firstparticle enters the lattice are assumed to be small enough to be ignored. The quantities T and T G without T satisfy T ≈ Lp s < Lp S < L − √ − p S (E2)and T G ≈ S X s =1 r s NQ ( p = p s ) > N − √ − p S , (E3)respectively. Therefore, T /T G reduces to T T G < L N . (E4)Under the proposition that N is large enough, we canassume L/ N << L/ N = 0 . 01 in the present study).In fact, observing the time series of the flows (199-steps central moving average) in Fig. 21, we find thatnearly the entire duration during transportation can beregarded to be in the steady state for large enough r s N .Note that we calculate the flows at time t by averagingnumber of moving particles per bond between t − t . Moreover, all the transportation times T ( T R , T H , and T M ) originally include T , so that this term disappearswhen they are subtracted from each other. Consequently, T (and therefore, the dependence of T G on the order ofeach group in the group sequence) can be assumed to beignorable. FIG. 21. (Color online) Simulation values of 199-stepscentral moving average of flow at time t with( N, α, β ) = (10 , , , S = 2 and S = 3, respectively,we set ( p , p ; r ) = (0 . , 1; 0 . p , p , p ; r , r , r ) = (0 . , . , . 8; 0 . , . , . Appendix F: Discussion of the sign of f ( α ) in Region 1 In this Appendix, we give a detailed derivation of Eq.(22) for Region 1, where α < − √ − p .8Eq. (21) gives f ( α ) f ( α ) = N ( p h − α ) α ( p h − α ) − S X s =1 r s N ( p s − α ) α ( p s − α )= Nα S X s =1 (cid:26) r s ( p h − α ) p h − α − r s ( p s − α ) p s − α (cid:27) = N ( α − p h − α ) Q Ss =1 ( p s − α ) C, (F1)where C = S X s =1 r s ( p h − p s ) Y k = s ( p k − α ) . (F2)The quantity C is calculated as follows: C = S X s =1 r s P St =1 r t /p t − p s ! Y k = s ( p k − α ) = 1 D S X s =1 r s S Y k =1 p k − p s S X t =1 r t Y k = t p k Y k = s ( p k − α ) = 1 D S X s =1 S X t =1 r s r t S Y k =1 p k − p s Y k = t p k Y k = s ( p k − α ) = 1 D S X s =1 X t = s r s r t ( p t − p s ) Y k = t p k Y k = s ( p k − α ) , (F3)where D = S X s =1 r s Y k = s p k . (F4)By regarding the sum of the term with ( s, t ) = ( x, y ) andthat with ( s, t ) = ( y, x ) as a new term for ∃ ( x, y ) ( x, y =1 , , ......, S, x < y ), we can rewrite Eq. (F3) as follows: C = 1 D S X s =1 X t 0, we obtain C > α − < C > 0, we finally obtain f ( α ) < . (F6) Appendix G: Continuity and differentiability of f ( α ) at each boundary In this Appendix, we briefly discuss the continuity anddifferentiability of f ( α ) at each boundary.Defining g ( x ) for 0 < x ≤ g ( x ) = p − x x ( p − x ) for 0 < x ≤ − √ − p, − √ − p for 1 − √ − p < x ≤ , (G1)where 0 < p ≤ 1, the following equations hold:lim x → q − g ( x ) = lim x → q +0 g ( x ) = 21 − √ − p (G2)and lim δ →− g ( x + δ ) − g ( x ) δ = lim δ → +0 g ( x + δ ) − g ( x ) δ = 0 , (G3)where q = 1 − √ − p . Therefore, g ( x ) is continuous anddifferentiable at x = q = 1 − √ − p , resulting in thecontinuity and differentiability of g ( x ) for 0 < x ≤ f ( α ) is represented as a linear sumof terms g ( α ), where p is substituted for p s or p h (0 (1 − p − p v +1 − + p s + p h − p s p h − > (1 − p − p h − + p s + p h − p s p h − p s (1 − p h ) > . (H3)We cannot specify the sign of df ( α ) /dα in this sub-region from Eqs. (H1), (H2), and (H3). However, nearthe boundary between Subregion 2–( u − 1) and 3– u , weobtain the following conditions:lim α → q h − (2 α − α − p h )= 2(1 − p − p h ) − (1 − p − p h ) − p h = 0 (H4)and lim α → q h − { ( α − + p s + p h − p s p h − } = (1 − p − p h − + p u − + p h − p u − p h − p u − (1 − p h ) > , (H5)where q h = 1 − √ − p h . Therefore, noting the ob-vious continuity of df ( α ) /dα for 1 − √ − p v ≤ α < − √ − p v +1 , the region of df ( α ) /dα > u − Appendix I: Discussion of the sign of df ( α ) /dα in Subregion 3– v In this Appendix, we give a proof on Eq. (32) in Sub-region 3– v , i.e., 1 − √ − p v − ≤ α < − √ − p v .From Eq. (31), df ( α ) /dα can be calculated as follows: df ( α ) dα ≈ S X s = v r s N ( p s − α + α ) α ( p s − α ) . (I1)For s = v, ......, S , the quantity p s − α + α satisfies p s − α + α > p s − − p − p v ) + (1 − p − p v ) = p s − p v > . (I2)From Eqs. (I1) and (I2), we finally obtain df ( α ) dα > . (I3) Appendix J: Discussion of the sign of f ( α ) in Region 4 In this Appendix, we give a detailed derivation of Eq.(36), where 1 − √ − p S ≤ α .From Eqs. (33) and (34), f ( α ) can be represented asfollows: T G ≈ N S X s =1 r s (1 + √ − p s ) p s = 2 N Q Ss =1 p s S X s =1 (cid:16) r s + r s p − p s (cid:17) Y t = s p t = 2 N Q Ss =1 p s ( S X s =1 r s Y t = s p t + S X s =1 r s p − p s Y t = s p t ) (J1)and T H ≈ N √ − p h p h = 2 N q − / P Ss =1 ( r s /p s )1 / P Ss =1 ( r s /p s )= 2 N Q Ss =1 p s ( S X s =1 r s Y t = s p t + vuuut S X s =1 r s Y t = s p t − S Y s =1 p s × S X s =1 r s Y t = s p t ) (J2)From Eqs. (35), (J1) and (J2), f ( α ) is given by f ( α ) ≈ N Q Ss =1 p s (vuut S X s =1 r s Y t = s p t ! − S Y s =1 p s × S X s =1 r s Y t = s p t ! − S X s =1 r s p − p s Y t = s p t !) = 2 N Q Ss =1 p s ( E − F ) , (J3)where E = vuuut S X s =1 r s Y t = s p t − S Y s =1 p s × S X s =1 r s Y t = s p t (J4)0and F = S X s =1 r s p − p s Y t = s p t . (J5)From Eqs. (J4) and (J5), E − F becomes E − F = S X s =1 r s Y t = s p t − S Y s =1 p s × S X s =1 r s Y t = s p t − S X s =1 r s p − p s Y t = s p t = S X s =1 r s Y t = s p t − S X s =1 X t = s r s r t Y k = s p k Y l = t p l − S Y s =1 p s × S X s =1 r s Y t = s p t − S X s =1 r s Y t = s p t + S X s =1 r s p s Y t = s p t − S X s =1 X t = s r s r t p − p s p − p t Y k = s p k Y l = t p l = S X s =1 (X t = s r s r t Y k = s p k Y l = t p l − r s (1 − r s ) p t Y k = s p k Y l = t p l − X t = s r s r t p − p s p − p t Y k = s p k Y l = t p l !) = S X s =1 X t = s " r s r t Y k = s p k Y l = t p l × ( − p t − p (1 − p s )(1 − p t ) ) . (J6)Here, regarding the sum of the term with ( s, t ) =( x, y ) and that with ( s, t ) = ( y, x ) as a new term for ∃ ( x, y ) ( x, y = 1 , , ......, S, x < y ), Eq. (J6) can be rewrit- ten as E − F = S X s =1 X t 2, both equations become more than quartic.From its definition of α cr , max , α cr , max can be writtenas α cr , max = max { α cr } , (K1)where { α cr } represents the set of α cr . Appendix L: Derivation of K ′ In this Appendix, we derive the approximate averagedminimal number of exchanges K ′ necessary to sort theparticles for two special cases: S = 2 and S = N . S = 2 First, for a general calculation of K ′ , τ G has to be fixedto be either of the two possible patterns. Once τ G is fixed, K ( τ G , τ R ) can be determined uniquely for all possible τ R .Without loss of generality, we can assume rN ≤ (1 − r ) N and τ G can be fixed as illustrated in the lower panel ofFig. 22.Suppose that for τ R , k (0 ≤ k ≤ rN ) particles ofspecies 1 are located in the Area 2, ( k particles of species2 are located in the Area 1, conversely) as described inthe lower of Fig. 22. Under this supposition, k -time ex-changes are necessary for sorting particles from τ R to τ G .1 TABLE VI. Explicit expressions for f ( α = α cr )=0Region No. Explicit expressions for f ( α = α cr ) = 02 N ( p h − α ) α cr ( p h − α cr ) − v X s =1 r s N − √ − p s − S X s = v +1 r s N ( p s − α ) α cr ( p s − α cr ) = 0 ∧ − √ − p v < α cr < − √ − p v +1 for 1 ≤ v ≤ u − − √ − p v < α cr < − √ − p h for v = u − 13 2 N − √ − p h − v − X s =1 r s N − √ − p s − S X s = v r s N ( p s − α ) α cr ( p s − α cr ) = 0 ∧ − √ − p h < α cr < − √ − p v for v = u − √ − p v − < α cr < − √ − p v for u + 1 ≤ v ≤ S r (cid:1) (cid:2) (cid:3) particles r (cid:1) (cid:2) (cid:3) particles (cid:1) (cid:1) (cid:3) particles (cid:3) particles (cid:1) (cid:1) (cid:1) particles (cid:1) (cid:1) r (cid:1) particles r (cid:1) particles (cid:1) (cid:2) Area 2 Area 1 FIG. 22. (Color online) Schematic illustration of τ G (upperpanel) and τ R (lower panel), where the red particles belongto species 1 and the green ones to species 2. In the upperpanel, we show one example from among all (cid:0) rNk (cid:1) × (cid:0) (1 − r ) Nk (cid:1) possible random sequences, whereas in the lower panel weshow one of the two possible group sequences. Considering that τ R satisfying this supposition possiblyhas (cid:0) rNk (cid:1) × (cid:0) (1 − r ) Nk (cid:1) sequences, a N = P ∀ τ R K ′ ( τ G , τ R )can be written as follows; a N = X ∀ τ R K ′ ( τ G , τ R )= rN X k =1 k (cid:18) rNk (cid:19)(cid:18) (1 − r ) Nk (cid:19) = rN X k =1 rN (cid:18) rN − k − (cid:19)(cid:18) (1 − r ) Nk (cid:19) = rN rN X k =1 ((cid:18) rNk (cid:19)(cid:18) (1 − r ) Nk (cid:19) − (cid:18) rN − k (cid:19)(cid:18) (1 − r ) Nk (cid:19)) . (L1)Using the Vandermonde convolution formula, Eq. (L1) can be rewritten as follows: a N = rN (cid:26)(cid:18) NrN (cid:19) − (cid:18) N − rN (cid:19)(cid:27) . (L2)Because the sequence τ R can take any of N ! / { ( rN )!((1 − r ) N )! } possible patterns with equal probability, we canfinally reduce K ′ to K ′ = ( rN )!((1 − r ) N )! N ! a N = r (1 − r ) N. (L3)Figure 23 compares the simulation (cir-cles) and theoretical (curves) values for various N ∈ { , , , , , } for S = 2.The simulations show very good agreement with ourexact analysis. r (cid:1) (cid:1) (cid:2) (cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:2)(cid:8)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:4)(cid:5)(cid:7)(cid:13)(cid:14)(cid:15)(cid:12)(cid:8)(cid:11)(cid:12)(cid:7)(cid:2)(cid:16)(cid:6)(cid:5)(cid:10)(cid:16)(cid:4)(cid:11)(cid:17)(cid:12)(cid:13) N = N = N = FIG. 23. (Color online) Simulation (circles) and theoretical(curve) values of K ′ as functions of r for various N ∈ { , , , , , } with S = 2.We obtained each of the simulation values by averaging over100 trials. S = N When S = N , τ G also has to be fixed as either ofthe two possible patterns—an ascending or a descend-ing sequence—for a general calculation of K ′ , as illus-trated in the upper panel of Fig. 24. Once τ G is fixed,2 K ′ ( τ G , τ R ) can be determined uniquely for all possible τ R . (cid:1) (cid:1) (cid:1) (cid:1) p (cid:1) p (cid:2) p (cid:3) p (cid:4)(cid:5)(cid:2) p (cid:4)(cid:5)(cid:1) p (cid:4) (cid:1) particles (cid:1) (cid:3) (cid:4) particles (cid:2) (cid:1) (cid:2) FIG. 24. (Color online) Schematic illustration of τ G (upperpanel) and τ G (lower panel) for the case S = N . In the lowerpanel, we show one example of all N × ( N − p < p < ... < p N for the ascendingsequence, whereas p > p > ... > p N for the descending one. If we regard the entire sequence as consisting of twoparts—the first (blue) particle and other ( N − 1) parti-cles, as described in the lower panel of Fig. 24—the sort-ing procedure can also be divided into two parts: sort-ing ( N − 1) particles plus the last exchange for the firstparticle. If the first particle corresponds to the particlewith hopping probability p l ( l = 1 , , ......, N ), and notingthat the sequence for the remaining ( N − 1) particles has( N − b N,l = P ∀ τ ′ R ,l K ( τ G , τ ′ R ,l ) as follows: b N,l = X ∀ τ ′ R ,l K ( τ G , τ ′ R ,l )= (cid:26) a N − for l = 1 ,a N − + ( N − l = 2 , , ......, N, (L4)where N > τ ′ R ,l represents the sequence for whichthe first particle is the particle with hopping probability p l . Note that the last sort is not necessary in the casewhere l = 1.Therefore, for N > 1, we can write a N = P ∀ τ R K ( τ G , τ R ): a N = X ∀ τ R K N ( τ G , τ R )= N X l =1 X ∀ τ ′ R ,l K N ( τ G , τ ′ R ,l )= N X l =1 b N,l = ( N − × ( N − N a N − (L5) Dividing both sides of Eq. (L5) by N !, we obtain c N = c N − + N − N = c + N X k =1 k − k , (L6)where c N = a N /N ! and N > 2. With the initial condi-tion c = a = 0, a N is finally reduced to a N = N ! N − N X k =1 k ! , (L7)which we note holds for the case N = 1.The sequence τ R can take N ! patterns with equal prob-ability, and therefore, K ′ is finally reduced to K ′ = a N N ! = N − N X k =1 k . (L8)Figure 25 compares the simulation (circles) and theo-retical (line) values for S = N . The simulations againshow a very good agreement with our exact analysis. N (cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:2)(cid:8)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:4)(cid:5)(cid:7)(cid:13)(cid:14)(cid:15)(cid:12)(cid:8)(cid:11)(cid:12)(cid:7)(cid:2)(cid:16)(cid:6)(cid:5)(cid:10)(cid:5)(cid:2)(cid:9)(cid:12) (cid:1) (cid:1) (cid:2) FIG. 25. (Color online) Simulation (circles) and theoretical(line) values of K N ′ as functions of N . We obtained each ofthe simulation values by averaging over 100 trials. [1] C. Appert-Rolland, S. Klein, M. Ebbinghaus, andL. Santen, in Proceedings of the Asia-Pacific Econo-physics Conference 2016 Big Data Analysis and Modelingtoward Super Smart Society (APEC-SSS2016) (2017) p.011001.[2] S. Croom, P. Romano, and M. Giannakis, Eur. J. Purch.Supply Manag. , 67 (2000).[3] T. Ezaki, D. Yanagisawa, and K. Nishinari, Physica A , 62 (2015).[4] Y. Taniguchi, R. Nishi, T. Ezaki, and K. Nishinari, Phys-ica A , 304 (2015).[5] M. Woelki, Phys. Rev. E , 062818 (2013).[6] T. Imai and K. Nishinari,Phys. Rev. E , 062818 (2015).[7] H. Yamamoto, D. Yanagisawa, and K. Nishinari, J. Stat.Mech. , 043204 (2017).[8] K. Christensen and Y. Sasaki, J. Artif. Soc. Soc. Simul. , 9 (2008).[9] N. Pelechano and A. Malkawi, Automat. Constr. , 377(2008).[10] D. Yanagisawa, A. Kimura, A. Tomoeda, R. Nishi,Y. Suma, K. Ohtsuka, and K. Nishinari, Phys. Rev. E , 036110 (2009).[11] T. Ezaki, D. Yanagisawa, K. Ohtsuka, and K. Nishinari,Physica A , 291 (2012).[12] J. E. Rothman, Nature , 55 (1994).[13] M. J. M¨uller, S. Klumpp, and R. Lipowsky, Proc. Natl.Acad. Sci. , 4609 (2008).[14] P. C. Bressloff and J. M. Newby, Rev. Mod. Phys. ,135 (2013).[15] C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin,Biopolymers , 1 (1968).[16] C. T. MacDonald and J. H. Gibbs, Biopolymers , 707(1969).[17] A. Parmeggiani, T. Franosch, and E. Frey, Phys. Rev.Lett. , 086601 (2003).[18] D. Chowdhury, A. Schadschneider, and K. Nishinari,Phys. Life Rev. , 318 (2005).[19] T. Chou, K. Mallick, and R. Zia, Rep. Prog. Phys. ,116601 (2011).[20] C. Appert-Rolland, M. Ebbinghaus, and L. Santen,Phys. Rep. , 1 (2015).[21] D. Helbing, Rev. Mod. Phys. , 1067 (2001).[22] H. Ito and K. Nishinari, Phys. Rev. E , 042813 (2014).[23] S. Tsuzuki, D. Yanagisawa, and K. Nishinari, Phys. Rev.E , 042117 (2018).[24] D. Yanagisawa, A. Tomoeda, R. Jiang, and K. Nishinari,JSIAM Lett. , 61 (2010).[25] C. Arita and D. Yanagisawa, J. Stat. Phys. , 829(2010).[26] C. Arita and A. Schadschneider, Math. Mod Meth. ,401 (2015).[27] J. de Gier and B. Nienhuis,Phys. Rev. E , 4899 (1999).[28] R. A. Blythe and M. R. Evans,J. Phys. A , R333 (2007).[29] B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier,J. Phys. A , 1493 (1993).[30] T. Minemura, K. Nishinari, and A. Schadschnei-der, in International Conference on Cellular Automata (Springer, 2010) pp. 593–599. [31] D. A. Adams, B. Schmittmann, and R. K. P. Zia,J. Stat. Mech. , P06009 (2008).[32] L. J. Cook, R. K. P. Zia, and B. Schmittmann, Phys.Rev. E , 031142 (2009).[33] L. J. Cook and R. K. P. Zia,J. Stat. Mech. , P02012 (2009).[34] L. J. Cook and R. K. P. Zia,J. Stat. Mech. , P05008 (2012).[35] O. Angel, J. Comb. Theory A , 625 (2006).[36] P. A. Ferrari, L. R. Fontes, and Y. Kohayakawa, J. Stat.Phys. , 1153 (1994).[37] C. Arita, J. Stat. Mech. , P12008 (2006).[38] K. Mallick, J. Phys. A , 5375 (1996).[39] B. Derrida and M. Evans, J. Phys. A , 4833 (1999).[40] A. Ayyer, J. L. Lebowitz, and E. R. Speer, J. Stat. Phys. , 1009 (2009).[41] V. Privman, Nonequilibrium statistical mechanics in onedimension (Cambridge University Press, 2005).[42] P. A. Ferrari, J. B. Martin, et al. , Ann. Probab. , 807(2007).[43] L. Cantini, J. Phys. A , 095001 (2008).[44] M. Uchiyama, Chaos Soliton Fract. , 398 (2008).[45] N. Crampe, K. Mallick, E. Ragoucy, and M. Vanicat, J.Phys. A , 175002 (2015).[46] B. Aneva, Eur. Phys. J. C , 403 (2003).[47] C. Arita, A. Kuniba, K. Sakai, and T. Sawabe, J. Phys.A , 345002 (2009).[48] E. Duchi and G. Schaeffer, J. Comb. Theory A , 1(2005).[49] M. Evans, J. Phys. A , 5669 (1997).[50] J. Krug and P. A. Ferrari, J. Phys. A , L465 (1996).[51] K. Mallick, S. Mallick, and N. Rajewsky, J. Phys. A ,8399 (1999).[52] N. Crampe, C. Finn, E. Ragoucy, and M. Vanicat, J.Phys. A , 375201 (2016).[53] C. Arita and K. Mallick, J. Phys. A , 085002 (2013).[54] S. Prolhac, M. R. Evans, and K. Mallick, J. Phys. A ,165004 (2009).[55] A. Ayyer and D. Roy, Sci. Rep. , 13555 (2017).[56] P. F. Arndt, T. Heinzel, and V. Rittenberg, J. Phys. A , L45 (1998).[57] M. Evans, Europhys. Lett. , 13 (1996).[58] T. Sasamoto, Phys. Rev. E , 4980 (2000).[59] M. Bengrine, A. Benyoussef, H. Ez-Zahraouy, J. Krug,M. Loulidi, and F. Mhirech, J. Phys. A , 2527 (1999).[60] E. Duchi and G. Schaeffer, Random Struct. Algor. ,434 (2008).[61] To be precise, if we regard vacant sites as particle 0, swap-ping between particle 0 and other kinds of particles (1,2, 3,......) occurs. However, swapping between particles 1,2, 3,...... does not occur.[62] This system can be simulated also with random updat-ing; however, we adopt parallel updating for the follow-ing three reasons. First, we consider that the results willbe mainly applied to real systems, where particles andagents move simultaneously. It is generally natural toadopt parallel updating in the contexts of pedestrian dy-namics or traffic flows [64, 65]. Second, it is more naturalto consider that time for swapping two particles is com-patible with a one-time step with parallel updating thanwith random-sequential updating. Finally, the simulation times can be suppressed drastically with parallel updat-ing compared to with random updating.[63] M. R. Evans, N. Rajewsky, and E. R. Speer, J. Stat.Phys. , 45 (1999). [64] K. Nagel and M. Schreckenberg, J. de Phys. I , 2221(1992).[65] C. Burstedde, K. Klauck, A. Schadschneider, and J. Zit-tartz, Physica A295 ∀ s ) and p h − α > F , thereby resulting in f ( α ) > . (J8) Appendix K: Specific conditions on α cr , max Here, we discuss the specific conditions on α cr , max .Table VI summarizes explicit expressions for f ( α = α cr ) = 0, where the upper (lower) expression holds inRegion 2 (Region 3). For S = 2, the lower expressionbecomes a quadratic equation in α cr . However, the upperexpression becomes a quartic equation that is too difficultto solve analytically those conditions. Note that for S >