Marco Schiró

Aditi Mitra

Transport Across an Impurity in One Dimensional Quantum Liquids Far From Equilibrium

We study the effect of a single impurity on the transport properties of a one dimensional quantum liquid highly excited away from its ground state by a sudden quench of the bulk interaction. In particular we compute the time dependent dc current to leading order in the impurity potential, using bosonization, and starting both from the limit of a uniform system, and from the limit of two decoupled semi-infinite systems. Our results reveal that the nonequilibrium excitation of bulk modes induced by the global quench has important effects on the conductor-insulator quantum phase transition, turning it into a crossover for any small quench amplitude, and destroying the exact duality between the conducting and insulating fixed points. In addition, the current displays a faster decay towards the steady-state as compared to the equilibrium case, a signature of quenchinduced decoherence.

I. INTRODUCTION

Transport phenomena in strongly correlated quantum systems are typically sensitive to interactions, inhomogeneities and low-dimensionality and often reveal intriguing and unexpected effects. A prominent example is the impact of a single potential barrier in a one dimensional interacting electron liquid. While the clean system would display, despite interactions, ideal conductance quantization in the presence of metallic leads [1][2][3][4] , the impurity at zero temperature is capable of turning the perfect conductor into an insulating link, giving rise to power law corrections to transport coefficients at finite temperature or finite voltage bias [5][6][7] . The physics of a single impurity in a Luttinger liquid has attracted enormous attention in the past decades [8][9][10] and still finds beautiful applications in systems such as the edges of a quantum spin Hall insulator 11,12 or in electronic quantum circuits [13][14][15] . Traditional condensed matter transport settings involve a system initially in thermal equilibrium which is perturbed away from it by the action of some external field conjugate to a conserved current. Within linear response theory one then obtains transport coefficients, which give information on the structure of the underlying groundstate and its low-energy excitations.

Recently, experimental advances in controlling and probing ultra-cold atomic gases has offered a new platform to study nonequilibrium time-dependent phenomena in a fully tunable setting 16,17 . First generation of experiments studied the dynamics induced by rapidly changing in time some system parameter [18][19][20] in an otherwise isolated system, a so called quantum quench. More recently, the experimental focus shifted towards realizing genuine transport experiments with cold atoms 21,22 . It is important to realize that these systems are very well isolated from their environment and hence intrinsically out of equilibrium, therefore the standard condensed matter idealizations do not necessarily apply.

Motivated by these experimental developments, in this paper we investigate the transport properties of nonequilibrium quantum many body states which are thermally isolated and highly excited above their ground state. While for generic, non-integrable and ergodic, quantum many-body systems one may expect the excitation energy to be effectively converted into temperature at long times 23,24 , and hence to induce decoherence, the situation is less clear at intermediate time scales where long-lived prethermal states may emerge [24][25][26][27][28][29][30][31] , especially parametrically close to an integrable point. The idea we pursue in this work is to use transport as a probe to unveil the structure of these prethermal excited states, to understand the relevant excitations and whether non-trivial quantum phenomena survive at these high-energies, and if so, in which form. It should be noted that there has been a recent spur of interest towards understanding similar dynamical quantum correlations in isolated many body systems excited after quantum quenches [32][33][34][35] . We have considered an example along these lines in a recent work on a time-dependent orthogonality catastrophe problem 36 where we introduced a novel dynamical Loschmidt echo, encoding the response of a highly excited state to a local perturbation. Here we study a similar problem from the point of view of transport. Specifically, in this paper we compute the timedependent current of a one dimensional system of spinless electrons described by the Luttinger model, which is excited by a sudden change of the two-particle interaction and a simultaneous switching-on of a local scattering potential. These two perturbations have very different effects, the former injecting extensive energy into the bulk modes, the latter creating a non-linear channel for local dissipation. We discuss the physics both in the limit of a uniform system, where the local potential induces a weak backscattering term, and in the limit of two decoupled semi-infinite systems where we switch on a local tunneling. In both cases, using bosonization, we formulate the problem in terms of a boundary sine-Gordon model in a nonequilibrium transient environment. Using perturbative approaches we discuss the fate of the conducting-insulating zero temperature quantum phase transition in the presence of a nonequilibrium bulk excitation, thus complementing and expanding previous results 37,38 .

The results for the steady state current reveal the emergence of a novel energy scale associated with a quench-induced decoherence effect, which turns the sharp equilibrium transition into a smooth crossover. Interestingly, the very same energy scale has been shown to play a key role in the transient orthogonality catastrophe problem, cutting off the renormalization group flow of the backscattering and resulting in an exponential decay of the Loschmidt echo 36 .

Our analysis, although perturbative in nature, ultimately suggests that the steady state impurity problem in a quenched Luttinger model has very different properties than in equilibrium at zero temperature. In particular the quench acts as a relevant perturbation on both sides of the equilibrium transition, even when the local non-linearity would be irrelevant at the tree level, generating an effective decoherence which drives the system away from the uniform and open-chain fixed points. While this behavior is reminiscent of an effective temperature, we will see that this analogy is only qualitative rather than quantitative, as the scaling of physical quantities in the steady state with respect to this emergent energy scale do not show signature of non-trivial powerlaws like the ones encountered at non-zero temperatures. At the same time we find that deviations from the asymptotic steady state regime, due to transient effects at finite time, display nonequilibrium power laws with characteristic Luttinger liquid exponents which may be a signature of a sort of Luttinger liquid universality out of equilibrium, as recently discussed 29,38,39 .

Note that the problem of an impurity in a quenched Luttinger model has been recently addressed in Refs. 37 and 38. In the former, the dynamics of the current in the tunneling regime was investigated after a global quench of the bulk interaction and a faster decay to the stationary state was found. In the latter, results for the temperature scaling of the conductance in the steady state after the quench have been obtained using different approaches such as bosonization and functional RG. Although we do not address temperature dependence of transport coefficients in this paper, and we work always at zero temperature, it is useful to comment on the relation between the results of Ref. 38 and the physical picture emerging from this work. We will do this in the discussion section.

Finally, it is worth mentioning that the question of nonequilibrium effects on local quantum criticality has been addressed recently also in the context of driven quantum systems 40 and specifically for single impurity in Luttinger liquids in the context of noisy driven Josephson junction circuits 41 or tunneling in biased quantum wires 42 using the out-of-equilibrium bosonization framework 43,44 . The picture emerging in these cases, where the nonequilibrium perturbation induces a decoherence mechanism which eventually cuts the quantum critical power-law behavior, is consistent with our results for the isolated quenched problem.

The paper is organized as follows. In section II we introduce the system, the nonequilibrium protocol, and a derivation of an expression for the time dependent current in terms of a bosonic Green's function. We then evaluate the current in the weak-backscattering limit in section II A. Such a perturbative approach eventually breaks down in a certain region of parameter space and this will lead us in section III to formulate the problem in the limit of two disconnected 1D systems coupled by a local tunneling term. We will compute perturbatively the tunneling current in section III A and thus obtain a complete picture of the problem.

II. TRANSPORT IN THE WEAK-BACKSCATTERING REGIME

We start by discussing the model and the nonequilibrium set up. We consider a one dimensional system of interacting spinless fermions described by the Tomonaga-Luttinger (TL) model. We assume the system to be initially (at time t ≤ 0) in the ground-state |Ψ 0 of the TL Hamiltonian

H 0 = H free + g 0 4 2 dx ρ 2 L (x) + ρ 2 R (x) +g 0 2 dxρ L (x)ρ R (x) (1)
where

H free = v F α=L,R s α dxψ † α i∂ x ψ α with s L = -s R = 1, ρ α = ψ †
α ψ α is the fermion density with given chirality α = L/R, and g 0 4 (g 0 2 ) are the strengths of the intra (inter)-branch scattering processes.

At time t = 0 we assume a sudden change (quench) in the value of the bulk interactions to g 2 , g 4 , so that the system will evolve with the Hamiltonian

H = H free + g 4 2 dx ρ 2 L (x) + ρ 2 R (x) +g 2 dxρ L (x)ρ R (x) (2)
This global quench injects extensive energy into the system and triggers a nonequilibrium occupation of the bulk modes. Dynamics arising simply due to this quench has been discussed extensively in the literature [45][46][47] . In addition to the bulk excitation, we assume that a static impurity potential V loc has also been suddenly switched on at time t = 0. This static impurity can induce an intra-branch (R, L) → (R, L) as well as inter-branch (R, L) → (L, R) scattering of fermions, that we write as the sum of two contributions

V loc = V fs ψ † L (0)ψ L (0) + L ↔ R + (3) +V bs ψ † L (0)ψ R (0) + h.c. . (4)
the first term representing the forward scattering while the latter the backward scattering. As a result the wave function at time t is |Ψ(t) = exp (-iH + t) |Ψ(0) , with

H + = H + V loc .
We employ bosonization to describe the system, thus introducing the bosonic fields φ(x) and θ(x) describing collective density and current excitations respectively,

ψ L/R (x) ∼ 1 √ 2πα e i(θ(x)±φ(x)) (5) ρ L/R (x) = -∂ x (φ(x) ± θ(x)) /2π (6)
The Hamiltonians before and after the quench in terms of these bosonic fields are

H 0 = u 0 2π dx K 0 (∂ x θ(x)) 2 + 1 K 0 (∂ x φ(x)) 2 (7)
and

H + = u 2π dx K (∂ x θ(x)) 2 + 1 K (∂ x φ(x)) 2 +V loc (8)
Within bosonization, the impurity potential in Eq. ( 3) can be written as

V loc = g fs ∂ x φ(0) + g bs cos 2φ(0) (9)
where the effective impurity couplings read respectively 48

g fs = - V fs π g bs = V bs πα (10)
The Luttinger parameter K 0 , K and the sound velocities u 0 , u are related to the Fermi velocity and the interaction parameters as

u 0 = v F (1 + g 0 4 /2πv F ) 2 -(g 0 2 /2πv F ) 2 (11)
K 0 = 1 + g 0 4 /2πv F -g 0 2 /2πv F 1 + g 0 4 /2πv F + g 0 2 /2πv F (12)
with similar relations holding for u, K as a function of g 2 , g 4 . In order to preserve Galilean invariance we choose the values of g 2 , g 4 such that u 0 K 0 = uK. This amounts to simply requiring that g 4 -g 2 = g 0 4 -g 0 2 . In order to probe transport through the system we will study linear response to a weak electric field E(x, t) applied after the quench,

E(x, t) = -∂ t A(x, t) (13)
with A(x, t) being the vector potential. In the presence of the electric field, the LL Hamiltonian is modified according to the minimal substitution which amounts to the shift

∂ x θ(x) → ∂ x θ(x) -eA(x, t) (14)
The current operator can be obtained as the functional derivative of H(A) with respect to the vector potential

J(x, t) = - δH δA(x, t) (15)
As usual, the current has two contributions, J(x, t) = J d (x, t) + J p (x, t), the diamagnetic (J d ) and the paramagnetic (J p ) one. The former is given by (hereafter = 1)

J d (x, t) = - e 2 uK π A(x, t) ≡ -DA(x, t) (16)
where we have introduced the diamagnetic term D = e 2 uK/π. The latter is

J p (x, t) = euK π ∂ x θ(x, t) = e π ∂ t φ(x, t) (17)
Finally, we notice that the total current can be equivalently obtained from the continuity equation

