Dynamical Gaussian state transfer with quantum error correcting architecture

Transferring a quantum state of a light field to a memory is of particular importance. However, this transfer is usually hampered because the memory system is subjected to some noise and this can limit the performance of the state transfer to a great extent. In this paper, we consider the transfer of a Gaussian state of light to a linear medium memory such as an opto-mechanical oscillator and propose a dynamical feedback controller that suppresses the noise in the memory system. To protect an unknown state, the feedback scheme employs the specific configuration of the quantum error correction; that is, a three-mode Gaussian state having appropriate syndromes is taken as the input. Correspondingly, the memory consists of three independent linear systems. The syndrome errors are estimated continuously in time through the measurement of the output field, and the results are then fed back to control the system. Because the input is Gaussian and the systems are all linear, it is possible to formulate the problem using the framework of the celebrated classical Kalman filtering and linear quadratic Gaussian control. A numerical simulation demonstrates the effectiveness of the control scheme.

In reality, however, the memory performance is essentially limited by environmental noise that inevitably occurs during the transfer process. In a real experiment several state-of-the-art techniques should be, at least implicitly, employed to suppress such noise. However, to the best of our knowledge there exists no systematic method for suppressing the noise. Therefore, control theory for quantum memory needs to be explored.
In working out this subject, a key fact is that due to the dynamical noise the information contained in the memory state is gradually erased, rather than that the noise suddenly vanishes it. Hence, if such an unwanted change of state can be monitored continuously in time, then the measurement result could be fed back to suppress the noise. In fact, in [2], the probe light field for state transfer is continuously measured at the terminal of the optical path to indirectly detect the error during the transfer process. In this paper, we propose a dynamical control scheme that uses such continuous measurement results not only for the error detection but further for * E-mail address: go@z8.keio.jp † E-mail address: yamamoto@appi.keio.ac.jp feedback control to obtain the better performance of the state transfer. Fortunately, measurement-based quantum feedback control theory is well developed [23,24,25] and has many practical uses [26,27,28,29,30,31,32]; hence, it will be also useful for our scheme.
Here we describe two specific features of the proposed system-controller configuration. The first one is that the memory is a linear medium system such as a collective atomic ensemble or an opto-mechanical oscillator. In addition, the input state to be transferred is assumed to be a Gaussian state of light [33], in which case the memory state also becomes Gaussian. Then, the feedback control scheme has the form of celebrated linear quadratic Gaussian (LQG) control with Kalman filtering [24,26,29,30,31]. The second feature is that the controller employes the schematic of quantum error correction (QEC) for the purpose of protecting an unknown input state. The basic idea of QEC is to encode an unknown state into an enlarged codespace so that an appropriate syndrome measurement yields sufficient information about the error, which can then be used to correct the error [34,35,36,37,38,39]. Use of this QEC strategy for quantum memory means that the input is a three-mode unknown Gaussian state having appropriate syndromes, and it is transferred to the memory that is correspondingly enlarged; during the transfer process, the output light fields are measured continuously in time, yielding the estimate of the syndrome errors that can be fed back to control the system.
The effectiveness of this control scheme is actually expected from the fact that the dynamical feedback control has successfully been applied to some QEC problems [40,41,42,43,44,45,46]. It should be noted that the problem considered in this paper is not a QEC problem itself; in QEC, the system state can be protected with the use of ancillary apparatus, whereas in our formalism, the initial state of the memory system cannot be protected, but rather the input state to be transferred to the system is what should be protected. This observation is vital in the sense that our scheme does not violate the no-go theorem for Gaussian QEC [47].
The paper is organized as follows. In Section II, as a preliminary, we describe a single-mode Gaussian state transfer. The dynamical control scheme is addressed in two sections. First, the state preparation and transfer processes are given in Section III. Section IV explains the error detection scheme (Kalman filter) and the feedback controller (LQG controller). Finally, in Section V, we numerically demonstrate the efficiency of our method by taking, as an example of the memory, an opto-mechanical oscillator manipulated under ultra-low temperature.

II. SINGLE-MODE GAUSSIAN STATE TRANSFER
FIG. 1: (Color online) Opto-mechanical oscillator as an example of the single-mode linear quantum memory. The right end-mirror of the Fabry-Perot cavity serves as the oscillator, whereâ andb are annihilation operators corresponding to the intra-cavity and oscillator modes, respectively. Adiabatically eliminating the cavity modeâ, Eq. (1) is obtained. F denotes a Faraday isolator that facilitates one-way coupling between the input and cavity modes.
As a preliminary, this section is devoted to describe a single-mode Gaussian state transfer from a light field to a memory, where the former is a coherent CW laser field and the latter is served by a linear system such as an opto-mechanical oscillator shown in Fig. 1 (see e.g., [14,17,48] for a detailed discussion). The input is a coherent state |α in , withâ in the annihilation operator of this coherent light field. The reflected output modê a out could be used to correct the errors occurring during the transfer process. In what follows we describe the dynamics of the memory system and show how much the input state |α in is transferred to the memory.
A general single-mode open linear dynamics is described by the following quantum Langevin equation: whereb denotes the system annihilation operator. Here, ν represents the coupling constant between the system and the input coherent light fieldâ in = α in +â 0 , where α in is the mean andâ 0 is the field annihilation operator. Moreover, we assume that the system couples with an unwanted single-mode noisy environment represented by the annihilation operatorξ with mean zero. Γ is the coupling strength.
Let us now take the white noise approximation on the outer field modesâ 0 andξ; i.e., â 0 (t)â † 0 (t ′ ) = δ(t − t ′ ) and ξ (t)ξ † (t ′ ) = (n + 1)δ(t − t ′ ), where n > 0 represents the strength of the noise. In the case of thermal noise, n is the averaged photon number. This approximation allows us to represent the dynamics ofb in terms of the Ito-type quantum stochastic differential equation (QSDE) [49,50]: whereÂ 0 (t) = t 0â 0 (s)ds andΞ(t) = t 0ξ (s)ds are the so-called quantum Wiener processes. Their infinitesimal increments satisfy the following quantum Ito-rule: In this QSDE representation the input mode is written by dÂ in = α in dt + dÂ 0 . Since the field state is the vacuum, we have dÂ in = α in dt. Also note that dΞ = 0. Now, let us see the steady state of the memory system. Using the above Ito rule we readily obtain the ordinary differential equations of b(t) and b(t) † b(t) , which give b (∞) = −2 √ να in /(ν + Γ) and 2i denote the dimensionless position and momentum operators of the system, respectively. This clearly illustrates that in the long time limit the system certainly acquires the information of the field input state |α in in the mean sense (note that ν and Γ are assumed to be known). However the fluctuation of the state could become much bigger than the vacuum fluctuation 1/2; in this case the state wrote down to the memory is far away from pure and has very small overlap with the input state.

III. PREPARATION AND TRANSFER PROCESSES OF THE INPUT STATE
In this paper, we study the system shown in Fig. 2. First, in the left box, the input state, which is a threemode Gaussian state of light, is prepared. This input state is transferred into the memory shown in the middle box in the figure, that correspondingly consists of three identical linear systems; in the figure, opto-mechanical oscillators are particularly depicted as the memory system. Finally, the output fields are measured, and the results are used for feedback control.
A notable feature of this transfer scheme is that the input state is allowed to contain some unknown parameters. More specifically, the input state is generated by embedding an unknown single-mode fieldÂ in into the three-mode light fields with the use of ancilla squeezed fields and beam splitters. This is no more than the schematic of QEC; actually,Â in corresponds to a source mode to be protected, and for this purpose it is encoded into a larger Hilbert space where we can construct appropriate syndromes that do not depend on the unknown source mode but only contain information about the errors. Those errors could be corrected by feedback control through the syndrome measurement depicted in the lower box of the figure. Note that the three-mode encoding allows us to detect only either the position or the momentum error, but the proposed scheme can be extended to a nine-mode code to detect errors acting on both the position and momentum operators [39].
In this section, we describe the above-mentioned state preparation and transfer processes, and then evaluate in detail the memory performance, when without any error correction.

