Traffic jams induced by rare switching events in two-lane transport

We investigate a model for driven exclusion processes where internal states are assigned to the particles. The latter account for diverse situations, ranging from spin states in spintronics to parallel lanes in intracellular or vehicular traffic. Introducing a coupling between the internal states by allowing particles to switch from one to another induces an intriguing polarization phenomenon. In a mesoscopic scaling, a rich stationary regime for the density profiles is discovered, with localized domain walls in the density profile of one of the internal states being feasible. We derive the shape of the density profiles as well as resulting phase diagrams analytically by a mean-field approximation and a continuum limit. Continuous as well as discontinuous lines of phase transition emerge, their intersections induce multi-critical behaviour.


Introduction
Non-equilibrium critical phenomena arise in a broad variety of systems, including nonequilibrium growth models [1], percolation-like processes [2], kinetic Ising models [3], diffusion limited chemical reactions [4], and driven diffusive systems [5]. The latter provide models for transport processes ranging from biological systems, like the motion of ribosomes along a m-RNA chain [6] or processive motors walking along cytoskeletal filaments [7,8], to vehicular traffic [9,10]. In this work, we focus on the steady-state properties of such one-dimensional transport models, for which the totally asymmetric simple exclusion process (TASEP) has emerged as a paradigm (for reviews see e.g. [11]- [13]). There, particles move unidirectionally from left to right on a one-dimensional lattice, interacting through on-site exclusion. The entrance/exit rates at the open left/right boundary control the system's behaviour; tuning them, one encounters different non-equilibrium phases for the particle densities [14].
Intense theoretical research has been devoted to the classification of such non-equilibrium phenomena. For example, within the context of reaction-diffusion systems, there is strong evidence that phase transitions from an active to an absorbing state can be characterized in terms of only a few universality classes, the most important being the one of directed percolation (DP) [15]. To search for novel critical behaviour, fruitful results have been obtained by coupling two reaction-diffusion systems [16,17], each undergoing the active to absorbing phase transition. Due to the coupling, the system exhibits a multi-critical point with unusual critical behaviour. Illustration of an exclusion model with two internal states, adopting the language of spin transport. Particles in states ↑ (↓) enter with rates α ↑ (α ↓ ), move unidirectionally to the right within the lattice, may flip at rate ω and leave the system at rates β ↑ (β ↓ ), always respecting Pauli's exclusion principle. 4 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT transport properties have in part been rationalized in terms of an effective single lane TASEP [26]- [28]. There, the case of strong coupling has been investigated: the time-scale of lane switching events is the same as of forward hopping. In our model, we explicitly want to ensure a competition between the boundary processes and the switching between the internal states. We therefore employ a mesoscopic scaling, i.e. we consider the case where the switching events are rare as compared to forward hopping. This is the situation encountered in intracellular traffic [7] where motors nearly exclusively remain on one lane and switch only very rarely. In the context of spin transport, it corresponds to the case where forward hopping occurs much faster than spin precession (weak external magnetic field).
The outline of the present paper is the following. In section 2, we introduce the model in the context of spin transport as well as two-lane traffic. Symmetries and currents are discussed, which play a key role in the following analysis. Section 3 describes in detail the mean-field approximation and the differential equations for the densities obtained therefrom through a continuum limit. The mesoscopic scaling is motivated and introduced, the details of the analytic solution for the spatial density profiles being condensed in appendix A. We obtain the generic form of the density profiles in section 4, and compare our analytic results to stochastic simulations. We find that they agree excellently, suggesting the exactness of our analytic approach in the limit of large systems. As our main result, we encounter the polarization phenomenon, where the density profiles in the stationary non-equilibrium state exhibit localized 'shocks'. Namely, the density of one spin state changes abruptly from low-density (LD) to highdensity (HD). The origin of this phenomenon is rationalized in terms of singularities in coupled differential equations. We partition the full parameter space into three distinct regions, and observe a delocalization transition. The methods to calculate the phase boundaries analytically are developed simultaneously. Section 5 presents details on the stochastic simulations which we have carried out to corroborate our analytic approach. The central result of this work is then addressed in section 6, where two-dimensional analytic phase diagrams are investigated. Our analytic approach identifies the phases where the polarization phenomenon occurs, as well as the continuous and discontinuous transitions that separate the phases. The nature of the transitions is explained by the injection/extraction limited current which is conserved along the track. As a second remarkable feature of the model, we uncover multi-critical points, i.e. points where two lines of phase boundaries intersect or the nature of a phase transition changes from discontinuous to a continuous one. Although multi-critical point are well-known in equilibrium statistical mechanics, a fundamental description for such a behaviour for systems driven far from equilibrium still constitutes a major challenge. A brief summary and outlook concludes this work.