∂ t ρ(x, t) + ∂ x J(x, t) = 0
by computing the time derivative of the density ρ(x, t) = -(e/π)∂ x φ(x, t)

∂ t ρ(x, t) = - ie π [H(A), ∂ x φ(x, t)]
using the LL Hamiltonian with the minimal substitution. The result of this calculation obviously recovers the expression for the current J(x, t) = J d (x, t) + J p (x, t).

We can now compute the average current within linear response theory, by doing perturbation theory in A(x, t). Since the diamagnetic part of the current J d (x, t) is already linear in the vector potential, we only need to take into account the paramagnetic contribution which gives

J p (x, t) = - ∞ 0 dt ′ dx ′ χ R JJ (xt; x ′ t ′ ) A(x ′ , t ′ ) ( 18
)
where we have defined the retarded current-current correlation function

χ R JJ (xt; x ′ t ′ ) = -iθ(t -t ′ ) [J p (x, t), J p (x ′ , t ′ )] = -i e 2 π 2 θ(t -t ′ ) [∂ t φ(x, t), ∂ t ′ φ(x ′ , t ′ )] (19)
and extended the time integral up to infinity. We stress here that the average in Eq.( 19) is taken over the initial density matrix (ground state of H 0 ) while the operators are evolved with the Hamiltonian H + . The resulting correlator is not the usual equilibrium one and in particular it is not time-translational invariant. This reflects the effect of the global and local quantum quenches that have been performed on the system at time t = 0. Finally we have for the current

J(x, t) = -DA(x, t) + J p (x, t) (20)
We now use an important identity relating the currentcurrent correlation (19) to the retarded Green's function of the field φ(x, t) defined as

G R (xt; x ′ t ′ ) = -iθ(t -t ′ ) [φ(x, t), φ(x ′ , t ′ )] (21)
We have

χ R JJ (xt; x ′ t ′ ) = e 2 π 2 ∂ t ′ ∂ t G R (xt; x ′ t ′ ) -2Dδ(x -x ′ )δ(t -t ′ ) (22)
Substituting this result into Eq. ( 20) and using Eq. ( 18), we find that the last term cancels the diamagnetic term exactly, and after an integration by parts we obtain

J(x, t) = - e 2 π 2 t 0 dt ′ dx ′ ∂ t G R (xt; x ′ t ′ ) E(x ′ t ′ ) (23)
This is in principle an exact result for the time dependent current through a finite-size LL and we can now further specify the electric field profile. We assume a sudden switching of the field, E(x, t) = θ(t)E(x) and in addition, since we are ultimately interested in the current due to the impurity (3), we assume that the potential drop occurs around x = 0, V (x) = V θ(-x), so that the electric field is effectively a delta-function

E(x, t > 0) = -∂ x V (x) = V δ(x).
Then, the current at x = 0, J(0, t) ≡ I(t) can be written only in terms of the exact retarded Green's function of the local field at the impurity site,

G R (t, t ′ ) = -iθ(t -t ′ ) [φ(t), φ(t ′ )] (24)
with φ(t) ≡ φ(0, t), and it reads

I(t) = - e 2 V π 2 t 0 dt ′ ∂ t G R (t, t ′ ) ( 25
)
This equation, which gives the time-dependent current in terms of an exact dynamical correlator of the local field, is the main result of this section. It generalizes to the timedependent quench problem the well known equilibrium result [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] . From this we can compute the current for a pure LL, using the expression for the bare retarded local

Green's function

G R 0 (t > t ′ ) = -K arctan Λ (t -t ′ ) ( 26
)
where Λ = u/α is an ultra-violet cut-off. We obtain

I 0 (t) = e 2 KV π 2 arctan Λt (27)
which approaches in the long time limit I 0ss /V = e 2 K/2π = e 2 K/h as expected. In the next section we will evaluate the corrections to the time dependent current due to the local potential. The forward scattering term in Eq. ( 9) turns out to not contribute to the current, as we discuss explicitly in appendix B, therefore in the next section we will focus our attention on the backward scattering term which is responsible for the interesting physical effects with V loc = g bs cos(2φ(0)).

A. Weak-Backscattering Corrections to Linear Conductance

We are now in the position to evaluate the weak backscattering correction to the conductance. We just need to evaluate the local Green's function for the field φ(t)

G ab (t, t ′ ) = -i φ a (t)φ b (t ′ ) (28)
where a, b = ± are Keldysh indices. We have defined the various Green's functions and the relations between them in appendix A. To lowest order in the backscattering g bs we obtain 49

δG ab (t, t ′ ) = - g 2 bs 2 dt 1 dt 2 cd (cd)F cd (t 1 , t 2 )Λ cd ab (t 1 , t 2 )
(29) where we define δG ab = G ab -G 0,ab ,

F cd (t 1 , t 2 ) = -i cos 2φ c (t 1 ) cos 2φ d (t 2 ) and Λ cd ab (t 1 , t 2 ) 4 = G 0,ac (t, 1)G 0,cb (1, t ′ ) + G 0,ad (t, 2)G 0,db (2, t ′ ) -G 0,ac (t, 1)G 0,db (2, t ′ ) -G 0,ad (t, 2)G 0,cb (1, t ′ ) ( 30
)
The retarded component is given by

G R (t, t ′ ) = - 1 2 ab b G ab (t, t ′ ) , (31)
Substituting the resulting expression into Eq. ( 25), and after some algebra we obtain for the current δI(t) = -KΛt

I(t) -I 0 (t) δI(t) = -C t 0 dt 1 ∂ t G R 0 (t, t 1 ) t1 0 dt 2 F R (t 1 , t 2 ) Φ(t 1 , t 2 ) (32) with C = 4e 2 V π 2 g
2 arctan Λt 2 - K 2 log 1 + Λ 2 t 2 1 1 + Λ 2 t 2 2
The retarded correlator entering the expression for the current is found to be

F R (t 1 > t 2 ) = - sin (2K arctan Λ(t 1 -t 2 )) 1 + Λ 2 (t 1 -t 2 ) 2 Kneq f (t 1 , t 2 )
(33) where we have introduced

f (t 1 , t 2 ) = 1 + 4Λ 2 t 2 1 1 + 4Λ 2 t 2 2 (1 + Λ 2 (t 1 + t 2 ) 2 )
Ktr/2 (34) and the exponents

K neq = K 0 2 1 + K 2 K 2 0 K tr = K 0 2 1 - K 2 K 2 0 ( 35
)
We can further simplify Eq. ( 32) by noticing that for

t, t 1 ≫ 1/Λ ∂ t G R 0 (t -t 1 ) ≃ -πKδ(t -t 1 )
Then we can write the correction to the time dependent current as

δI(t) = -C t 0 dτ sin [2K arctan Λτ ] (1 + Λ 2 τ 2 ) Kneq f (t, t -τ ) Φ(t, t -τ ) (36) with C = πKC = 4e 2 KV g 2
bs /π. In the next section we will discuss the behavior of this quantity as a function of time and quench amplitude.

B. Discussion

We start by discussing the equilibrium zero temperature case corresponding to K 0 = K. In this case we have K tr = 0 and K neq = K and the integral in Eq. ( 36) simplifies. We obtain for the transient current at long times (Λt ≫ 1)

δI(t)/V ∼ 1 t 2(K-1) (37)
i.e., for K > 1 the weak-backscattering correction to the current vanishes in the long time limit as a power law, the junction is perfectly conducting and the backscattering is irrelevant, this is the well known result from Kane and Fisher 6 . As soon as the quench amplitude is non-zero, K = K 0 , we find a number of interesting differences. We first focus on the long time steady state value of the backscattering current, δI ss = I ss -I 0ss , this is

δI ss /V = -4πK 2 g 2 bs G 0 ∞ 0 dτ τ sin (2K arctan Λτ ) (1 + Λ 2 τ 2 ) Kneq (38) G 0 = e 2 /2π = e 2 /
h the quantum of conductance, which we set to one from now on, G 0 ≡ 1. We first notice that the integral is well defined as long as K neq > 1, as expected from perturbative RG calculations 36 , since only in this regime the local backscattering is an irrelevant perturbation for the steady state. We plot in figure 1 the region of parameters where this condition is satisfied, and our calculations are therefore controlled. In figure 2 we plot the backscattering correction to the steady state current for different values of K in region I, as a function of K 0 . Quite interestingly, the correction to the steady state current is different from zero at any finite quench amplitude K = K 0 , in other words the systems deviates from the perfect conduction limit. For a small quench amplitude one finds that the steady state correction to the conductance is proportional to (

K -K 0 ) 2 , δI ss /V ∼ -g 2 bs (K -K 0 ) 2 (39)
We now discuss the transient behavior and the approach to the steady state. To this extent we evaluate numerically the integral in Eq. ( 36) and plot the result in figure 3. We find that the current decays to the steady state in a power law fashion,

δI(t) -δI ss ∼ 1 t α (40)
with an exponent α(K, K 0 ) whose dependence on the quench parameters is shown in the bottom panel of figure 3, as a function of K 0 at fixed K = 1.25 in the region I. We notice that the exponent α behaves nonmonotonically with K 0 and reaches a minimum value for K = K 0 , i.e., the decay of the current is stronger in the presence of a finite quench amplitude. We can get an analytical understanding of this by looking at the integral expression for δI(t), Eqn. (36). Using the fact that for large t we have

f (t, t -τ ) ∼ 1 + O(1/t) while Φ(t, t -τ ) ∼ τ + O(1/t), we obtain α = 2(K neq -1).
We notice that this is indeed confirmed by the numerical evaluation of the integral, although for larger values of the bulk quench, deviations start to appear that we attribute to the finite time resolution of the numerical integration.

We have analyzed only a few values of K while tuning K 0 , yet we expect this behavior to hold throughout the region I in figure 1. We however expect the perturbative expansion to eventually break down as we approach the regime K neq = 1 and enter region II, where the backscattering becomes a relevant perturbation, i.e., its strength grows under RG. In equilibrium, a complementary approach to explore this strong coupling regime amounts to starting from the limit of two decoupled semi-infinite systems and switching on a local tunneling. In the next section we will suitably generalize this approach to the quench case and discuss to what extent it can be used to extract useful information about the strong coupling regime.

III. TRANSPORT IN THE TUNNELING REGIME

We now consider a different quench protocol which will allow us to access the strong coupling regime where the back-scattering is relevant. Let us suppose at t = 0 we have two disconnected and identical semi-infinite TL models describing 1D interacting spinless fermions. By definition, each lead contains left and right moving fermions interacting in the bulk through intra-branch (g 0 4 ) and inter-branch (g 0 2 ) scattering processes, as in the previous case. A major difference however arises due to the presence of a sharp edge in each lead, which imposes open boundary conditions for the fermionic field. As a result, and differently from the translational invariant case of the previous section, the two fermionic species are not independent of each other. This suggests an equivalent and more convenient representation of each lead i = 1, 2 in terms of a single chiral fermionic field ψ i (x) defined on an infinite system 50 (i.e., obeying periodic boundary conditions). In terms of this field, the initial Hamiltonian in each lead becomes

H 0i = H free + g 0 4 2 dx ρ i (x)ρ i (x) + g 0 2 2 dxρ i (x) ρ i (-x) (41)
where