A. State preparation
For discrete variable systems, a series of CNOT gates is often used for encoding. A possible continuous-variable analogue is realized in linear optics circuit shown in the left box of Fig. 2; first, the information source modeÂ in is mixed via a 1 : 2 beam splitter with an ancilla squeezed vacuum field denoted byÂ 2 , and one of the outputs is further combined with the second ancilla squeezed vacuum fieldÂ 3 . This combination of beam splitters is called the tritter [39].
Let us describe the above encoding process in detail. As in the single-mode case discussed in Sec. II, we write the source mode as dÂ in = α in dt + dÂ 1 , where α in and A 1 are the mean amplitude and the quantum fluctuation, respectively. In particular, we assume implying that the momentum element of the mean is known to be zero. As mentioned above, in this case, the three-mode encoding is sufficient to detect the position errors. The quantum fluctuationÂ j (j = 1, 2, 3) is a stochastic process satisfying dÂ j = 0 and the following quantum Ito rule: The fluctuation parameters N 1 and M 1 of the source mode are unknown. For the ancilla modesÂ 2 andÂ 3 , we set which satisfy N j (N j + 1) = M 2 j ; that is, the ancilla modes are the identical pure squeezed vacuum fields with squeezing parameter µ. Let us collect the field quadratures in a single vector aŝ 2i. Then, through the tritter,Ŵ 1 becomes dŴ ′ 1 = βdt + T dŴ 1 , where T is the orthogonal matrix corresponding to the combination of the above-mentioned 1 : 2 and half beam splitters: Using the Ito rule, it is found that the position quadra- This quantity is essentially equivalent to the power spectrum density of (Q ′ i −Q ′ j ) 2 at the center frequency of the career laser field [51]. We then find that Eq. (4) becomes zero for every i, j when taking the limit µ → −∞. This means that, in the Schrodinger picture, this ideal input state generated through the tritter lives in the codespace spanned by {|Q, Q, Q } and is of the following GHZ-like form : where ψ(Q) is a Gaussian wave function of the source state. The tilde means that the state is the ideal one. Next let us focus on the following quantity: which leads to IfÂ 1 ,Â 2 , andÂ 3 are all coherent fields, we have dP fd /dt = 7.5, which thus can be interpreted as the classical limit. This means that a non-classical input state is generated when dP fd /dt < 7.5, which is now equivalent to that µ < 0 and M 1 ≥ 0. Moreover, it is known that a symmetric state such as the one taken here is entangled if dP fd /dt < 6 [52,53]. For instance setting a coherent source state, i.e., N 1 = M 1 = 0, and taking the ideal limit µ → −∞, we obtain dP fd /dt → 4.5; hence in this case the GHZ-like state (5) is indeed entangled.

B. Transfer process
Let us next describe the transfer process of the modê W ′ 1 to the memory served by the three identical linear systems. The dynamics of the memory is given by the combination of Eq. (2) with the input field replaced by dŴ ′ 1 ; that is, the vector of quadratures, 2i, satisfies the following QSDE: whereŴ = (Ŵ 1 ,Ŵ 2 ) ⊤ and The fluctuation vectorŴ 1 is given by Eq. (3) and is the vector of noise processes to which the linear systems are subjected; henceΞ j (j = 1, 2, 3) satisfies Again, n > 0 represents the noise strength. Finally, u ∈ R 6 represents the control input for the memory. Note here it is assumed that both the quadratures of each linear system can be controlled. In what follows of this section we set u = 0. The state transfer can be evaluated in terms of the mean vector and the covariance matrix, as discussed in Sec. II. First, the mean vector x = ( q 1 , . . . , . That is, the system state becomes Gaussian with mean vector parallel to β, implying that the memory certainly acquires information of the input state in the mean sense (note that ν and Γ are assumed to be known). Next let us consider the covariance matrix of the system: Using the Ito rule, the time evolution of V is found in the following form called the Lyapunov differential equation: where Note Λ 2 = Λ 3 by assumption. Now, to explicitly evaluate the state transfer, let us take the source state to be coherent, in which case M 1 = N 1 = 0; then the steady solution of Eq. (8) is obtained as Note V 2 = V 3 . Then we have the explicit form of the fidelity between the input state |ψ in and the steady state of the memory with its mean appropriately displaced: where V in = T ΛT ⊤ is the covariance matrix of the input state |ψ in . Note that Eq. (5) is the ideal limit of |ψ in .
Eq. (9) shows that, when Γ = 0, we have F = 1 without respect to the value of µ. This means that the input state is perfectly transferred into the memory when it is not coupled to the noisy environment. But in reality the system must be subjected to some noise; QEC is expected to overcome this issue. Motivated by some continuous variable QEC protocols found in the literature [36,37,38,39], let us take the operatorq i −q j (i, j = 1, 2, 3) as a syndrome. Actually, this has desirable properties shown as follows. First, it can be seen from the structure of the system dynamics (7) that the syndrome does not contain the unknown parameters α in , N 1 , and M 1 . Second, in the ideal limit µ → −∞, the variance of each syndrome is given by which takes zero only when the system does not couple to the noisy environment. This means that, if in the small interval [t, t + dt) only one of the three linear systems is subjected to the noise, simultaneous "measurement" of the syndrome operators can detect in which system that error has occurred and how much it is. (The reason why the double quotation is taken is that the syndrome operator cannot be directly measured in our setting, but it needs to be estimated; see the next section.) Because of these two properties, in the limit µ → −∞, the operatorŝ q i −q j certainly plays a role of a syndrome. Here it should be noted that F takes the maximum value when µ = 0; that is, without error correction the encoding process merely degrades the performance of the state transfer [54]. Hence we really need the correction process. Before closing this section, let us evaluate the entanglement of the system state, using essentially the same quantity as Eq. (6). Again in the case N 1 = M 1 = 0, we have (4.5 + 3e µ )ν ν + Γ + (7.5 + 15n)Γ ν + Γ .
Hence, together with ∆(q i −q j ) 2 = 0 in this case, the memory state becomes an entangled GHZ-like state. But in the realistic situation with Γ > 0 such entanglement can vanish. Particularly in the ideal case µ → −∞ the sufficient condition for entanglement, P sys < 6, leads to n < 0.1(ν/Γ) − 0.1; that is, in the case of thermal noise, the averaged photon number of the environment field must be less than about one order of magnitude below the S/N rate ν/Γ.

IV. SYNDROME FILTER AND DYNAMICAL FEEDBACK CONTROL
In the previous section it was found thatq i −q j serves as the syndrome through which the error can be detected. Since this observable cannot be measured directly in our setting, it should be appropriately estimated through indirect measurement. In this section we first describe a continuous-time estimator, i.e., the filter, for the syndromes. We then present the dynamical feedback controller that is based on the syndrome filter.

A. Syndrome Filter
Let us focus on the following two kinds of operator vectors:ŝ HereB 1 := Z 1 T ⊤ andB 2 := Z 2 T ⊤ are isometric matrices with If only one of the three linear systems is subjected to the noise,ŝ 2 can detect in which system that error has occurred. Hence,ŝ 2 serves as the syndrome operator. The first element ofŝ 1 is used to evaluate entanglement of the memory state. Our task is to estimate the operatorŝ 1 orŝ 2 through certain indirect measurement. The lower box in Fig. 2 depicts an optical configuration that achieves this goal for the case ofŝ 1 ; that is, we construct the tritter acting on the output fields that are reflected at the field-system coupler. The observables to be measured by homodyne detectors with appropriate LO phase are then given by Hence measuringŶ 1 actually implies the indirect measurement ofŝ 1 . Note this contains the source fluctuation dP 1 , hence in this case only the mean value of the input state, β, can be kept unknown. On the other hand, for the case ofŝ 2 , two of the three output fields are measured; Y 2 does not contain both the mean and covariance of the source state, thus we can estimate the syndromes without respect to the unknown parameters. Note here that the position squeezing of the ancilla fields reduces the fluctuation dQ 2 and dQ 3 ; in particular taking the ideal limit µ → −∞ we have dŶ 2 = √ 2νŝ 2 dt, thus this is no more than the syndrome measurement.
We now describe the filter ofŝ i , which is constructed with the measurement results ofŶ i (i = 1, 2). Let us write Eq. (13) or (14) in the following general form: For the case of estimatingŝ 1 , the matrices correspond to ). For the simple notation, we do not put the index i onŶ , C, and D. The continuous measurement ofŶ enables us to construct the filter for the dynamics (7): where Z implies Z 1 or Z 2 . The set of above equations is called the quantum Kalman filter [24,26,29,30,31]. The 6-dimensional c-number vector π(x t ) = E(x t |Y t ) represents the quantum conditional expectation ofx t conditioned on the set of measurement data Y t = {y s | 0 ≤ s ≤ t} with y s the measurement result ofŶ s . Note π(x t ) is the least mean square estimate ofx t . V c is the conditional covariance matrix V c = π(∆x∆x ⊤ +(∆x∆x ⊤ ) ⊤ )/2 with ∆x =x − π(x), and it obeys the following Riccati differential equation: where K is given by Eq. (17). Note this is not a stochastic process. The filtering equation for the syndromeŝ 1 orŝ 2 can be directly obtained from Eq. (16) as follows: Due to the coupling to the noisy environment, the system state must escape from the codespace spanned by {|q, q, q }, but this error can be detected by estimating the syndrome operators. Therefore, it is expected that a feedback controller minimizing the estimated value of the syndromes corrects that error. Note that the controller must be a dynamical one so that it can deal with the dynamical noise; hence, the minimization should be carried out through the whole transfer process. Fortunately, because of the linearity of both the dynamics and the output equation, this requirement is satisfied when employing the LQG control strategy [24,26,29,30,31]. That is, we can construct an optimal control law u * that minimizes the following quadratic-type cost function: where R = R ⊤ > 0 represents the penalty for the control input. Corresponding to the syndrome operatorsŝ 1 or s 2 , the weighting matrix Q is set to: It is immediately seen that Therefore, for both cases the optimal control tries to minimize the syndrome errors. In the case ofŝ 1 , the optimal control further takes into account the effect of entanglement and tries to keep the GHZ-like form of the state. The LQG control theory gives the explicit form of the (stationary) optimal controller: where P is the solution to the following algebraic Riccati equation:Ã Let us here assume that R = rI 6 > 0. Then, corresponding toŝ 1 andŝ 2 , the matrix P is respectively obtained in the following form: where We here give two brief remarks. (i) The first one is about the structure of the optimal feedback control input u * . In the case ofŝ 1 , it is given by The first element can be regarded as a linearization to the following sign function (λ = f 2 /3r): λ sgn(π(q 2 −q 1 )) + sgn(π(q 3 −q 1 )) , where sgn(x) = +1 if x > 0 and −1 otherwise. Hence, in the small interval [t, t + dt), if the position shift ofq 1 due to the thermal noise is larger than the others, the controller adds the largest inverse shift −2λdt on dq 1 so as to cancel out the error. The other elements of u * have similar meanings. It should be maintained that the optimal controller obtained in the LQG control setup has this kind of natural structure found in the usual QEC scheme. (ii) The second remark is related to the nogo theorem [47] mentioned in Sec. I. The QEC problem corresponding to our setting is that the system's initial state is to be protected by feedback control. However, this goal is never accomplished, because, in general, a stable Kalman filter forgets the initial values of π(x 0 ) and V c (0); that is, the initial Gaussian state with mean π(x 0 ) and covariance V c (0) cannot be protected. This is a way to understand the no-go theorem for Gaussian QEC from the dynamical control viewpoint.