The model
In this section, we describe our model in terms of spin transport as well as two-lane traffic. Though we will preferentially use the language of spins in the subsequent sections, the two-lane interpretation is of no lesser interest, and straightforwardly obtained. Furthermore, we introduce two symmetries which are manifest on the level of the dynamical rules.

Dynamical rules
We consider hopping transport on a one-dimensional lattice, composed of L sites, with open boundaries, see figure 1. Particles possess internal states, which we restrict to two different α II as well as exiting rates, β I resp. β II .
kinds; adopting a spin notation, they are referred to as spin-up (↑) and spin-down (↓). They enter at the left boundary at rates α ↑ resp. α ↓ , and move unidirectionally from left to the right through the lattice. The time-scale is fixed by putting the rate for these hopping events to unity. Within the bulk, particles may also flip their spin state, from spin-up to spin-down and back, at rate ω. Finally, having reached the right boundary, particles may exit the system at rates β ↑ resp. β ↓ , depending on their spin state. We allow all of these processes only under the constraint of Pauli's exclusion principle, meaning that every lattice site may at most be occupied by one particle of a given state. Spin-up and spin-down thus may simultaneously occupy the same site, however two particles with identical spin polarization cannot share a lattice site. In summary, our dynamical rules are the following: i. at site i = 1 (left boundary), particles with spin-up (spin-down) may enter at rate α ↑ (α ↓ ), ii. at site i = L (right boundary), particles with spin-up (spin-down) leave the lattice at rate β ↑ (β ↓ ), iii. particles may hop at unit rate from site i − 1 to the neighbouring site i for i ∈ {2, . . . , L}, i.e. within bulk, iv. within bulk, particles can flip their spin state with rate ω, i.e. spin-up turns into spin-down and vice versa, always respecting Pauli's exclusion principle. Processes (i)-(iii) constitute the TASEP for the two different states separately, while rule (iv) induces a coupling between them. Indeed, when the spin-flip rate ω vanishes, we recover the trivial situation of two independent TASEPs, while we will show that a proper treatment of ω through a mesoscopic scaling induces nontrivial effects.

Two-lane interpretation
Having introduced our model in the language of semi-classical spin transport, where Pauli's exclusion principle is respected while phase coherence completely ignored, we now want to show that it also describes transport with site exclusion on two parallel lanes. As schematically drawn in figure 2, we consider two parallel lanes, each consisting of L sites, labelled as upper lane (I) and lower lane (II). They are identified with the internal states of the particles considered before: a particle with spin-up (spin-down) now corresponds to a particle on lane I (lane II). The processes (i) and (ii) describe entering of particles at lane I (II) at rate α I ≡ α ↑ (α II ≡ α ↓ ) and exiting of lane I (II) at rate β I ≡ β ↑ (β II ≡ β ↓ ). Due to (iii), particles hop unidirectionally to the right on each individual lane; at rate ω, they may switch from lane I to II and back. Pauli's exclusion principle translates into simple site exclusion: all the above processes are allowed under the constraint of admitting at most one particle per site. Again, we clearly observe that it is process (iv) that couples two TASEPs, namely the ones on each individual lane, to each other.

Symmetries
Already on the level of the dynamical rules (i)-(iv) presented above, two symmetries are manifest that will prove helpful in the analysis of the system's behaviour. We refer to the absence of particles with certain state as holes with the opposite respective state 2 . Considering their motion, we observe that the dynamics of the holes is governed by the identical rules (i)-(iv), with 'left' and 'right' interchanged, i.e. with a discrete transformation of sites i ↔ L − i as well as rates α ↑,↓ ↔ β ↓,↑ . The system thus exhibits a particle-hole symmetry. Even more intuitively, the two states behave qualitatively identical. Indeed, the system remains invariant upon changing spin-up to spin-down states and vice versa with a simultaneous interchange of α ↑ ↔ α ↓ and β ↑ ↔ β ↓ , constituting a spin symmetry (in terms of the two-lane interpretation, it translates into a lane symmetry). When analysing the system's behaviour in the five-dimensional phase space, constituted of the entrance and exit rates α ↑,↓ , β ↑,↓ and ω, these symmetries allow to connect different regions in phase space, and along the way to simplify the discussion.

