Quantum learning of coherent states

We develop a quantum learning scheme for binary discrimination of coherent states of light. This is a problem of technological relevance for the reading of information stored in a digital memory. In our setting, a coherent light source is used to illuminate a memory cell and retrieve its encoded bit by determining the quantum state of the reflected signal. We consider a situation where the amplitude of the states produced by the source is not fully known, but instead this information is encoded in a large training set comprising many copies of the same coherent state. We show that an optimal global measurement, performed jointly over the signal and the training set, provides higher successful identification rates than any learning strategy based on first estimating the unknown amplitude by means of Gaussian measurements on the training set, followed by an adaptive discrimination procedure on the signal. By considering a simplified variant of the problem, we argue that this is the case even for non-Gaussian estimation measurements. Our results show that, even in absence of entanglement, collective quantum measurements yield an enhancement in the readout of classical information, which is particularly relevant in the operating regime of low-energy signals.


Introduction
Programmable processors are expected to automate information processing tasks, lessening human intervention by adapting their functioning according to some input program. This adjustment, that is, the process of extraction and assimilation of relevant information to perform efficiently some task, is often called learning, borrowing a word most naturally linked to living beings. Machine learning is a well-established and interdisciplinary research field, broadly fitting within the umbrella of Cybernetics, that seeks to endow machines with this sort of ability, rendering them able to "learn" from past experiences, perform pattern recognition and identification in scrambled data, and ultimately self-regulate [1,2]. Algorithms featuring learning capabilities have numerous practical applications, including speech and text recognition, image analysis, and data mining.
Whereas conventional machine learning theory implicitly assumes the training set to be made of classical data, a more recent variation, which can be referred to as quantum machine learning, focuses on the exploration and optimisation of training with fundamentally quantum objects. Quantum learning [3], as an area of strong foundational and technological interest, has recently raised great attention. Particularly, the use of programmable quantum processors has been investigated to address machine learning tasks such as pattern matching [4], binary classification [5,6,7,8], feedback-adaptive quantum measurements [9], learning of unitary transformations [10], 'probably approximately correct' learning [11], and unsupervised clustering [12]. Quantum learning algorithms provide not only performance improvements over some classical learning problems, but they naturally have also a wider range of applicability. Quantum learning has also strong links with quantum control theory [13], and is thus becoming an increasingly significant element of the theoretical and experimental quantum information processing toolbox.
In this paper, we investigate a quantum learning scheme for the task of discriminating between two coherent states. Coherent states stand out for their relevance in quantum optical communication theory [14,15,16], quantum information processing implementations with light, atomic ensembles, and interfaces thereof [17,18], and quantum optical process tomography [19]. Lasers are widely used in current telecommunication systems, and the transmission of information can be theoretically modelled in terms of bits encoded in the amplitude or phase modulation of a laser beam. The basic task of distinguishing two coherent states in an optimal way is thus of great importance, since lower chances of misidentification translate into higher transfer rates between the sender and the receiver.
The discrimination of coherent states has been considered, so far, within two main approaches, namely minimum-error and unambiguous discrimination, although the former is more developed. Generally, a logical bit can be encoded in two possible coherent states |α and |−α , via a phase shift, or in the states |0 and |2α , via amplitude modulation. Both encoding schemes are equivalent, since one can move from one to the other by applying Weyl's displacement operatorD(α) to both states 2 . In the minimum-error approach, the theoretical minimum for the probability of error is given by the Helstrom formula for discriminating two pure states [20]. A variety of implementations have been devised to achieve this task, e.g., the Kennedy receiver [21], based on photon counting; the Dolinar receiver [22], a modification of the Kennedy receiver with real-time quantum feedback; and the homodyne receiver 3 . Concerning the unambiguous approach to the discrimination problem, results include the unambiguous discrimination between two known coherent states [24,25], and its programmable version, i.e., when the information about the amplitude α enters the discrimination device in a quantum form [26,27,28].
The goal of this paper is to explore the fundamental task of discriminating between two coherent states with minimum error, when the available information about their amplitudes is incomplete. The simplest instance of such problem is a partial knowledge situation: the discrimination between the (known) vacuum state, |0 , and some coherent state, |α , where the value of α is not provided beforehand in the classical sense, but instead encoded in a number n of auxiliary modes in the state |α ⊗n . Such discrimination scheme can be cast as a learning protocol with two steps: a first training stage where the auxiliary modes (the training set) are measured to obtain an estimate of α, followed by a discrimination measurement based on this estimate. We then investigate whether this two-step learning procedure matches the performance of the most general quantum protocol, namely a global discrimination measurement that acts jointly over the auxiliary modes and the state to be identified.
Before proceeding with the derivation of our results and in order to motivate further the problem investigated in this paper, let us define the specifics of the setting in the context of a quantum-enhanced readout of classically-stored information.
Imagine a classical memory register modelled by an array of cells, where each cell contains a reflective medium with two possible reflectivities r 0 and r 1 . To read the information stored in the register, one shines light into one of the cells and analyses its reflection. The task essentially consists in discriminating the two possible states of the reflected signal, which depend on the reflectivity of the medium and thus encode the logical bit stored in the cell. In a seminal paper on quantum reading [29], the author takes advantage of ancillary modes to prepare an initial entangled state between those and the signal. The reflected signal is sent together with the ancillae to a detector, where a joint discrimination measurement is performed. A purely quantum resource-entanglement-is thus introduced, enhancing the probability of a successful identification of the encoded bit 4 . This model has been later extended to the use of error correcting codes, thus defining the notion of quantum reading capacity [30] also studied in the presence of various optical limitations [31]. The idea of using nonclassical light to improve the performance of classical information tasks can be traced back to precursory works on quantum illumination [32,33], where the presence of a low-reflectivity object in a bright thermal-noise bath is detected with higher accuracy when entangled light is sent to illuminate the target region. For more recent theoretical and experimental developments in optical quantum imaging, illumination and reading, including studies on the role of nonclassical correlations beyond entanglement, refer e.g. to Refs. [34,35,36,37,38,39,40,41,42,43,44].
In this paper we consider a reading scenario with an imperfect coherent light source and no initial entanglement involved. The proposed scheme is as follows (see Fig. 1). Wê 3 While the latter is the simplest procedure, it does not achieve optimality. However, for weak coherent states (|α| 2 0.4), it yields an error probability very close to the optimal value P e , and it is optimal among all Gaussian measurements [23]. In fact, among the three mentioned, only the Dolinar receiver is globally optimal. 4 In particular, in [29] a two-mode squeezed vacuum state is shown to outperform any classical light, in the regime of few photons and high reflectivity memories. Figure 1. A quantum reading scheme that uses a coherent signal |α , produced by a transmitter, to illuminate a cell of a register that stores a bit of information. A receiver extracts this bit by distinguishing between the two possible states of the reflected signal, |0 and |α , assisted by n auxiliary modes sent directly by the transmitter. model an ideal classical memory by a register made of cells that contain either a transparent medium (r 0 = 0) or a highly reflective one (r 1 = 1). A reader, comprised by a transmitter and a receiver, extracts the information of each cell. The transmitter is a source that produces coherent states of a certain amplitude α. The value of α is not known with certainty due, for instance, to imperfections in the source, but it can be statistically localised in a Gaussian distribution around some (known) α 0 . A signal state |α is sent towards a cell of the register and, if it contains the transparent medium, it goes through; if it hits the highly reflective medium, it is reflected back to the receiver in an unperturbed form. This means that we have two possibilities at the entrance of the receiver upon arrival of the signal: either nothing arrives-and we represent this situation as the vacuum state |0 -or the reflected signal bounces back-which we denote by the same signal state |α . To aid in the discrimination of the signal, we alleviate the effects of the uncertainty in α by considering that n auxiliary modes are produced by the transmitter in the global state |α ⊗n and sent directly to the receiver. The receiver then performs measurements over the signal and the auxiliary modes and outputs a binary result, corresponding with some probability to the bit stored in the irradiated cell.
We set ourselves to answer the following questions: (i) which is the optimal (unrestricted) measurement, in terms of the error probability, that the receiver can perform? and (ii) is a joint measurement, performed over the signal together with the auxiliary modes, necessary to achieve optimality? To accomplish the set task, we first obtain the optimal minimum-error probability considering collective measurements (Section 2). Then, we contrast the result with that of the standard estimate-and-discriminate (E&D) strategy, consisting in first estimating α by measuring the auxiliary modes, and then using the acquired information to determine the signal state by a discrimination measurement tuned to distinguish the vacuum state |0 from a coherent state with the estimated amplitude (Section 3). In order to compare the performance of the two strategies, we focus on the asymptotic limit of large n. The natural figure of merit is the excess risk, defined as the excess asymptotic average error per discrimination when α is perfectly known. We show that a collective measurement provides a lower excess risk than any Gaussian E&D strategy, and we conjecture (and provide strong evidence) that this is the case for all local strategies (Section 4). We conclude with a summary and discussion of our results (Section 5), while some technical derivations and proofs are deferred to Appendices.