H free = v F dxψ † i i∂ x ψ i , while ρ i (x) = ψ † i (x)ψ i (x)
, and one should notice the non-local g 0 2 interaction that results from the single chiral field representation.

We prepare the system in the ground state |Ψ 0 of H 0 = i H 0i and then, for t > 0, we evolve the system with a different Hamiltonian H

H = i=1,2 H i + i=1,2 µ i N i + H T ( 42
)
where we have (a) quenched the bulk interactions,

H i = H free + g 4 2 dx ρ i (x)ρ i (x) + g 2 2 dxρ i (x) ρ i (-x) (43)
(b) switched on a tunneling term coupling the two wires at the edge

H T = λ 0 ψ † 1 (0)ψ 2 (0) + h.c. ( 44
)
and finally (c) switched on a bias voltage µ 1,2 = ±eV /2 that couples to the charge imbalance. We are interested in computing the time dependent current

I(t) = Ψ(t)| Î|Ψ(t) = Ψ 0 | e iHt Î e -iHt |Ψ 0 ( 45
)
where the current operator is

Î = ie[H, (N 1 -N 2 )] = -2ieλ 0 ψ † 1 (0)ψ 2 (0) -h.c. (46)
In order to proceed, it is useful to perform a unitary transformation to eliminate the bias voltage in H, which amounts to going to the rotating frame defined by an operator Ω † (t) = exp (i i µ i N i t). Inserting Ω into the expression for the current (45) we obtain

I(t) = Ψ(t)| Î(t)| Ψ(t) (47)
where the new state

| Ψ(t) = Ω † (t)|Ψ(t) (48)
satisfies the Schrodinger equation

i∂ t | Ψ(t) = H(t)| Ψ(t) (49)
with a rotated-frame Hamiltonian

H(t) = i∂ t Ω † Ω + Ω † HΩ ≡ i H i + H T (t) (50)
While the bias voltage has disappeared, the tunneling has become explicitly time-dependent

H T (t) = Ω † (t)H T Ω(t) = λ 0 e ieV t ψ † 1 (0)ψ 2 (0) + h.c. (51)
where we have used the fact that

Ω † (t)ψ i (x)Ω(t) = e -iµit ψ i (x) (52)
The current operator in (45) also acquires an explicit time dependence

Î(t) = Ω † (t) ÎΩ(t) = -2ieλ 0 e ieV t ψ † 1 (0)ψ 2 (0) -h.c. (53)
The time dependent state | Ψ(t) can be written in the interaction picture O I = e iHit Oe -iHit with respect to the unperturbed Hamiltonian of the uncoupled wires as

| Ψ(t) = e -i i Hit T e -i t 0 dt ′ HT,I (t ′ ) |Ψ 0 (54)
We can now evaluate the time dependent current to the lowest order contribution in the tunneling, which gives

I(t) = -i t 0 dt 1 Ψ 0 | ÎI (t), H T,I (t 1 ) |Ψ 0 ( 55
)
where all the operators are explicitly time dependent as in previous equations and evolved with the Hamiltonian of the uncoupled wires, after the global quench.

It is now useful to use bosonization, in its openboundary formulation 50 to proceed with the evaluation of the tunneling current. We first write the fermionic operator in each lead in terms of a single chiral bosonic mode

ψ i (x) = 1 √ α e iϕi(x) (56)
as well as the electron density in the lead i as

ρ i (x) = ∂ x ϕ i (x)/2π (57)
The Hamiltonian before and after the global quench of the interactions (Eqns 41 and 43) become now quadratic in this bosonic field and can be easily diagonalized with a Bogoliubov rotation (see appendix C). For example we have for H 0i

H 0i = u 0 4 1 K 0 + K 0 dx 2π (∂ x ϕ i (x)) 2 + - u 0 4 1 K 0 -K 0 dx 2π ∂ x ϕ i (x)∂ x ϕ i (-x) (58)
The resulting harmonic theory is again described in terms of two parameters, the sound velocity u 0 , u and the Luttinger parameter K 0 , K, which by construction are identical in the two leads, and whose expressions in terms of the coupling constant g 0 2 , g 0 4 and g 2 , g 4 are the same as in the translational invariant case (see equations (11)). We rewrite them here for convenience

u 0 = v F (1 + g 0 4 /2πv F ) 2 -(g 0 2 /2πv F ) 2 (59)
K 0 = 1 + g 0 4 /2πv F -g 0 2 /2πv F 1 + g 0 4 /2πv F + g 0 2 /2πv F (60)
with similar relations holding for u, K as functions of g 2 , g 4 . In addition we need to write the tunneling and current operator in the bosonic language. To this extent it is useful to introduce the combinations

ϕ ± (x) = ϕ 1 (x) ± ϕ 2 (x) √ 2 (61)
and to notice that only the ϕ -combination enters both the tunneling and the current. Indeed we have disregarded the Klein factors which can be shown to be unimportant. Thus,

H T (t) = 2λ cos √ 2ϕ -(t) -eV t (62) Î(t) = -4eλ sin √ 2ϕ -(t) -eV t (63)
where we have defined ϕ -(t) ≡ ϕ -(0, t) and λ = λ 0 /α.

A. Tunneling Correction To Current

In order to proceed further, we plug Eqns. (62,63) in the expression for the tunneling current to obtain, after some simple algebra,

I(t) = -8 eλ 2 t 0 dτ sin (eV τ ) F R (t, t -τ ) ( 64
)
where F R (t, t ′ ) is the retarded component of the correlator In the white area (region I) corresponding to K dual neq > 1, perturbation theory is well behaved, while the shaded area (region II) corresponds to the regime where the tunneling is a relevant perturbation in the steady state and the perturbative expansion breaks down. The dashed line is a reference only to the validity of the weak backscattering expansion previously discussed.

F ab (t, t ′ ) = -i cos √ 2ϕ - a (t) cos √ 2ϕ - b (t ′ ) =
= - i 2 e i √ 2ϕ - a (t) e -i √ 2ϕ - b (t ′ ) (65)
It is convenient to express this correlator in terms of the fields ϕ i (t) in the two semi-infinite leads. We have

F ab (t, t ′ ) = - i 2 e iϕ1a(t) e -iϕ 1b (t ′ ) e iϕ 2b (t ′ ) e -iϕ1a(t)
where we have used the fact that in the absence of tunneling the two leads are decoupled. To compute this correlator we need therefore the local Green's function of the field ϕ i after a quench of the bulk interaction parameter in a semi-infinite chain, which we give in appendix C. The result for the retarded component,

F R = -ab b F ab /2 is found to be F R (t 1 > t 2 ) = - sin 2 K arctan Λ(t 1 -t 2 ) 1 + Λ 2 (t 1 -t 2 ) 2 K dual neq f (t 1 , t 2 ) (66)
where we have introduced

f (t 1 , t 2 ) = 1 + 4Λ 2 t 2 1 1 + 4Λ 2 t 2 2 (1 + Λ 2 (t 1 + t 2 ) 2 ) 2 K dual tr /2
(67)

with the exponents

K dual neq = 1 2K 0 1 + K 2 0 K 2 K dual tr = 1 2K 0 1 - K 2 0 K 2 (
68) We stress that Eq. ( 64) is perturbative in the tunneling but contains the bias voltage eV to all orders. In fact in the absence of the bulk quench, K = K 0 , we can recover from this the well known result 6 for the nonlinear I-V characteristic of the wire, I(V ) ∼ V 2/K-1 .

In the following we will instead consider mostly the low bias regime, for a finite quench amplitude and to this extent we will evaluate the current in the low bias regime and equivalently the non-linear differential conductance at zero bias. Note that the latter I/V | V =0 = 0 in the absence of quench. some simple manipulations, we can write the time dependent current at small bias voltage as

I(t) = 8e 2 V λ 2 t 0 dτ τ sin 2 K arctan Λτ (1 + Λ 2 τ 2 ) K dual neq f (t, t -τ ) (69)
We notice the result we have obtained for the current in the tunneling regime is dual to the one in the weakbackscattering case in the sense that the two currents are related by a transformation of the Luttinger parameters, before and after the quench, into

K, K 0 -→ 1 K , 1 K 0 (70)
This is the natural generalization of the equilibrium duality [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] . Yet, as one can immediately see from the expressions (35) for K neq , K tr and (68) for K dual neq , K dual tr , the exponents controlling the decay of the correlation in the tunneling and weak-backscattering regime are not related by this simple duality, i.e.,

K dual neq = 1/K neq (71)
K dual tr = 1/K tr (72)
unless of course K = K 0 . This will have important consequences as we are going to discuss in the next section.

A simple way to understand the origin of this result is to notice that while for the Hamiltonian after the quench the duality still holds, the initial condition of the problem (the ground state of the Luttinger model with interaction parameter K 0 ) is not dual when written in the basis of the natural eigenmodes of the system after the quench (Luttinger model with interaction parameter K). The fact that the Luttinger model is integrable and therefore never fully loses memory of this initial condition implies the breakdown of duality in the stationary state.

B. Discussion

As in the weak-backscattering case, we start discussing the equilibrium zero temperature result corresponding to K 0 = K. In this case we have K tr = 0 and K neq = K, and the integral in Eq. (69) simplifies. For the transient current at long times (Λt ≫ 1) we obtain,

I(t)/V ∼ 1 t 1/K-1 (73)
i.e. for K < 1 the correction to the current vanishes in the long time limit as a power law, the junction is perfectly insulating and the tunneling is irrelevant 6 . We now consider the finite quench case, K = K 0 , and first focus on the long time steady state value of the tunneling current which reads

I ss /V = 8e 2 λ 2 ∞ 0 dτ τ sin 2 K arctan Λτ (1 + Λ 2 τ 2 ) K dual neq (74)
We notice that the integral is well defined as long as K dual neq > 1, a condition which is satisfied in the region of parameters plotted in figure 4 (region I). Interestingly, this region is not simply the dual of the region I for the weak-backscattering case, (light shaded line inside region I in figure 4) as it would be in the equilibrium case. This is a consequence of Eqns (71,72). Therefore we conclude that the quenched impurity model does not display a full duality, contrarily to other nonequilibrium realizations such as the driven noisy one 41 . We will come back to this point in the next section. Let us now discuss the behavior of the steady state current. Quite interestingly, the steady state current is different from zero at any finite quench amplitude K = K 0 , in other words the system deviates from the perfect insulating limit. In figure 5 we plot the steady state value of the tunneling current at fixed K and as a function of K 0 . As in the weak-backscattering limit we find that for small quench amplitudes the steady state current is

I ss /V ∼ λ 2 (K -K 0 ) 2 (75)
We now discuss the transient behavior and the approach to the steady state. We evaluate numerically the integral in eq. ( 69) and plot the result in figure 6. We find that the current decays to the steady state as a power law

I(t) -I ss ∼ 1 t α (76)
with an exponent α(K, K 0 ) whose dependence on the quench parameters is shown in figure 6. As in the weak- backscattering case previously discussed, we find that the decay of the current is in general faster for a finite quench amplitude, in accordance with the analysis of Ref. 37.

In particular the exponent α reaches its minimum for K = K 0 and behaves as α ∼ 2(K dual neq -1).

IV. DISCUSSION