Mean-field equations, currents and the continuum limit
In this section, we shall make use of the dynamical rules introduced above to set up a quantitative description for the densities and currents in the system. Within a mean-field approximation, their time evolution is expressed through one-point functions only, namely the average occupations of a lattice site. Such mean-field approximations have been successfully applied to a variety of driven diffusive systems, see e.g. [12]. We focus on the properties of the non-equilibrium steady state, which results from boundary processes (entering and exiting events) as well as bulk ones (hopping and spin-flip events). Both types of processes compete if their time-scales are comparable; we ensure this condition by introducing a mesoscopic scaling for the spin flip rate ω. Our focus is on the limit of large system sizes L, which is expected to single out distinct phases. To solve the resulting equations for the densities and currents, a continuum limit is then justified, and it suffices to consider the leading order in the small parameter, namely the ratio of the lattice constant to system size. Such a mesoscopic scaling has been already successfully used in [29,30] in the context of TASEP coupled to Langmuir dynamics.

Mean field approximation and currents
Let n ↑ i (t) resp. n ↓ i (t) be the fluctuating occupation number of site i for spin-up resp. spin-down state, i.e. n ↑,↓ i (t) = 1 if this site is occupied at time t by a particle with the specified spin state and n ↑,↓ i (t) = 0 otherwise. Performing ensemble averages, the expected occupation, denoted by ρ ↑ i (t) 7

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
and ρ ↓ i (t), is obtained. Within a mean-field approximation, higher order correlations between the occupation numbers are neglected, i.e. we impose the factorization approximation Equations of motion for the densities can by obtained via balance equations. The timechange of the density at a certain site is related to appropriate currents. The spatially varying spin current j ↑ i (t) quantifies the rate at which particles of spin state ↑ at site i − 1 hop to the neighbouring site i. Within the mean-field approximation, equation (1), the current is expressed in terms of densities as and similarly for the current j ↓ i (t). The sum yields the total particle current . Due to the spin-flip process (iv), there also exists a leakage current j ↑↓ i (t) from spin-up state to spin-down state. Within mean-field and similarly for the leakage current j ↓↑ i (t) from spin-down to spin-up state. Now, for i ∈ {2, . . . , L − 1} we can use balance equations to obtain the time evolution of the densities, This constitutes an exact relation. Together with the mean field approximation for the currents, equations (2) and (3), one obtains a set of closed equations for the local densities d dt At the boundaries of the track, the corresponding expressions involve also the entrance and exit events, which are again treated in the spirit of a mean-field approach Due to the spin symmetry, i.e. interchanging ↑ and ↓, an analogous set of equations holds for the time evolution of the density of particles with spin-down state.
In the stationary state, the densities ρ ↑(↓) i (t) do not depend on time t, such that the time derivatives in equations (5)-(7) vanish. Therefrom, we immediately derive the spatial conservation of the particle current: indeed, summing equation (4) with the corresponding equation for the density of spin-down states yields

8
DEUTSCHE PHYSIKALISCHE GESELLSCHAFT such that the particle current does not depend on the spatial position i. Note that this does not apply to the individual spin currents, they do have a spatial dependence arising from the leakage currents.
In a qualitative discussion, let us now anticipate the effects that arise from the non-conserved individual spin currents as well as from the conserved particle current. The latter has its analogy in TASEP, where the particle current is spatially conserved as well. It leads to two distinct regions in the parameter space: one where the current is determined by the left boundary, and the other where it is controlled by the right one. Both regions are connected by the discrete particle-hole symmetry. Thus, in general, discontinuous phase transitions arise when crossing the border from one region to the other. In our model, we will find similar behaviour: the particle current is either determined by the left or by the right boundary. Again, both regions are connected by the discrete particle-hole symmetry, such that we expect discontinuous phase transitions at the border between both. Except for a small, particular region in the parameter space, this behaviour is captured quantitatively by the mean-field approach and the subsequent analysis, which is further corroborated by stochastic simulations. The phenomena linked to the particular region will be presented elsewhere [31].
On the other hand, the non-conserved spin currents may be compared to the current in TASEP coupled to Langmuir kinetics, see [29,30]. Due to attachment and detachment processes, the in-lane current is only weakly conserved, allowing for a novel phenomenon, namely phase separation into a LD and a HD region separated by a localized domain wall. The transitions to this phase are continuous considering the domain wall position x w as the order parameter. In our model, an analogous but even more intriguing phase will appear as well, with continuous transitions being possible.