V. SIMULATION
In this section, a numerical simulation is provided to demonstrate effectiveness of the proposed control scheme. As an example of the memory we here take an optomechanical oscillator shown in FIG. 1. The oscillator must be subjected to thermal noise even at ultra-low temperature; in this case the noise strength n corresponds to the averaged photon number n = (e ωm/kB T − 1) −1 , where k B , T , and ω m denote the Boltzmann constant, the temperature, and the career frequency of the thermal channel, respectively.
The parameters are taken as follows. The oscillator couples with the input field with strength ν/2π = 30 kHz, while for the thermal channel it is Γ/2π = 1 Hz. The thermal channel is with frequecy ω m /2π = 10 MHz and with temperature T = 4 K, which leads to n = 8.8 × 10 3 . The mean value of the input state is taken as α in = −230 Hz 1/2 , which leads to − √ νβ = [100, 0, 100, 0, 100, 0]. Note that |α in | 2 is the mean photon number per unit time.

A. Coherent source state
First let us consider the case where the source state is taken as a coherent state |α in . This means that the input fluctuation is set to N 1 = M 1 = 0, implying that, in addition to the syndrome operators,p 1 +p 2 +p 3 can be estimated. Hence, we can now apply the feedback control scheme based on π(ŝ 1 ). Figure 3 shows how much the encoding and the feedback control improve the fidelity between the three-mode input state and the oscillators' steady state. They are quantified as follows. (i) Encoding strength is directly evaluated by the squeezing parameter µ; actually, as indicated by Eq. (4), making −µ bigger means that the input state is going to approximate the GHZ-like state (5). (ii) The LQG optimal control input u * can take possibly a bigger value by tuning the penalty parameter r smaller (Recall that we have set R = rI 6 in Eq. (21).) Hence, −log 2 r is interpreted as the control strength. (iii) The fidelity is given by the same form as Eq. (9) with V ∞ now replaced by the covariance matrix of the ensemble average of the controlled system state, say V ′ ∞ ; that is, F = 1/ det(V ′ ∞ + V in ). Note that the controlled system variablex t stochastically changes in time according to the QSDE (7) with the input u * a function of the estimate π(ŝ 1 ), i.e., Eq. (23). Hence V ′ ∞ corresponds to the first 6 × 6 block matrix of the 9 × 9 covariance matrix of (x ⊤ , π(ŝ 1 ) ⊤ ) ⊤ that satisfies the controlled Lyapunov equation (A1); see Appendix for the detailed calculation.
Let us now discuss the performance of the LQG optimal feedback control. It is observed from Fig. 3 that bigger control strength indicates bigger fidelity for all value of µ. This means that the feedback control can always reduce the excess fluctuation brought from the thermal environment, without respect to how much the source state is encoded. Regarding the encoding strength, however, care should be taken in its choice. Actually, while larger squeezing enables us to detect a bigger error signal which can induce more efficient feedback control, at the same time the memory state becomes more sensitive to the thermal noise. In other words, too much squeezing makes the input state fragile to the noise; consequently there exists an optimal value of the squeezing parameter, and it is about µ * = −0.4. Note this value can be reached within the current technology. In this case, the amount of improvement of the fidelity via the feedback control compared to the case without encoding and control is about 0.05; but this is not a big improvement. One reason to this limitation is that the error taken in this paper is the worst one in the sense that each system is subjected to the thermal noise for all time; this kind of error is not fully tractable via the standard QEC protocol. Therefore, it is expected that the proposed control scheme could show possibly much better performance of the state transfer against some weak noise such as simple displacement acting only one of the oscillators. and (b) π(q1). For both plots, the less (red) and larger (blue) fluctuating lines correspond to the cases where the optimal feedback control is performed or not, respectively.
Next, Fig. 4 shows the time evolutions of (a) the syndrome π(q 2 −q 3 ) and (b) the first oscillator's position π(q 1 ). For both plots, the less (red) and larger (blue) fluctuating lines correspond to the cases where the optimal feedback control is performed or not, respectively. The error bars show the standard deviation of the estimation error, which is calculated from the solution to the Riccati equation (19). The parameters are taken as µ = −0.4 and r = 10 −9 , which were shown in Fig. 3 to be the optimal values attaining the maximum fidelity.
Note that, because α in is unknown,q 1 cannot be exactly estimated; but we here plot the exact value of π(q 1 ) just for demonstration. Figure 4 (a) demonstrates that the fluctuation of the syndrome is certainly suppressed by the feedback control; this is the result that should be expected, since the controller is designed to minimize the syndrome errors. On the other hand, it is not straightforwardly expected that controllingq 1 would actually work well, but Fig. 4 (b) demonstrates that the feedback control based on the syndrome estimation is fairly effective for correcting the error acting on each position of the oscillator.

B. Squeezed source state
Here we are concerned with the situation where the source state is squeezed. The squeezing parameter is taken as µ 1 = log(2M 1 + 2N 1 + 1) with M 1 real, and N 1 and M 1 are both unknown in addition to α in . We then have to use the filter forŝ 2 , which does not contain N 1 and M 1 . But for comparison we further consider a squeezed source state with known covariance and unknown mean value; in this case the Kalman filter forŝ 1 can be used for feedback control. This comparison will reveal how much the additional information (i.e., the estimate ofp 1 +p 2 +p 3 ) improves the control performance at the expense of limiting the class of input states. Figure 5 shows the fidelity between the input state and the controlled oscillators' steady state, versus the squeezing parameters µ and µ 1 . The control penalty is r = 10 −9 . In the figure the upper and lower surfaces correspond to the filter forŝ 1 andŝ 2 , respectively. First, a notable fact observed from the lower surface is that, when aiming to transfer the input state with the completely unknown source state, the maximum fidelity is attained when µ 1 = 0; that is, squeezing the source state always decreases the fidelity. This is consistent with the standard understanding that an unknown squeezed state is in general fragile to the thermal noise, because it randomly rotates the phase of the state, which as a result brings that state into a mixed state. On the other hand, the upper surface in the figure shows that, when transferring the input state with known fluctuation, which means thatŝ 1 can be used for estimation and control, squeezing the source state improves the fidelity. Note that the improvement is observed if µ 1 > 0 and µ < 0; that is, the source state is now momentum-squeezed while the ancilla states are position-squeezed. Actually in this case the input state more closely approximates the GHZ-like state, which is an eigenstate ofp 1 +p 2 +p 3 . Therefore, the additional information π(p 1 +p 2 +p 3 ) should certainly improves the control performance of the state transfer. However, as in the previous case, too much squeezing degrades the fidelity, hence µ 1 must be optimized.

VI. CONCLUSION
In this paper we proposed a new feedback control scheme for state transfer, which has the form of QEC.
In particular, due to the Gaussianity of the input state and the linearity of the memory dynamics, the celebrated Kalman filtering and LQG feedback control were employed. We have considered a specific system, but the proposed control strategy is applicable to any linear quantum system, as long as the output process can be properly defined. Although it was shown in the example that the controller is not a very effective one, it is expected that it can show much bigger improvement for some linear systems having relatively weak noise. Also possible combination with a feedforward control scheme [21,22] could be helpful to attain better performance for Gaussian state transfer.

(A1)
Now, the covariance matrix of the controlled memory state, V ′ , corresponds to the first 6 × 6 block matrix of V z . In this case, the steady solution V ′ ∞ can be explicitly obtained as follows: Note that as r → 0 (i.e., cheap control), we have V ′ ∞ → V c (∞), which is the steady solution to Eq. (19). That is, in this limit, the quantum fluctuation of the controlled memory state is maximally reduced down to the level of estimation error.