Let us summarize the results of previous sections and discuss their physical consequences. We started from the limit of a uniform system and found that for K neq > 1, corresponding to region I of fig. 1, perturbation theory in the backscattering is well behaved. While this means that in equilibrium, K → K 0 , the conductance approaches the unitary limit, a bulk quench K = K 0 gives rise to a finite deviation from this perfect conducting regime. This is consistent with the RG analysis we developed in Ref. 36. There we showed that although for K neq > 1 the backscattering is an irrelevant perturbation for the steady state, it nevertheless gives a sizeable effect, generating an effective temperature for the local degree of freedom, which explains, at least qualitatively, the deviation from the unitary limit. When K neq < 1, the backscattering becomes relevant, and perturbative approaches in the back-scattering breaks down, and one enters the strong coupling regime.

A possible approach to grasp the behavior of the system in this limit of strong backscattering is to start from the strong coupling fixed point, assuming the effective temperature is not enough to cut the growth of the backscattering under RG flow, and perturb this system of two decoupled wires with the sudden switching on of a local tunneling. We have considered this regime in section III and found that for K dual neq > 1, corresponding to region I in fig. 4, perturbation theory in the tunneling is well behaved. While this could naively suggest the strong coupling fixed point remains stable as in equilibrium, we have found that a finite tunneling current appears in the steady state for any finite quench amplitude K = K 0 . This result is again reminiscent of a thermal behavior and we may speculate that an effective temperature would indeed be generated under RG by the local tunneling, in analogy with the backscattering case. Taken together, these results suggest that a global bulk quench has a dramatic effect on the problem, smearing out the sharp distinction between the conducting and the insulating phase that exists in equilibrium at zero temperature. To further appreciate the role played by the quench amplitude in the problem, it is useful to look again at the steady state current in both regimes. In figure 7 we plot the weak-tunneling and weakbackscattering correction to the current as a function of K and at fixed K -K 0 . At zero quench amplitude this would give the usual sharp transition from an ideal insulator for K < 1 to a perfect conductor for K > 1. However as we clearly see, any small finite quench amplitude turns this into a smooth crossover. Finally, it is worth noticing that, as we mentioned earlier, the two regions K neq > 1 and K dual neq > 1 are not dual to each other. This means that there exist regions of the parameter space (K, K 0 ) where neither the backscattering nor the tunneling are relevant, i.e., perturbation theory is well behaved. We interpret this result as a further signature of the effective temperature behavior of the problem.

I ss /V K -K 0 = 0 K -K 0 = 0.01 K -K 0 =0.02 K -K 0 = 0.03 K -K 0 = 0.04 K -K 0 = 0.
These results for transport offer therefore another example of the quench-induced decoherence mechanism that we had previously identified in the orthogonality catastrophe problem 36 . As in that case, the relevant energy scale, in the limit of small quenches, is proportional to the quench amplitude squared

γ ⋆ ∼ g 2 (K -K 0 ) 2 (77)
where g in the equation above indicates any source of local non-linearity, either due to the backscattering potential in the weak impurity limit, or due to the tunneling term in the limit of two decoupled chains. While the qualitative behavior is suggestive of an effective thermal behavior, it is important to stress that such an analogy does not fully carry over. In particular, the characteristic power law structure exhibited by the finite temperature current in equilibrium problem 5,6,[START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] , on both sides of the transition, appears to be washed out by nonequilibrium effects. Indeed, as we have discussed earlier, the scaling of the current in the steady state does not show signature of non-trivial power laws, even in the small quench amplitude regime where transport is essentially set by the quench-induced decoherence scale γ ⋆ , with very weak additional dependence from the Luttinger parameters. This suggests that the analogy between this scale and an effective temperature, meant here as an infrared cutoff, cannot be pushed too far.

Finally it is worth comparing our results for the steady state impurity physics with the results and the predictions of Ref. 38. Here the authors computed the zero frequency charge susceptibility in a uniform LL after an interaction quench starting from a free system. For repulsive interactions they found it to diverge at the backscattering wave vector and concluded that, as in equilibrium, even a weak single impurity strongly disturbs the homogeneous LL or, in an RG language, the perfect chain fixed point is unstable. Then, they computed the local density of states close to an open boundary and found that it saturates to a non-zero value, with power-law corrections at low frequency. From the scaling of these corrections they conclude that the steady state analog of the open-chain fixed point is stable, with a power-law finite temperature conductance, which eventually crosses over to a non-zero value at low enough temperature. Based on their results the authors argue that the fixed point structure of a single impurity in a nonequilibrium steady state LL is similar to the one in equilibrium 38 .

As we have seen, our transport results combined with our analysis of the transient orthogonality catastrophe problem 36 , suggest a rather different picture and highlight the importance of inelastic effects for the steady state impurity physics after the quench. In particular in the back-scattering case we found that, due to the global quench, a local non-linearity generates an effective decoherence in the problem, even in the regime where it would nominally be irrelevant at the one loop level. This quench-induced decoherence energy scale reduces the value of the zero temperature conductance away from the perfect uniform limit and, more importantly, competes with the growth of the backscattering as one enters the strong coupling phase, corresponding to K < 1 for the case K 0 = 1 considered in Ref. 38. The result of this competition is clearly a subtle issue to establish firmly. Our results for the Loschmidt echo 36 suggested that well within the putative strong coupling phase this new energy scale acts as an effective infrared cut off, resulting in an exponential decay of the echo. This result would suggest that, quite differently from the equilibrium case, a weak single impurity does not disturb substantially the uniform system, i.e., one should not trust the perturbative break down for K neq < 1.

To further confirm this intuition, we have considered the opposite limit of an open chain, perturbed by an interaction quench and by a sudden switching on of a local tunneling. For any finite bulk quench amplitude we find a non-vanishing zero temperature conductance essentially set by the quench-induced decoherence scale. This result suggests that, differently from the equilibrium case, one never reaches the open-chain fixed point when the bulk has been quenched because the weak tunneling is always effective in generating a non-zero conductance. We may interpret this in an RG picture for the tunneling problem, dual to what we discussed in Ref. 36 for the backscattering case, where the local non-linearity generates an effective decoherence that cuts off the RG flow. While a nonzero tunneling conductance was also found in Ref. 38, the physics behind this effect was not discussed. We now fully clarify its origin as a quench-induced decoherence phenomenon.

V. CONCLUSIONS

In this paper we have studied the problem of transport through a localized impurity in a Luttinger model far from equilibrium due to a bulk interaction quench. This model in its ground state has a familiar and well studied zero temperature transition between a conducting and an insulating phase, depending on the value of the bulk interaction parameter K, and has an associated finitetemperature and finite-voltage correction which shows universal power law behavior.

Here we have presented perturbative calculations in the strength of the impurity potential starting both from the limit of a uniform liquid and from the one of two decoupled semi-infinite systems. In the former case a standard bosonization approach gives rise to a local backscattering term (described by a boundary sine-Gordon problem in a time-dependent nonequilibrium bath) that we treat to leading order in perturbation theory. In the latter, we use open boundary bosonization to treat the problem and obtain a dual formulation again in terms of a boundary sine-Gordon model where the non-linearity comes from a local tunneling between the isolated wires. Our results quite generically reveal that the nonequilibrium excitation of bulk modes induced by the global quench has important and peculiar effects on the conducting-insulating transition, which gets smeared out into a crossover for any small quench amplitude K -K 0 = 0. In addition, the dynamics of the current displays a faster decay towards the steady state as compared to the equilibrium zero temperature case.

All together this suggests that the global quench effectively induces a decoherence mechanism for the local degrees of freedom and we have identified an energy scale γ ⋆ associated with this. While this behavior is qualitatively similar to a finite effective temperature, the decoherence energy scales enters the problem differently than a temperature. In particular the steady state current in both limits is set by the decoherence scale but does not show power laws in this energy scale.

An interesting question concerns the generalities of these results beyond the quenched Luttinger model. We may speculate that quench-induced decoherence is a generic phenomenon that occurs in other interacting quantum impurity problems coupled to environments which are non-thermal at long times but rather flow to a Generalized Gibbs ensemble steady state. We leave the investigation of this intriguing question to future work. where H free = v F α=L,R s α L 0 dxψ † α i∂ x ψ α with s L = -s R = 1. The crucial observation is that imposing open boundary conditions introduces a constraint between right and left moving fermions which are no longer independent, but rather satisfy 50

ψ L (-x) = -ψ R (x) (C2)
and is a direct consequence of the existence of a single Fermi point. We can therefore fold the line into (-L, L) and use Eq. (C2) to express the Hamiltonian only in terms of right movers

H 0 = H free + g 0 4 2 dx ρ R (x)ρ R (x)+ g 0 2 2 dxρ R (x)ρ R (-x) (C3
) where H free = v F dxψ † R i∂ x ψ R and one has to notice the non-local g 0 2 interaction that results from the single chiral field representation. The right moving fermions now satisfy periodic boundary conditions and can be bosonized as usual in terms of a single bosonic field

ψ(x) = 1 √ α e iϕ(x) (C4)
in terms of which the Hamiltonian reads

H 0 = v F 2 1 + g 0 4 2πv F dx 2π (∂ x ϕ i (x)) 2 + - g 0 2 4π dx 2π ∂ x ϕ i (x)∂ x ϕ i (-x) (C5)
To diagonalize this problem it is useful to decompose the field into normal modes ϕ(x) = q>0 π qL e -aq/2 b q e iqx + b † q e -iqx (C6) in terms of which the Hamiltonian becomes

H 0 = q>0 A 0 q b † q b q -B 0 q b † q b † q + b q b q (C7)
where

A 0 = v F 1 + g 0 4 /2πv F (C8) B 0 = g 0 2 4π (C9)
We can diagonalize the Hamiltonian via a Bogoliubov rotation parameterized by the angle α 0 .

b q = cosh α 0 η q -sinh α 0 η † q (C10)

b † q = -sinh α 0 η q + cosh α 0 η † q (C11)

The Hamiltonian when written in terms of the fields η q is, H 0 = q v 0 q η † q η q (C12) = -i K 0 dq q e -αq coth( ε q 2T ) cos ε q (t -t ′ ) 1 + K 2 0 K 2 + cos ε q (t + t ′ ) 1 -

K 2 0 K 2 (C32)
The overall result is that the local Green's functions for the boundary field in the OBC case can be obtained from the bulk one by sending

K 0 , K → 2 K 0 , 2 K (C33)
a result which naturally generalizes the well know equilibrium result. A simple interpretation [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] is that the field in the OBC case corresponds to the θ field in the translational invariant case (the field φ being frozen due to the boundary condition [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] ), which is known to have correlators with K → 1/K. On top of this, the field is now living in a semi-infinite space, which explains the extra factor of 2. Finally, let us compute the following correlator which will be useful in evaluating the tunneling current f ab (t 1 , t 2 ) = e iγϕa(t1) e -iγϕ b (t2) (C34)

As before, this can be written only in terms of the local Green's function of the field ϕ a (t)