Mesoscopic scaling and the continuum limit
3.2.1. Mesoscopic scaling. Phases and corresponding phase transitions are expected to emerge in the limit of large system size, L → ∞, which therefore constitutes the focus of this work. We expect interesting phase behaviour to arise from the coupling of spin-up and spin-down states via spin-flip events, in addition to the entrance and exit processes. Clearly, if spin-flips occur on a fast time-scale, comparable to the hopping events, the spin degree of freedom is relaxed, such that the system's behaviour is effectively the one of a TASEP. Previous work on related two-lane models [27,26] focused on the physics in that situation. In this work, we want to highlight the dynamical regime where coupling through spin-flips is present, however not sufficiently strong to relax the system's internal degree of freedom. In other words, we consider physical situations where spin-flips occur on the same time-scale as the entrance/exit processes. Defining the gross spin-flip rate = ωL yields a measure of how often a particle flips its spin state while traversing the system. To ensure competition between spin-flips with boundary processes, a mesoscopic scaling of the rate ω is employed by keeping fixed, of the same order as the entrance/exit rates, when the number of lattice sites becomes large L → ∞.

Continuum limit and first order approximation.
The total length of the lattice will be fixed to unity and one may define consistently the lattice constant = 1/L. In the limit of large systems → 0, a continuum limit is anticipated. We introduce continuous functions first order in the lattice constant, the difference equations (5)-(7) turn into differential equations. Observing that ω = is already of order , we find that the zeroth order of equation (5) vanishes, and the first order in yields Similarly, the same manipulations for ρ ↓ yield The expansion of equations (6) and (7) in powers of , yields in zeroth order which impose boundary conditions. Since two boundary conditions are enough to specify a solution of the coupled first-order differential equations, the system is apparently overdetermined. Of course, the full analytic solution, i.e. where all orders in are incorporated, will be only piecewise given by the first-order approximation, equations (10)- (12). Between these branches, the solution will depend on higher-orders of , therefore, these intermediate regions scale with order and higher. They vanish in the limit of large systems, → 0, yielding domain walls or boundary layers. Let us explain the latter terms. At the position of a domain wall, situated in bulk, the density changes its value discontinuously, from one of a LD region to one of a HD. Boundary layers are pinned to the boundaries of the system. There as well, the density changes discontinuously: from a value that is given by the corresponding boundary condition to that of a LD or HD region which is imposed by the opposite boundary.

Symmetries and currents revisited.
In the following, we reflect important properties of the system, symmetries and currents, on the level of the first-order approximation, equations (10)- (12). The explicit solution of the latter can be found in appendix A.
The particle-hole symmetry, already inferred from the dynamical rules, now takes the form Interchanging ↑ and ↓ in the densities as well as the in and outgoing rates yields the spin symmetry, The individual spin currents as well as the particle current have been anticipated to provide further understanding of the system's behaviour. In the continuum limit the zeroth order of the spin currents is found to be The terms on the right-hand side, arising from the spin-flip process (iv), are seen to violate the spatial conservation of the spin currents. However, due to the mesoscopic scaling of the spin flip rate ω, the leakage currents between the spin states are only weak, see equation (3), locally tending to zero when → 0, such that the spin currents vary continuously in space. This finding imposes a condition for the transition from one branch of first-order solution to another, as described above: such a transition is only allowed when the corresponding spin currents are continuous at the transition point, thus singling out distinct positions for a possible transition.
Finally, summing the two equations in equation (15) yields the spatial conservation of the particle current: ∂ x J = 0.

Partition of the parameter space and the generic density behaviour
The parameter space of our model, spanned by the five rates α ↑,↓ , β ↑,↓ , and , is of high dimensionality. However, in this section, we show that it can be decomposed into only three basic distinct regions: the maximal-current (MC) region as well as the injection-limited (IN) and the extraction-limited one (EX). While trivial phase behaviour occurs in the MC region, our focus is on the IN and EX region (connected by particle-hole symmetry), where a striking polarization phenomenon occurs. The generic phase behaviour in these regions is derived, exhibiting this effect.