Collective strategy
The global state that arrives at the receiver can be expressed as either [α] ⊗n ⊗ [0] or [α] ⊗n ⊗[α], where the shorthand notation [ · ] ≡ | · · | will be used throughout the paper. For simplicity, we take equal a priori probabilities of occurrence of each state. We will always consider the signal state to be that of the last mode, and all the previous modes will be the auxiliary ones. First of all, note that the information carried by the auxiliary modes can be conveniently "concentrated" into a single mode by means of a sequence of unbalanced beam splitters 5 . The action of a beam splitter over a pair of coherent states |α ⊗ |β yields where T is the transmissivity of the beam splitter, R is its reflectivity, and T + R = 1. A balanced beam splitter (T = R = 1/2) acting on the first two auxiliary modes thus returns |α ⊗ |α −→ | √ 2α ⊗ |0 . Since the beam splitter preserves the tensor product structure of the two modes, one can treat separately the first output mode and use it as input in a second beam splitter, together with the next auxiliary mode. By choosing appropriately the values of T and R, the transformation | √ 2α ⊗ |α −→ | √ 3α ⊗ |0 can be achieved. Applying this process sequentially over the n auxiliary modes, we perform the transformation Note that this is a deterministic process, and that no information is lost, for it is contained completely in the complex parameter α. This operation allows us to effectively deal with only two modes. The two possible global states entering the receiver hence become [ . The parameter α is not known with certainty. This lack of information can be embedded into the analysis by considering averaged global states over the possible values of α, where the choice of the prior probability distribution accounts for the prior knowledge that we might already have. One readily sees that a flat prior distribution for α, representing a limiting situation of complete ignorance, is not reasonable in this particular setting. On the one hand, such prior would yield divergent average states of infinite energy, since the phase space is infinite. On the other hand, in a real situation it is not reasonable at all to assume that all amplitudes α are equally probable 6 . The usual procedure in these cases is to consider that a small number of auxiliary modes is used to make a rough estimation of α, such that our prior becomes a Gaussian probability distribution centred at α 0 , whose width goes as ∼ 1/ √ n 7 . Under these considerations, we express the true amplitude α as where the parameter u follows the Gaussian distribution 5 See, e.g., Section III A in [45] for details. 6 Nonetheless, for finite dimensional systems, assuming a uniform prior distribution can be better justified and very useful; see, e.g., [46]. 7 Since we are interested in comparing the asymptotic performance of discrimination strategies in the limit of large n, the number of modes used for the rough estimation is negligible, i.e.,ñ = n 1− . Then, it can be shown that α belongs to a neighbourhood of size n −1/2+ centred at α 0 , with probability converging to one (this is shown, though in a classical statistical context, in [47]). Moreover, this happens to be true for any model of i.i.d. quantum states ρ (regardless their dimensionality), hence the analysis of the asymptotic behaviour of any estimation model of this sort can be restricted to a local Gaussian model, centred at a fixed state ρ 0 . This is known as local asymptotic normality [48,49,50,51].
To avoid divergences, we have introduced the free parameter µ as a temporal energy cutoff that defines the width of G(u). After obtaining expressions for the excess risks in the asymptotic regime of large n, we will remove the cut-off dependence by taking the limit µ → ∞.
Exploiting the prior information acquired through the rough estimation, that is using Eqs. (3) and (4), we compute the average global states arriving at the receiver The optimal measurement to determine the state of the signal is the Helstrom measurement for the discrimination of the states σ 1 and σ 2 [20], that yields the average minimum-error probability 8 where M 1 = tr √ M † M denotes the trace norm of the operator M. The technical difficulty in computing P opt e (n) resides in the fact that σ 1 − σ 2 is an infinite-dimensional full-rank matrix, hence its trace norm does not have a computable analytic expression for arbitrary finite n. Despite this, one can still resort to analytical methods in the asymptotic regime n → ∞ by treating the states perturbatively.
To ease this calculation, we first apply the displacement operator to the states σ 1 and σ 2 , whereD 1 (D 2 ) acts on the first (second) mode, and we obtain the displaced global states Since both states have been displaced by the same amount, the trace norm does not change, i.e., σ 0 − σ 1 1 = σ 0 −σ 1 1 . Eq. (9) directly yields where and {|k } is the Fock basis. Note that, as a result of the average, the first mode in Eq. (11) corresponds to a thermal state with average photon number µ 2 . Note also that the n-dependence is entirely inσ 2 . In the limit n → ∞, we can expand the second mode ofσ 2 as 8 Note that, sensu stricto, the dependence of P opt e (n) on the localisation parameter α 0 should be made explicit. Keep in mind that, in general, all quantities computed in this paper will depend on α 0 . Thus for the sake of notation clarity, we omit it hereafter when no confusion arises.
Then, up to order 1/n its asymptotic expansion gives Inserting Eq. (13) into Eq. (10) and computing the corresponding averages of each term in the expansion, we obtain a state of the form We can now use Eqs. (11) and (14) to compute the trace norm σ 1 −σ 2 1 in the asymptotic regime of large n, up to order 1/n, by applying perturbation theory. The explicit form of the terms in Eq. (14), as well as the details of the computation of the trace norm, are given in Appendix A. Here we just show the result: the average minimum-error probability P opt e (n), defined in Eq. (7), can be written in the asymptotic limit as P opt e ≡ P opt e (n → ∞) where Λ (2) ± is given by Eq. (A. 19).

Excess risk
The figure of merit that we use to assess the performance of our protocol is the excess risk, that we have defined as the difference between the asymptotic average error probability P opt e and the average error probability for the optimal strategy when α is perfectly known. As we said at the beginning of the section, the true value of α is α 0 + u/ √ n for a particular realisation, thus knowing u equates knowing α. The minimum-error probability for the discrimination between the known states |0 and α 0 + u/ √ n , P * e (u, n), averaged over the Gaussian distribution G(u), takes the form To compute this integral we do a series expansion of the overlap in the limit n → ∞ and integrate the resulting terms (see Appendix D). After some algebra we obtain P * e ≡ P * e (n → ∞) where The excess risk is then given by Eqs. (15) and (17) as Finally, we remove the cut-off imposed at the beginning by taking the limit µ → ∞ and we obtain Note that the excess risk only depends on the module of α 0 , i.e., on the average distance between |α and |0 . The excess risk is thus phase-invariant, as it should. Eq. (20) is the first piece of information we need to address the main question posed at the beginning, namely whether the optimal performance of the collective strategy is achievable by an estimate-and-discriminate (E&D) strategy. We now move on towards the second piece.

Estimate & Discriminate strategy
An alternative-and more restrictive-strategy to determine the state of the signal consists in the natural combination of two fundamental tasks: state estimation, and state discrimination of known states. In such an E&D strategy, all auxiliary modes are used to estimate the unknown amplitude α. Then, the obtained information is used to tune a discrimination measurement over the signal that distinguishes the vacuum state from a coherent state with the estimated amplitude. In this Section we find the optimal E&D strategy based on Gaussian measurements and compute its excess risk R E&D . Then, we compare the result with that of the optimal collective strategy R opt .
The most general Gaussian measurement that one can use to estimate the state of the auxiliary mode | √ nα is a generalised heterodyne measurement, represented by a positive operator-valued measure (POVM) with elements i.e., projectors onto pure Gaussian states with amplitudeβ and squeezing r along the direction φ. The outcome of such heterodyne measurementβ = √ nβ produces an estimate for √ nα, hence β stands for an estimate of α 9 . Upon obtainingβ, the prior information that we have about α gets updated according to Bayes' rule, so that now the signal state can be either [0] or some state ρ(β). The form of this second hypothesis is given by where p(α|β) encodes the posterior information that we have acquired via the heterodyne measurement. It represents the conditional probability of the state of the auxiliary mode being √ nα , given that we obtained the outcomeβ. Bayes' rule dictates where p(β|α) is given by (see Appendix B) p(α) is the prior information of α before the heterodyne measurement, and is the total probability of giving the estimate β.
The error probability of the E&D strategy, averaged over all possible estimates β, is then Note that the estimate β depends ultimately on the number n of auxiliary modes, hence the explicit dependence in the left-hand side of Eq. (26). We are interested in the asymptotic expression of Eq. (26), so let us now focus on the n → ∞ scenario. Recall that an initial rough estimation of α permits the localisation of the prior p(α) around a central point α 0 , such that α ≈ α 0 +u/ √ n, where u is distributed according to G(u), defined in Eq. (4). Consequently, the estimate β will also be localised around the same point, i.e., β ≈ α 0 + v/ √ n, v ∈ C. As a result, we can effectively shift from amplitudes α and β to a local Gaussian model around α 0 , parameterised by u and v. According to this new model, we make the following transformations: where, for simplicity, we have assumed α 0 to be real. Note that this can be done without loss of generality. Note also that, by the symmetry of the problem, this assumption implies φ = 0. The shifting to the local model transforms the trace norm in Eq. (26) as where To compute the explicit expression of ρ(v) we proceed as in the collective strategy. That is, we expand [u/ √ n] in the limit n → ∞ up to order 1/n, as in Eq. (13), and we compute the trace norm using perturbation theory (see Appendix C for details). The result allows us to express the asymptotic average error probability of the E&D strategy as where ∆ E&D is given by Eq. (C.5).

Excess risk
The excess risk associated to the E&D strategy is generally expressed as where P * e is the error probability for known α, given in Eq. (17), and P E&D e is the result from the previous section, i.e., Eq. (33). The full analytical expression for R E&D (r) is given in Eq. (C.6). Note that we have to take the limit µ → ∞ in the excess risk, as we did for the collective case. Note also that all the expressions calculated so far explicitly depend on the squeezing parameter r (apart from α 0 ). This parameter stands for the squeezing of the generalised heterodyne measurement in Eq. (21), which we have left unfixed on purpose. As a result, we now define, through the squeezing r, the optimal heterodyne measurement over the auxiliary mode to be that which yields the lowest excess risk (34), i.e., To find the optimal r, we look at the parameter estimation theory of Gaussian models (see, e.g., [51]). In a generic two-dimensional Gaussian shift model, the optimal measurement for the estimation of a parameter θ = (q, p) is a generalised heterodyne measurement 10 of the type (21). Such measurement yields a quadratic risk of the form where p(θ) is some probability distribution,θ is an estimator of θ, and G is a two-dimensional matrix. One can always switch to the coordinates system in which G is diagonal, G = diag(g q , g p ), to write It can be shown [51] that the optimal squeezing of the estimation measurement, i.e., that for which the quadratic risk Rˆθ is minimal, is given by We can then simply compare Eq. (37) with Eq. (34) to deduce the values of g q and g p for our case. By doing so, we obtain that the optimal squeezing reads where Eq. (39) tells us that the optimal squeezing r is a function of α 0 that takes negative values, and asymptotically approaches zero when α 0 is large (see Fig. 2). This means that the optimal estimation measurement over the auxiliary mode is comprised by projectors onto coherent states, antisqueezed along the line between α 0 and the origin (which represents the vacuum) in phase space. In other words, the estimation is tailored to have better resolution along that axis because of the subsequent discrimination of the signal state. This makes sense: since the error probability in the discrimination depends primarily on the distance between the hypotheses, it is more important to estimate this distance more accurately rather than along the orthogonal direction. For large amplitudes, the estimation converges to a (standard) heterodyne measurement with no squeezing. As α 0 approaches 0 the states of the signal become more and more indistinguishable, and the projectors of the heterodyne measurement approach infinitely squeezed coherent states, thus converging to a homodyne measurement. Inserting Eq. (39) into Eq. (35) we finally obtain the expression of R E&D as a function of α 0 , which we can now compare with the excess risk for the collective strategy R opt , given in Eq. (20). We plot both functions in Fig. 3. For small amplitudes, say in the range 0.3 α 0 1.5, there is a noticeable difference in the performance of the two strategies, reaching more than a factor two at some points. We also observe that the gap closes for large amplitudes α 0 → ∞; this behaviour is expected, since the problem becomes essentially classical when the energy of the signal is sufficiently large. Interestingly, very weak amplitudes α 0 → 0 also render the two strategies almost equivalent.

General estimation measurements
We have showed that a local strategy based on the estimation of the auxiliary state via a generalised heterodyne measurement, followed by the corresponding discrimination measurement on the signal mode, performs worse than the most general (collective) strategy. However, the considered E&D procedure does not encompass all local strategies. The heterodyne measurement, although with some nonzero squeezing, still detects the phase space around α 0 in a Gaussian way, i.e., up to second moments. In principle, one might expect that a more general measurement that produces a non-Gaussian probability distribution for the estimate β might perform better in terms of the excess risk, and even possibly match the optimal performance, closing the gap between the curves in Fig. 3. Here we show that the  observed difference in performance between the collective and the local strategy is not due to lack of generality of the latter. We do so by considering a simplified although nontrivial version of the problem that allows us to obtain a fully general solution.
The intuitive reason why one could think, at first, that a non-Gaussian probability distribution for β might give an advantage is the following. Imagine that α is further restricted to be on the positive real axis. Then, the true α is either to the left of α 0 or to the right, depending on the sign of the local parameter u. In the former case, α is closer to the vacuum, so the error in discriminating between them two is larger than for the states on the other side. One would then expect that an ideal strategy should better estimate the negative values of the parameter u, compared to the positive ones. Gaussian measurements like the heterodyne detection do not contemplate this situation, as they are translationally invariant, and that might be the reason behind the gap in Fig. 3.
To test this, we design the following simple example. Since the required methods are a straightforward extension of the ones used in the previous sections, we only sketch the procedure without showing any explicit calculation. Imagine now that the true value of α is not Gaussian distributed around α 0 , but it can only take the two values α = α 0 ± 1/ √ n, representing the states that are closer to the vacuum and further away. Having only two possibilities for α allows us to solve analytically the most general local strategy, since estimating the auxiliary state becomes a discrimination problem between the states √ nα 0 + 1 and √ nα 0 − 1 . The measurement that distinguishes the two possibilities is a two-outcome POVM E = {[e + ], [e − ]} 11 . We use the displacement operator (8) to shift to the local model around α 0 , such that the state of the auxiliary mode is now either |1 or |−1 . Note that, without loss of generality, one can confine the POVM vectors to the (Bloch) plane spanned by |1 and |−1 , so that |e + and |e − are real linear combinations of these. Indeed, any component orthogonal to this plane would give no aid to distinguish the hypotheses. This allows us to express the probabilities of correctly identifying each state as where χ = 1|−1 = e −2 , and the overlap c completely parametrises the measurement E. If the optimal estimation measurement is indeed asymmetric, it should happen that p + < p − , i.e., that the probability of a correct identification is greater for the state |−1 than for |1 . From now on we proceed as for the Gaussian E&D strategy. We first compute the posterior state of the signal mode according to Bayes' rule. Then, we compute the optimal error probability in the discrimination of [−α 0 ] and the posterior state, which is a combination of [1/ √ n] and [−1/ √ n], weighted by the corresponding posterior probabilities. The cdependence is carried by these probabilities. Going to the asymptotic limit n → ∞, applying perturbation theory to compute the trace norm, and averaging the result over the two possible outcomes in the discrimination of the signal state, we finally obtain the asymptotic average error probability for the local strategy as a function of c. The asymptotic average error probability for the optimal collective strategy in this simple case is obtained exactly along the same lines as shown in Section 2, and the one for known states is given by the asymptotic expansion of Eq. (16), substituting the average over G(u) appropriately. Now we can compute the excess risk for the local and collective strategy, and optimise the local one over c. As already advanced at the beginning, the optimal solution yields c * = 1 + χ + 1 − χ /2, and therefore p + = p − . That is, the POVM E is symmetric with respect to the vectors |1 and |−1 , hence both hypotheses receive the same treatment by the measurement in charge of determining the state of the auxiliary mode. Moreover, the gap between the excess risk of both strategies remains. This result leads us to conjecture that the optimal collective strategy performs better than any local strategy.

Discussion
In this paper we have proposed a learning scheme for coherent states of light. We have presented it in the context of a quantum-enhanced readout of classically-stored binary information, following a recent research line initiated in [29]. The reading of information, encoded in the state of a signal that comes reflected by a memory cell, is achieved by measuring the signal and deciding its state to be either the vacuum state or some coherent state of unknown amplitude. The effect of this uncertainty is palliated by supplying a large number of auxiliary modes in the same coherent state. We have presented two strategies that make different uses of this (quantum) side information to determine the state of the signal: a collective strategy, consisting in measuring all modes at once and making the binary decision, and a local (E&D) strategy, based on first estimating-learning-the unknown amplitude, then using the acquired knowledge to tune a discrimination measurement over the signal. We have showed that the former outperforms any E&D strategy that uses a Gaussian estimation measurement over the auxiliary modes. Furthermore, we conjecture that this is indeed the case for any (even possibly non-Gaussian) local strategy, based on evidence obtained within a simplified version of the original setting that allowed us to consider completely general measurements.
Previous works on quantum reading rely on the use of specific preparations of nonclassical-namely, entangled-states of light to improve the reading performance of a classical memory [29,34,35,36].
Our results indicate that, when there exists some uncertainty in the states produced by the source (and, consequently, the possibility of preparing a specific entangled signal state is highly diminished), alternative quantum resources-namely, collective measurements-still enhance the reading of classical information using uncorrelated, classical coherent light. It is worth mentioning that there are precedents of quantum phenomena of this sort providing enhancements for statistical problems involving coherent states. As an example, in the context of estimation of product coherent states, the optimal measure-and-prepare strategy on identical copies of |α can be achieved by local operations and classical communication (according to the fidelity criterion), but bipartite product states |α |α * require entangled measurements [52].
On a final note, the quantum enhancement found here is relevant in the regime of lowenergy signals 12 (small coherent amplitudes). This is in accordance to the advantage regime provided by nonclassical light sources, as discussed in other works [29,35,37]. A low energy readout of memories is, in fact, of very practical interest. While, mathematically, the success probability of any readout protocol could be arbitrarily increased by sending signals with diverging energy, there are many situations where this is highly discouraged. For instance, the readout of photosensitive organic memories requires a high level of control over the amount of energy irradiated per cell. In those situations, the use of signals with very low energy benefits from quantum-enhanced performance, whereas highly energetic classical light could easily damage the memory.
as the result of taking the limit n → ∞ up to order 1/n in Eq. (10). Computing the arising averages (see Appendix D), the terms in Eq. (A.2) take the explicit form We now apply perturbation theory to compute the trace norm σ 1 −σ 2 1 in the asymptotic limit n → ∞, up to order 1/n, using Eqs. (A.1) and (A.2). We start by expressing the trace norm as where A =σ 1 −σ (0) 2 , B = −σ (1) 2 , C = −σ (2) 2 , and γ j is the jth eigenvalue of Γ, which admits an expansion of the type γ j = γ (0) j + γ (1) j / √ n + γ (2) j /n. The matrix Γ belongs to the Hilbert space H ∞ ⊗ H 3 , i.e., the first mode is described by the infinite dimensional space generated by the Fock basis, and the second mode by the three-dimensional space spanned by the linearly independent vectors {|−α 0 , |0 , |1 } (we will see that the contribution of |2 vanishes, hence it is not necessary to consider a fourth dimension). Writing the eigenvalue equation associated to γ j and separating the expansion orders, we obtain the set of equations where ψ j is the eigenvector associated to γ j , which also admits the expansion ψ j = ψ (0) j + ψ (1) j / √ n + ψ (2) j /n. Eq. (A.7) tells us that γ (0) j is an eigenvalue of A with associated eigenvector ψ (0) j . We multiply (A.8) and (A.9) by ψ (0) j to obtain Note that Eq. (A.11) assumes that there is no degeneracy in the spectrum of Γ at zero order (as we will see, this is indeed the case). From the structure of A we can deduce that the form of its eigenvector ψ (0) j is where we have replaced the index j by the pair of indices i, ε. The index i represents the Fock state |i in the first mode, and the vectors |v ε are eigenvectors of [−α 0 ] − [0] and form a basis of H 3 in the second mode. Every eigenvalue of Γ is now labelled by the pair of indices i, ε, where i = 0, . . . , ∞ and ε = +, −, 0: the second mode in A has a positive, a negative, and a zero eigenvalue, to which we associate eigenvectors |v + , |v − and |v 0 , respectively. It is straightforward to see that the first two are 14) The third eigenvector |v 0 is orthogonal to the subspace spanned by |−α 0 and |0 , and corresponds to the eigenvalue γ (0) i,0 = 0 13 . This eigenvector only plays a role through the overlap 1|v 0 , which arises in Eqs. (A.10) and (A.11). We thus do not need its explicit form, but it will suffice to express 1|v 0 in terms of known overlaps.  (2) i,ε as γ (2) i Now that we have computed the eigenvalues of Γ, we are finally in condition to evaluate the sum in the right-hand side of Eq. (A.6). Incorporating the relevant eigenvalues, given by Eqs. (A.14) and (A.15), it reads 13 Note that the zero-order eigenvalues γ (0) i,ε are nondegenerate, hence Eq. (A.11) presents no divergence problems. where Appendix B. Conditional probability p(β|α), Eq. (24) Given two arbitrary Gaussian states ρ A , ρ B , the trace of their product is where V A and V B are their covariance matrices and δ is the difference of their displacement vectors. For the states ρ A ≡ [ √ nα] and ρ B ≡ Eβ, we have where α = a 1 + ia 2 ,β =b 1 + ib 2 , r is the squeezing parameter, and φ indicates the direction of squeezing in the phase space. In terms of α andβ, Eq. (B.1) reads