extfieldsofstudy:
Computer Science
Mathematics

text:
Information theory of massively parallel probe storage channels

Motivated by the concept of probe storage, we study the problem of information retrieval using a large array of N nano-mechanical probes, N ~ 4000. At the nanometer scale it is impossible to avoid errors in the positioning of the array, thus all signals retrieved by the probes of the array at a given sampling moment are affected by the same amount of random position jitter. Therefore a massively parallel probe storage device is an example of a noisy communication channel with long range correlations between channel outputs due to the global positioning errors. We find that these correlations have a profound effect on the channel's properties. For example, it turns out that the channel's information capacity does approach 1 bit per probe in the limit of high signal-to-noise ratio, but the rate of the approach is only polynomial in the channel noise strength. Moreover, any error correction code with block size N>>1 such that codewords correspond to the instantaneous outputs of the all probes in the array exhibits an error floor independently of the code rate. We illustrate this phenomenon explicitly using Reed-Solomon codes the performance of which is easy to simulate numerically. We also discuss capacity-achieving error correction codes for the global jitter channel and their complexity.

I. INTRODUCTION
T HE invention of atomic force microscopy in 1986 [1] opened the possibility of storing information at nanometer scales resulting in proposals for achieving aerial information densities of tens or even hundreds of Terabits per square inch. The basic idea is that information can be stored by altering certain features of the storage medium at the scale of nanometers. These changes can be then sensed and information retrieved by nanoscale probes similar to those used in atomic force microscopy. For example, binary information can be stored in crystalline dots created in amorphous media or amorphous dots in crystalline media and retrieved using electric probes (phase change storage, [2]). Alternatively, the information can be stored using indentations or, more recently, variable length grooves made in polymer media and retrieved using thermoelectric probes (thermo-mechanical probe storage, [3] 1. The layout of the 'Millipede' -a thermomechanical probe storage device. Courtesy of IBM -Zurich, [7] The whole concept of storing and retrieving information using nano-scale probes became known as probe storage, see [4], Chapter 4, for a comprehensive review of the current state of the field. Perhaps the most widely known concept of the probe storage device is the IBM's 'Millipede' [8], [5] -an array of thermo-electric probes with sharp tips used to create and sense indentations in the polymer media. The layout of the Millipede is shown in Fig. 1, the basic principle of thermomechanical reading and writing is explained in Fig. 2. It has been demonstrated that data can be retrieved with low bit error rate per probe (about 10 −4 ) at densities of up to 2 Terabit per square inch [6].
The particular feature of the Millipede shared by all existing concepts of probe storage is the presence of a large (∼ 2 12 ) number of probes reading and writing the information in parallel. This feature makes probe storage very different from more traditional storage devices such as magnetic or optical disks or flash memory. Each probe reads and writes information in its own field. The array of probes moves as a whole to allow each probe to explore its field. For aerial densities of information of several terabit per square inch, the array has to be moved by a distance of the order of 10 nanometers from one set of the sampling points to the next and repositioned with a sub-nanometer precision for writing or reading. Inevitably, positioning errors affect every single probe in the array. We refer to the combined effect of errors in the positioning of the array at both the reading and writing stages as global positioning jitter or just global jitter. The aim of this paper is to investigate the performance of error correcting codes operating on the output of the probe array subject to highly correlated Fig. 2. The principle of thermomechanical reading and writing. Top left picture, writing bit 'zero': cold probe pressed against the polymer surface leaves no mark in the media. Top right picture, writing bit 'one': the probe heated above the polymer's melting temperature and pressed against the polymer surface leaves an indentation in the media shown in the bottom left picture. Bottom right pictures, reading bits 'one' and 'zero': the probe inserted in the indentation is cooler than the probe pressed against the flat surface. These temperature variations can be captured using thermo-resistive sensors. Courtesy of IBM Zurich, [8] disturbances due to global jitter. We will investigate both the performance of special codes (Reed-Solomon) and address the general question of existence of good error correction codes for the global jitter channel by calculating channel capacity and studying Gallager's random coding bound.
Despite the fairly abstract tools used in the paper, the purpose of our investigation is to answer a very practical question: given a 'good' communication channel (e. g. the single probe's channel with bit error rate 10 −4 ), how good is a communication system (e.g. a Millipede) consisting of thousands of good channels subjected to correlated error events?
Our system-level analysis allows one to understand advantages and limitations of complicated communications devices without actually building one. In particular, it turns out that information theoretical performance limits derived in this paper have crucial system-level implications for the design of error correction codes for probe storage.
The rest of the paper is organised as follows. In Section II we introduce a simple model for the probe storage channel, the global jitter channel, which accounts for the effects of the probe array's position jitter and calculate the probability distribution function of signal amplitude for Gaussian isolated pulse response and Gaussian statistics of jitter. In Section III we introduce a low complexity signal detection scheme for the global jitter channel and prove its optimality in the limit of large array sizes. In Section IV we calculate block error rate for non-interleaved Reed-Solomon codes applied to the probe storage channel and analyze their error floor behaviour. In Section V we calculate Shannon's capacity of the global jitter channel and show that it approaches 1 bit according to a power law in the limit of large signal-to-noise ratio. In Section VI we calculate the average block error rate for non-interleaved codes sampled from Gallager's ensemble (the random coding bound or RCB) and show that it exhibits an error floor behaviour as a function of SNR. In Section VII we show that there exists no non-interleaved codes with a positive rate which can be used for error-free retrieval of information using large ar rays of probes in the presence of global jitter. In Section VIII we discuss minimal requirements for capacity-achieving codes for global jitter channels. We conclude our work with Section IX which contains the summary of the results of our investigation. For completeness, we present the derivations of more technical results obtained in the paper in the Appendix.

II. CHANNEL MODEL
We consider the system consisting of an array of N probes reading/writing in parallel. We assume that channel coding has been used (an RLL code for example) and that the symbol pitch is large enough so that inter-symbol interference can be ignored. The sampled readback signal at the k-th probe at the t-th moment in time is modelled as: where a (k) t ∈ {0, 1} is the bit written to the medium by the k-th probe at time t, J t is the global positioning error (jitter) at time t, p(J) is the channel impulse response and {σn (t) k } is a sequence of random variables modelling the combined effect of electronics and media noise.
Experiments carried out at IBM Zurich [9] have confirmed that for thermo-mechanical storage media the combined electronics/media noise is well modelled by Gaussian random variables. Thus we assume that {n (k) t } are independent identically distributed Gaussian variables with mean zero and unit variance. Parameter σ is the standard deviation of the resulting additive white Gaussian noise (AWGN).
It has also been demonstrated experimentally that when using a Millipede-like positioning control loop the random jitter J t is also well modeled by a mean-zero Gaussian random variable [9]. Note that this observation applies to any probe storage device, for instance a system based on phase change media, as long as a similar positioning system is used. Let σ J be the standard deviation of Gaussian jitter.
Most of the qualitative results reported below do not depend on the detailed assumptions about the statistics of jitter or the precise shape of the impulse response. For quantitative analysis and numerical simulations of the global jitter channel we will use the Gaussian impulse response: where W is a parameter related to pulse width. In what follows it will be more convenient to work directly with the random variable p(J) ∈ (0, 1) which measures signal amplitude degradation due to position jitter J. A calculation leads to the following answer for the probability density function of p: Finally we define the signal-to-noise ratio (SNR) of the global jitter channel measured in decibels as follows:

III. OPTIMAL CHANNEL DETECTOR
The global jitter channel is characterised by strong equal time correlations between outputs of all N -probes in the array. In the limit N → ∞ it is possible to exploit these correlations and derive an asymptotically optimal low-complexity detection scheme.
To motivate the rigorous argument given below, let us ask the following question: what is the optimal channel detector conditional on the knowledge of the value of jitter at time t? Conditional on the known amplitude value p t ≡ p(J t ) = p, channel outputs r (k) t 's are independent. Consequently, the optimal maximum a posteriori (MAP) detector [11] is simply the collection of N independent optimal threshold detectors for AWGN channel: The problem with the above 'Genie-assisted' detector is that the value of the amplitude p is a priori unknown. However the value of p can be reliably estimated from the string r if N is sufficiently large. The key observation concerns the sample average of all N received signals at time t: Assume that channel inputs a (k) t are independent and uniformly distributed. Then by the strong law of large numbers [10], both sample averages in (6) converge almost surely to their expected values as N → ∞: Substituting into (6) we find in the large N limit it is actually possible to compute the previously unknown distortion p t = p(J t ) as follows: Feeding this estimate of p to N independent threshold detectors (5) we obtain estimatesâ 1 ,â 2 , . . . ,â N of data bits which should become optimal in the limit N → ∞.
The heuristic argument given above leads to the following asymptotically optimal detection algorithm for the global jitter channel: 1) Input: Channel output r (1) t , r (2) t , . . . r (N ) t at time t. 2) Estimate the threshold: 3) Perform bit-by-bit detection: 4) Output: The estimate of channel inputâ We will now prove that detector (10) is indeed optimal in the limit N → ∞. Firstly, note that detector (5) is an optimal MAP detector which infers the most probable channel input conditional on r t , p t . Therefore its bit error rate is smaller or equal than bit error of any detector which infers channel input conditional on channel input r t only. So we can establish the asymptotic optimality of (10) by proving that its bit error rate approaches the bit error rate of (5) in the limit N → ∞: A long but straightforward calculation presented in Appendix A shows that Taking the large-N limit of both sides of the above inequality we arrive at (11). Thus the asymptotic optimality of (10) is established. We will refer to detector (10) as the LLN -detector and detector (5) as Genie − detector.
The performance of the LLN-detector for large but finite N is compared with the performance of Genie-detector in Fig. 3 for Gaussian pulse (2) and various values of jitter strength σ J . The size of the probe array is N = 1000, the pulse width is W = 0.5 and the simulation has been run for 10 5 consecutive readings of the entire probe array. As a reference the performance of the channel with no jitter which corresponds to the classical binary-input additive white Gaussian noise (AWGN) channel has been included. It is observed that as jitter strength is increased the BER performance degrades significantly but for all cases the performance of the LLN-detector is virtually indistinguishable from Genie-assisted detector. It is interesting to note that the gap between BER of Genie-detector and the LLN-detector increases with SNR, see Fig. 3. This observation is consistent with the bound (12): The right hand side of the bound grows if the block size N is kept fixed but the additive noise strength σ is reduced. The intuition behind the observed divergence is very simple: for smaller noise one needs a more precise threshold estimate to stay near the performance of the Genie-detector. The precision of the threshold estimate depends on the rate of convergence of the sample sum and scales as 1/ √ N in accordance with Central Limit Theorem (see [10] for a review). It is interesting to compare the LLN detector to the simplest detection scheme proposed for Millipede, see. [12]. The Millipede detector consists of N independent MAP detectors (one detector per probe). The k-th detector estimates the data bit a The maximum-likelihood detector of the single probe output modeled by (1) is a simple threshold detector: the most likely bit given the output of the k-th probe,â (k) t is given by: where r 0 is the optimal threshold -a number between 0 and 1 which can be either measured experimentally or computed theoretically for a given the channel model by solving the maximum likelihood equation see [13] for more details.
In Figure 4 the BER of the Millipede detection scheme is compared with that of the LLN-detector for two values of jitter strength σ J . The result is quite striking: the optimal detector for the global jitter channel outperforms the set of N independent threshold detector by over a decibel for BER∼ 10 −4 .
Finally, let us analyse the complexity of the LLN-detector. The complexity of adding N fixed precision numbers in (9) scales as N log(N ), the complexity of the detection step (10) is O(N ) so the overall detection complexity is O(N log(N )), the detection complexity per detected bit is O(log(N )).
Having constructed the optimal detection scheme we can investigate the performance of error correction codes for the global jitter channel starting with Reed-Solomon codes which featured in the original Millipede proposal. The typical probe storage array size N is 64 × 64 = 4096, which is close to the sector size of the previous generation of hard disk drives. 1 The simplest error correction coding (ECC) scheme for probe storage follows the example of hard drives: for the latter data is encoded sector-by-sector, for the former ECC is applied independently to K-bit strings of data I 1 , I 2 , I 3 . . . to be recorded on the media by the probe array at times t 1 , t 2 , t 3 , . . .. We refer to the described application of error correction coding as non-interleaved meaning that channel outputs corresponding to different moments of sampling time cannot belong to the same code block.
The non-interleaved ECC block size is equal to the number of probes N , the ECC rate is R = K/N . During the reading stage, the single-time outputâ of the channel detector is fed into the ECC decoder resulting in an estimatê I t of K recorded bits at each sampling time t.
We start our study of error correction for global jitter channel with classical Reed-Solomon (RS) codes, see [14] for review. These are the (N s , K s ) symbol block codes over the Galois field GF (2 n ). Here N s = N/n is the number of RS symbols per block, K s = K/n is the number of information symbols. Symbols are represented by n-bit binary strings. The maximal block size of the code is N max = 2 n − 1 symbols. An RS code with any block size N s ≤ N max can be constructed by treating missing (N max − N s ) symbols as zeros (this operation is called shortening). A Reed-Solomon code with rate R = K s /N s can correct up to incorrectly detected symbols, which makes it a maximal distance separable (MDS) code.
The fact that the error event of a RS code depends on the number of incorrectly detected symbols rather than bits, makes it a very good code for channels dominated by relatively short bursts of noise. It is therefore not very well suited for the global jitter channel. However we will see below that conclusions drawn from analyzing the performance of RS codes in the presence of global jitter can be applied to any non-interleaved block ECC.
We start with presenting results of numerical simulations. Fig. 5 illustrates the performance of rate-0.8 RS codes with symbol size n = 10 bits. Three codes are considered: with block size N s = 1023 (the maximal block size) and two shortened RS codes -with N s = 511 and N s = 255. These block sizes correspond to hypothetical probe arrays with 2550, 5110 and 10230 probes. To model global jitter we use the channel (1) with jitter strength σ j = 0.2, and pulse width W = 0.5. For each of the codes we measure the probability of sector error or sector error rate (SER) as a function of signal-to-noise ratio (SNR). SER is measured as the number of times the number of symbol errors in the received string exceeded RS threshold N s 1−R 2 . The total number of strings is 10 6 to ensure that the total error count used to estimate SER is at least 10 2 for the highest SNR point.
The first striking feature of the curves presented in Fig. 5 is that there no discernible performance loss using the shortened codes. This is in stark contrast from the known behaviour of RS codes for the AWGN channel where Therefore, the numerical evidence suggests that the probability of RS codes applied to the global jitter channel does not depend on the code's block size.
The second unusual feature of SER vs. SNR curves shown in Fig. 5 is the exponential law of the decay of SER with SNR -all the curves look like straight lines on the semilogarithmic plot. The exponential rather than the 'waterfall' shape of SER curves is often referred to as an 'error floor'. Its appearance is normally attributed to the use of suboptimal decoding algorithms such as belief propagation, rather than channel properties. Here we are driven to a conclusion that non-interleaved RS codes decoded with an optimal maximum likelihood hard input algorithm exhibit an error floor, the position of which is independent of the block size.
As it turns out, the breakdown of Reed-Solomon codes for channels with global jitter can be understood analytically using the machinery of large deviations, [16].
To achieve this, we need to introduce some notations. Let be the maximal fraction of correctable symbols for the code at hand. Let the random variable ξ k take values {0, 1} with probabilities {e 0 , e 1 } respectively. Here e 1 is the probability that the k-th symbol is detected incorrectly. Due to statistical homogeneity of the channel model (1) the probability of symbol error does not depend on the symbol index k. The event ξ k = 1 corresponds to the k-th symbol being detected erroneously and ξ k = 0 the k-th symbol being detected correctly. Conditionally on the event p t = p the probability of RS sector error is simply given by a multinomial formula [15]. Unfortunately, the multinomial formula is not very useful for quantitative analysis in the region of low SER. Instead, we are going to use simple asymptotic expressions for RS sector error rate based on Cramer's theory [16]. The event of sector error corresponds to the fraction of incorrectly detected symbols exceeding RS threshold (14). Hence, the probability of a block being decoded incorrectly can be written: An application of Chernoff's bound [10] results in the following upper bound on Pr(SE) for any λ > 0: where E[• | p t = p] stands for p t -conditional expectation value. Recall that conditional on the value of signal amplitude p t the global jitter channel is memoryless. Therefore symbol error events ξ k are conditionally independent and identically distributed and the bound in (16) can be re-written as: The bound in (17) holds for any λ > 0 and thus we can choose the tightest bound possible by minimising over all λ's: Function I(τ, p) is known in the theory of large deviations as the rate function and is given by: The question remains: is the bound in (18) tight enough to be useful? The answer is provided by an application of Cramer's Theorem [16] which states that provided E [ξ | p t = p] < τ (which is certainly true in the limit of low sector error-rates) the bound (18) is tight in the limit of large block size: An explicit expression for the rate function (19) can be found by solving the critical point equation: The function differentiated in (21) is convex with the unique point of global minimum given by: where e 1 (p) = E(ξ | p t = p) is the conditional symbol error rate, e 0 (p) = 1 − e 1 (p). If the condition of Cramer's theorem E [ξ | p t = p] = e 1 (p) < τ is satisfied then λ c > 0 and substituting back into (19) we find the rate function can be expressed as a Kullback-Leibler divergence [17]: Recall that for any two stochastic vectors P and Q, On the other hand if e 1 (p) ≥ τ then λ c ≤ 0 and due to convexity the minimum of (19) is achieved at λ = 0 resulting in the trivial bound I(τ, p) = 0. Let us summarise our findings so far: where Now we can derive an upper-bound for the performance of non-interleaved RS codes in the global jitter channel using the upper bound (25, 26) on the conditional SER: where dµ P (p) = ρ P (p)dp is the probability measure of jitterdependent signal amplitude with density (3). The range of integration can now be split around the critical value of signal degradation p c which is the unique solution of the equation Then using the upper bound (25) for the conditional probability of RS error we find: Where F τ (p) is expressed in terms of Kullback-Leibler divergence: Note that F τ (p) has a unique non-degenerate minimum at p c where it takes the value F (p c ) = 0. The proof of this fact follows easily from Gibbs inequality [17]. Therefore, we can apply the Laplace formula [18] to the second integral in equation (28) to derive the large-N s asymptotic of the probability of sector error: Note that the resulting expression (29) tends to zero as 1 Ns in the limit N s → ∞. Therefore what we have discovered is that for the global jitter channel in the limit N s → ∞ the probability of sector error is upper-bounded as follows: We therefore confirmed theoretically that for a fixed level of noise and in the limit of large sector size, the probability of sector error for information encoded with a Reed-Solomon code and transmitted over the global jitter channel does not depend on the code's block size. This conclusion is in perfect Using (30) we can also explain the shape of the SER vs SNR curves thus confirming the appearance of the error floor analytically. Moreover, we will be able to determine the position of the error floor as a function of code rate R = 1−2τ .
Up to this point our considerations did not depend on the specific shape of the impulse response or on jitter statistics. From now on we will assume the Gaussian impulse response given by equation (2) and Gaussian position jitter. Then it follows from (30) that where J c = W ln(p −1 c ) is the critical value of jitter that causes signal degradation p c . Bounding the integral in equation (31) with elementary functions we arrive at the following: If we further assume that a high-rate Reed-Solomon is used so that τ ≪ 1 then it follows that the symbol error rate conditional on p c , e 1 (p c ) = τ is also much less that 1. In this limit we can approximate symbol error rate as follows: where f (p c ) is the bit error rate conditional on p c which admits the following upper bound: Note that this upper bound is tight in the limit of large signal to noise ratios -the region relevant for studying the error floor. Therefore we can estimate the critical point p c byp c ≥ p c as follows:p Substituting (35) into equation (32) we arrive as the following bound on sector error rate for the global jitter channel with the exponential impulse response for high-rate Reed-Solomon codes: This is a disastrous result: the probability of sector error decays as power law in σ with an exponent that does not depend on code rate or block size! The algebraic dependence of SER on the noise strength implies the exponential dependence of the probability of sector error on signal-to-noise ratio which explains straight lines on the semi-logarithmic SER-SNR plot in Fig. 5. Moreover, we observe that the position of the error floor (determined by the pre-factor in (36)) depends very weakly (logarithmically) on the code rate R = 1 − 2τ .
In Figure 6 we compare numerical simulations of RS codes of various rates with expression (31) which is valid beyond the high rate approximation used to derive (36). We find a good agreement with our theoretical prediction: the rate of the exponential decay does not depend on either the block size or the code rate, the position of the error floor changes slowly with the code rate.
We conclude that non-interleaved Reed-Solomon codes applied to channel (1) exhibit an error floor: the probability of sector error decays exponentially as a function of SN R. Moreover, the rate of the exponential decay does not depend on either the code rate or the block size. The position of the error floor varies as a logarithm of the code rate. Non-interleaved Reed-Solomon codes are thus not suitable for use in a probe storage system that suffers from global positioning errors.
The above discussion suggests that the performance of Reed-Solomon codes can be improved with interleaving: by spreading the codewords over multiple time samples, the effect of occasional strong jitter leading to p < p c (τ ) can be mitigated. We will discuss the complexity of this solution in Section VIII.
But first, we will address the following foundational question: do good codes for a global jitter channel exist in principle?
V. SHANNON CAPACITY OF THE GLOBAL JITTER CHANNEL.
Shannon's capacity C is one of the most important measures of quality for any communication channel. According to Claude Shannon's 1948 Channel Coding Theorem [19], there exists no error correction code of rate R which would achieve an arbitrarily small probability of error for a channel with capacity smaller than code rate, C < R. The 'positive' part of Shannon's theorem states that for any R < C and for any ǫ > 0 then there exists a block code C with block size M (ǫ) and rate less than or equal to R and a decoding algorithm such that the maximal probability of block error is less than ǫ. Thus by computing Shannon's capacity for the global jitter channel we will establish an upper bound on the rate of good error correction codes for this channel.
The definition of capacity rests on the notion of mutual information, see [17] for a review of fundamental notions of information theory. The mutual information between a discrete channel input a taking values in Ω a and a continuous output ensemble r taking values in Ω r is defined as follows: is Radon-Nykodim derivative of the conditional probability measure Pr R|A with respect to the marginal probability measure Pr R . If the probability densities ρ R|A (r | a) and ρ R (r) of the probability measures Pr R|A and Pr R exist, the Radon-Nykodim derivative is simply the ratio of densities: Recall that for the global jitter channel, a is an N -bit data string and r is a string of N real valued signals generated by the probe array at a given sample time.
Shannon capacity is defined as the maximal mutual information over all input probability distributions per bit of input: As it is easy to see from the definition, 0 ≤ C ≤ 1 (bit). In particular, C = 1 for noiseless channels. No information can be communicated over the channel with C = 0 in finite time.
For a complicated channel where the maximisation over Pr A is difficult to perform it is common to study a weaker form of capacity: where Pr i.u.d. is the probability distribution for which channel inputs are chosen independently and with equal probability. It is clear that Therefore if we find that C i.u.d. is close to 1, then C is also close to one and we can be certain that there is a high rate ECC scheme which will achieve low probability of error for the global jitter channel. From this point onwards we will be only concerned with C i.u.d. and assume that Pr A (a) = 1 2 N , i. e. that all 2 N inputs sequences a = (a 1 , . . . , a N ) are sampled independently according to a uniform distribution.
The calculation of channel capacity is a notoriously difficult problem. Analytically it can be evaluated for the very simplest channels only, such as binary symmetric, binary erasure or AWGN channels [17], Chapter II. There exist efficient numerical algorithms for calculating capacity for channels with rapidly decaying correlations such as ISI channels, [20]. Unfortunately, the numerical evaluation of capacity is not an option for channels with long range correlations such as the global jitter channel: due to the strong correlations between all signals received at the same time, the calculation of the N -dimensional integral in the right hand side of (37) cannot be reduced to a set of low dimensional problems which makes the numerical evaluation of capacity extremely inefficient.
Fortunately, it turns out that capacity of channel (1) can be calculated asymptotically in the limit of large array size N >> 1. The simplification which allows capacity calculation in the large-N limit is easy to understand: on the one hand, the received signals r are independent conditionally on the value of the amplitude p = p(J). Conditionally on p, channel capacity is given by the well known expression for binary AWGN channels. On the other hand, the value of p(J) can be extracted from the string of N >> 1 received signals with relative accuracy of the order of 1/ √ N due to the law of large numbers.
This simple argument suggests the following answer for capacity of the global jitter channel: where C AW GN (p) is the capacity of the AWGN channel with fixed signal amplitude p given by: where f = p/σ, and σ is the standard deviation of additive white noise, see [17], Chapter II, for more details. In what follows we will also use the asymptotic expansion of C AW GN valid in the limit of weak noise, f >> 1: As it turns out, the derivation of (40) is very simple and relies on some basic properties of mutual information which are both fundamental and intuitively obvious. Firstly, we notice that the mutual information for AWGN channel is maximized for the uniform distribution of inputs. Then a simple rearrangement of terms in (37) gives the following: I(R; A) − N E P (C AW GN (p)) = I(P ; R) − I(P ; (A, R)), (43) where I(P ; R) is the mutual information between signal amplitude degraded by jitter and the received signal, I(P ; (A, R)) is the mutual information between signal amplitude and the joint ensemble of channel input and output. Clearly, I(P ; (A, R)) ≥ I(P ; R).
(44) (The information we learn about P from observing A and R must be greater or equal to the information about P contained in R.) In case the above argument fails to convince a rigorousminded reader, here is the proof based on Jensen inequality [10]: Ωp×Ωr d Pr(P, R) = log 1 = 0.
Using (44), relation (43) leads to the following inequality: Intuitively, the validity of the above claim is fairly obvious: imagine for example that random variable P is represented by m-bit numbers. Then the mutual information I(P ; (A, R)) cannot exceed m bits and the limit in the right hand side (45) is trivially zero. For the proof of (45) in full generality, the reader is referred to Appendix B. The problem of computing capacity of the global jitter channel for large values of N is solved in principle: the right hand side of (40) is a finite-dimensional integral which depends on the probability distribution of p(J). In particular, it is well suited for numerical study. In Figure 7 the capacity of the probe storage channel suffering a global Gaussian jitter is shown for various values of jitter strength σ J . For σ J = 0 the channel is equivalent to the binary AWGN channel with signal amplitude equal to one. As seen from the plot, this channel has the highest capacity. For a fixed SNR, the capacity of global jitter decreases as σ J increases. For example at the SNR point corresponding to C AW GN = 0.9 bits, the capacity of global jitter channel with σ J = 0.3 is C i.u.d = 0.62 bits, which rules out the use of high rate linear error correction codes for this channel. The good news is that capacity seems to approach 1 bit per channel symbol in the limit of high SNR even in the presence of global jitter. According to Shannon's theorem this means that by reducing channel noise one can read and write information reliably using large parallel probe arrays at a small redundancy cost. Moreover, the global jitter channel can actually have a larger capacity than a collection of independent probes subject to individual jitter distortions of the same strength: In Figure 8 the capacity of the global jitter channel is compared against the reference channel where each probes suffers an independent Gaussian jitter distortion of the same strength for which the i.u.d capacity can be easily computed numerically. We observe that at low SNR the reference channel has the same capacity than the global jitter channel but at high SNR the global jitter channel actually has a larger capacity! We can explain this phenomenon by the fact that strong correlations between position jitters across the whole array can be used to extract extra information about the hidden parameters of the channel. This extra information can be used for example to build better signal detection algorithms, as discussed in Section III.
To complete our investigation of C i.u.d. we still need to: (i) confirm the numerical observation that capacity approaches 1 bit in the limit of low channel noise; (ii) estimate the corre-sponding rate of convergence. We will solve these problems for the exponential impulse response assuming weak additive noise and weak jitter, i. e. σ << 1, σ J << W. (46) While both of these conditions are reasonable in the context of applications, the weakness of jitter leads to a significant technical simplification of the argument given below. Our starting point is the following bound on the capacity of the AWGN channel: for some positive constants U and L. Eq. (47) results from the straightforward yet tedious analysis of (41). Averaging the above inequality over p we obtain the following bound for the capacity of the global jitter channel: where L ′ and U ′ are positive constants, Γ(x) is the Γ-function, As the right hand side of the lower bound (48) approaches 1 bit in the low noise limit σ = 0, we conclude that which confirms the results of numerical simulations. The upper bound (49) shows that the convergence to one is not faster than the power law σ W 2 2σ 2 J , which is much slower than the exponential convergence of C AW GN (p) ∼ 1 − Ae −p 2 /8σ 2 to one -only a very 'clean' global jitter channel will have capacity close to the capacity of AWGN channel.
Finally, comparing the lower capacity bound with the upper bound we see that with logarithmic precision the convergence of C i.u.d. for large N is indeed given by the power law: Therefore, the capacity approaches 1 bit in exactly the same way as the Reed-Solomon sector error rate vanishes in the limit of zero additive noise, see (36). The power law approach of the channel characteristics to their noiseless limiting values seems to be a feature of the global jitter channel. We will encounter the law (50) again in the next Section, when we compute the Gallager's coding bound for the global jitter channel.