Effective rates
The entrance and exit rates as well as the carrying capacity of the bulk impose restrictions on the particle current. For example, the capacity of the bulk limits the individual spin currents j ↑(↓) to maximal values of 1/4. The latter occurs at a density of 1/2, as seen from the previous result To illustrate the influence of the injection and extraction rates, we first consider an 'open' right boundary i.e. β ↑ = β ↓ = 1. Particles then leave the system unhindered, such that only the entrance rates may limit the particle current. Provided one of these rates, say α ↑ , exceeds the value 1/2, the current of the corresponding state (↑) is limited by the capacity of the bulk to a value of 1/4 in the vicinity of the left boundary. A boundary layer thus forms in the density profile of spin-up state at the left boundary, connecting the value of the injection rate α ↑ to the value 1/2. Up to this boundary layer, the density profile ρ ↑ (x) is identical to the one where α ↑ takes a value of 1/2, cf figure 3. Similar reasoning holds for the extraction rates β ↑(↓) . They as well behave effectively as 1/2 when exceeding this value. To treat these findings properly, we introduce the effective rates that only the capacity of the bulk and the entrance rates limit the spin currents. The injection rate α ↑ > 1 2 effectively acts as 1 2 . The analytic predictions correspond to the solid lines, the results from stochastic simulations for L = 10 000 are indicated by the wiggly line. With increasing spatial position, the densities approach a common value ρ e . The parameters used are α ↑ = 0.7, α ↓ = 0.15 and = 0.5.
The system's bulk behaviour will only depend on them, and, in particular, remain unaffected when a rate is varied at values exceeding 1/2.

IN, EX, and MC region
Equipped with these results, in the case of an 'open' right boundary, the spin currents in the vicinity of the left boundary are given by resulting in a particle current J IN imposed by the injection rates: . The analogous relations, with the injection and extraction rates interchanged, hold for the case of an 'open' left boundary, α ↑ = α ↓ = 1. The particle current is then controlled by the right boundary: In general, depending on which imposes the stronger restriction, either the left or the right boundary limits the particle current: J min(J IN , J EX ). Indeed, J = min(J IN , J EX ) holds except for an anomalous situation, where the current is lower than this value 3 . Depending on which of both cases applies, two complementary regions in phase space are distinguished: Since they are connected by discrete particle-hole symmetry, we expect discontinuous phase transitions across the border between both, to be referred as IN-EX boundary.
Right at the IN-EX boundary, the system exhibits coexistence of LD and HD phases, separated by domain walls. Interestingly, this phase coexistence emerges on both lanes (states), which may be seen as follows. Recall that a domain wall concatenates a region of low and another of HD. However, while the densities exhibit a discontinuity, the spin currents must be continuous. In other words, the spin currents, and therefore the particle currents, it turns out that there, they do indeed form, and are delocalized. We refer to our forthcoming publication [31] for a detailed discussion of this phenomenon. Away from the IN-EX boundary, it follows that at most on one lane (state) a domain wall may appear. When both entrance rates α ↑ , α ↓ as well as both exit rates β ↑ , β ↓ exceed the value 1/2, the particle current is limited by neither boundary, but only through the carrying capacity of the bulk, restricting it to twice the maximal value 1/4 of the individual spin currents: J = 1/2. The latter situation therefore constitutes the maximal current region.

The generic state of the densities
As we have seen in the previous section, particularly simple density profiles emerge in the MC region. There, up to boundary layers, the density profiles remain constant at a value 1/2 for each spin state. Another special region in parameter space is the IN-EX boundary, characterized by the simultaneous presence of domain walls in both spin states, as we discuss elsewhere [31].
Away from these regions, the generic situation for the density profiles is illustrated in figure 4. Here, we have considered parameters belonging to the IN region; the behaviour in the EX region follows from particle-hole symmetry. A domain wall emerges for one spin state and a boundary layer for the other one. For specificity, we consider a domain wall for the spin-up state, the other situation is obtained from spin symmetry. The density profiles ρ ↑(↓) l close to the left boundary are given by the solution of the first-order differential equations (10) and (11) IN region). Therefore, the densities satisfy right boundary conditions which are given by ρ ↑ r (x = 1) = 1 − β ↑ eff ; and ρ ↓ r (x = 1) is found from the conservation of the particle current: At some point x w in bulk, the left and right solutions have to be concatenated by a domain wall for spin-up. To determine the position x w of this domain wall, we use the continuity of the spin currents; see figure 4(b) 4 . This continuity condition singles out a distinct spatial position for the domain wall: denote by ρ ↑ l (x w ) the value of the density to the left of x w , and ρ ↑ r (x w ) the value to the right.
for the domain wall position 5 . From the conservation of the particle current J, it follows that the density ρ ↓ is continuous at the position x w . When considering the internal states as actual spins, the appearance of a domain wall in the density profile of one of the spin states results in a spontaneous polarization phenomenon. Indeed, while both the density of spin-up and spin-down remain at comparable low values in the vicinity of the left boundary, this situation changes upon crossing the point x w . There, the density of spin-up jumps to a high value, while the density of spin-down remains at a low value, resulting in a polarization in this region.
Comparing the generic phase behaviour to the one of TASEP, we observe that the IN region can be seen as the analogue to the LD region there: within both, a LD phase accompanied by a boundary layer at the right boundary arises. Following these lines, the EX region has its analogue in the HD region, while the MC region is straightforwardly generalized from the one of TASEP. Furthermore, the delocalization transition across the IN-EX boundary is similar to the appearance of a delocalized domain wall at the coexistence line in TASEP.