Appendix C. Trace norm for the E&D strategy
For assessing the performance of the E&D strategy, we want to obtain the error probability in discriminating the state [0] and the posterior state ρ(β), resulting from a heterodyne estimation of the state of the auxiliary mode that provides the estimate β. Under a local Gaussian model around α 0 parametrised by the complex variables u and v, these states transform into [−α 0 ] and ρ(v), respectively, where the second is given by and where p(u|v) is given by Eq. (30). The error probability is determined by the trace norm To compute it, we first series expand ρ(v) in the limit n → ∞, up to order 1/n. We name the appearing integrals of u, u * , |u| 2 , u 2 , and (u * ) 2 over the probability distribution p(u|v) as I 1 , I * 1 , I 2 , I 3 , and I * 3 , respectively. This allows us to write the trace norm as and λ κ is the κth eigenvalue of Φ, which admits the perturbative expansion λ κ = λ (0) κ +λ (1) κ / √ n+ λ (2) κ /n, just as its associated eigenvector ϕ κ = ϕ (0) κ + ϕ (1) κ / √ n + ϕ (2) κ /n. Up to order 1/n, the matrix Φ has effective dimension 4 since it belongs to the space spanned by the set of linearly independent vectors {|−α 0 , |0 , |1 , |2 }. Hence the index κ has in this case four possible values, i.e., κ = +, −, 3, 4. The zero-order eigenvalues λ (0) κ , which correspond to the eigenvalues of the rank-2 matrix A , are . Their associated eigenvectors are |ϕ (0) κ = |v κ , where |v ± is given by Eq. (A.13), and, by definition, v κ |−α 0 = v κ |0 = 0 for κ = 3, 4. From analogous expressions to Eqs. (A.10) and (A.11) we can write the first and second-order eigenvalues as The needed overlaps for computing λ (1) κ and λ (2) κ are given by Eqs. (A.16), (A.17), and The expressions for the overlaps (C.1) and (C.2) actually depend on the dimension of the space that we are considering (four in this case), and they are not unique: there are infinitely many possible orientations of the orthogonal pair of vectors {|v 3 , |v 4 } such that both of them are orthogonal to the plane formed by {|−α 0 , |0 }, which is the only requirement we have. However, one can verify that this choice does not affect the trace norm Φ 1 , thus we are free to choose the particular orientation that, in addition, verifies v 3 |2 = 0, yielding the simple expressions (C.1) and (C.2).
Finally, we write down the trace norm as which we use now to obtain the asymptotic expression for the average error probability, defined in Eq. (26). Recall Eq. (29) and note that we have to average Eq. (C.3) over the probability distribution p(v). Regarding this average, it is worth taking into account the following considerations. First, the v-dependence of the eigenvalues comes from I 1 , I 2 , I 3 , and its complex conjugates. The integrals needed are given in the last part of Appendix D. Second, the integration yields λ (1) κ = 0 and hence the order 1/ √ n term vanishes, as it should. And third, the second-order eigenvalues λ (2) 3 and λ (2) 4 are v-independent and positive, so we can ignore the absolute values in Eq. (C.3). Putting all together, we can express the asymptotic average error probability of the E&D strategy as P E&D e ≡ P E&D e (n → ∞)

Appendix D. Gaussian integrals
At many points in this paper, we integrate complex-valued functions over the complex plane, weighted by the bidimensional Gaussian probability distribution G(u). This section gathers the integrals that we need. Recall that G(u) is defined as G(u) = 1 πµ 2 e −|u| 2 /µ 2 , u ∈ C . Expressing u either in polar or Cartesian coordinates in the complex plane, i.e., u = re iθ = u 1 + iu 2 , one can readily check that G(u) is normalised: Note that the result of averaging a coherent state over G(u) is nothing else than a thermal state with average photon number µ 2 .