f (

FIG. 1 .

 1 FIG.1. Parameter space for the weak-backscattering problem. In the white area (region I) corresponding to Kneq > 1, perturbation theory is well behaved, while the shaded area (region II) corresponds to the regime where backscattering is a relevant perturbation in the steady state, and the perturbative expansion breaks down. The dashed line is a reference only to the validity of the strong coupling (weak-tunneling) expansion (see further below).

2 FIG. 2 .

 22 FIG.2. Backscattering corrections to steady state current (top panel) as a function of K0 and for different values of K. A finite quench amplitude restores transport with deviations from K0 = K (arrows) vanishing quadratically for small quenches, δIss/V ∼ -g 2 bs (K -K0) 2 (see bottom panel, log-log scale). Parameters: Λ = 1, g bs = 0.1

2 FIG. 3 .

 23 FIG.3. Power law decay in time of the weak backscattering correction to the current, for fixed K = 1.25 and different values of K0. We notice the exponent α depends nonmonotonically on K0 and takes a minimum value in the equilibrium case K = K0 (see bottom panel).

FIG. 4 .

 4 FIG.4. Parameter space for the tunneling problem. In the white area (region I) corresponding to K dual neq > 1, perturbation theory is well behaved, while the shaded area (region II) corresponds to the regime where the tunneling is a relevant perturbation in the steady state and the perturbative expansion breaks down. The dashed line is a reference only to the validity of the weak backscattering expansion previously discussed.

2 FIG. 5 .

 25 FIG. 5. Steady state tunneling current as a function of K0and for different values of K (top panel). We see that a nonzero quench amplitude results in a finite tunneling current, vanishing quadratically as K0 → K (arrows), Iss/V ∼ λ 2 (K -K0) 2 (see bottom panel, log-log scale). Parameters: Λ = 1, λ = 0.1

FIG. 6 .

 6 FIG.6. Power law decay in time of the weak tunneling correction to the current, for fixed K = 0.9 and different values of K0. We notice the exponent α depends non-monotonically on K0 and takes a minimum value in the equilibrium case K = K0 (see bottom panel).

  Steady state current versus Luttinger parameter K0 at fixed values of the quench amplitude K -K0, in the tunneling regime (left panel, for λ = 0.1) and in the weak backscattering limit (right panel, for g bs = 0.1). We see that the sharp transition between insulator and conductor driven by K0 gets smeared out by the global quench.

		0.9996
	05	0.999
		0.9984
		0.9978
	1 1.2 1.4 1.6 1.8 2
	K 0
	FIG. 7.

  t 1 , t 2 ) ab = exp -γ 2 2 Φ ab (t 1 , t 2 ) (C35)where we have defined Φ ab (t 1 , t 2 ) as the following combination of local Green's functionsΦ ab (t 1 , t 2 ) = , t 1 ) -iG K ϕϕ (t 1 , t 2 ) + +i aG A ϕϕ (t 1 , t 2 ) + bG R ϕϕ (t 1 , t 2 )Using the above results we can writef ab (t 1 , t 2 ) = f (t 1 , t 2 ) exp [-iΥ ab (t 1 , t 2 )] 1 + Λ 2 (t 1 -t 2 ) (1 + Λ 2 (t 1 + t 2 ) 2 ) 2

	i 2	G K ϕϕ (t 1 , t 1 ) +	i 2	G K ϕϕ (t 1 2	neq Kdual	(C36)
	where we have defined
		Kdual neq/tr =	γ 2 4K 0	1 ±	K 2 0 K 2	(C37)
	as well as
	f (t 1 , t 2 ) =	1 + 4Λ 2 t 2 1	1 + 4Λ 2 t 2 2	Kdual tr	/2	(C38)
	Υ ab (t 1 , t 2 ) =	γ 2 2	aG A

ϕϕ (t 1 , t 2 ) + bG R ϕϕ (t 1 , t 2 ) (C39)

We thank Ehud Altman, Volker Meden, Achim Rosch and Hubert Saleur for helpful discussions. AM was supported by nsf-dmr 1303177.

Appendix A: Green's Functions -Definitions and Identities

We define the contour ordered Green's function for a real bosonic field φ(t) in the α, β = ± Keldysh basis

where we have, following the convention of Ref. 49

We then define the retarded, advanced and Keldysh components as

and find the following useful relation

from which the following relations follow

Appendix B: Forward Scattering Contribution to the Current

In this appendix we show that the forward scattering does not contribute to the time dependent current. To see this we follow standard steps 51,[START_REF] Gogolin | Bosonization and Strongly Correlated Systems[END_REF] and introduce first the even and odd combinations of the LL fields θ, φ, defined as

In terms of these fields the bulk LL Hamiltonian decouples in the two channels

It is particularly convenient to write the local scattering in terms of these fields, due to their well defined properties under inversion. We have

i.e. the forward scattering couples to the antisymmetric mode while the backward term only to the symmetric one. Now we recall Eq. ( 25) of the main text, which shows that the time-dependent current I(t) is expressed only in terms of the Green's function of the local field φ(0), which turns out to be related only to the symmetric combination φ(0) = K/2Φ s (0), from which we conclude that the current does not depend on the forward scattering potential.

Appendix C: Quenches in the Open-Boundary Tomonaga Luttinger Model

In this appendix we briefly discuss the quench problem in a TL model with open boundary conditions, and compute in particular the local bosonic correlator relevant for obtaining the tunneling current. We first recall how to diagonalize the initial Hamiltonian H 0 , as the same strategy will be used for the Hamiltonian after the quench, and this will highlight the major differences with the bulk case. We follow the treatment in Ref. 50 to which we refer the reader to for further details.

Let us start with the problem of spinless fermions ψ(x) defined on a line [0, L] with open boundary conditions (OBC) with the Hamiltonian

for an angle α 0 given by

The sound velocity reads

while the angle α 0 is defined as

We can now proceed to discuss the problem of an interaction quench of the bulk interactions g 0 2 , g 0 4 → g 2 , g 4 . This amounts to suddenly change A 0 , B 0 → A, B in the bosonic Hamiltonian. We can repeat the previous steps and introduce operators γ † q , γ q

giving a quadratic Hamiltonian when written in terms of the fields γ q H = q vq γ † q γ q (C18)

for an angle α given by

Using the above transformations, we can easily express the time evolved bosonic operators in terms of those diagonalizing the initial Hamiltonian

where (assuming as usual a Galilean invariance preserving quench such that v 0 K 0 = vK) f q (t) = cosh α 0 cos ε q t -i cosh (2α -α 0 ) sin ε q t g q (t) = sinh α 0 cos ε q t + i sinh (2α -α 0 ) sin ε q t

In particular we can write down the time-evolved local bosonic field ϕ(x = 0, t) as

where

from which all the local Green's function can be evaluated. We start with the retarded component

Inserting the time-dependent field expansion (C23), we obtain

Similarly for the Keldysh component

we obtain

e -aq q coth ε q 2T Re B q (t) B * q (t ′ ) (C28) After some algebra and using the fact that ReB q (t)B * q (t ′ ) = 1 K 0 cos ε q t cos ε q t ′ + + K 0 K 2 sin ε q t sin ε q t ′ (C29)

Im B q (t) B * q (t ′ ) = -

we obtain the following results for the local Green's functions

e -aq q sin ε q (t -t ′ ) = = -2 K arctan (Λ (t -t ′ )) and G K ϕϕ (t, t ′ ) = -iπ K 0 L q e -αq q coth( ε q 2T ) cos ε q (t -t ′ ) 1 + K 2 0 K 2 + cos ε q (t + t ′ ) 1 -


est destinée au dépôt

  Elena M. Bennett, Cecile H. Albert, Aaron Ball, Jeffrey Cardille, Karine Dancose, Sylvestre Delmotte, Andrew Gonzalez, Hsin Hui-Huang, Martin J. Lechowicz, Katie Liss, et al. To cite this version: Elena M. Bennett, Cecile H. Albert, Aaron Ball, Jeffrey Cardille, Karine Dancose, et al.. The Montérégie Connection: Understanding How Ecosystems Can Provide Resilience to the Risk of Ecosystem Service Change. Atlas of Ecosystem Services: Drivers, Risks, and Societal Responses, Springer International Publishing, 2019, 978-3-319-96228-3. hal-02010820


H Berry
email: hugues.berry@inria.fr

Nonequilibrium phase transition in scattered cell communities coupled by auto/paracrine-like signalling

Keywords: Signal transmission, Autocrine relay, Stochastic models, Critical phenomena, Directed percolation

Auto/paracrine cell-to-cell communications via diffusive messengers can be coupled to a positive feedback loop in which cell stimulation by a messenger results in the production of new messengers. This yields a potential mechanism for relay transmission of the emitted message. This paper investigates the influence of noise on this mutual coupling of the cells with their environment, using numerical simulations of a stochastic minimal model. The results demonstrate that the deterministic (mean-field) approximation of this stochastic process fails short of predicting its behaviour because of the presence of strong noise-induced fluctuations. Instead, the behaviour of the model can be explained by the occurrence of a nonequilibrium phase transition, which is found to be in the universality class of directed percolation. This provides a theoretical framework to understand signal transmission in these stochastic systems.

Introduction

Complex behaviours in cell communities such as self-organization and emergent phenomena may result from coupling of the cells with an environment they dynamically modify. For instance, cells often respond to molecules in their environment via intracellular signalling pathways that eventually result in altered concentration of the very extracellular molecular species that triggered the pathway. A well-known example is auto/paracrine signalling. In this paradigm, cells emit a peptidic factor (e.g. EGF) that diffuses in the extracellular space until it reaches a neighbouring cell (paracrine signaling) or the source cell that emitted it (autocrine signaling) [START_REF] Wiley | Computational modeling of the EGFreceptor systems: a paradigm for systems biology[END_REF]. Stimulation by the diffusive factor may in turn trigger intracellular signalling cascades (e.g. the MAPK pathway) that eventually lead to the release of new diffusive factor molecules in the environment (positive feedback loop) (Pribyl et al., 2003a;[START_REF] Freeman | Feedback control of intercellular signalling in development[END_REF].

Strikingly, this process can be considered as a rather generic communication pattern, that is not restricted to members of the EGF family or cytokines, but is also encountered in cAMP-mediated communications between Dictyostelium discoideum cells; acylated homoserine lactone-mediated quorum sensing in Vibrio fischeri [START_REF] James | Luminescence control in the marine bacterium Vibrio fischeri: An analysis of the dynamics of lux regulation[END_REF] and other bacteria (deKievit & Iglewski, 2000); or even airborne virus spreading -which has been implied in the dissemination of the foot-and-mouth disease [START_REF] Gloster | Re-assessing the likelihood of airborne spread of foot-and-mouth disease at the start of the 1967-1968 UK foot-and-mouth disease epidemic[END_REF] or the influenza virus flu [START_REF] Hammond | Impact of atmospheric dispersion and transport of viral aerosols on the epidemiology of influenza[END_REF] for instance-.

Although dissimilar, these examples share common features. First, the messenger travels at random, through diffusional or nondirectional transport. Moreover, encounter with the messenger molecule alters cell functioning in such a way that new messengers are ultimately excreted. Finally, stimulation by the messenger molecules is usually followed by a refractory or lag period during which the cell does not respond to new encounters with messenger molecules.

Broadly speaking, this process can be thought of as implementing relay broadcasting: "if a message is received, relay it to one of your nearest neighbours".

Hence, though the spatial range of a single messenger molecule may be limited [START_REF] Shvartsman | Spatial range of autocrine signalling: modelling and computational analysis[END_REF], the presence of the feedback loop may allow messages to be transmitted over long distances, very much like the spreading of epidemics.

For instance, in a detailed deterministic continuous biophysical model, Pribyl et al. (2003) showed that this kind of relay autocrine signalling may give rise to messenger travelling waves that spread over the entire cell population.

However, like many processes in cell biology, this mechanism comes with inherent noise or stochasticity at several levels. First, because of the diffusive nature of the messenger movements, the target cell of an emitted messenger is random, i.e. cannot be precisely specified. Secondly, because it relies on intrinsically stochastic biochemical reactions, the triggering of an intracellular signalling pathway upon cell-messenger molecule interaction is probabilistic. Finally, at every moment, a messenger molecule can be removed from the system, either by extracellular proteolysis, or by scavenging in the extracellular space.

In recent years, several papers have presented biophysical models for auto/paracrine signalling. However, to the best of our knowledge, these studies consisted either in deterministic models of auto/paracrine relay transmission (Pribyl et al., 2003a;2003b) or in stochastic models for messenger-cell interactions, but without feedback relay loops [START_REF] Batsilas | Stochastic model of autocrine and paracrine signals in cell culture assays[END_REF][START_REF] Shvartsman | Spatial range of autocrine signalling: modelling and computational analysis[END_REF].

H BERRY

The goal of the present work is thus to study the influence of noise in auto/paracrine-like relay broadcasting systems. In particular, it investigates the collective behaviour exhibited by the mutual coupling between cells and their environment, and how messages can be transmitted in stochastic conditions. In this framework, we focused on the basic mechanisms underlying these generic processes, and voluntarily restricted our model to the elementary ingredients identified above. The results demonstrate that classical deterministic modelling fails to predict the behaviour observed in stochastic simulations because of large fluctuations due to the intrinsic noise. Moreover, we show that the behaviour of this stochastic system is due to the occurrence of a nonequilibrium phase transition and present evidence that this nonequilibrium phase transition is in the universality class of directed percolation.

The model

As already mentioned above, we are interested here in the most fundamental processes implied in auto/paracrine-like relay signalling. Hence, we study a "minimal" model that incorporates only the few basic ingredients enumerated above (messenger diffusion, feedback-positive loop for messenger relay, refractory phase). Furthermore, the model is intrinsically stochastic regarding cell stimulation by messenger molecules, as well as messenger molecule survival in the extracellular space. For the same reasons, we study a one-dimensional version of the model, i.e. messenger molecules and cells evolve in a one-dimensional lattice of linear size L with periodic boundaries (i.e. a circle made of L equidistant lattice sites).

a) Messenger molecules

Messenger molecules undergo isotropic diffusion on the lattice via non-interacting nearest-neighbour random walks without excluded volume conditions (i.e. several messengers can occupy a same site at a given time). Every messenger can also be removed from the lattice with rate λ per time step.

b) Cells