VI. THE RANDOM CODING BOUND FOR THE GLOBAL
JITTER CHANNEL.
According to the results of the previous Section, the capacity of the global jitter is positive and even approaches 1 bit in the limit of low noise. Therefore, it follows from Shannon's theorem that there are error correction codes with rate up to capacity C which can be used to transmit information over the channel with vanishingly small probability of error. Yet as we have seen in Section IV, non-interleaved Reed-Solomon codes exhibit an error floor for any block size and any rate.
Therefore, we are led to the following question: are there any error correction block codes (linear or non-linear) which yield a small probability of block error if the code block coincides with the instantaneous output of the probe array?
Let C be a code with rate R and block size N (bits) which is also the number of probes in the array. Let Pr(SE | C) be the probability of block decoding error for this code.
This probability is difficult to compute for any non-trivial code. However, we can use the idea of Gallager [21] and calculate the average probability Pr(SE | R, N ) of block error under the maximum likelihood decoding, where the average is taken over all random codes with a given rate R and block size N . Given Pr(SE | R, N ) there must exist a code C 0 with rate R and block size N such that Therefore, if we find that Pr(SE | R, N ) is sufficiently low (e. g. 10 −12 ) for a desired code rate R, we will know that there are codes which can be used to correct errors for information communicated over the global jitter channel.
To perform the calculation of the random coding bound we need to define the space of random codes and a probability measure on this space. Following Gallager, we will construct a random binary code by picking K = N R binary codewords from the N -dimensional binary space {0, 1} N independently and uniformly (Gallager ensemble).
Using the sum rule, Pr(SE | C) can be re-written as follows: where Pr(SE | C, p) is the probability of block decoding error for the binary AWGN channel with a fixed signal amplitude p. Therefore, where Pr(Se | p, R, N ) is the random coding bound for the binary AWGN channel with a fixed signal amplitude p.
The following set of results can be easily extracted from the original Gallager's paper [21]: where E(R, p) is the error exponent for the binary AWGN channel. In what follows we do not need an explicit expression for the error exponent, but its basic properties listed below will be important to us: Using the notion of the error exponent for AWGN channel, we can re-write the random coding bound for the global jitter channel as follows: Let p c be the unique solution to Then (57) can be rewritten as: The crucial observation is that for p ≤ p c , C AW GN (p) ≤ R. As a result, for p ≤ p c , E(R, p) ≡ 0 due to (53). Therefore the average probability of decoding error for the global jitter channel is bounded: The second integral in the right hand side of (60) can be evaluated for N >> 1 using Laplace method: due to (54, 55) the main contribution to the integral comes from a small neighbourhood of p c . It follows from (55,56) that as a function of p, E(R, p) has a non-degenerate critical point at p = p c . A calculation based on the above points yields: It is remarkable that in the limit N → ∞ the probability of block decoding error is bounded by a function independent of block size: Therefore, the random coding bound for the global jitter channel does not vanish in the limit of the large block size! In Figure 9 the random coding bound (61) is shown for various code rates for a global Gaussian positioning error with strength σ J = 0.2. The random coding bound for a channel where each probe in an array of N = 1000 suffers an independent Gaussian positioning error of the same strength is also given for the same code rates. It is observed that the probability of decoding error (SER) is vastly worse for the case of global positioning errors. For the independent channel the SER decay with a waterfall (super-exponential) shape whereas for the global jitter channel the SER decays exponentially with SNR, i. e. exhibits an error floor.
The presence of the error floor in the random coding bound can be demonstrated analytically for Gaussian impulse response and Gaussian jitter in the limit of high SNR and high Where J c = W ln(1/p c ) is the positive value of jitter that results in signal loss p. Recall from (2) that W is parameter related to pulse width. Thus the average probability of decoding error is bounded by: Furthermore in the limit of weak noise σ → 0 and high rate R → 1, it is possible to use the asymptotic expansion for AWGN capacity C AW GN (p) given by (42) to derive an expression for p c to the logarithmic precision: Thus for high-rate codes the probability of decoding error can be approximately bounded: where C(σ J , w, R) is the σ-independent constant: We conclude the average probability of decoding error in the limit of weak noise and high code rate is asymptotically independent of the code block size, exhibits an error floor with an exponent independent of the code rate R and an amplitude with depends on the rate via log(1−R). An identical behaviour has been observed for Reed-Solomon codes in Section IV. Similar conclusions concerning the random coding bound can be reached for low code rates as well, see Fig. 10 where the bound (61) is shown as a function of rate R for a fixed SNR. Note that the average probability of block error approaches zero only in the limit of zero code rate. The fact that the random coding bound exhibits an error floor behavior identical to that of RS code suggests to us that all non-interleaved error correction codes suffer identical performance degradation due to global jitter. This suggestion is confirmed in the following Section.