Phases and phase boundaries
In the generic situation of figure 4, the density of spin-down is in a homogeneous LD state, while for spin-up, a LD and a HD region coexist. We refer to the latter as the LD-HD IN phase, as the phase separation arises within the IN region, to be contrasted from a LD-HD EX phase which may arise within the EX region. Clearly, the LD-HD IN phase is only present if the position x w of the domain wall lies within bulk. Tuning the system's parameter, it may leave the system through the left or right boundary, resulting in a homogeneous phase. Indeed, x w = 1 marks the transition between the LD-HD IN phase and the pure LD state, while at x w = 0 the density changes from the LD-HD IN to a homogeneous HD state. Regarding the domain wall position x w as an order parameter, these transitions are continuous. Implicit analytic expressions for these phase boundaries, derived in the following, are obtained from the first-order approximation, equations (10) and (11).
Spin symmetry yields the analogous situation with a domain wall appearing in the density profile of spin-down, while particle-hole symmetry maps it to the EX region, where a pure HD phase arises for one of the spins. Discontinuous transitions accompanied by delocalized domain walls appear at the submanifold of the IN-EX boundary (see [31] for a detailed discussion).
The phase boundaries may be computed from the condition x w = 0 and x w = 1 in the situation of figure 4. Consider first the case of x w = 0. There, the density profiles are fully given by the first-order approximation ρ ↑(↓) r satisfying the boundary conditions at the right. The condition (18) translates to which yields an additional constraint Again, the latter is a constraint on the parameters and defines the hyper-surface in the IN region where x w = 1 is found, being the phase boundary between the LD-HD IN and the homogeneous LD phase. The conditions (19) and (20) yield implicit equations for the phase boundaries. The phase diagram is thus determined up to solving algebraic equations, which may be achieved numerically. Further insight concerning the phase boundaries is possible and may be obtained analytically, which we discuss next.
Firstly, we note that in the case of equal injection rates, α ↑ = α ↓ , the density profiles in the vicinity of the left boundary are constant. If in addition α ↑ = α ↓ = β ↑ < 1/2, we observe from equation (20) that a domain wall at x w = 1 emerges. Therefore, this set of parameters always lies on the phase boundary x w = 1, independent of the value of .
Secondly, we investigate the phase boundary determined by x w = 0. Comparing with figure 4, we observe that the first-order approximation ρ ↑ r for the density of spin-up may reach the value 1 2 at a point which is denoted by x 1/2 : ρ ↑ r (x 1/2 ) = 1 2 . This point corresponds to a branching point of the first-order solution. Increasing , the value of x 1/2 increases as well. The domain wall in the density of spin-up can only emerge at a value x w x 1/2 . At most, x w = x 1/2 , in which case a domain wall with infinitesimal small height arises. For the phase boundary specified by x w = 0, this implies that it only exists as long as x 1/2 0. The case x w = x 1/2 = 0 corresponds to a domain wall of infinitesimal height, which is only feasible if α ↑ eff = 1 2 . Now, for given rates α ↑ eff = 1 2 , α ↓ , β ↑ , the condition x 1/2 = 0 yields a critical rate * (α ↓ , β ↑ ), depending on the rates α ↓ , β ↑ . The situation x w = 0 can only emerge for rates * (α ↓ , β ↑ ). Varying the rates α ↑ , α ↓ and β ↑ , the critical rate * (α ↓ , β ↑ ) changes as well. In appendix A, we show that its largest value occurs at α ↓ = β ↑ = 0. They yield the rate C ≡ * (α ↓ = β ↑ = 0), which is calculated to be The critical * (α ↓ , β ↑ ) are lying in the interval between 0 and C : * (α ↓ , β ↑ ) ∈ [0, C ], and all values in this interval in fact occur. The rate C defines a scale in the spin-flip rate : For C , the phase boundary determined by x w = 0 exists, while disappearing for > C . Thirdly, we study the form of the phase boundaries for large , meaning C . In this case, the phase boundary specified by x w = 0 is no longer present. Furthermore, it turns out that in this situation, the densities close to the left boundary quickly approximate a common value ρ e . The latter is found from conservation of the particle current: 2ρ e (1 − ρ e ) = J. We now consider the implications for the phase boundary determined by This condition specifies the phase boundary x w = 1, asymptotically for large . It constitutes a simple quadratic equation in the in and outgoing rates, independent of β ↓ , and contains the set

