S Cronenberger

H Boukari

D Ferrand

J Cibert

D Scalbert

Long-range spin jump diffusion revealed by dynamic light scattering

Long-range spin jump diffusion revealed by dynamic light scattering

S Cronenberger, H Boukari, D Ferrand, Joel Cibert,

I. INTRODUCTION

Donor electron spins in semiconductors are being actively studied as they hold great promise in the context of quantum technologies [1][2][3]. Such spins are among the most coherent quantum systems in the solid state and can be efficiently exploited by using advanced semiconductor technology [4,5]. Long-distance displacement of individual electron spins is another key ingredient for scalability in a spin-based semiconductor quantum circuit. This is one option to perform long-range interaction between distant spin qubits [6,7]. It has been proposed to move the electron spin coherently over large distances in a well-defined tunnel coupled array of empty donors or dots. In practice, this has been realized in quantum-dot arrays containing up to four electrostatic dots in GaAs/GaAlAs heterostructures, with size not exceeding a few micrometers [8][9][10].

In this paper, we demonstrate spontaneous spin jumps between donors separated by a few micrometers in a many-body interacting disordered spin system formed by isolated donors. Although the jumps are, in this situation, uncontrolled, understanding the underlying mechanisms will open new prospects for coherent long-distance displacement of spins in semiconductors and help to determine to what extent unwanted jumps can be suppressed.

Experimentally, we reveal the jump-diffusion regime by passively probing the stochastic spin fluctuations at thermal equilibrium with off-resonant light. This noise-based approach, issuing from atomic physics [11][12][13] and quantum optics [14][15][16], has been also actively developed in semiconductors since 2005 [17,18]. In these materials, homodyne and heterodyne detection opened the route to nearly perturbationfree measurements [19][20][21]. Different theoretical proposals have pointed out the interest in probing spatial spin correlations both in semiconductors [22][23][24] and in cold fermionic gases [25]. Recently, we experimentally realized spatiotemporal spin noise spectroscopy applied to semiconductors (ST-SNS) and obtained the wave-vector-resolved dynamical spin response function of electrons [26,27]. Here, we further extend the accessible range of wave vectors q in the spirit of dynamic light scattering. This allows us to go beyond the q values which have been explored so far [28][29][30][31][32][33]. This way, we clearly reveal a deviation from the classical q 2 diffusion law and observe a plateau at large q; these findings are the signature of the jump-diffusion regime.