VII. THE NON-EXISTENCE OF NON-INTERLEAVED ERROR
CORRECTION CODES FOR THE GLOBAL JITTER CHANNEL.
The existence of an N -independent error floor in the average probability of block error has the following simple explanation: provided the strength of additive noise is positive and R > 0, there is an N -independent critical value of jitter p c (R) beyond which the conditional channel capacity C AW GN (p) is smaller than the code rate R. Therefore, according to the negative part of Shannon's theorem, information transmission without errors is impossible with N -independent probability Pr(p < p c (R)).
The above consideration can be turned into a rigorous argument which shows that non-interleaved encoding of the large array's outputs cannot ensure error free information retrieval no matter what error correction code with positive rate is used: Let C(R, N ) be a block code with rate R and block size N . Clearly, Recall that p c is the unique solution to R = C AW GN (p). Notice that for every p in the region of integration, R > C AW GN (p). By Fano's inequality [22], Therefore, Therefore, the probability of block error is bounded below by a constant asymptotically independent of the block size. Integrating (69) by parts it is easy to show that where Expression (70) proves that the probability of block error is bounded away from zero by an N -independent constant for any non-interleaved error correction code with a positive rate R.
Therefore, error-free information storage using large arr ays of probes is impossible without interleaving error correction codes between multiple array outputs.
It is worth noting that there is a counterpart of Fano's inequality for the global jitter channel -a non-trivial fact given strong correlations between all probe channels within the array. Using the fact that lim N →∞ C i.u.d. = E p (C AW GN (p)) which we established in Section V we can derive the following weaker bound from the bound (70): which shows the impossibility of error-free information transmission over the global jitter channel using non-interleaved codes with R > lim N →∞ C i.u.d. . Finally, let us show that any non-interleaved code used to encode information transmitted over Gaussian global jitter exhibits an error floor identical to Reed-Solomon error floor (32, 36) or the error floor in the random coding bound (63, 65). Note that the capacity of AWGN channel C AW GN (p) is a function of p/σ only, see (41). Substituting expression (3) for the probability measure dµ P (p) for Gaussian jitter into the integral in the right hand side of (69) and changing the integration variable p = σx we arrive at the following result: where γ = W 2 2σ 2 J and x c is the unique positive solution to the equation Notice that x c σ = p c < 1. The following inequality is valid provided the additive noise is weak: If log 1 x c σ > 1/2 : and for any x : 0 < x ≤ x c , 1 Using this inequality in (72) we conclude the following: where We conclude that any non-interleaved error correction code with rate R and a large block size used to transmit information over the global jitter channel exhibits an exponential error floor: log Pr(SE) is a linear function of signal-to-noise ratio with exponent The position of the error floor is determined by the function I(R), which is a slow function of code rate, but we do not study it here. Note that the upper bound on the error floor we derived using elementary tools only (Fano's inequality) coincides with the error floor observed directly for Reed-Solomon codes and the random coding bound for the global jitter channel.