Stochastic simulations
To confirm our analytic findings from the previous section, we have performed stochastic simulations. The dynamical rules (i)-(iv) described in section 2.1 were implemented using random sequential updating. In our simulations, we have performed averages over typically 10 5 time steps, with 10 × L steps of updating between successive ones. Finite size scaling singles out the analytic solution in the limit of large system sizes, as exemplified in figures 3 and 4. For all simulations, we have checked that the analytic predictions are recovered upon approaching the mesoscopic limit. We attribute the apparent exactness of our analytic approach in part to the exact current density relation in the steady state of the TASEP [34]. The additional coupling of the two TASEPs in our model is only weak: the local exchange between the two states vanishes in the limit of large system sizes. Correlations between them are washed out, and mean-field is recovered.
The observed exactness of the analytic density profiles within the mesoscopic limit implies that our analytic approach yields exact phase diagrams as well. The latter are the subject of the subsequent section.

Two-dimensional phase diagrams
In this section, we discuss the phase behaviour on two-dimensional cuts in the whole fivedimensional parameter space. Already the simplified situation of equal injection rates, α ↑ = α ↓ , yields interesting behaviour. There as well as in the general case, we investigate the role of the spin-flip rate by discussing the situation of small and large values of .

Equal injection rates
For simplicity, we start our discussion of the phase diagram with equal injection rates, α ↑ = α ↓ . Then, the spin polarization phenomenon, depicted in figure 4, becomes even more striking. Starting from equal densities at the left boundary, and hence zero polarization, spin polarization suddenly switches on at the domain wall position x w . The particular location of x w is not triggered by a cue on the track, but tuned through the model parameters.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
The phase transitions from LD to the LD-HD IN arising in the IN region take a remarkably simple form. Their location is found from x w = 1, and is determined by equation (20) (if phase coexistence arises for spin-up). Since ρ ↑ (x) = ρ ↓ (x) = α = constant for x < x w , equation (20) turns into α = β ↑ . The latter transition line intersects the IN-EX boundary, given by J IN = J EX , at β ↑ = β ↓ = α, i.e. at the point where all entrance and exit rates coincide. At this multi-critical point A, a continuous line intersects a discontinuous one. The same transition in the density of spin-down state is, from similar arguments, located at α = β ↓ , and also coincides with the IN-EX boundary in A. Neither the multi-critical point A nor these phase boundaries depend on the magnitude of the gross spin flip rate . Therefore, qualitatively tuning the system's state is possible only upon changing the injection or extraction rates. The other phase transitions within the IN region, namely from the HD to the LD-HD IN phase, are more involved. The analytic solution given by (A.12) and (A.13) has to be considered together with the condition (19) for the transition. However, at the end of section 4.4, we have found that these transitions (determined by x w = 0) disappear for sufficiently large > C . For larger values of α, a localized domain wall emerges for spin-up (implying a LD-HD EX phase), and a pure HD phase for spin-down. If α is further increased, the domain wall in the spin-up density profile leaves the system through the left boundary (at x w = 0), and pure HD phases remain for both spin states.
While we have found the transitions within the IN region by simple expressions in the previous subsection, the ones emerging in the EX region are more complex and involve the full analytic solutions (A.12) and (A.13). Their most notable feature is that the width of the corresponding coexistence phase decreases with increasing spin-flip rate , until it finally vanishes in the limit → ∞. This may be seen by considering the analogue of equation (22) in the EX region, which describes the phase boundary as it is asymptotically approached when → ∞: it coincides with the IN-EX boundary. occurs in the EX region. In figure 6, we show resulting phase diagrams for the spin-up (left panel) and spin-down (right panel), resp. The additional transition lines intersect the IN-EX boundary (bold) at additional multi-critical points B IN and B EX . Also, they partly substitute the IN-EX boundary as a phase boundary: across some parts of the latter, phase transitions do not arise. This behaviour reflects the decoupling of the two states for decreasing spin-flip rate . Indeed, for → 0, the states become more and more decoupled, such that the IN-EX boundary, involving the combined entrance and exit rates of both states, loses its significance.

Multi-critical points.
Although the shapes of most of the transition lines appearing in the phase diagrams shown in figure 6 are quite involved, they also exhibit simple behaviour. Pairwise, namely one line from a transition in spin-up and another from a related transition in spin-down states, they intersect the IN-EX boundary in the same multi-critical point. This