Besides messenger molecules, the lattice also contains n immobile point cells that are regularly scattered over the lattice. Each cell occupies a single lattice site, so that, as n<L, the distance between two nearest-neighbour cells is several lattice sites. Each cell i=1…n is associated with a dynamical variable φ i (t) that takes integer values between 0 and K-1 and represents the state of cell i at simulation time step t. Upon stimulation by a messenger molecule, a cell enters a refractory phase which lasts K-1 time steps and during which the cell is not responsive to messenger molecules (see figure 1). The state φ i = 0 is the quiescent state. The transition from φ i = 0 to φ i = 1 represents the stimulation step and occurs with probability ω (per time step) if a messenger molecule is present on cell i lattice site. Transition from the state φ i = j to φ i = j+1 then occurs with probability 1 (at each time step) for j = 1...K-1. After K steps, the cell is back in the quiescent state i.e. goes from φ i = K-1 to φ i = 0 with probability 1. Finally, the transition from φ i = K p to φ i = K p +1 (K p ∈ 1…K-1) is accompanied by the release of a new messenger on cell i lattice site.

c) Initialization

The initial state of the system corresponds to "full lattice" initial conditions: one messenger molecule is initially positioned on each lattice site and the state of each cell φ i (0), ∀i=1…n, is randomly chosen from a uniform distribution between 0 and K-1.

d) Simulations

Numerical simulations basically implement the above listed rules (Fig. 1 , sections 2.a. & 2.b). For clarity, we present below the translation of these rules in the point of view of the simulation algorithm.

After initialization, the system is updated as follows at each simulation step t>0:

• The state φ i of each cell i is first updated. If a messenger molecule is present on cell i lattice site and if i is quiescent (φ i (t-1) = 0), its state is updated to φ i (t) = 1 with probability ω. Note that, in the presence of m messengers on the same cell site, the effective stimulation probability is

P(φ i = 0 → φ i = 1) = 1-(1 -ω) m . Cells in the refractory phase (φ i (t-1)>0) are updated according to φ i (t) = φ i (t-1)+1 or φ i (t) = 0 if φ i (t-1) = K-1.
Finally, new messengers are created at each cell site for which the updated state variable φ i (t) = K p +1.

• The messenger molecules are then updated according to two substeps.

Each messenger is first independently and simultaneously moved to one of its two nearest-neighbour sites (chosen at random). Each molecule is then removed from the lattice with uniform probability λ.

e) Parameter values

Unless otherwise stated, standard parameter values used in the present paper were: L=500×10 3 sites, n = 50×10 3 cells, K = 20 and K p = 10. For each parameter set, the results presented are averages over 10 different realizations of the initial conditions and of the simulation run. The system behaviour was mainly studied through variations of the stimulation probability ω and the messenger removal probability λ.

General behaviour a) Mean-field analysis

To predict the behaviour of the model, a first approach neglects the inherently stochastic nature of the model and assumes that the fluctuations of the number of messenger molecules on the lattice are not significant. This so-called "mean-field" approach fundamentally approximates the model as a deterministic process, using differential equations that are actually similar to the mass-action laws used in classical (bio)chemical kinetics (see [START_REF] Berry | Nonequilibrium phase transition in a self-activated biological network[END_REF] for further details). Under these assumptions, the system is predicted to asymptotically reach a steady state (i.e. a stable fixed point). The messenger density ρ (= m/L where m is the total number of messenger molecules on the lattice) in the steady state is given by:

     ≥         - - < = c c c ss if K if ω ω ω ω ω ω ρ 1 1 1 1 0 (1) with n L c λ ω = (2)
Hence, mean-field predictions state that, if ω<ω c , messenger molecules are eventually cleared from the lattice (while every cell eventually becomes quiescent).

In a way, this corresponds to the death of the system and will be referred to as the "absorbing phase". Conversely, if ω≥ω c , messenger molecules persist indefinitely on the lattice. This phase will be referred to as the "active phase".

When ω is close to ω c , a Taylor series expansion at first order of (1/ω c -1/ω) is (ωω c )/ω c 2 . Thus Eq. ( 1) can be rewritten as: It must however be emphasized that the mean-field analysis above is valid only if the spatial fluctuations of messenger molecules are negligible. Yet simulations of the model such as those presented in Fig. 2 show that at every time step, the distribution of messenger molecules is rather inhomogeneous in space, with "white" zones (devoid of messengers) of various sizes and "black" regions where messenger density is high. This is a fundamental expression of the inherently stochastic nature of the model. It indicates that spatial fluctuations in the system are high, and could invalidate some of the predictions made by the deterministic mean-field analysis. To quantify and characterize the behaviour of the system with its natural stochasticity and fluctuations, intensive simulations of the model are needed. The following of the paper is devoted to these analyses.

H BERRY

Simulation results

a) Dynamics of the messenger density

Figure 3 presents the time evolution for the density of messenger molecules ρ in the active (Fig. 3A) or absorbing (Fig. 3B) phases. In the active phase, ρ decreases according to three main regimes. After an initial slow decay regime (up to ≈20 time steps), ρ decreases abruptly for ≈10 3 time steps. Finally, at long times, ρ dynamics reaches a regime where it fluctuates around a constant average value.

Fundamentally, this stationary fluctuating regime is the stochastic equivalent of the steady state predicted by the deterministic mean-field analysis (Eq. 1). The average value of the fluctuations at long times in the simulations will thus as well be noted ρ ss . Note however that, while both simulations and mean-field analysis indicate a stationary state in the active phase at long times, the predicted values are different.

For instance, with the parameters of Fig. 3, Eq. 1 predicts ρ ss = 0.276 (for ω = 0.7) and 0.292 (for ω = 0.9) while simulations yield 0.131 and 0.186, respectively. This is a first evidence for failure of the mean-field analysis due to the high fluctuations in the model.

In the absorbing phase (Fig. 3B), the behaviour at long times is quite contrasted.

Instead of reaching a stationary regime, messenger density vanishes very fast.

Indeed, this curve is plotted in log-linear coordinates so that the observed linear behaviour indicates an exponentially fast decay.

Hence, in the active phase at long times, ρ settles onto a fluctuating stationary regime, while in the absorbing phase, messengers vanish according to an exponential decay. The difference between these two behaviours can be used to identify the boundary between the two phases. Figure 4A shows the dynamics of ρ for values of ω that are close to the transition between active and absorbing phases.

For the smallest values of ω (the bottom curves in the figure), the curves bends down. Note that the curves are plotted here in log-log coordinates so that this downward curvature reflects the exponential decay characteristic of the absorbing phase. Conversely, the curves obtained for the largest values of ω (the upper curves in fig. 4A) inflect upwards and tend towards the stationary state, which is a characteristic of the active phase. At the boundary between the two phases, i.e. for the critical value ω c , the decay of ρ at long times is thus expected to be linear in log-log coordinates (neither curving upwards nor downwards), which corresponds to the power-law:

( ) α ρ - ∝ t t ,
for ω = ω c , and t→∞, (4) Hence, the value of ω that yields a power-law behaviour can be used as an estimate for ω c . Using this principle, Figure 4B shows that an estimate of ω c = 0.6230 is obtained. Furthermore the estimated value for the exponent of the corresponding power-law in this case is α ≈ 0.165. It must be noted that the prediction of the H BERRY mean-field analysis fails here. In the case of Fig. 4, Eq. 2 predicts ω c = 0.150, which is less than one fourth of the value obtained in the simulations. Averaging over these simulations, an improved estimate for the exponent α (Eq. 4) is obtained: α = 0.159(7) (the number between parentheses indicates the estimated uncertainty on the last digit). However, a fundamental difference between the phase transition observed here and classical phase transitions lies in its nonequilibrium nature. Indeed, most of the classical phase transitions (including the ferromagnetic transition) are equilibrium phenomena, i.e. they obey the so-called "detailed balance" condition [START_REF] Hinrichsen | Non-equilibrium phase transitions[END_REF]: the transition rate (or more precisely fluxes) from any state i to any state j is equal to that of the reverse transition rate, from j to i. This condition implies the presence of real equilibrium states and is a fundamental property of these systems.

In particular, it allows applying usual thermodynamic tools and concepts. In the system studied here, the absorbing state induces a violation of the detailed balance condition: this state can be reached with some probability from the active state, but cannot be escaped. In turn, this prevents the application of several of the classical thermodynamic concepts. The general properties of these kinds of phase transitions, called nonequilibrium phase transitions, are thus still a mater of debate in the physics community, but several points are well established (for reviews see [START_REF] Hinrichsen | Non-equilibrium phase transitions[END_REF][START_REF] Odor | Universality classes in nonequilibrium lattice systems[END_REF]or Hinrichsen, 2000a).