VIII. ON GOOD CODES FOR THE GLOBAL JITTER
CHANNEL.
In the two previous sections we established the impossibility of error-free information retrieval for large parallel arr ays of probes subject to global jitter using non-interleaved error correction codes of any positive code rate.
We also established in Section V that the capacity C i.u.d. of global jitter channel is positive and approaches 1 bit as the additive noise strength σ goes to 0. Therefore, by Shannon's theorem there must exist families of error correction codes with rates 0 < R < C i.u.d. which ensure vanishingly small probability of block error in the limit of large block sizes.
Can we say anything about the structure of these codes? We understand that non-interleaved codes cannot perform well on the global jitter channel due to rare strong jitter fluctuations leading to small effective AWGN amplitude p σ. Due to this fluctuations N bits of information get lost regardless of error correction code used.
To avoid this we must spread information over many time slices by encoding blocks of array outputs. Therefore good error correction codes for the global jitter channel must be interleaved.
Let us estimate the depth of interleaving and the block size of the corresponding codes. Let r t∈Z be the time series of array outputs. Recall that each output is an N -dimensional vector, where N is the number of probes in the array. For a feedbackloop based positioning system used in Millipede, jitter random variables J t are correlated in time. Let L be the correlation length measured in the number of sampling periods. Our previous discussion on the influence of correlations on the performance of error correction codes can be summed up as follows: the probability of decoding error of the block code of rate R < C i.u.d. and block size B approaches zero as the number of independent groups of samples B N L goes to infinity. In other words, if N L strongly correlated samples are treated as a single symbol, we get the usual statement of Shannon's theorem for memoryless channels: the probability of block error goes to zero as the number of symbols per block goes to infinity.
The minimal block size of high rate error correction codes used in modern storage devices which achieve the probability of error of the order of 10 −10 is of the order of 10 3 bits. Clearly, this is the lower bound on the interleaving depth. Therefore we can estimate the block size of the interleaved error correction code for probe storage as Thus for a probe storage device with N ∼ 10 3 probes in the array the error correction block size should be of the order B 10 6 bits .
The codes of similar sizes (5 · 10 5 bits) are already used in optical storage (Blu-Ray disks), however we thought it unusual to have discovered the need for large block sizes for a nonremovable media storage device.
IX. CONCLUSION In this paper we introduced and analyzed a simple model of the massively parallel probe storage channel suffering from global positioning errors. This model can be viewed as an example of communication channel consisting of N >> 1 parallel sub-channels all of which are strongly correlated via a common distortion event.
We solved the problem of optimal signal detection for the global jitter channel. Namely we found a detection algorithm of O(N logN ) complexity such that the estimated sequence converges bitwise to the maximum a posteriori estimate in the limit N → ∞. This is not an entirely trivial result, as in general, the complexity of optimal detection grows exponentially with memory length (for example, the optimal MAP detector for the channel with inter-symbol interference of length I and AWG noise has 2 I states). For the channel at hand the detection algorithm simplifies due to statistical independence of the channel outputs conditional on the value of jitter. This value can be found by analyzing all channel outputs using the law of large numbers.
We analyzed the performance of Reed-Solomon codes applied to the global jitter channel without interleaving channel inputs corresponding to different moments of time. We discovered both numerically and theoretically that in the limit of large probe array, any non-interleaved Reed-Solomon code will exhibit an error floor the position of which is independent of the code's block size N and only weakly dependent on the code rate. This is a surprising result, as the phenomenon of the error floor is usually associated with sub-optimality of the decoding algorithm (e. g. belief propagation) rather than specifics of the channel. For the case of global jitter the origin of the error floor can be traced back to rare instances of strong jitter J such that the signal amplitude p(J) becomes of the same order as channel noise σ.
Motivated by these findings we addressed the following question: are there any good error correction codes for prob e storage channels? The answer turned out to be two-fold: firstly, we calculated the capacity C i.u.d. of the probe storage channel and discovered that it does approach one bit per channel symbol in the limit of weak channel noise σ, albeit very slowly: 1 − C i.u.d. approaches zero as σ γ , where the exponent γ is equal to the error floor exponent for Reed-Solomon codes. However, we also found that Gallager's random coding bound for channel encoding which does not interleave between different time slices exhibits exactly the same error floor as Reed-Solomon codes! Moreover, an application of Fano inequality allowed us to prove that any non-interleaved error correction code applied to a massively parallel probe storage channel exhibits an error floor behaviour which is at least as bad as the Reed-Solomon error floor.
With the benefit of hindsight, the appearance of the universal error floor in Reed-Solomon block error rate, Shannon capacity, Gallager's random coding bound and Fano's low bound on sector error rate for an arbitrary fixed code can be related to the statistics of large jitter fluctuations: if p(J) << σ then regardless of the code used regions of confuse-ability (a hypersphere of radius σ) for any two codewords intersect and all information and parity is lost! In this sense the global jitter channel for weak channel noise can be viewed as an effective block-wise erasure channel: either N bits are detected with no errors, or all N bits are lost. For such a channel, it is clearly impossible to improve performance by increasing the number of probes in the array if the encoding block size is kept equal to N . In order to reach Shannon's limit for the global jitter channel it is necessary to spread information between many outputs of the probe array thus mitigating the effects of occasional strong jitter. We estimate the necessary block size of good error correction codes for probe storage channel to be of the order of 10 6 -an unusually large number for a non-removable storage device.
The mathematical results reported in the paper have severe practical implication for probe storage. Either the positioning system must be accurate enough so that global jitter is made very weak or we must interleave the error correction code in the time direction. The former option poses a huge engineering challenge: for indentation sizes of the order of 10 nm, our results suggest that the precision of the positioning system must be significantly better than 1 nm for the effects of global jitter to become insignificant. The latter option also presents significant implementation difficulties especially since jitter is known to be strongly correlated in the time direction and as a result we must interleave very deeply to overcome its effects.
Due to the universal nature of information theoretic performance limits we discovered, we see no easy 'engineering' way around the problem of global jitter. For example one can attempt to re-read the information recovered by the array of probes following a strong jitter event. This is similar to the way off-track errors are dealt with in magnetic hard drives. Our results mean however that the throughput of such a system will be severely degraded due to a large number of necessary re-reads.
From a more theoretical point of view, we developed a general approach for solving the problem of performance We now minimise with respect to ǫ to obtain the tightest bound possible. The critical value of ǫ is given by: Substituting ǫ min into (84) we arrive at the bound (12). APPENDIX B THE CALCULATION OF lim N →∞ 1 N I(P ; (A, R)). We will prove that lim N →∞ 1 N I(P ; (A, R)) = 0 not just for Gaussian jitter, but any distribution of p ∈ (0, 1) with finite differential entropy: By definition, the mutual information between signal strength P and the channel's input and output (A, R) is: Using the independence of signal strength P and channel input A the above expression can be re-written as follows: For Gaussian additive channel noise, the fraction under the sign of the logarithm takes the form F (P, A, R) ≡ d Pr(R | A, P ) d Pr(R | A) Conditionally on P , random variables are X and Y are equal to sums of independent identically distributed random variables and Therefore, conditionally on P , X and Y converge strongly to their respective expectation values. This motivates the introduction of new random variables strongly converging to zero: Using X N , Y N , the problem of calculating lim N →∞ 1 N I(P ; (A, R)) can be formulated as follows: compute where F (P, X, Y ) = Let us fix ǫ > 0. A calculation which essentially repeats the derivation of Chebyshev inequality shows that Now let us re-write L using the partition of unity as the sum of two terms: where N >ǫ logF (P, X N , Y N ), (99) Under the condition that the differential entropy (the logarithmic moment) of p(J) exists, it is easy to show that Similarly, where the last line is obtained using Schwarz inequality. Using (96) and the fact that E(2|X N | + |Y N |) 2 = 1 N 1 2 + p 2 + 2σ 2 2 we find that The expectation value in the above expression always exists as the random variable P is bounded. Combining (100) and (101) we find that for any ǫ > 0 and any N , for some ǫ, N -independent constant C. Choosing ǫ = N −1/2 and taking the limit N → ∞, we find that