Jump diffusion is a stochastic process that involves both discrete jumps and continuous motion. It finds applications in many different fields ranging from physics [34,[START_REF]Diffusion in Condensed Matter[END_REF] to epidemiology [START_REF] Brauer | Mathematical Models in Population Biology and Epidemiology[END_REF], biology [37,38], and risk management in financial markets [39,[START_REF] Cont | Financial Modelling with Jump Processes[END_REF]. In condensed matter physics, jump-diffusion processes control, e.g., the atomic diffusion in crystals and in some liquids, such as water and lead [41,42]. Because jump diffusion normally occurs at the atomic scale, it has been extensively studied by means of quasielastic neutron scattering [43,[START_REF] Zaliznyak | Modern Techniques for Characterizing Magnetic Materials[END_REF]. We rely on this large body of work to develop below a model which gives the first complete description of spatiotemporal spin dynamics in an ensemble of disordered localized spins.

Among the presently investigated materials it has been suggested that wide-gap II-VI semiconductors combine most of the advantages of solid-state systems [START_REF] Greilich | [END_REF][46][47]: high natural abundance of isotopes with zero nuclear spins, which also permits isotopic purification and reduction of decoherence of electron spins due to hyperfine interaction [48]; a direct band gap allowing for an efficient spin-photon entanglement; and large binding energies which confer high thermal stability to donor-bound electrons. Long spin coherence T 2 times up to 50 μs, a prerequisite for quantum computing, have been measured, for example, in ZnO [47] and could be further extended with isotope purification. In fluorine-doped ZnSe, inhomogeneous spin dephasing times up to 33 ns were observed below 40 K [START_REF] Greilich | [END_REF].

We demonstrate jump diffusion in an insulating n-CdTe sample by resolving effective spin jumps over distances ∼ 2.7 μm. ST-SNS alone is not sufficient for this task because high spatial frequencies are cut, due to the finite resolution of the imaging system [26]. This hurdle is bypassed if the dynamically scattered light is collected in an arbitrary direction. Here, we combine three detection schemes: classical SNS (q = 0 in forward scattering), ST-SNS (q 3 μm -1 ), and backscattering SNS (q b = 4π n/λ 50 μm -1 , with n being the refractive index and λ being the wavelength). High spatial frequencies up to q b become accessible, corresponding to a spatial resolution ∼ λ/10.

In the next section we present our model for the wavevector-resolved spin noise spectra. There, we derive a phenomenological equation, which incorporates the structure factor necessary to account for the light scattering in an arbitrary direction, and the different spin decoherence mechanisms. Section III is devoted to the experimental results and their analysis with the model developed in Sec. II. Finally, we conclude in Sec. IV.

II. WAVE-VECTOR-RESOLVED SPIN NOISE SPECTRA IN THE RANDOM-WALK APPROACH

Our calculation is inspired by previous theories [41,42,[49][50][51][52]. In the same spirit, we compute the spin polarization decay within a random-walk approach and express it as an infinite sum over the number of spin jumps between localized sites. Our approach, unlike theories aimed at calculating the spin resonance or spin noise line shapes, introduces the optical phase factor due to light scattering at a finite angle.

A. Phase coherence decay

The decay of phase coherence of electron spin rotation is calculated in the presence of jump diffusion taking into account the jump length distribution.

The random-walk description

We first recall the main lines of the random-walk model for the decay of phase coherence of rotating spins which perform jumps in an ideal crystal lattice [50]. The spins precess with the Larmor frequency ω 0 under the influence of an external magnetic field. The phase coherence of the rotation of different spins is gradually diminished through the random orientation of nuclear fields at different sites. In the presence of spin jump diffusion, the individual spins sense more or less averaged nuclear fields, depending on the jump rate. The decay of the polarization is reduced compared with the case of no diffusion; that is, motional narrowing occurs. The basic feature of the model is a separation of the polarization decay P(t ) into the contribution of a fixed number of generalized jumps,

P(t ) = ∞ l=0 R l (t ).
(1) R l (t ) is the contribution to the polarization decay of all spins which have made l jumps until the time t. The initial condition is P(t = 0) = 1. In particular, we have

R 0 (t ) = P 0 (t ) exp(-t/τ c ), (2)
where τ c is the residence time at each site and exp(-t/τ c ) is the probability that the spin has performed no jump until t. P 0 (t ) = exp(-1 2 σ 2 t 2 ) is the decay of the phase coherence due to the Gaussian distribution (σ 2 second moment) of spin precession frequencies at different localized sites and given by the ensemble average over many spins. If we assume no correlations between the precession frequencies at different sites, the ensemble average can be performed at each jump, and the decay of phase coherence at each step can be calculated as a function of R 0 (t ); we have, after one jump,

R 1 (t ) = t 0 dt R 0 (t -t ) 1 τ c R 0 (t ), (3)
and after l jumps

R l (t ) = t 0 dt R 0 (t -t ) 1 τ c R l-1 (t ), (4)
and the series in Eq. ( 1) is summed up by the integral equation

P(t ) = R 0 (t ) + 1 τ c t 0 dt R 0 (t -t )P(t ). (5)
The integral equation can be solved by Laplace transformation

P(s) = R0 (s) ∞ l=0 1 τ c R0 (s) l = R0 (s) 1 -1 τ c R0 (s) , ( 6
)
where P(s) and R0 (s) are given by

P(s) = ∞ 0 dt exp(-st )P(t ), (7)
R0 (s) = ∞ 0 dt exp(-st )R 0 (t ) ( 8 ) = √ π √ 2σ exp (s + 1 τ c ) 2 2σ 2 erfc s + 1 τ c √ 2σ . ( 9
)
R0 (s) can be conveniently expressed with the Faddeeva function, for which efficient implementations exist for numerical calculations,

W (z) = e -z 2 erfc(-iz) (10) so that R0 (s) = √ π √ 2σ W i(s + 1 τ c ) √ 2σ . ( 11
)
P(s) is related to the spin noise spectrum S(ω),

S(ω) = 1 π Re{ P[s = i(ω -ω 0 )]}. ( 12
)
An angular rotation R so occurs from the jump event when spin-orbit coupling is present. In the case where all angular rotations due to spin orbit are the same at each spin jump, Eq. ( 6) becomes [51]

P(s) = R0 (s) ∞ l=0 R l so 1 τ c R0 (s) l (13) = R0 (s) 1 -R so 1 τ c R0 (s) , ( 14
)
where R so = 1 -4 3 sin 2 ( γ 2 ) for the spatially averaged spinorbit rotation.

One can incorporate an intrasite relaxation independent of hopping by multiplying the right-hand side of Eq. ( 2) by exp(-ν s t ), with ν s being the intrasite spin relaxation rate. Thus, with an intrasite relaxation, R0 (s) becomes R0 (s + ν s ).

Optical phase factor

In the following we will introduce the effect of a spatial phase factor which appears when the spin system is probed at scattered angles, i.e., when the signal is detected in a direction different from that of the incident light. This kind of configuration of noise measurement is usually called dynamic light scattering (DLS) measurement. The spatial phase factor comes from "interferences" between the incident light and the scattered light; it is then equal to the relative optical phase between the incident and scattered fields. A particle at a given position r i scatters the incident plane wave, initially propagating in the direction u 0 , in all directions. When detected far from the sample at a position r d , the scattered light has an additional phase factor e -iq•r i , which is equal to 1 when detected in the direction of the incident light, where q = q 0q d , with q 0 and q d being the incident and scattered wave vectors, respectively. q = 4π λ sin( θ 2 ), with θ = acos( u 0 •(r d -r i ) |u 0 ||(r d -r i )| ). Thus, when a particle at a given time t and at a given position r i jumps to another position at another time, its phase changes; then the phase coherence diminishes at each jump. This results in an additional broadening of the line described by Eqs. (6) or (13), which contains the spatial jump structure.

The Kehr-Richter-Honig model described in Sec. II A 1 has been developed assuming a uniform lattice, but it is valid (except for the inclusion of the spin-orbit coupling) for jumps with different lengths while the system is probed in transmission. We revisit the steps of the random-walk description leading to Eq. ( 6) taking into account the spatial distribution of jumps and the optical spatial phase factor appearing in a DLS measurement. Let us consider that the jumps occur between a given site and N others. A spin at a position r i can jump to N other positions r i + d ν , where d ν (ν = 1, 2, . . . , N) is a set of N jump vectors connecting the site with the others. The jump rate from a site r i to a site r i + d ν is ν . We assume that the spatial structure of the jumps is the same whatever the site position, i.e., a spin initially at r i after one jump will occupy N different sites at r i + d ν depending on the jump rates and after a second jump it will occupy ∼N 2 /2 sites at positions r i + d ν + d ν (ν, ν = 1, 2, . . . , N) (the factor 1/2 accounts for the fact that the same site can be reached by two different paths). At a site a spin spends a time τ c with 1/τ c = N ν=1 ν , and the time required to jump from site to site is neglected. For zero jump, Eq. ( 2) remains unchanged, but after one jump we have to take into account the phase variation δφ = q • d ν due to the spatial displacement d ν of the spin; then Eq. ( 3) becomes

R 1 (t ) = t 0 dt R 0 (t -t ) N ν=1 ν e -iq•d ν R 0 (t ) (15) = τ c N ν=1 ν e -iq•d ν t 0 dt R 0 (t -t ) 1 τ c R 0 (t ) (16)
and the Laplace transform of the series in Eq. ( 1) is written as

P(s) = R0 (s) ∞ l=0 τ c N ν=1 ν e -iq•d ν 1 τ c R0 (s) l (17) = R0 (s) 1 -[1 -τ c (q)] 1 τ c R0 (s) , ( 18
)
with R0 (s) given by Eq. ( 8) and (q) being the Fourier transform of the transition rate matrix

i,i , (q) = i,i i,i e -iq•(r i -r i ) = N ν=1 ν (1 -e -iq•d ν ). i,i
satisfy the master equation of the conditional probability P(r i t | 00) of finding a particle at site r i at time t when it was at site 0 at time 0,

d dt P(r i t | 00) = i i,i P(r i t | 00), (19)
with

i,i = ⎧ ⎪ ⎨ ⎪ ⎩ ν ν if r i = r i -ν if r i -r i = d ν 0 otherwise.
The master equation has a convolution form which takes a more practical form in Fourier space,

d dt P(q, t ) = -(q)P(q, t ), (20)
where P(q, t ) = i e iq.r i P(r i t | 00) = e -(q)t . The form of (q) has been calculated for various situations such as atom or molecule diffusion in solids or in liquids. One can name the well-known models of Chudley-Elliott (CE) [41], Singwi-Sjölander (SS) [42], and Hall-Ross (HR) [53], largely used in neutron scattering theories and experiments in liquids. These models are called jump-diffusion models where the particle is in its equilibrium position for a mean time τ c and then diffuses by continuous motion for a mean time τ d . When τ c τ d , the diffusion takes place as successive jumps. The different models have been developed for different jump length distributions. The CE model assumes jumps of fixed distance but in random directions. In the SS model the particle has a probability e -t/τ d to remain in diffusive motion; this leads to an exponential jump length distribution as given in Fig. 1 [54]. The HR model deals with a Gaussian jump length distribution. In Fig. 1 the forms of (q) and of the spatial jump distribution ρ(r) are given for the different models. One can note that, whatever the model, in the small-q limit we recover the q 2 law for continuous diffusion for (q) and in the large-q limit all models yield (q) = 1 τ c .

FIG. 1. The Fourier transform of the transition rate matrix (q) for a hypercubic arrangement of sites (chain, square lattice, simple cubic lattice, etc.) and for different models of jump diffusion with a distribution of jump lengths ρ(r). l2 is the mean-square displacement, which is related to the characteristic jump length l 0 usually used in the models.

B. Interplay of hyperfine interaction and spatially resolved spin jump distribution

Motional narrowing suppression: Probing no jump

The interplay of the hyperfine interaction and carrier hopping in spin dynamics was studied in works on muon diffusion [50] (approach used in this paper) and in molecular systems [51,55,56]. More recently, an analytical expression of the spin noise spectrum has been derived for an ensemble of localized electrons, taking into account nuclear effects and electron hopping between the localized sites [52]. We have not seen reports of the interplay of the hyperfine interaction and spin jump in a configuration of DLS measurement. We will discuss some behaviors suggested by the wavevector-dependent spin noise spectrum S 0 (ω, q) deduced from Eqs. ( 18) and ( 12),

S 0 (ω, q) = 1 π Re R0 (s) 1 -[1 -τ c (q)] 1 τ c R0 (s) , ( 21
)
with s = i(ω -ω 0 ). First, in a transmission measurement with q = 0, (0) = 0, and Eq. ( 21) reduces to the real part of Eq. ( 6), which reproduces the Gaussian decay of the polarization for large σ τ c due to hyperfine interaction and motional narrowing inducing exponential decay for small σ τ c [50]. For signals detected in a direction different from the incident light, for intermediate σ τ c , and for not-too-large wave vectors, the spin noise line broadens in a peculiar way. As known in motional narrowing, when σ τ c decreases, the line shape changes from Gaussian to Lorentzian in such a way that the Gaussian part decreases and the Lorentzian part increases. At intermediate σ τ c , when the line shape is partly Gaussian and partly Lorentzian, an increase in q broadens the spin noise line with an increase in the Gaussian part rather than the Lorentzian part. This kind of behavior is usually encountered with an increase in σ τ c because the spins become more and more localized. The fact that we observe a similar behavior for the line shape when q increases suggests that we look at spins which experience some but not all jumps of the distribution. The detection selects some diffusion paths; thus the spins seem more localized. Moreover, when q has a value q c such that (q c ) = 1/τ c , Eq. ( 21) reduces to

S 0 (ω, q c ) = 1 π Re{ R0 (s)}. ( 22
)
Then for large wave vectors, in the case of continuous jump length distributions or for periodic values of q in the case of jump on uniform lattices, the spin motion is suppressed by the detection configuration, i.e., the spatial optical phase averages out the spin jump distribution, and the spin dynamics before the first jump is measured. Thus no motional narrowing occurs when σ τ c decreases, and the Gaussian decay of the polarization is accelerated by the probability that the spin leaves its equilibrium position.

Spin-orbit coupling

The inclusion of angular rotation due to spin-orbit coupling is quite simple if the diffusive motion is a direct hopping between nearest neighbors on uniform lattices because the hopping rotation is independent of the time (it appears only when the hop occurs) and it is independent of space since the rotation is the same at each hop. Thus the inclusion of a polarization decay due to spin-orbit coupling R so in Eq. ( 15) is obvious, and we get

R 1 (t ) = R so t 0 dt R 0 (t -t ) N ν=1 ν e -iq•d ν R 0 (t ) = R so τ c N ν=1 ν e -iq•d ν t 0 dt R 0 (t -t ) 1 τ c R 0 (t ); (23)
the Laplace transform of the series in Eq. ( 1) is written as

P(s) = R0 (s) 1 -R so (1 -τ c (ul ) (q)) 1 τ c R0 (s) , ( 24
)
with (ul ) (q) being the diffusion broadening for uniform lattices, (ul ) (q) = 1 τ c {1 -1 3 [cos(q x a) + cos(q y a) + cos(q z a)]}, for a simple cubic lattice of spacing a. It is interesting to note that if (ul ) (q) = 1/τ c , then P(s) reduces to τ c R0 (s), which corresponds to the situation explained in the previous section. The spin-orbit rotation is suppressed because the correlations are lost (by the spatial averaging by the optical spatial phase) before the spin reaches its second equilibrium position (its equilibrium position after one jump). One can expect similar behavior at large q in the case of continuous jump length distributions, i.e., spin-orbit-induced broadening vanishes when (q) = 1/τ c . Therefore the hopping-induced rotation will be space dependent leading to a wave-vector dependency, which would be convoluted with 1 -τ c (q) as suggested by the term N ν=1 R so,ν ν e -iq•d ν in the equation below taking into account the jump length dependency of R so in the expression of R 1 (t ),

R 1 (t ) = t 0 dt R 0 (t -t ) N ν=1 R so,ν ν e -iq•d ν R 0 (t ) = τ c N ν=1 R so,ν ν e -iq•d ν t 0 dt R 0 (t -t ) 1 τ c R 0 (t ), (25)
with R so,ν being a function of d ν which depends on the mechanism of the diffusive motion leading to the jumps from one site to another. The Laplace transform of the series in Eq. ( 1) becomes

P(s) = R0 (s) 1 -τ c N ν=1 R so,ν ν e -iq•d ν 1 τ c R0 (s) . ( 26
)
To evaluate the term

N ν=1 R so,ν ν e -iq•d ν , we first must know the form of N ν=1 ν e -iq•d ν . We know that τ c N ν=1 ν e -iq•d ν = 1 -τ c (q)
; in the CE model the sum ν is replaced by an integral over space and orientational averaging. It is possible to determine (q) whatever the jump length distribution from CE (q), the diffusion broadening obtained with the CE model, through the following relation:

(q) = ∞ 0 CE (q)ρ(r)dr = ∞ 0 1 τ c (1 -sin(qr) qr )ρ(
d ν →r (27)
When q → ∞, the left and right terms in Eq. ( 27) tend to 0, and P(s) reduces to τ c R0 (s) as expected. For small q, while sin(qr) ∼ qr, Eq. ( 27) gives Rso , the mean value of R so . For any q, ∞ 0 sin(qr) qr R so (r)ρ(r)dr = Rso (q) * [1 -τ c (q)], with * being the convolution product, which can be replaced by Rso [1 -τ c (q)] when R so (r) is a smooth function compared with the jump length distribution.

Let us discuss the case where each jump between localized sites consists of numerous fast (compared with τ c ) and short-distance (compared with the characteristic jump length l 0 ) elementary jumps in random directions. One can assume, for example, fast hops between nearest neighbors separated by a distance r c . At each hop the polarization decay due to spin-orbit coupling is given by R so = 1 -1 3 γ 2 , if the angular rotation γ is small. Thus, when a spin at a site r i jumps to another site r i + d ν , it will experience N ν = (d ν /r c ) 2 uncorrelated rotations by an angle γ . Thus one can basically write 2 , and thus

R so,ν = 1 -1 3 γ 2 ( d ν r c ) 2 or R so (r) = 1 -1 3 γ 2 ( r r c )
Rso = 1 - 1 3 γ 2 l2 r 2 c , ( 28
)
where l2 = ∞ 0 r 2 ρ(r)dr is the mean-square displacement. If one introduces the spin-orbit length L so , which characterizes the jump distance at which the spin orientation is lost, we may FIG. 2. Half-width at half maximum (HWHM) of the real part of Eq. ( 27) in reduced units for Gaussian (thin lines) and exponential (thick lines) jump length distribution, as a function of q with = √ l2 . Solid lines, Eq. ( 27) with R so (r) = 1 -1 3 ( r Lso ) 2 . Dashed lines, Eq. ( 27) assumed to be equal to Rso [1 -τ c (q)].

write γ = r c /L so [57], so that, in fine,

Rso = 1 - 1 3 l2 L 2 so . ( 29
)
We note that this expression is identical to the one expected for a single jump between distant donors separated by an average mean-square distance l2 .

Figure 2 shows the broadening of the line deduced from the real part of Eq. ( 26) calculated for Gaussian and exponential jump length distributions, for different values of the spin-orbit length and using Eq. ( 27) or its averaged form. For L so 5 l2 the relation Rso [1 -τ c (q)] = Rso (q) * [1 -τ c (q)] is almost true, while for L so 2 l2 it is difficult in an experiment to differentiate between the convolution or the averaged forms.

Keeping the averaged form and using Eqs. ( 12) and ( 26), the q-resolved spin noise spectra S 0 (ω, q) can be expressed as

S 0 (ω, q) = 1 π Re R(ω) 1 -R so F (q) 1 τ c R(ω) , R(ω) = 1 σ π 2 W ω -ω 0 + i(ν s + 1/τ c ) √ 2σ ,
F (q) = 1 N N ν=1 exp(-iq • d ν ). ( 30
)
It contains the one-site response function R(ω) = R0 (s + ν s ), and we have introduced the structure factor F (q) = 1τ c (q). FIG. 3. Sketch of an electron spin trajectory along a chain of trapped states (green circles) starting from an initial state "0." Only the first five states are represented. We consider a spin initially in a bright fringe. The fringes represent the spin visibility, which depends on the relative phase between the incident (wave vector k i ) and the scattered light (wave vector k s ) (see Ref. [26]). The spin is represented as a green arrow and precesses in a local nuclear field (dashed arrow). The black curve represents the spin trajectory along fast diffusion channels. From left to right, three spin noise detection schemes used in this paper are represented. (a) Forward-scattering (FWD) geometry (classical SNS): Spin visibility is maximum and constant along the trajectory. All terms in Eq. ( 32) equally contribute to the signal. (b) ST-SNS: The spin visibility is slowly modulated. The terms contribute less and less as the number of jumps increases. (c) Backscattering (Back) geometry: When the modulation period of the spin visibility is smaller than the distance between traps, only the initial state "0" contributes to the signal (for random trap positions).

Time of flight

Next, we introduce the finite resolution in q via the convolution [26]

S(ω, q) = ∞ -∞ d 2 q S 0 (ω, q ) Î2 0 (q -),
where Î0 (q) = e -q 2 /2Q 2 /( √ π Q) is the wave-vector distribution of the probe beam, with Q = √ 8 ln(2)/w, where w is the laser spot size on the sample. For a distribution that is not too large, one finds an analytical expression (see Appendix A)

S(ω, q) = 1 π Re R(ω) 1 -R so G(q) 1 τ c R(ω) , G(q) = F (q) + 1 4 Q 2 F (2) (q). ( 32
)
F (2) (q) is the second derivative (Laplacian in q space) of the structure factor F (q). The corresponding term in Eq. ( 32) can be interpreted as the polarization decay due to time of flight (TOF) of the spin fluctuation across the probed area [58]. In effect the associated broadening is inversely proportional to the time for the spin fluctuation to leave the sampling area.

Equations ( 30) and (32) are only valid at large external field satisfying ω 0 σ . To obtain the noise spectra in zero field, we use for R(ω) the one-site response calculated by Glazov {the term τ ω A(τ ω ) in Eqs. ( 6)-( 9) of Ref. [52]}.

Interpretation

Equation (30) has a simple interpretation. R(ω) is the response of spins during their residence time on one site. At each step of the random walk, the spin coherence decays by a constant factor, which adds, to the decoherence during the time spent in the trap described by 1 τ c R(ω), the spin rotation induced by SOC during the jump described by 1 -R so , and the structure factor F (q). Hence Eq. ( 30) simply results from the summation of a geometric series, where the term of order p corresponds to p jumps.

Figure 3 graphically represents how the different terms of the geometric series contribute to the signal depending on the different detection schemes we use. Only one spin trajectory is represented, and we disregard TOF effects (Q = 0) and relaxation by SOC (R so = 1). In backscattering, where q is much larger than the inverse of the jump distance, the structure factor F (q) averages to zero, and S(ω, q) ∝ Re[R(ω)]. The spectrum corresponds to the one-site response, in agreement with Fig. 3(c). Under an applied magnetic field, it follows an exact Voigt profile with σ as Gaussian parameter and γ = 1/τ c + ν s as Lorentzian parameter.

In classical SNS [q = 0, Fig. 3(a)], F (0) = 1, so that all terms contribute. In other words, one follows the spin all along the trajectory, and the signal results from averaging over an ensemble of such trajectories. S(ω, q) remains very close to a Voigt profile, but the Gaussian parameter is smaller than σ : The linewidth decreases by motional narrowing of the local nuclear fields. In ST-SNS [Fig. 3(b)] a finite number of terms contribute; in other words, the trajectory fades away. This adds a broadening mechanism linked to motion, quadratic at small q (q 1), from which the spin diffusion coefficient D s can be deduced. Finally, in the presence of SOC a competition sets in between broadening due to relaxation by SOC and motional narrowing.

III. SPIN NOISE EXPERIMENTS

We studied an Al-doped CdTe layer grown by molecular beam epitaxy on a Cd 0.96 Zn 0.04 Te(001) substrate, mechanically polished on its back face. The donor density, below the Mott transition, was estimated from secondary ion mass spectroscopy to be ∼2 × 10 16 cm -3 , corresponding to an average interdonor distance of ∼50 nm. The layer thickness is 10 μm. In a previous study we have observed a satellite line in the zero-field spin noise spectra of this sample (M3460), assigned to the electron spin precession in the local nuclear field ∼1.2 mT [27]. A similar effect was also observed in high-purity GaAs [59].

Details of the experimental setup can be found in Refs. [26,27]. Here, the laser was detuned 9.3 meV below the donor-bound exciton at E = 1.5926 eV in a region of low absorption. Spin noise spectra with and without applied magnetic field were obtained simultaneously by repeated FIG. 4. Series of spin noise spectra for three different temperatures in forward scattering and backscattering (self-homodyne detection), in zero field (left) and in a transverse magnetic field (right). The curves are calculated with Eq. (32). Red (blue) curves are calculated at q = 0 (q = ∞). At 4 K and B = 0 we highlight the difference between forward scattering and backscattering by plotting the curves calculated for both configurations (same color code).

interleaved acquisition in a single experimental run. The applied field 16 mT was chosen much larger than the local nuclear field.

Figure 4 shows typical spin noise spectra, in forward scattering and backscattering (intermediate spectra using ST-SNS are not shown). A good fit was obtained for the spectra with field using Voigt profiles, for all values of q, and the halflinewidth at half maximum ω 1 2 (q) is shown in Fig. 5(a) [60]. The strong dependence on q is evidence of the role of spin motion. Figure 5(a) shows that ω 1 2 (q) increases quadratically at small q (see also Ref. [26]): This is the usual continuous diffusion regime described by a diffusion coefficient D s with a diffusion damping rate γ diff = D s q 2 . However, ω 1 2 (q) exhibits a plateau at large q. This behavior suggests that a characteristic length in the micrometer range characterizes the motion. Below we will interpret as the mean hopping distance and show that indeed the regime q 1 is reached for the backscattering spectra. Under such conditions [Fig. 3(c)], the signal is related to the one-site response: It can be characterized by an inhomogeneous broadening {for instance, the Gaussian distribution of precession frequencies in the trapped states p(ω) = (2πσ 2 ) -1/2 exp[-(ω -ω 0 ) 2 /2σ 2 ] due to hyperfine interaction (HFI)}, an intrasite spin relaxation rate ν s (ISR), and the spin residence time τ c (hopping probability 1/τ c ). The near-Gaussian profile at 4 K suggests that the linewidth is inhomogeneous and a good estimate of σ . The evolution towards a broader Lorentzian profile in backscattering (q = ∞) as the temperature is increased shows that ν s + 1/τ c increases. The evolution towards a narrower Lorentzian profile in forward scattering (q = 0) is interpreted as the result of motional narrowing, with an expected linewidth ∼ σ 2 τ c . Finally, the broad Lorentzian profile at higher temperature is interpreted as due to the increased role of time of flight (TOF; ∼ A quantitative check using Eq. ( 32) is shown in Figs. 4 and5. We obtain a simultaneous agreement between all experimental and theoretical spectra (including those at finite q) for a Gaussian distribution with = 2.7 μm, and for σ = 1.12 × 10 8 s -1 and L so = 5 μm (see Appendix C for details of the fitting procedure).

The value w = 16.8 μm is estimated from the measurement of the spot size. The main parameter which changes with temperature is τ c . In order to reproduce the line shape of the satellite at q = 0, we have to assume an additional broadening mechanism at 4 K. Although the origin of this broadening is still unclear, we take it into account by assuming 32) (parameters in the text). In (a), ω 1 2 (q) is calculated for a Gaussian (solid lines) and an exponential distribution of jump lengths (dashed lines). In (b) the thin solid black line is calculated for q = ∞; the other black lines are calculated for q = 0 with only HFI included (dotted line), HFI + SOC (dashed line), HFI + SOC + TOF (thick solid line). The asymptotic behaviors are also indicated for σ τ c → 0 and σ τ c → ∞. The rightmost point in (b) was obtained at P = 0.25 mW. FIG. 6. HWHM ω 1/2 of the spin noise line as a function of the wave-vector distribution Q LO of the LO beam, experiments (symbols: triangles and circles, self-homodyne and homodyne detection, respectively) and models (lines). The black solid line shows the calculation with Eq. ( 31) including the optical transfer function of the collection long-working-distance objective (Mitutoyo ×10), the dash-dotted black line shows the calculation with Eq. ( 31), and the dotted black line shows the calculation with Eq. ( 32). The calculations are done for a Gaussian jump length distribution, with Q 0 = 115 000 m -1 , τ c = 11 ns, σ = 0.112 ns -1 , l = 2.7 μm, Rso = 0.903, and q = 0 (except for the dashed black where q = ∞). Gray solid and dash-dotted lines show the calculation with Eq. ( 31) (q = 0) with and without the optical transfer function of the imaging system, respectively, and for Brownian motion with a diffusion coefficient

D s = l 2 /6τ c . a small ISR [ν s = (2π ) × 4.3 MHz].
Remarkably the satellite at zero field is better resolved in backscattering (see spectra at 4 K in Fig. 4). This difference is also well reproduced by our calculation. Indeed, in this situation only the initial state "0" not affected by spin motion is observed [Fig. 3(c)]. In other words, the one-site response is being probed. Finally, once τ c and are determined, one deduces D s = 2 /6τ c , which varies in the range 0.44-14 cm 2 /s from 4 to 18 K.

The time-of-flight broadening appears in the expression of the spin noise spectrum in the same way as the spin-orbit coupling; it is then impossible to discriminate between these two mechanisms without knowing D s and the spatial extent of the beam. A way to differentiate between both mechanisms would be to measure the spin noise linewidth as a function of the spot size on the sample. However, focusing more or less may, for example, induce a temperature change and then a change in the correlation time. A better way is to use the homodyne detection technique in order to keep unchanged (and as large as possible) the spot size on the sample and to select, with a local oscillator (LO), the probed area in the detection plane. In this case, in the expression of Î2 0 (q), Q 2 must be changed by 1 2 (Q 2 0 + Q 2 LO ), where Q 0 and Q LO are linked to the widths of the wave-vector distributions of the beam probing the sample and the local oscillator, respectively. Calculations and measurements are depicted in Fig. 6. In the experiment a large probe beam (Q 0 = 115 000 m -1 ) illuminates the sample, while a LO beam with a varying Q LO probes the correlations in the sample's image. Triangles are for self-homodyne measurements (without the LO beam) in transmission (forward) and reflection (back), while open circles correspond to homodyne measurement (with LO beam) in transmission. Calculations are made taking into account or not the transfer function of the imaging system. When optical transfer function h(q) is considered, S(ω, q) must be replaced by S(ω, q)h(q); this case is never considered in the approximate solution given by Eq. (32).

We found that the broadening due to the TOF agrees well with the calculated one. This is important to deduce the value of ω 1 2 (q) corrected for the TOF [Fig. 5(b)].

IV. CONCLUSION

To conclude, our experiments have uncovered a spin-jumpdiffusion regime in an n-doped CdTe layer, which manifests as a plateau at large wave vectors for the broadening of the q-resolved spin noise spectra. We have proposed phenomenological equations, which simultaneously fit all spin noise spectra in zero and nonzero magnetic field, for all q vectors and temperatures, with a fixed set of parameters σ , L so , and . Besides, τ c decreases with temperature and is mainly determined from the width of the noise spectra in backscattering. Remarkably, this determination of τ c is consistent with the observation or not of the satellite line. It turns out that the parameters which govern the spin relaxation by HFI and SOC in CdTe, σ and L so , respectively, are quite similar to the values in GaAs (σ = 1.29 × 10 8 s -1 [59], L so 5 μm [57]). Systematic measurements for different doping levels will be necessary for a quantitative comparison with GaAs. However, we have strong indications that the faster relaxation in CdTe is not due to material parameters, but rather to the specific long-distance spin jump diffusion that we have uncovered. The inhomogeneities of the system may plausibly play an important role in this phenomenon [61][62][63].

Our experimental value of L so is close to the theoretical value L so 3.37E b a 4 b /γ c = 4.1 μm [57], where E b = 13 meV is the donor binding energy, a b = 5.3 nm is its Bohr radius, and γ c = 8.5 eV Å 3 is the Dresselhaus coefficient taken from theoretical estimations [64]. σ can be estimated from the Bohr radius and the hyperfine constants of Cd and Te. Unfortunately, reliable values for the hyperfine constants are missing (see discussion in Ref. [27]). The longest spin relaxation time, corrected for a small TOF effect, that we measure in the regime of motional narrowing is T * 2 = 15 ns (see Appendix B for the precise definition of T * 2 and how it is determined).

The long-range spin hopping uncovered here in CdTe has important consequences for the electron spatiotemporal spin dynamics and hence for qubit operations based on donors: It should be sought in other materials which are under focus for quantum technologies [3,[START_REF] Greilich | [END_REF]47], such as GaAs, ZnSe, ZnO, and, if possible, Si and Ge. Our approach could also answer fundamental questions which underpin such operations, among which is the predicted random spin rotation during a hop between two donors, which arises due to the uncertainty in the quantum trajectory [61].

Let us consider the spin noise line whose spectrum at any wave vector, in the presence of hyperfine interaction, spinorbit coupling, and time of flight, is given by Eq. (A5). Under an external magnetic field applied transverse to the incident light, R(ω) identifies with R0 (s + ν s ), R0 (s) being given by Eq. (8). For the two limiting cases σ τ c 1 and σ τ c 1, one can easily deduce the line shape, which will be Gaussian in the former case and Lorentzian in the latter case. The damping parameter , defined as the inverse of the time during which the temporal signal has decayed by 1/e, is a known parameter for both cases. equals the damping rate γ for a Lorentzianshaped line and the square root of half the variance for a Gaussian-shaped line. Thus having acquaintance with the line shape, knowing the linewidth, one can deduce the damping parameter. This is true only for short and long τ c . Indeed, there is no analytical expression for the damping parameter corresponding to Eq. (A5) for any σ τ c , neither for the width of the line, nor for the relation between these two parameters. It is then difficult to estimate correctly the spin relaxation time from the noise power spectrum. A good approximation of the spin relaxation time is half the amplitude of the spin noise line when its surface is normalized to unity. This relaxation time, called T * 2 and introduced by Bloembergen et al. [65], is defined by

T * 2 = 1 2 g(ν) max dνg(ν) = π g(ω) max dωg(ω) , ( B1
)
where g(ν) [g(ω)] is the spectrum in the frequency (pulsation) ν (ω) space and g(ν) max [g(ω) max ] is its maximum amplitude. (T * 2 ) -1 equals the damping rate for a Lorentzian-shaped line and overestimates by an amount of less than 13% the damping parameter for a Gaussian-shaped line, so that (T * 2 ) -1 ≈ and one can call it the damping parameter. T * 2 is thus a good estimate of the spin relaxation time; moreover, it can be estimated experimentally without any assumption as to the line shape. A joint determination of (T * 2 ) -1 and the half-width at half maximum ω 1/2 gives information about the line shape. If (T * 2 ) -1 / ω 1/2 = 1, the line is Lorentzian, and if (T * 2 ) -1 / ω 1/2 = 0.678, the line is Gaussian. Definition (B1) applied to Eq. (A5) leads to the following expression for the damping parameter: 2 ) -1 deduced from Eqs. (A5) and (B2) for q = 0 (gray solid line) and q = ∞ (black solid line). Right: Half-width at half maximum ω 1/2 from experiments (points with error bars) and from Eq. (A5) (solid lines) [black and gray stand for backscattering SNS (q = ∞) and classical SNS (q = 0), respectively].

(T * 2,q ) -1 = R -1 (ω 0 ) - η(q) τ c , ( B2
with η(q) = Rso G(q). Let us discuss this quite interesting formulation for the case q = 0 and η(0) = 1. In this case, Eq. (B2) expresses the line broadening due to hyperfine interaction only, leading to motional narrowing with decreasing τ c . R -1 0 (ω 0 ) is the broadening of a Voigt profile which mixes the inhomogeneous broadening due to hyperfine interaction and the homogeneous broadening due to the probability 1/τ c that the spin leaves a localized site. This homogeneous broadening does not contribute in a classical measurement (q = 0); it is then subtracted by the second term in Eq. (B2), so that only the competition between the inhomogeneous and the homogeneous parts remains. At short τ c this is quite visible because the line broadening is proportional to σ σ γ , with γ = 1/τ c being the damping rate of the homogeneous part. When η < 1, the other sources of broadening add simply as homogeneous terms.

When looking at the damping parameter in the two limiting cases q = 0 and q = ∞, they differ only by the quantity η(0)/τ c . Then a measure of the difference (T

* 2,∞ ) -1 -(T * 2,0 ) -1 is a direct determination of η(0)/τ c . Plotting (T * 2,q ) -1 as a function of (T * 2,∞ ) -1 -(T * 2
,0 ) -1 will place the data automatically at the right abscissa position getting rid of one fitting parameter, the correlation time τ c .

In Fig. 7 we plot (T * 2 ) -1 and ω 1/2 extracted manually (without using any fitting function) from the noise power spectrum in measurements made on the n-CdTe sample described in the main text, under a transverse magnetic field of 16 mT, for different temperatures and for two different detection schemes [backscattering SNS, (T * 2,b ) -1 from noise of the reflected beam; forward scattering SNS, (T * 2, f ) -1 from noise of the transmitted beam]. On the left-hand side, the experimentally determined (T * 2,b/ f ) -1 are plotted as a function of their differences, while (T * 2,q ) -1 calculated with Eq. (B2) for q = 0 and q = ∞ is plotted as a function of η/τ c , with η = η(q = 0). A good agreement between theory and experiment is obtained for σ = 1.12 × 10 8 s -1 , η = 0.87, and ν s = 0. On the right-hand side, the experimentally determined half-width at half maximum is plotted as a function of η/[(T * 2,b ) -1 -(T * 2, f ) -1 ], while the numerical solution of the equation S( ω 1/2 , q) = S(ω 0 , q)/2 (with σ , η, and ν s determined previously) for q = 0 and q = ∞ is plotted as a function of τ c . From the ratio (T * 2 ) -1 / ω 1/2 , one can deduce that the line is Lorentzian at short τ c (high temperature T > 16 K) and almost Gaussian at long τ c (low temperature T = 4 K).

Finally, inspired by Eq. (B2), one can attempt to give an analytical expression for the damping parameter of Eq. ( 5) using the damping parameter of R 0 (t ) = exp(-1 2 σ 2 t 2 ) exp(-t/τ c ). Let us call σ,τ c the damping parameter of R 0 (t ), such that R 0 ( -1 σ,τ c ) = 1/e. Following the structure of Eq. (B2), we propose

= σ,2τ c - 1 2τ c . ( B3
)
In Fig. 8 we compare Eq. (B2) [for η(q = 0) = 1] and Eq. (B3) with extracted from Fig. 3 of Ref. [50], deduced from P(t ) obtained by inversion of the Laplace transform of Eq. ( 6). We find a maximum deviation of 1.7% for σ τ c ≈ √ π between Eq. (B3) and extracted from Fig. 3 of Ref. [50].

When considering other sources of homogeneous broadening, then η(q = 0) < 1, so that Eq. (B3) is written as

= σ,2τ c - 2η(q) -1 2τ c . ( B4
)

APPENDIX C: DETERMINATION OF THE SPIN RELAXATION PARAMETERS

In order to determine the spin relaxation parameters σ , L so , τ c , and , we make use of the results presented in Appendixes A and B, together with the q dependence of ω 1 2 [Fig. 5(a)]. As we have seen in Appendix B, using the spin relaxation rates is very convenient for this task because it allows for the direct determination of the parameter η(q)/τ c . However, the experimental determination of (T * 2,q ) -1 requires a high enough signal-to-noise ratio in order to have a good estimate of the wings of the line, and of the background level. This is a bit more difficult to estimate in the ST-SNS configuration, where the optical transfer function of the collection objective filters out the signal as q increases. Therefore we could determine precisely only (T * 2, f ) -1 and (T * 2,b ) -1 (see Fig. 7). This is why in Fig. 5 we plot only ω 1 2 . However, from Fig. 5

T -1 2, f = T -1 2,b - η τ c . ( C2
)
The difference between back and forward spin relaxation rates already fixes the value of the parameter η/τ c , where η = R so [1 -(Q ) 2 /6]. At long τ c (low temperatures) these two relaxation rates converge towards the common value √ 2 π σ . This approximately determines σ . At short τ c (high temperatures), (T * 2,b ) -1 tends towards 1/τ c . Thus one can roughly deduce both η and τ c . Then the fine-tuning of the parameters is made by seeking for the best agreement for all temperatures between experiment and theory (see Fig. 7). We finally obtain σ = 1.12 × 10 8 s -1 and η = 0.87. τ c is thus automatically determined at all temperatures.

We are left with the determination of L so and . = 2.7 μm is deduced from the variation of ω 1 2 with q [Fig. 5(a)]. Then from the measured width of the spot size w = 16.8 μm, we get Q = 140 200 m -1 ; hence R so = 0.903, and finally L so = 5 μm.

  r)dr, with ρ(r) being the jump length distribution and for ∞ 0 ρ(r)dr = 1. Then we get τ c N ν=1 R so,ν ν e -iq•d ν ≡ ∞ 0 sin(qr) qr R so (r)ρ(r)dr.

2 w 2

 22 τ c ) and spin-orbit coupling (SOC; ∼ 2 L 2 so τ c ): Both rates increase as τ c decreases. Those trends are schematized in Fig. 5(b).

FIG. 5 . HWHM ω 1 2

 51 FIG. 5. HWHM ω 1 2 vs q (a) and vs τ c (b). Experimental values (symbols) are measured at 16 mT (probe laser power P = 2 mW). Data are obtained by self-homodyne detection either in forward scattering (open circles) or backscattering (dark circles), and by ST-SNS (gray circles). All curves are calculated using Eq. (32) (parameters in the text). In (a), ω 1

)FIG. 7 .

 7 FIG. 7. Left: (T * 2,b/ f ) -1 vs (T * 2,b ) -1 -(T * 2, f ) -1 deduced with Eq. (B1) from the measured spin noise power spectrum for two different detection schemes: backscattering SNS [(T * 2,b ) -1 ] and classical SNS [(T * 2, f ) -1 ] (black and gray points with error bars, respectively). (T *2 ) -1 deduced from Eqs. (A5) and (B2) for q = 0 (gray solid line) and q = ∞ (black solid line). Right: Half-width at half maximum ω 1/2 from experiments (points with error bars) and from Eq. (A5) (solid lines) [black and gray stand for backscattering SNS (q = ∞) and classical SNS (q = 0), respectively].

FIG. 8 .

 8 FIG. 8. Comparison of the damping constant /σ for the Markovian random-walk model (symbols; extracted from Fig. 3 of Ref. [50]), Eq. (B2) (dashed black line), and Eq. (B3) (solid black line).

  one can see that T * 2,b ≡ T * 2,∞ , while T * 2, f ≡ T * 2,0 is always true. Hence with Eq. (B2) we can writeT -1 2,b = R -1 (ω 0 ),(C1)

ACKNOWLEDGMENTS

S.C. and D.S. wish to thank C. Lhenoret for technical assistance.

APPENDIX A: TIME OF FLIGHT

The spin may diffuse outside of the probed area leading to an additional broadening [58] that we will call time-of-flight broadening because it is inversely proportional to the time for the spin to leave the sampling area. It can be taken into account by a convolution of the q-dependent spectrum with the wave-vector distribution of the probe beam intensity [26],

with Î0 (q) being the space Fourier transform of I 0 (r), the intensity profile of the probe beam. It is possible to give a simple analytical expression of integral in Eq. (A1) by using the Taylor series expansion of F (q ) at the measured q,

∂q y

] q =q e y . Then, if Î0 (q) is symmetrical, the convolution integral (A1) reduces to an integral,

If the wave-vector distribution of the beam is not too large or if F (q) varies a little in the width of this distribution such that [∇ q F (q)] • q + 1 2 ∇ q {[∇ q F (q)] • q } • q F (q), the integration can be put directly in the denominator,

Linear in q and mixed second-order partial derivative terms are equal to 0 after integration. For a Gaussian beam of the form Î2 0 (q) = e -q 2 /Q 2 /(π Q 2 ) the integral reduces to 1 4 Q 2 F (2) (q), and S(ω, q) becomes

with

is then the broadening due to the probability of the spin leaving the sampling area. In a transmission experiment with q = 0, F (0) = 1, and γ TOF = Q 2 l2 6τ c = Q 2 D s , with D s being the spin diffusion coefficient, whatever the jump length distribution.


Juraj Szavits-Nossan

Luca Ciandrini

M Carmen Romano

Deciphering mRNA Sequence Determinants of Protein Production Rate

Keywords: numbers: 87.16.aj, 87.10.Mn, 05.60.-k

One of the greatest challenges in biophysical models of translation is to identify coding sequences features that affect the rate of translation and therefore the overall protein production in the cell. We propose an analytic method to solve a translation model based on the inhomogeneous totally asymmetric simple exclusion process, which allows us to unveil simple design principles of nucleotide sequences determining protein production rates. Our solution shows an excellent agreement when compared to numerical genome-wide simulations of S. cerevisiae transcript sequences and predicts that the first 10 codons, which is the ribosome footprint length on the mRNA, together with the value of the initiation rate, are the main determinants of protein production rate under physiological conditions. Finally, we interpret the obtained analytic results based on the evolutionary role of codons' choice for regulating translation rates and ribosome densities.

Translation is one of the major steps in protein biosynthesis. During this process, the nucleotide sequence of a messenger RNA (mRNA) is translated into a functional protein. Each nucleotide triplet, called codon, codes for a specific amino acid, the proteins' building block. There is experimental evidence that the rate at which a certain mRNA is translated depends on its specific codon sequence [1][2][3][4], especially in the case of eukaryotes. Identifying sequence features that determine protein production rate, also commonly referred to as translation rate or efficiency, is a fundamental open question in molecular biology [1,3].

Translation is performed by molecular motors called ribosomes, which move unidirectionally along the mRNA. The amino acids are delivered to the ribosome by molecules called transfer RNAs (tRNAs), which are specific to the codon and the amino acid they deliver (Fig. 1). The dwelling time of a ribosome on a specific codon depends primarily on the abundance of the corresponding tRNA [5,6]. Understanding how codon sequences determine protein production rates can potentially unlock many synthetic applications [7,8].

The standard biophysical model of translation is known as the totally asymmetric simple exclusion process (TASEP), which captures the concurrent motion of ribosomes on the mRNA [9,10]. In this model ribosomes progress along the mRNA codon by codon, provided that the codon a ribosome moves onto is not occupied by another ribosome. The protein production rate FIG. 1: Sketch of the mRNA translation process involving initiation (a), elongation consisting of tRNA delivery (b) and translocation (c), followed by termination (d). We emphasise that this is an oversimplified scheme of the process and that actual ribosomes cover = 10 codons. At each elongation step, the ribosome receives an amino acid from a tRNA that matches the codon occupied by the A-site of the ribosome (the "reading" site). After the amino acid is added to the growing polypeptide, the ribosome translocates one codon forward, and the process is repeated.

can then be identified as the ribosomal current of this arXiv:1708.07678v2 [physics.bio-ph] 8 Mar 2018 driven lattice gas. Due to the net current, TASEP is not in equilibrium and its steady state is only known for a few special cases [11,12]. Unfortunately, most biologically relevant variants of TASEP can be studied only numerically [13][14][15], and efficient methods of exploring a large number of parameters are lacking. The steady state of the TASEP with non-uniform hopping rates is a long outstanding problem in nonequilibrium statistical physics [15,16]. An analytic prediction of protein production rates is also needed in order to interpret recently developed ribosome profiling experiments that are capable of monitoring ribosome positions along the mRNA [17]. Such analytical prediction is missing in previous models of ribosome dynamics implicitly or explicitly based on the TASEP [1,2,19,20,22], and it is key to decipher sequence determinants of protein production rates.

In this Letter, we develop a versatile analytic method to solve TASEP-based models of translation. Our analytic approach, integrated with simulations and experimental data, allows us to efficiently identify the main features of mRNA codon sequence that determine the rate of protein production. The genome-wide comparison of our analytic predictions with numerical simulations shows an excellent agreement for the model organism S. cerevisiae (baker's yeast).

Stochastic model of mRNA translation. We focus on the model for translation introduced in [1,23,24]. The mRNA is represented by a one-dimensional lattice consisting of L discrete sites (codons), where site 1 designates the start codon. Ribosomes are represented by particles that occupy = 10 lattice sites, which is the ribosome footprint length measured in ribosome profiling experiments [17]. We identify the position 1 ≤ i ≤ L of a ribosome with the position of its A-site, which is located d = 5 lattice sites from the trailing end of the ribosome (Fig. 1a). A ribosome that waits for a tRNA at position i is labelled by 1 i and a ribosome that has already received the correct tRNA and is ready to move is labelled by 2 i . The set of labels of all translating ribosomes on the lattice is called a configuration C of the system. For example, C = 1 1 2 13 denotes a lattice with 2 ribosomes, one at site 1 waiting for a tRNA and another one at site 13 that has already received the correct tRNA.

Ribosomes initiate translation at rate α by binding to the mRNA so that their A-site is at the start codon, provided the sites 1, . . . , -d + 1 are empty (Fig. 1a). A ribosome at site i makes the transition 1 i → 2 i at rate k i dependent on tRNA abundances (Fig. 1b) and moves one site forward at rate γ (Fig. 1c). Due to steric interactions between particles, the particle at site i can move only if there is no ribosome at site i + . Termination occurs when the particle at site L receives the last amino acid, releases the final protein and detaches from the mRNA, which we integrate into a single step occurring at rate β (Fig. 1d).

Our main goal is to compute the rate of protein production as a function of the parameters of the model, which are α, β, γ and k i for each of the L codons used. We assume that translation takes place under steadystate conditions, so that the ribosomal current is constant along the mRNA and is equal to the rate of protein production, which in turn is equal to rate at which ribosomes load onto the mRNA and initiate translation. To this end we define the codon occupation number τ i to be equal to 1 if the i-th codon is occupied by an A-site, and 0 if otherwise. By definition, the exact steady-state ribosomal current J then reads

J = α C i=1 (1 -τ i (C)) P (C), (1)
where τ i (C) denotes the i-th codon occupation number for configuration C, P (C) denotes the steady-state probability that the lattice is in configuration C, and the summation goes over all configurations C. Other quantities of interest that we compute are the local and total particle densities ρ i = τ i and ρ = (1/L) L i=1 ρ i , respectively. Series expansion method for computing P (C). In order to find P (C) one has to solve the steady-state master equation M P = 0, where P is a column vector whose N elements are the steady-state probabilities P (C) of being in configuration C, and N denotes the total number of configurations. The transition rate matrix M C,C is given by W C →C for C = C and -e(C) for C = C , where W C →C is the transition rate from C to C and e(C) = C =C W C→C is the total exit rate from C. The exact solution of the master equation can be formally written as

P (C) = detM (p,p) N q=1 detM (q,q) , ( 2
)
where p is the position of configuration C in the column vector P and detM (p,p) is a determinant of the matrix obtained by removing p-th row and p-th column from M (see Supplemental Material for details). Unfortunately, calculating this determinant is feasible only for unrealistically small system sizes.

To circumvent this problem, we exploit the fact that detM (p,p) is a multivariate polynomial in the variables α, k 1 , . . . , k L , β and γ [25]. In the biological literature it is often assumed that the initiation rate α is a major limiting step of the translation process, mainly determined by the presence of secondary structures [26,27]. For this reason we use α as an expansion parameter and we assume that α γ, β, k 1 , . . . , k L . By collecting terms with the same power of α, we can rewrite P (C) as an univariate polynomial f (C) in the variable α with unknown coefficients f n (C)

P (C) = f (C) C f (C) , f (C) = K(C) n=0 f n (C)α n , (3)
where the coefficients f n (C) depend on the transition rates k 1 , . . . , k L , γ and β (in order to ease the notation, we leave out the explicit dependence on those parameters). Since we expect that f (C) can be well approximated by the first few terms, the value of K(C) in Eq. (3 is irrelevant in our study.

In order to find the unknown coefficients f n (C), we insert Eq. ( 3) into the master equation M P = 0, collect all the terms with the same n-th power of α and equate their sum to zero. For a given power n, the resulting equation is similar to the original master equation in which P (C) is replaced by f n (C) unless the coefficient multiplying P (C) is α, in which case P (C) is replaced by f n-1 (C). Starting with n = 0, we note the equations for the coefficients f 0 (C) have the same form as the original master equation, but with α = 0, leading to the trivial solution f 0 (C) = 0 if C = ∅. Since the terms f (C) are not normalized, we have the freedom to choose any value for f 0 (∅), which we set to 1.

For n ≥ 1, the equation for f n (C) reads

e 0 (C)f n (C) + f n-1 (C) C I C,C = C I C ,C f n-1 (C ) + C (1 -I C ,C )W C →C f n (C ) , (4)
where A key observation in our analysis is that f n (C) = 0 whenever the number of particles in C is larger than n. This follows from the result f 0 (C) = 1 (0) if C = ∅ (C = ∅) in conjunction with the hierarchical structure of Eq. ( 4), which connects configurations differing in the number of particles by no more than one (a proof for n = 1 is presented in the Supplemental Material). This allows us to write Eq. ( 4) taking into account configurations with only one particle and discarding all the others, which yields

I C,C = 1 if W C→C = α
f 1 (1 1 ) = 1 k 1 f 0 (∅), f 1 (1 L ) = β k L f 1 (2 L ) (5a) f 1 (1 i ) = γ k i f 1 (2 i-1 ), 2 ≤ i ≤ L, (5b)
f 1 (2 i ) = k i γ f 1 (1 i ), 1 ≤ i ≤ L -1. ( 5c
)
The solution to Eqs. ( 5) is given by

f 1 (1 i ) = 1 k i , f 1 (2 i ) = 1 γ i = 1, . . . , L -1 1 β i = L (6)
The equations for f 2 (C) involving configurations with one and two particles are presented in the Supplemental Material.

Once we determine the coefficients f n (C) up to a desired order n, we can compute the steady-state average of any observable O(C) by inserting P (C) from Eq. ( 3) and expanding O around α = 0,

O(C) = K n=0 C O(C)f n (C)α n K n=0 C f n (C)α n = ∞ n=0 c n α n . (7)
For example, the first three coefficients c 0 , c 1 and c 2 are given by

c 0 = a 0 b 0 , c 1 = a 1 -c 0 b 1 b 0 , c 2 = a 2 -c 1 b 1 -c 0 b 2 b 0 , (8)
where a n and b n are defined as

a n = C O(C)f n (C), b n = C f n (C). ( 9
)
Notice that the expansion in Eq. ( 7) is slightly different for the current J due to an extra α in Eq. ( 1) and is given by J = ∞ n=0 c n α n+1 . Using the expressions for f n (C) and f n (C) computed earlier yields

J = α 1 - i=1 1 k i + 1 γ i α + O(α 2 ) , (10)
ρ = 1 L L i=1 1 k i + 1 γ i α + O(α 2 ), (11)
ρ i = 1 k i + 1 γ i α + O(α 2 ), (12)
where γ i = γ+(β-γ)δ i,L . These equations constitute our main result. Equation (10) shows that if the initiation rate α is small compared to k 1 , . . . , k L and γ, than the protein production rate J depends predominately on the initiation rate, along with the translocation and elongation rates of the first 10 codons, corresponding to the ribosome footprint . The importance of the first 10 codons is a direct consequence of the excluded volume interactions: any ribosome already present in that region will prevent a new ribosome from binding the mRNA. It is also important to emphasize that (10)-( 12) are exact series expansions around α = 0; the approximation is made only when the series is truncated. Independent Particle Approximation (IPA). Interestingly, the excluded volume interactions between particles have no effect on the first two terms in the series expansion of J. This motivates us to ask how the expansion in Eq. (10) would look like if we assumed that all particles are independent, i.e. not experiencing any exclusion interaction. In our model, the IPA amounts to replacing P (C) with

P IPA (C) = 1 Z L N (C) j=1 w 1 (θ X(j) ), (13)
where N (C) is the number of particles in a configuration C, θ is one of the two particle states 1 and 2, X(i) is the position of the i-th particle on the lattice and Z L is the normalization constant. The weights w 1 (1 i ) and w 1 (2 i ) for i = 1, . . . , L are obtained by solving the master equation for a single particle and are given by α/k i and α/γ i , respectively. The corresponding expressions for the current J and local density ρ i read

J IPA = α i=1 (1 + p i ) , ρ IPA i = p i 1 + p i , (14)
where p i = α(1/k i + 1/γ i ). We will use these results later in order to determine the importance of ribosome collisions in real genetic sequences. We also note that the IPA in ( 13) provides a good approximation to f n (C) for n = N (C) when the particles in C are far apart from each other.

Application to mRNA translation in yeast. We now apply our results to the transcriptome of S. cerevisiae using realistic model parameters. The values of α in the range 0.005 -4.2 s -1 with the median value of 0.09 s -1 have been previously estimated in Ref. [1] using genomewide experimental values of the ribosomal density. We assume that the rates k i are mainly proportional to the gene copy number of tRNAs delivering the corresponding amino acid [28]; the rates are normalized so that the average codon translation rate is equal to the experimental value of 10 codons/s [3]; the estimates of all elongation rates along with the distribution of α is presented in the Supplemental Material. The translocation rate γ is fixed to γ = 35 codons/s [29], and termination is assumed to be fast and comparable to translocation [3], β ≈ γ, so that γ i = γ ∀i in Eqs. ( 10)- (12).

In total, we analyzed 5836 gene sequences; for each gene we calculated J, ρ and ρ i for 1 ≤ i ≤ L up to and including the second order of the perturbative expansion at the corresponding physiological value of α. The results were then compared to the exact values obtained numerically with stochastic simulations using the Gillespie algorithm [31] by calculating the percent error .

For the protein production rate J, with the zeroth order of the perturbative expansion (predicting that J = α) we obtain an error of < 5% for only 11% of the genes. Remarkably, that percentage jumps to 80.7% when the first-order coefficients in Eq. ( 6) are taken into account (Fig. 2). Including the second-order coefficients (computed numerically from Eq. ( 4)) does not significantly improve results, due to a large value of α 0.15 s -1 in about 20% of the genes. Since the coefficients in Eq. (7) typically alternate in sign, truncating the series will ultimately lead to a wrong result when the value of α is large enough. On the other hand, the IPA does not suffer from this problem and leads to an error < 5% in 94% of genes, whereby only 1% of genes have > 20%. The success of the IPA also suggests that ribosome collisions and traffic jams have a minor effect on the rate of translation, which is in accordance with recent experimental evidence [26,32]. This is also apparent from the density profile ρ i , which is very well approximated by the linear approximation in Eq. ( 12), even for larger values of α (Fig. 2 b andc).

Identifying determinants of mRNA translation. Our analytic prediction allows us to decompose the contributions from initiation and elongation to the rate of translation J, thereby addressing a long-standing question about main determinants of protein production rate. Remarkably, the expressions for the current obtained with both the first-order and the independent particle approximation involve only the first = 10 codons. This result therefore strongly indicates that, together with the initiation rate, the first 10 codons of the mRNA are the key determinants of the protein production rate by preventing a new ribosome from binding which effectively decreases the initiation rate.

If we assume that the cell maximizes the rate of protein production, given the above result, we would ex- pect to find a signature in the genome for selecting efficient fast codons at the beginning of each gene. We test this hypothesis by computing the ribosomal current for the fastest (J F ) and the slowest (J S ) set of first synonymous codons for each gene, using the IPA. A score η J = (J -J S )/(J F -J S ) is then assigned to each gene, which is 1 (0) when the sequence corresponds to the fastest (slowest) codon sequence. On the other hand, one might assume that the cell not only tries to maximize protein production rates, but at the same time it tries to minimize the ribosome density ρ on mRNAs. This assumption is motivated by the fact that ribosomes are limiting [22,33] and highly costly in terms of cellular energy resources [34], and therefore ribosome queues are to be avoided. Hence, we also compute η ρ = (ρ-ρ S )/(ρ F -ρ S ) for each gene, where ρ F and ρ S denote the ribosome density for the fastest and slowest set of synonymous codons, respectively.

Figure 3 shows the histogram of η J (yellow) and η ρ (blue) computed for 5836 genes of S.cerevisiae. Both histograms show an average of 0.7 suggesting the selection of fast codons near the start codon to maximize J, as well as an overall selection of fast codons along the mRNA to minimize ρ. However, the width of the distribution of η ρ is substantially smaller than the one of η J . This might indicate that the optimization of protein production rate is strongly dependent on the particular gene, since different proteins are needed at different concentrations. In contrast, the minimization of the number of ribosomes on mRNAs could be a more general constraint.

Conclusions. We have presented an analytic method that allows us to quantify sequence determinants of protein production rates, using a model for translation that is based on an inhomogeneous exclusion process. Our results demonstrate that the rate of protein production is largely determined by the initiation rate and the elongation rates of the first 10 codons (the ribosome footprint length on the mRNA), which control how fast ribosomes load onto the mRNA, whereby ribosome collisions and queues have a minor effect under physiological conditions.

ARGUMENT FOR f (n) (C) = 0 WHENEVER THE NUMBER OF PARTICLES IN A CONFIGURATION C IS LARGER THAN n
As noted in the main text, a key observation in our analysis is that f (n) (C) = 0 whenever the number of particles in C is larger than n. We present here the argument for the case of n = 1; a similar argument applies to higher orders. For n = 1, using Eq. ( 5) in the main text we obtain

f (1) (C) = 1 e(C) C I C →C δ C ,∅ + C (1 -I C →C )W C →C f (1) (C ) -δ C,∅ . (S10)
Consider C to be one full lattice configuration with all N = (L -1 + )/ particles in state 1. The first term on the right hand side of Eq. (S10) is zero because there is no direct transition from the empty configuration ∅ to C.

The second term on the right hand side is also zero, because C can only be accessed through the transition rate α.

Finally, the third term on the right hand side is clearly zero. Hence, f (1) (C) = 0. Now, it is easily seen that all full lattice configurations C that can be accessed only from C will fulfil f (1) (C ) = 0. Iterating further this procedure to the configurations that can be accessed from those, it becomes clear that all full lattice configurations C have f (1) (C) = 0. This argument can be further extended to all configurations with N -1, N -2, N -3, ..., 2 particles, resulting in f (1) (C) = 0 for all of them. The situation changes when we consider configurations with only 1 particle on the lattice, since then there is a direct transition from the empty configuration ∅ to one of those configurations, leading to f (1) (C) = 0. The rest of the coefficients f (1) (C) for 1-particle configurations will depend on each other, and hence, they do not vanish. Therefore, we can conclude that only 1-particle configurations have non-vanishing f (1) (C). For n = 1, this key observation allows us to write Eq. ( 5) taking into account configurations with only one particle and discarding all the others, which yields Eqs. (6a)-(6c) in the main text.

EQUATIONS FOR THE SECOND-ORDER COEFFICIENTS f (2) (C)
As stated in the main text, all second-order coefficients f (2) (C) whereby a configuration C has more than two particles are equal to zero. The equations for the non-zero coefficients f (2)(C) are presented below, where we use the the notation γ j = γ for 1 ≤ j ≤ L -1 and γ L = β.

First, we look at configurations with one particle at site i = 1 and the other at site j = i + l, . . . , L:

f (2) (1 1 1 j ) = 1 k 1 + k j f (1) (1 j ) + (1 -δ 1+l,j )γ j-1 f (2) (1 1 2 j-1 ) (S11a) f (2) (1 1 2 j ) = 1 k 1 + γ j f (1) (2 j ) + k j f (2) (1 1 1 j ) (S11b) f (2) (2 1 1 j ) = 1 (1 -δ 1+l,j )γ 1 + k j k 1 f (2) (1 1 1 j ) + (1 -δ 1+l,j )γ j-1 f (2) (2 1 2 j-1 ) (S11c) f (2) (2 1 2 j ) = 1 (1 -δ 1+l,j )γ 1 + γ j k 1 f (2) (1 1 2 j ) + k j f (2) (2 1 1 j ) . ( S11d
)
where we used the Kronecker delta to account for the excluded volume interaction. Next, we look at configurations with particles at sites i = 2, . . . , L -l and j = i + l, . . . , L.

f (2) (

1 i 1 j ) = 1 k i + k j γ i-1 f (2) (2 i-1 1 j ) + (1 -δ i+l,j )γ j-1 f (2) (1 i 2 j-1 ) (S12a) f (2) (1 i 2 j ) = 1 k i + γ j γ i-1 f (2) (2 i-1 2 j ) + k j f (2) (1 i 1 j ) (S12b) f (2) (2 i 1 j ) = 1 (1 -δ i+l,j )γ i + k j k i f (2) (1 i 1 j ) + (1 -δ i+l,j )γ j-1 f (2) (2 i 2 j-1 ) (S12c) f (2) (2 i 2 j ) = 1 (1 -δ i+l,j )γ i + γ j k i f (2) (1 i 2 j ) + k j f (2) (2 i 1 j ) . (S12d)
Finally, the equations for configurations with only one particle are given by

f (2) (1 1 ) = 1 k 1 γ L f (2) (1 1 2 L ) + f (1) (∅) , (S13a)
f (2) (2 1 ) = k 1 γ 1 f (2) (1 1 ) + γ L γ 1 f (2) (2 1 2 L ), (S13b)
f (2) (1 i ) = 1 k i γ i-1 f (2) (2 i-1 ) + θ(L -l + 1 -i)γ L f (2) (1 i 2 L ) -f (1) (1 i ) , i = 2, . . . , L (S13c)
f (2) (2 i ) = 1 γ i k i f (2) (1 i ) + θ(L -l + 1 -i)γ L f (2) (2 i 2 L ) -f (1) (2 i ) , i = 2, . . . , L (S13d)
where θ(n) = 0 for n < 0 and 1 for n ≥ 0.

The equations (S11)-(S13) can be easily solved numerically by iteration. We first solve Eqs. (S11) for i = 1, j = l+1. We then iterate the recursion relation in Eqs. (S11) for fixed i = 1 and j = 2 + l, . . . , L. We then solve equations (S12) for fixed i = 2 and j = l + 3, . . . , L and repeat this procedure until i = L -l and j = L. Finally, we calculate the coefficients f (2) (1 i ) and f (2) (2 i ) for the single-particle configurations using f (2) (1 i 2 L ) and f (2) (2 i 2 L ) that we solved in the previous steps. We note that the only unknown that we cannot determine is f (1) (∅), but it turns out that all terms containing f (1) (∅) will cancel out later when we calculate the coefficients in the series expansions ( 11)-( 13) in the main text.

ESTIMATES OF THE CODON ELONGATION RATES

The total translation rate ω i of the codon i is:

1 k i + 1 γ = 1 ω i , (S14)
where the translocation rate γ is set to 35s -1 as mentioned in the main text. The codon elongation rate k i , which represents the average arrival and recognition time of the cognate tRNA, is determined by following the procedure introduced in [1], which we report here.

For each of the 41 tRNAs types j we consider their gene copy number (GCN), which allows us to provide a first estimate of the rate k j :

k j = r GCN j 41 j=1 GCN j , (S15)
where GCN j is the gene copy number of the tRNA of type j with j = 1, . . . , 41, and r is a proportionality constant. These rates were then adjusted to take into account experimental evidence suggesting that the translation rates of codons using the G-U wobble are reduced by 39% compared to their G-C counterparts; analogously, codons using the wobble I-C and codons using the wobble I-A are reduced by 36% relative to their I-U counterparts [2]. To calculate the proportionality constant r, we used the experimental value of 10 codons/s for the average codon translation rate ω i defined as

ω i = 61 i=1 k i γ k i + γ n i n (S16)
where n i /n is the relative abundance of all different codon types in the cell, n i is total number of codons in the cell of exactly type i, and n = 61 i=1 n i (the three STOP codons are excluded) [1]. Table I summarises the resulting elongation rates k i for 61 codons.

The estimates provided are valid in physiological conditions, when amino acids are not limiting. In case of amino acid starvation the rates of the codons affected should be modified (the elongation rates depend on the abundance of tRNAs charged with the correct amino acid).

DISTRIBUTION OF THE TRANSLATION INITIATION RATES

The translation initiation rates α were estimated for each gene in Ref. [1] by matching the average ribosome density predicted by the model with the experimental value of the ribosome density obtained in Ref. [3] for 5836 genes that we analyzed in the main text are in the range 0.005 -4.2 s -1 with the median value of 0.09 s -1 .

The histogram of the rates is presented in Figure S1, showing that the initiation rate α for the majority of genes is indeed much smaller than the smallest elongation rate, which has the value of 1.02 codons/s for the CUU codon (see Table I). The list of all translation initiation rates can be found in the Supplemental Material of Ref. [1].

  and is 0 otherwise, and e 0 (C) = C (1-I C,C )W C→C is the total exit rate from C excluding the rate α.

FIG. 2 :

 2 FIG. 2: (a): Histogram of the percent error measuring the discrepancy between the protein production rate obtained by stochastic simulations and the perturbative expansion including the independent-particle approximation (IPA). (b) and (c): Ribosomal density profile obtained by stochastic simulations (black) compared to Eq. (12) (red) for two values of the initiation rate α, one close to the median value ≈ 0.09 s -1 (gene YAL045C, b) and the other that is close to the 90th percentile value ≈ 0.23 s -1 (gene YGL034C, c).

FIG. 3 :

 3 FIG. 3: Histogram of ηJ (yellow) and ηρ (blue) for the S. cerevisiae genome.

  FIG.S1: Histogram of the estimated initiation rates for 5836 genes from the S. cerevisiae genome. The 95th percentile value of α is ≈ 0.3, which is much smaller that the smallest elongation rate, which is 1.02 codons/s for the CUU codon (vertical dashed line). Only 15 genes have α larger than mini{ki} = 1.02 codons/s.

  . The resulting values

	tRNA anti-codon codon GCN ki [1/s]	tRNA anti-codon codon GCN ki [1/s]
	Ala1 Ala1 Ala2 Ala2 Arg1 Arg2 Arg2 Arg2 Arg3 Arg4 Asn Asn Asp Asp Cys Cys Gln1 Gln2 Glu3 Glu4 Gly1 Gly1 Gly2 Gly3 His His Ile1 Ile2 Ile2 Leu1 Leu1	IGC IGC UGC UGC CCU ICG ICG ICG UCU CCG GUU GUU GUC GUC GCA GCA UUG CUG UUC CUC GCC GCC UCC CCC GUG GUG UAU IAU IAU UAG UAG	GCU 11.00 18.34 GCC 11.00 11.74 GCA 5.00 8.34 GCG 5.00 8.34 AGG 1.00 1.67 CGU 6.00 10.01 CGC 6.00 6.40 CGA 6.00 6.40 AGA 11.00 18.34 CGG 1.00 1.67 AAU 10.00 10.17 AAC 10.00 16.68 GAU 15.00 15.26 GAC 15.00 25.01 UGU 4.00 4.07 UGC 4.00 6.67 CAA 9.00 15.01 CAG 1.00 1.67 GAA 14.00 23.35 GAG 2.00 3.34 GGU 16.00 16.28 GGC 16.00 26.68 GGA 3.00 5.00 GGG 2.00 3.34 CAU 7.00 7.12 CAC 7.00 11.67 AUA 2.00 3.34 AUU 13.00 21.68 AUC 13.00 13.87 CUA 3.00 5.00 CUG 3.00 5.00	Leu3 Leu4 Leu5 Leu5 Lys1 Lys2 Met Phe Phe Pro1 Pro1 Pro2 Pro2 Ser2 Ser2 Ser3 Ser3 Ser4 Ser5 Thr1 Thr1 Thr2 Thr3 Trp Tyr Tyr Val1 Val1 Val2 Val2b	CAA UAA GAG GAG CUU UUU CAU GAA GAA UGG UGG IGG IGG IGA IGA GCU GCU UGA CGA IGU IGU CGU UGU CCA GIA GIA IAC IAC UAC CAC	UUG 10.00 16.68 UUA 7.00 11.67 CUU 1.00 1.02 CUC 1.00 1.67 AAG 14.00 23.35 AAA 7.00 11.67 AUG 5.00 8.34 UUU 10.00 10.17 UUC 10.00 16.68 CCA 10.00 16.68 CCG 10.00 16.68 CCU 2.00 3.34 CCC 2.00 2.13 UCU 11.00 18.34 UCC 11.00 11.74 AGU 4.00 4.07 AGC 4.00 6.67 UCA 3.00 5.00 UCG 1.00 1.67 ACU 11.00 18.34 ACC 11.00 11.74 ACG 1.00 1.67 ACA 4.00 6.67 UGG 6.00 10.01 UAU 8.00 8.14 UAC 8.00 13.34 GUU 14.00 23.35 GUC 14.00 14.94 GUA 2.00 3.34 GUG 2.00 3.34

TABLE I :

 I Elongation rates ki considering supply (gene copy number of tRNAs) and wobble base-pairing.

* Electronic address: jszavits@staffmail.ed.ac.uk

Acknowledgments. MCR and LC contributed equally to this work. JSN was supported by the Leverhulme Trust Early Career Fellowship. MCR was supported by the Biotechnology and Biological Sciences Research Council (BBSRC) BB/N017161/1 and the Scottish Universities Life Sciences Alliance. LC thanks the CNRS for having granted him a demi-délégation (2017-18).

The exact solution, Eq. ( 3) in the main text, pertains to any ergodic Markov jump process with a finite number of states. In order to derive the exact solution, let us write again the steady-state master equation in a matrix form

where P is a column vector whose elements are steady-state probabilities P i , the N states are indexed by i = 1 . . . , N and the N × N matrix M is given by

where W ij is a transition rate from a state i to j. We note from Eq. (S2) that the sum of all elements in each column of M is zero

Consequently, the sum all row vectors of M is zero, which means that the vectors are linearly dependent and thus the determinant of M is equal to zero, detM = 0. (S4)

Combining Eq. (S4) with the Laplace expansion of a determinant yields

Here C ij is a cofactor of M defined as C ij = (-1) i+j detM (i,j) , where the matrix M (i,j) is obtained from M by removing i-th row and j-th column . Inserting Eq. (S2) into (S5) gives detM =

from which we conclude that

for any i and j, i.e. the cofactor C ij does not depend on the state i. Inserting Eq. (S7) back into Eq. (S4) gives

which is precisely the starting steady-state master equation, Eq. (S1). We thus conclude that

N j=1 detM (j,j) . (S9)


Siyuan Niu
email: siyuan.niu@lirmm.fr

Aida Todri-Sanial

Analyzing crosstalk error in the NISQ era

Keywords: Quantum computing, Crosstalk, Error mitigation

Noisy Intermediate-Scale Quantum (NISQ) hardware has unavoidable noises, and crosstalk error is a significant error source. When multiple quantum operations are executed simultaneously, the quantum state can be corrupted due to the crosstalk between gates during simultaneous operations, decreasing the circuit fidelity. In this work, we first report on several protocols for characterizing crosstalk. Then, we discuss different crosstalk mitigation methods from the hardware and software perspectives. Finally, we perform crosstalk injection experiments on the IBM quantum device and demonstrate the fidelity improvement with the crosstalk mitigation method.

I. INTRODUCTION

Quantum Computing is expected to solve specific classical intractable problems. Several companies such as IBM, Google, and Rigetti have released their quantum chips with 65, 72, and 31 qubits, respectively. However, these chips are qualified as Noisy Intermediate-Scale Quantum (NISQ) hardware, which have less than one hundred qubits and suffer from unavoidable noises.

Crosstalk is one of the major noise sources not only in superconducting but also in trapped-ion devices. It can corrupt the qubit state when multiple quantum operations are executed simultaneously. Crosstalk has a significant impact on quantum gate error. For instance, [START_REF] Murali | Software mitigation of crosstalk on noisy intermediate-scale quantum computers[END_REF] shows an increase of CNOT error up to 11 times caused by crosstalk. Different protocols were proposed in [START_REF] Radoslaw | Quantum process tomography of a universal entangling gate implemented with josephson phase qubits[END_REF]- [START_REF] Jay M Gambetta | Characterization of addressability by simultaneous randomized benchmarking[END_REF] to detect and characterize crosstalk in quantum devices. After assessing crosstalk, there are two types of methods to mitigating it. The first one is based on hardware strategies, such as tunable coupling [START_REF] Mundada | Suppression of qubit crosstalk in a tunable coupling superconducting circuit[END_REF], or frequency allocation [START_REF] Li | Towards efficient superconducting quantum processor architecture design[END_REF]. From the software perspective, crosstalk has been mitigated by changing simultaneous CNOT operations with high crosstalk to execute separately while trading off the increase of the decoherence time as in [START_REF] Murali | Software mitigation of crosstalk on noisy intermediate-scale quantum computers[END_REF].

In this work, we first introduce how to characterize a quantum device's crosstalk using protocol Simultaneous Randomized Benchmarking (SRB) [START_REF] Jay M Gambetta | Characterization of addressability by simultaneous randomized benchmarking[END_REF]. Second, we present various state-of-the-art crosstalk mitigation methods. Third, we inject crosstalk on IBM quantum devices to show its impact on output fidelity. Finally, we evaluate the crosstalk mitigation method using several benchmarks to demonstrate the fidelity improvement.

This work is funded by the QuantUM Initiative of the Region Occitanie, France, University of Montpellier, France, and IBM Montpellier, France.

II. METHODS

A. Crosstalk characterization using SRB

To characterize the crosstalk effect of a quantum device, we choose the most commonly used protocol -Simultaneous Randomized Benchmarking (SRB) [START_REF] Jay M Gambetta | Characterization of addressability by simultaneous randomized benchmarking[END_REF] based on its ability to quantify the impact of parallel instructions, such as the crosstalk effect of one quantum gate to another when they are executed at the same time. To characterize the crosstalk effect between gate g i and g j , we first perform Randomized Benchmarking (RB) on both gates separately and obtain their independent error rate E(g i ) and E(g j ). Then, applying SRB on both gates yields the correlated error rate E(g i |g j ) and E(g j |g i ) for simultaneous executions. If there exists crosstalk between them, the relation between independent and correlated errors should comply with

E(g i |g j ) > E(g i ) or E(g j |g i ) > E(g j ).
Previous work [START_REF] Murali | Software mitigation of crosstalk on noisy intermediate-scale quantum computers[END_REF] has demonstrated that crosstalk is significant due to simultaneous CNOT executions, hence, in this work, we only focus on characterizing the crosstalk effect between CNOT pairs. We use the ratio of correlated error to independent error r(g i |g j ) = E(g i |g j )/E(g i ) as the indicator of the crosstalk effect of the CNOT pairs. We choose IBM Q 7 Casablanca as an example to show the crosstalk effect on this device. As CNOT errors vary over each calibration, we characterize the crosstalk effect twice respectively on each possible simultaneous CNOT pair of IBM Q Casablanca.

The results of average crosstalk data are presented in Fig. 1. The crosstalk effect remains stable regardless of the variation of gate errors across days. For each experiment, there are 150 Cliffords in the sequence and each sequence is repeated five times due to the expensive runtime cost. The crosstalk effect is not as severe in this device as in other devices such as IBM Q 20 Poughkeepsie, where the gate error grows up to 11 times caused by crosstalk [START_REF] Murali | Software mitigation of crosstalk on noisy intermediate-scale quantum computers[END_REF]. However, the error rate is still amplified up to 3 times in IBM Q Casablanca. Sometimes ratios are less than one, which is probably due to the limits of crosstalk metric in SRB protocol [START_REF] David C Mckay | Correlated randomized benchmarking[END_REF].

B. Crosstalk mitigation

After characterizing the crosstalk effect, the next challenge is how to mitigate it. The hardware-oriented crosstalk mitigation strategies have several physical constraints. For example, the tunable coupling [START_REF] Mundada | Suppression of qubit crosstalk in a tunable coupling superconducting circuit[END_REF] is only applicable for tunable coupling superconducting devices, whereas the frequency allocation [START_REF] Li | Towards efficient superconducting quantum processor architecture design[END_REF] Fig. 2: Crosstalk injection of CSWAP circuit. We introduce two other qubits (q 3 , q 4 ) and apply CNOT operations on them. Barriers are inserted to make sure the CNOTs are executed simultaneously. Note that, crosstalk injection are performed on CNOT pairs with strong crossalk error.

is mainly designed for fixed-frequency transmon qubit devices. Despite hardware improvement approaches to reduce crosstalk, there are also software methods to address crosstalk. Here, we report on hardware-agnostic software crosstalk mitigation approaches.

An intuitive approach for software crosstalk mitigation is to make simultaneous CNOTs execute serially to avoid simultaneous high crosstalk CNOTs. However, serial instructions can increase the circuit depth and cause decoherence errors. Therefore, a scheduler is required to trade-off between them. The state of the art [START_REF] Murali | Software mitigation of crosstalk on noisy intermediate-scale quantum computers[END_REF] proposes a crosstalk-adaptive scheduler (labeled as XtalkSched) based on SMT optimization that reschedules a quantum intermediate representation (IR), taking into account both crosstalk and decoherence error. It inserts barriers between the simultaneous CNOTs with a significant crosstalk effect to avoid parallel instructions while accounting for the qubit coherence time.

III. EVALUATION A. Methodology

Benchmarks. We first use a 3-qubit CSWAP gate circuit and inject a different number of simultaneous CNOTs to show the impact of crosstalk on the output fidelity of a quantum circuit. An example of demonstrating crosstalk injection by enabling two simultaneous CNOTs is shown in Fig. 2. Then, we examine benchmarks collected from the previous works [START_REF] Li | Qasmbench: A low-level qasm benchmark suite for nisq evaluation and simulation[END_REF], [START_REF] Murali | Software mitigation of crosstalk on noisy intermediate-scale quantum computers[END_REF], including SWAP circuits, Bell state circuit, etc., to show the performance of XtalkSched. For each benchmark, we select its mapping to ensure it includes at least one pair of CNOTs with a high crosstalk effect.

Algorithm configuration. For crosstalk characterization of IBM Q Casablanca, we set the crosstalk threshold to 2, and w used in XtalkSched is set to 0.5 to trade-off circuit depth and crosstalk mitigation. The Qiskit version is 0.23.6.

Comparison. We compare XtalkSched [START_REF] Murali | Software mitigation of crosstalk on noisy intermediate-scale quantum computers[END_REF] with ParSched, which is the current scheduler used in Qiskit to make instructions execute in parallel without considering the crosstalk.

B. Experimental results

The experimental results of crosstalk injection are shown in Fig. 3. The probability of obtaining the right output state is reduced with the increase of simultaneous CNOTs due to the injected crosstalk error. In the worst case, the fidelity is decreased by 33.3% compared to the non-crosstalk case.

The comparison between XtalkSched with ParSched in terms of output fidelity and circuit depth is shown in Fig. 4. The fidelity is improved by 9.4%, whereas the circuit depth increased by 39.1%. Although the fidelity is enhanced because of XtalkSched, the circuit depth increase is not negligible. This raises the need for more effective and practical crosstalk mitigation methods.

IV. CONCLUSION

Crosstalk is a non-negligible error source in quantum computers. This work examines the crosstalk characterization protocol based on simultaneous randomized benchmarking, and we report on the crosstalk effects on the IBM quantum [START_REF] Radoslaw | Quantum process tomography of a universal entangling gate implemented with josephson phase qubits[END_REF][START_REF] Li | Qasmbench: A low-level qasm benchmark suite for nisq evaluation and simulation[END_REF] represents the SWAP circuit which aims to realize a connection between q 0 and q 4 through SWAP operations.

device. We also evaluate the crosstalk-adaptive scheduler and illustrate its impact on output fidelity and circuit depth.

Fig. 1 :

 1 Fig.1:The SRB result of IBM Q 7 Casablanca. The number represents the ratio of correlated error to independent error. Note that the SRB experiment is performed on CNOTs that can be executed in parallel, which means they do not share the same qubits. There are no SRB experiments on blank spots due to the sharing qubits.

Fig. 3 :

 3 Fig. 3: The output state probability distribution of CSWAP gate circuit. The right output state should be '101'. We inject a different number of simultaneous CNOTs to this circuit. The baseline circuit has zero simultaneous CNOT.

Fig. 4 :

 4 Fig. 4: (a) Fidelity. (b) Circuit depth. Note that[START_REF] Radoslaw | Quantum process tomography of a universal entangling gate implemented with josephson phase qubits[END_REF][START_REF] Li | Qasmbench: A low-level qasm benchmark suite for nisq evaluation and simulation[END_REF] represents the SWAP circuit which aims to realize a connection between q 0 and q 4 through SWAP operations.