Like their equilibrium counterparts, nonequilibrium phase transitions exhibit generic power-law behaviours close to the critical threshold. Eq. 4 is a first example of such a power-law, but we will encounter more of them in the 1 Note that the apparent discontinuity close to the critical value in Fig. 6 is a numerical artefact due to the so-called "critical slowing down". Indeed, as the critical threshold is approached, the fluctuating stationary regime starts at increasingly long simulation times. As this regime has to be reached to estimate the corresponding value of ρ ss , the number of time steps in the simulations must be substantially increased as the critical threshold is approached. Hence the apparent discontinuity of Fig. 6 only reflects computational limitations hindering measurements with ω values closer to the threshold.

H BERRY

following. Several quantities found in the power-laws depend on details of the simulation and the system considered: they are called "non-universal". For instance, the value of the critical threshold ω c is known to be non-universal, and thus depends on the probability λ for instance, as seen in Fig. 5A (just like the value of the percolation threshold depends on the type of percolation model considered, i.e. bond or link percolation). However, close to the critical threshold, several quantities (called universal quantities) are largely insensitive to microscopic details of the system (as in equilibrium phase transitions). For instance, the value of the exponent α (the critical exponent of power-law Eq. 4)

does not vary when λ changes (Fig. 5B).

One great success of the theory of equilibrium phase transitions is the explanation that transitions arising in different systems can share the same set of critical exponents. This phenomenon is called "universality". For example, the critical exponents at the liquid-gas critical point are independent of the chemical composition of the liquid. Furthermore, several models for transport in porous media (Sahimi, 1994) share the same critical exponent values as models of sol-gel transitions [START_REF] Adam | Sol-gel transition[END_REF]: they are in the same universality class (that of percolation, in this case). This means that the dominant processes implied in these systems are explainable by the same fundamental mechanisms.

Nonequilibrium phase transitions may as well be classified into several universality classes. Among them, the universality class of Directed Percolation (DP) has proven the most robust. Apparently very different systems, ranging from population dynamics (Lipowski & Lipowska, 2000), epidemic spreading [START_REF] Dammer | Epidemic spreading with immunization and mutations[END_REF], forest fires [START_REF] Albano | Critical behaviour of a forest fire model with immune trees[END_REF], biological evolution [START_REF] Ferreira | Nonequilibrium phase transition in a model for the origin of life[END_REF], Ca 2+ signalling [START_REF] Timofeeva | Directed percolation in a two-dimensional stochastic fire-diffusefire model[END_REF][START_REF] Bär | Discrete Stochastic Modeling of Calcium Channel Dynamics[END_REF] to morphology dynamics of growing surfaces (Hinrichsen, 2000a), belong to the DP universality class. This indicates that, beyond their diversity, all these systems rely fundamentally on the same ground processes. However, several other universality classes for nonequilibrium phase transitions have been uncovered, such as the parity conserving universality class [START_REF] Zhong | University class of two-offspring branching-annihilating random walks[END_REF], or the conserved threshold transfer process universality class [START_REF] Lübeck | Universal finite-size scaling behavior and universal dynamical scaling behavior of absorbing phase transitions with a conserved field[END_REF].

Hence, determining the universality class of a nonequilibrium phase transition is an important step towards the understanding of its fundamental mechanisms.

Practically, this consists in estimating the value of the critical exponents of the system. This task is carried out in the following sections.

c) Universality class of the phase transition

Besides Eq.4, an important power-law for nonequilibrium phase transitions relates the average value of the density of messengers in the fluctuating stationary state, ρ ss , to the distance to the critical threshold, ω-ω c [START_REF] Hinrichsen | Non-equilibrium phase transitions[END_REF]:

( ) β ω ω ρ c ss - ∝ , for ω → ω c in the active phase (5)
The corresponding critical exponent β is a universal quantity. Actually Eq. 5 is identical to Eq. 3 with β = 1, which means that the mean-field analysis predicts β = 1. Figure 6 presents simulation results showing how ρ ss varies as a function of the distance to the threshold for various values of λ. These curves show roughly parallel straight lines at small distances, indicating a power-law behaviour with a critical exponent that does not depend on λ. Averaging over the 6 conditions of the figure, we obtain an estimate of β = 0.29(1). Here again, note that the value predicted by the mean-field analysis is more than threefold the measured one. Figure 8 shows the results obtained for λ = 0.0150 (ie Fig 4A). Using the estimated value of α = 0.159 (see above), the best collapse of the curves (at long times) was H BERRY found for ν || ≈ 1.7. Similar treatments of the data for other values of λ yielded comparable quantities (not shown). On average, the estimate for the critical exponent is ν || = 1.66(9).

A third important critical exponent is the so-called dynamic exponent z. Roughly speaking, this exponent relates the behaviour of the system to the size of the spatial domain (i.e. the number of sites L in the lattice in our case) [START_REF] Hinrichsen | Non-equilibrium phase transitions[END_REF]. It is important to realize that the various power-laws mentioned above are strictly valid for infinite spatial domains only, which is not accessible with computer simulations. The lattice size used in the simulations presented up to this point (L = 5×10 5 ) was large enough so that the effects of this finite size could be neglected. In fact, for simulation durations of 10 6 time steps, we have found that finite size effects are negligible as soon as L ≥ 5×10 3 sites. Smaller L values however evidence drastic effects. Figure 9A shows the dynamics of ρ for small lattices (L ≤ 10 3 ) at criticality (ω = ω c ) for λ = 0.0150 and constant cell density (i.e. the number of cell is n = 0.10×L). In this case, the system has always a certain probability to jump into the absorbing phase (the corresponding curves bend downwards) even though the parameters correspond to the critical threshold (of course, this probability decreases when lattice size increases, and becomes very low for large but finite systems). Theoretical arguments [START_REF] Hinrichsen | Non-equilibrium phase transitions[END_REF] predict that in this case, plots of ρ(t) × t α as a function of t/L z should yield a collapse of the curves shown Fig. 9A on a single curve. As with ν || above, this gives a method for estimating the dynamic exponent z. Fig. 9B shows the best collapse, obtained with z = 1.55. The results of Fig. 9 can be replicated using other values for λ, and yield comparable estimates for z (not shown). The average value obtained is z = 1.54(5).

Table 1 summarizes the values obtained for the four critical exponents estimated in this study. It also indicates an estimate for the critical exponent related to the spatial fluctuations, ν ⊥ . This exponent was not measured directly but deduced from the definition z = ν || /ν ⊥ [START_REF] Hinrichsen | Non-equilibrium phase transitions[END_REF]. The table also indicates the corresponding values for the directed percolation (DP) universality class in one space dimension. Clearly, for each exponent, the obtained estimates are in good agreement with those of DP. This demonstrates that the nonequilibrium phase transition observed in the present model belongs to the directed percolation universality class. [START_REF] Hinrichsen | Non-equilibrium phase transitions[END_REF]). Note that the value of z for the present model is not a measure, but is deduced from z = ν || /ν ⊥ . Numbers in parentheses indicate uncertainty in the last digit.

H BERRY

Message Broadcasting

The characterization of the phase transition allows a better understanding and prediction of several aspects of the model, such as the characteristics of message broadcasting. Let us imagine that the lattice is in the absorbing (dead) phase, i.e.

devoid of any messenger molecule, with all cells quiescent. At time t 0 , a single messenger molecule is produced by a cell at a random position on the lattice (see corresponding simulations in Figure 10). Will this embryo of a message persist over time? Will it finally be transmitted to every cell in the lattice? In physicists' terms, this raises the question of the correlation length in the system. It turns out that, unlike their equilibrium counterparts, nonequilibrium phase transitions possess two independent correlation lengths: a spatial correlation length ξ ⊥ , and a temporal correlation length ξ || (Hinrichsen, 2000a;[START_REF] Hinrichsen | Non-equilibrium phase transitions[END_REF]. These lengths are directly related to the spreading of the message seed. In the absorbing phase, the message survives for a certain amount of time and tends to slowly expand in space (Fig . 10A). The theory of nonequilibrium phase transition tells us that the average survival time is given by (or at least proportional to) ξ || ∝ |ω -ω c | -ν|| and the maximal expansion of the message in space by ξ ⊥ ∝ |ω - (Hinrichsen, 2000a). Similarly, in the active phase, message spreads as a cone (Fig. 10B) whose opening angle is given by ξ ⊥ / ξ || . Using these theoretical concepts, it becomes easy to understand that, as the critical threshold is approached from below in the simulations, the message spreading from the initial seed survives for increasing times and reaches a growing fraction of the cells. At the critical threshold, the effect of a single seed is even eventually transmitted (relayed) all along the lattice, whatever the physical distance to the initial seed site (both ξ || and ξ ⊥ →∞). As ω increases further (above the threshold) and enters the active phase, the opening angle of the spread cone grows as |ω -ω c | ν||-ν⊥ . Thus, as ν || -ν ⊥ = 0.637 > 0 in DP, the time needed for the message to invade the whole lattice decreases with increasing ω. Hence identifying the phase transition allows predicting both qualitatively and quantitatively the behaviours shown by the simulations.

ω c | -ν⊥

Discussion