The general case
Having focused on the physically particularly enlightening case of equal entering rates in the previous subsection, we now turn to the general case. To illustrate our findings, we show phase diagrams depending on the injection and extraction rates for spin-up states, α ↑ and β ↑ . Similar behaviour as for equal entrance rates is observed. The multi-critical point A now splits up into two distinct points A IN and A EX . The transition from LD to the LD-HD IN phase in the IN region asymptotically takes the form of equation (22), and the one from HD to the LD-HD EX phase in the EX region is obtained by particle-hole symmetry. All phase boundaries, including the IN-EX boundary, are thus given by simple quadratic expressions.
Phase diagrams with different topologies that can emerge are exhibited in figures 7 and 8. As in the previous subsection, we show the phases of spin-up (spin-down) states on the left (right) panels. The phase boundaries between the LD and the LD-HD IN  We now discuss the influence of the spin-flip rate on the continuous transition lines for spin-up. In section 4.4 the manifold defined by α ↑ = β ↑ = α ↓ eff was found to be a submanifold of the phase boundary specified by x w = 1 in the IN region. Independent of , the point α ↑ = β ↑ = α ↓ eff , denoted by N IN , thus lies on the boundary between the LD and the LD-HD IN phase (determined by x w = 1). For large , this boundary approaches the one given by equation (22).
Regarding the transition from the HD to the LD-HD IN within the IN region (determined by x w = 0), section 4.4 revealed that for increasing it leaves the IN region at a critical transfer rate * (α ↓ , β ↑ ). In the limit → 0, the densities ρ ↑ (x) and ρ ↓ (x) approach constant values, and both the curve x w = 1 as x w = 0 for spin-up in the IN region approach the line β ↑ = α ↑ for α ↑ 1 2 . The phase in the upper right quadrant in the phase diagram converges to the MC phase, such that in this limit, the case of two uncoupled TASEPs is recovered.

Conclusions
We have presented a detailed study of an exclusion process with internal states recently introduced in [19]. The TASEP has been generalized by assigning two internal states to the particles. Pauli's exclusion principle allows double occupation only for particles in different internal states. Occasional switches from one internal state to the other induce a coupling between the transport processes of the separate states. Such a dynamics encompasses diverse situations, ranging from vehicular traffic on multiple lanes to molecular motors walking on intracellular tracks and future spintronics devices.
We have elaborated on the properties of the emerging non-equilibrium steady state focusing on density and current profiles. In a mesoscopic scaling of the switching rate between the internal states, nontrivial phenomena emerge. A localized domain wall in the density profile of one of the internal states induces a spontaneous polarization effect when viewing the internal states as spins. We provide an explanation based on the weakly conserved currents of the individual states and the current-density relations. A quantitative analytic description within a mean-field approximation and a continuum limit has been developed and solutions for the density and current profiles have been presented. A comparison with stochastic simulations revealed that our analytic approach becomes exact in the limit of large system sizes. We have attributed this remarkable finding to the exact current-density relation in the TASEP, supplemented by the locally weak coupling of the two TASEPs appearing in our model: ω → 0 in the limit of large system sizes. Local correlations between the two internal states are thus obliterated, as particles hop forward on a much faster time-scale than they switch their internal state.
Furthermore, the parameter regions that allow for the formation of a localized domain wall have been considered. Analytic phase diagrams for various scenarios, in particular the case of equal entrance rates, have been derived. The phase diagrams have been found to exhibit a rich structure, with continuous as well as discontinuous non-equilibrium phase transitions. The discontinuous one originates in the conserved particle current, which is either limited by injection or extraction of particles. At the discontinuous transition between both regimes, delocalized domain walls emerge in the density profiles of both internal states. Multi-critical points appear at the intersections of different transition lines organizing the topology of the phase diagrams. Two classes of multi-critical points are identified, one of them arises only for sufficiently small gross spin-flip rate < C . The value C , calculated analytically, provides a natural scale for the rate .
It would be of interest to see which of the described phenomena qualitatively remain when generalizing the model to include more than two internal states. Indeed, within the context of molecular motors walking on microtubuli [7], between 12 and 14 parallel lanes are relevant. Also, the internal states might differ in the sense of different switching rates from one to another [28] and the built-in asymmetry may result in different phases. In the context of intracellular transport it appears worthwhile to investigate the consequences of a coupling to a bulk reservoir, cf [29,30,35]; in particular, to study the interplay of domain wall formation induced by attachment and detachment processes as well as rare switching events.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
Summing them we find constitutes a first integral. Remember that such that J is given by the total current: This equation suggests the following parameterization: The derivative reads which leads to the differential equation and I is an constant of integration.