The idea that the stochasticity inherent to cellular processes has a strong influence on elementary behaviours such as toggle switch properties (Tian & Burrage, 2006;[START_REF] Vilar | Modeling network dynamics: the lac operon, a case study[END_REF], clock oscillations [START_REF] Suel | An excitable gene regulatory circuit induces transient cellular differentiation[END_REF] or signal transduction [START_REF] Bhalla | Signaling in small subcellular volumes. I. Stochastic and diffusion effects on individual pathways[END_REF], has recently emerged as a new paradigm in cell biology (for a review, see [START_REF] Rao | Control, exploitation and tolerance of intracellular noise[END_REF]. Furthermore, the role played by noise and stochasticity at the cell or organism level is beginning to be realized (for recent reviews see [START_REF] Samoilov | From fluctuations to phenotypes: the physiology of noise[END_REF][START_REF] Raser | Noise in gene expression: origins, consequences, and control[END_REF].

The results of the present work add auto/paracrine-like relay signalling to the list of basic cell processes that may be highly influenced by noise-induced fluctuations.

Our main objective was to study the influence of noise in a model for stochastic communications where cells respond to a diffusive environment that they themselves have produced through a positive feedback loop. Using numerical simulations, we showed that the deterministic (mean field) approximation of this process fails short of predicting the behaviour of the system. Instead, the latter could be understood as a nonequilibrium phase transition that was evidenced to be in the universality class of directed percolation.

That the model belongs to the directed percolation universality class is actually not totally unexpected. Indeed, DP has been conjectured [START_REF] Grassberger | On phase transitions in Schlögl's second model[END_REF] to be the universality class of models exhibiting a phase transition from a fluctuating active phase to a unique absorbing phase; whose dynamic rules involve only short-range processes; and that have no unconventional attributes such as additional symmetries (the actual conditions are in fact a bit more restrictive, see Hinrichsen, 2000a). The present model verifies these conditions thus supporting this conjecture.

Modelling studies of auto/paracrine relay transmission have rarely been addressed so far. In two studies, Pribyl and coworkers (Pribyl et al. 2003a;2003b) proposed a mechanistic model that presents similarities with the model studied here. Their model takes messenger diffusion, reversible binding to cell surface receptors and a positive feedback from ligand binding to new ligand release. The main differences with our model are that the feedback mechanism in their model is slightly more H BERRY detailed (from a biochemical perspective) and most notably, their model is purely deterministic. Hence, from a stochastic point of view, this model could be considered as a mean-field limit. These authors show that their model supports travelling wave solutions that spread throughout the cell layer. In these waves, the front connects a steady state with nonzero messenger density to a steady state devoid of messengers. These behaviours are at first sight reminiscent of the spreading waves encountered in directed percolation, for instance for "seed" initial conditions (see section 5), and of the active and adsorbing stationary states.

Whether this model is related to directed percolation is however not trivial and would necessitate further works. For instance, the spreading waves in DP are not travelling waves (i.e. the shape of the wave front changes during its propagation).

Starting from the discrete version of their model (Pribyl et al., 2003b), biologically realistic stochastic processes could be introduced at several points. Possible phase transitions and universality class of the resulting model could then be assessed through the study of the time evolution of the probability P(t) that a spreading waves survives up to time t, that scales in the subcritical regime as P(t) ~ t -δ with δ=0.1595 or 0.451 for DP in d=1 or 2 dimensions, for instance.

Biological interpretation of the present work must be handled with great care.

Indeed, the aim of the present work was to study the basic mechanisms at play in such kinds of diffusive communications, so that we voluntarily restricted the model to a set of a few elementary ingredients. Of course, this limits the biological realism of the model, and raises the question of the persistence of the observed behaviours in more detailed and biologically plausible models. While this question can only be answered by further modelling and simulation works, we give here some remarks related to this point.

First of all, the universality class of a phase transition is, by definition, not changed by details of the model. And directed percolation, in particular, has proven a very robust universality class to this respect. Hence, it seems unlikely that minor modifications in the model (such as changing the lattice geometry or using offlattice conditions, using stochastic refractory periods for the cells…) would modify drastically the observed behaviour. For instance, preliminary investigations showed that changing the values of the cycle parameters (K and K p ) modifies the value of the critical threshold but not the transition itself. Actually, shorter cycle lengths K (not shown) decrease the duration of the initial phase i.e. the non power-law regime of ρ dynamics (the first 10 3 time steps in Fig. 4B).

An important parameter is the Euclidean dimension of the lattice. We choose here to study a one-dimensional model because 1-D models are much cheaper in computation time, thus allowing intensive simulations. However, biological realism would rather impose two-(or three-)dimensional spaces. Actually a twodimensional version of the model presented here has previously been studied in [START_REF] Berry | Nonequilibrium phase transition in a self-activated biological network[END_REF]. However, because simulations in this previous paper were not refined enough, the universality class could not be determined unambiguously.

Hence, future work will focus on intensive simulations of the model in two and three space dimensions.

Furthermore, another significant assumption in the present work is the regular spacing of the cells on the lattice. However, preliminary simulations (not shown)

indicate that random cell positioning does not modify the occurrence of the nonequilibrium phase transition, but has the deleterious effect that the density of messenger molecules decays as a power-law in the absorbing phase (and not exponentially, as observed for regular spacing)2 . This invalidates the estimation of the critical threshold based on the criterion of regime change between subcritical and critical conditions, that was used in Fig. 4, for instance. While more computationally expensive, other methods based on finite-size effects for instance [START_REF] Ortega | Critical behavior of a probabilistic cellular automaton describing a biological system[END_REF][START_REF] Binder | Monte Carlo Simulation in Statistical Physics[END_REF] can however be applied in this case and will be used in future studies of the model with random cell positioning.

A surprising point concerning directed percolation is that in spite of its amazing robustness in models/simulations, this critical behaviour has still not been evidenced in experiments (see Hinrichsen, 2000b). Hence, the possibility that autocrine relays may provide experimental evidence of DP should be considered seriously. Possible experimental setups could for instance consist in cell and tissue culture assays, where a confluent (or subconfluent) layer of autocrine cells is covered by a liquid medium, in which the soluble messenger is secreted and diffuses. Alternatively, such a situation is naturally found in Drosophilia eggs (see Pribyl 2003b). The evolution of messenger concentration in the liquid or of related reporters could be monitored in situ, for instance using fluorescence. Such in situ recordings of spatio-temporal evolutions are already accessible to cell biologists (see e.g. [START_REF] Nikolic | Role of boundary conditions in an experimental model of epithelial wound healing[END_REF]. Now, messenger death probability (λ in the present model) could be modulated through the addition of various levels of proteinases (degrading the messenger) in the liquid medium. Likewise, the stimulation probability ω could be modulated through the addition of various levels of antagonists of the messenger receptor, or intracellular inhibitors of the signalling pathway supporting the positive feedback. Finally, the A431 carcinoma cell line could be a good candidate for these studies [START_REF] Dent | Radiation-induced release of transforming growth factor alpha activates the epidermal growth factor receptor and mitogen-activated protein kinase pathway in carcinoma cells, leading to increased proliferation and protection from radiation-induced cell death[END_REF].

Figure 1 :

 1 Figure 1: Overview of the cell activation cycle. In the sequence, the cell state φi is indicated by a bold number below the cell. Transition probabilities (at each time step) between two cell states are indicated by italicized numbers below corresponding arrows. Intervening messenger molecules are symbolized by white triangles. During the refractory phase of the cycle (light grey background), the cell is not sensible to messenger molecules.

  Figure2presents simulation results illustrating the occurrence of the absorbing and active phases. In this figure, the x-axis is the 1-D lattice. A black dot is drawn at each lattice site containing a messenger molecule. The evolution of the system in time is plotted along the y-axis. When ω is small (Fig.2A), messenger molecules persist for circa 100 time steps, before undergoing massive extinction. Eventually, after ≈ 600 time steps, all the messengers have been removed from the lattice, and the system freezes into the absorbing (dead) phase. Conversely, for large ω values (Fig.2B), the messenger molecules survive for the entire time span of the simulation (active phase). Albeit the system is constantly changing, visual inspection of the figure indicates that messenger density tends to a roughly constant value.

Figure 2 :

 2 Figure 2: Simulations of the model illustrating the occurrence of an absorbing (A) and an active (B) phase. Parameters were: L=400 sites, n = 40 cells, K = 20, Kp = 10, λ=0.0150 and (A) ω=0.100 or (B) ω=0.6235

Figure 3 :

 3 Figure 3: Dynamics of the messenger density ρ(t). Simulation results for "Full lattice" initial conditions in the active (A) or absorbing (B) phases. In (A) the dashed lines indicate the corresponding average value in the fluctuating stationary state ρ ss . In (B) dashed lines represent fits with exponentially decreasing functions. The values of ω are indicated beside the curves; λ=0.0150; L=500×10 3 sites, n = 50×10 3 cells, K = 20 and K p = 10. Results are from a single simulation run.

Figure 4 :

 4 Figure 4: Dynamics of the messenger density ρ(t) close to the boundary between active and absorbing phases. (A) Each curve corresponds to a different value of ω, i.e., from bottom to top, ω =0.6220, 0.6225, 0.6230, 0.6235, 0.6240, 0.6260, 0.6270, 0.6280 and 0.6290. (B) Time-variation of ρ when ω is set to its estimated critical value ω = ω c = 0.6230(5). (The number between parentheses denotes the estimated incertitude on the last digit). The dashed curve indicates a power-law decay with exponent -0.165. Parameters were λ=0.0150 and other parameters from the standard set (c.f. 2.e).

  The critical behaviour for various values of the messenger death (removal) probability λ is shown in Figure5. Simulations of the model show that the critical value ω c increases exponentially fast with λ (Fig.5A, open circles). Here again, the prediction of the mean-field analysis (Fig.5A, dashed line) clearly underestimates the actual value of the critical threshold. The dynamics of ρ for various values of λ at criticality (i.e. setting ω to its critical value ω c at the corresponding value of λ) are plotted Fig.5B. Albeit the final power-law regimes start at longer times when λ decreases, all the curves obtained eventually reach a power-law regime with similar value of the exponent α (i.e. parallel straight lines in log-log coordinates).

Figure 5 :

 5 Figure 5: Critical behaviour for various values of the messenger death probability λ. (A) Estimates of the critical stimulation probability ω c as a function of λ (open circles). The dashed line indicates prediction from the mean-field analysis (Eq. 2). The full line is an exponential fit. (B) Dynamics of the messenger density ρ(t) for various values of λ at the corresponding critical value ω c . λ values are (from top to bottom): 0.0050; 0.0060; 0.0100; 0.0125;0.0150 and 0.0165. The dashed line indicates a power-law decay with exponent 0.159. All other parameters according to the standard values (section 2.e).

Figure 6 :

 6 Figure 6: Variation of the messenger molecule density in the (fluctuating) stationary state, ρ ss , as a function of the stimulation probability ω. Open circles show the results from the simulations of the models (the full line is a guide for the eyes). Parameters were λ=0.0150, and standard values for the others. Error bars show ± 1 s.d.

Figure 7 :

 7 Figure 7: Scaling of the messenger density in the fluctuating stationary states with the distance to the threshold for different values of the messenger death probability λ. The dashed line indicates a power law with exponent 0.29. All other parameters are set to their standard values. Error bars show ± 1 s.d.

Figure 8 :

 8 Figure 8: Scaling test for the value of the critical exponent ν || . The data of Fig.4A are plotted here as ρ(t) × t α as a function of t ×(ωω c ) ν|| with α=0.159 . The best collapse of the curves is obtained for ν || =1.70. Parameters: λ=0.0150 and all other parameters according to their standard values.

Figure 9 :

 9 Figure 9: Finite size scaling test at criticality (ω = ω c ) and λ=0.0150. (A) Dynamics of the messenger density ρ(t) at criticality for different sizes of the lattice L = 100; 200; 500;800; and 1,000 (from bottom to top). Corresponding values of the number of cells n = L/10 (constant cell density). (Β) Scaling test for the value of the critical exponent z. The data of (A) are replotted here as ρ(t) × t α as a function of t/L z with α=0.159 . The best collapse of the curves is obtained for z=1.55. Parameters: λ=0.0150 and all other parameters according to their standard values.

Figure 10 :

 10 Figure 10: Simulations of the model illustrating the spreading of a unique initial messenger "seed" in the absorbing phase (A, ω =0.50 ) or in the active phase (B, ω = 0.95). Parameters were: L=400 sites, n = 40 cells, K = 20, Kp = 10, λ=0.0150. Each cells were initiated in the quiescent state (φ i (0) = 0;∀ i = 1…n).

Table 1 :

 1 Comparison of the estimates obtained for the critical exponents studied in the present model and those of Directed Percolation (DP) in one space dimension (values are from

	Model	α	β	ν ||	z	ν ⊥
	Present 0.159(7)	0.29(1)	1.66(9)	1.54(5)	1.08
	DP	0.159464(6)	0.276486(8)	1.733847(6)	1.580745(10) 1.096854(4)

Note that this phenomenon has already been encountered in the literature; see[START_REF] Szabo | Generalized contact process on random environments[END_REF] for instance.