Optimal control of non-Markovian dynamics in a single-mode cavity strongly coupled to an inhomogeneously broadened spin ensemble Dmitry O. Krimer§,1,∗ Benedikt Hartl§,1 Florian Mintert,2 and Stefan Rotter1 1Institute for Theoretical Physics, Vienna University of Technology (TU Wien), Wiedner Hauptstraße 8-10/136, 1040 Vienna, Austria 2Department of Physics, Imperial College, SW7 2AZ London, United Kingdom Ensembles of quantum mechanical spins offer a promising platform for quantum memories, but proper functionality requires accurate control of unavoidable system imperfections. We present an efficient control scheme based on a Volterra equation, which relies solely on very weak control pulses and allows us to overcome the challenge posed by macroscopically large Hilbert spaces. The viability of our approach is demonstrated in terms of explicit storage and readout sequences for a 7 spin ensemble strongly coupled to a single-mode cavity. 1 0 2 In the last decade, we have witnessed tremendous the inversion of atomic detunings during the temporal n progress in the implementation of elementary operations evolution. Mostofthetechniquesdevelopedforthispur- a for quantum information processing. Single qubit gates posearebasedonphoton-echotypeapproachesincavity J can be realized with fidelities reaching 1−10−6 [1], and or cavity-less setups, such as those dealing with spin- 2 also two-qubit gates can be implemented in a variety refocusing [14, 15], with atomic frequency combs (AFC) 2 of systems [2, 3]. With all these elements at hand, it [12, 16–20] or with electromagnetically induced trans- ] is nowadays possible to implement quantum algorithms parency (EIT) [21]. Traditionally, these architectures h on architectures with a few (on the order of five) qubits operateintheopticalregionandrequireadditionalhigh- p - [4]. Implementingquantumlogicsonlargerarchitectures, intensity control fields. The resulting large number of nt however, will most likely require a separation between excitations is prone to spoil the delicate quantum infor- a quantum processing units and quantum memory units, mationthatisencodedinstateswithextremelylownum- u where qubits in the former units admit fast gate opera- bers of excitations. It would therefore be much better to q tions and the qubits in the latter units offer long coher- work with low-intensity control fields, which, however, [ ence times. have the other problem to become easily correlated with 1 Since extended coherence times naturally imply weak the quantum memory. For the identification of control v interactions with other degrees of freedom, the suffi- strategies, this implies that one may no longer treat the 4 ciently fast swapping of quantum information between manydifferentmemoryspinsasindependentobjects,but 2 2 processing and memory units is a challenging task. The thatthe(macroscopically)largeensembleneedstobede- 6 most promising route to overcome slow swapping is the scribed by a quantum many-body state. This makes any 0 encoding of quantum information as a collective excita- description of dynamics and an identification of control . 1 tion in a large ensemble composed of many (N) con- strategies a seemingly hopeless task. 0 stituents, since this increases the swapping speed by a √ In this Letter we develop a very efficient optimization 7 factor of N. Among promising realizations of such en- technique based on a Volterra integral equation, which 1 sembles those based on spins, atoms, ions or molecules : allows us to write information into a spin ensemble cou- v areofparticularinterest[5–11]. Inmanycases, however, pled to a single cavity mode by means of optimized mi- i system imperfections result in broadening effects giving X crowave pulses and to retrieve it at some later time in rise to rapid dephasing of ensemble constituents - a re- r the form of well separated cavity responses. In contrast a strictionthatlimitsthecoherencetimesofsuchcollective to established echo techniques our scheme only involves quantum memories. low-intensity signals and therefore diminishes the influ- As a result, various protocols to ensure the controlled ence of noise caused by writing and reading pulses. The andreversibletemporaldynamicsinthepresenceofinho- efficiency of our approach is demonstrated in conjunc- mogeneous broadening were the subject of many studies tion with a spectral hole-burning technique [22–24] that during the last decade. One of the proposed techniques allows us to reach storage times going far beyond the in this context is the so-called controlled reversible inho- dephasing time of the inhomogeneously broadened en- mogeneous broadening (CRIB) approach [12, 13], which semble. is based on a rather subtle preparation method and on To be specific, we consider an ensemble of spins strongly coupled to a single-mode cavity via magnetic or electricdipoleinteractionassketchedinFig.S1. Alltyp- §Theseauthorscontributedequallytothiswork. icalparametervalues arechosenhereinaccordancewith ∗ [email protected] the recent experiment [25], the dynamics of which can 2 Starting from a polarized state with all spins in their groundstate, eachofthelogicalstates|0(cid:105)and|1(cid:105)canbe prepared by driving the cavity with the corresponding write pulses, η(W)(t) and η(W)(t) during the write sec- |0(cid:105) |1(cid:105) tion (i). In the strong-coupling regime the information is then transferred from the cavity to the spin ensemble. During the delay section (ii) the information is subject to dephasing by the inhomogeneous ensemble broaden- ing and the external drive is optimized here to reduce the cavity amplitude during the entire delay (to prevent the information in the spin ensemble from leaking back FIG. S1. (Color online.) Schematics of a single-mode cavity to the cavity prematurely). In the readout section (iii) characterizedbyafrequencyωc andalossrateκ, coupledto we switch on the readout pulse η(R)(t), which transfers an ensemble of two level atoms (violet spheres) with individ- the information from the spin ensemble back to the cav- ual transition frequencies ω and a loss rate γ (cid:28)κ. Colored curves designate optimizedkinput and (non-overlapping) out- ity, resulting in the cavity amplitude A(R)(t) or A(R)(t) |0(cid:105) |1(cid:105) put signals. respectively. (Forthesakeofsimplicitywedonotexplic- itly specify the delay pulse but impose on the readout pulse η(R)(t) a constraint such that the cavity responses beexcellentlydescribedintermsoftheTavis-Cummings aremaximallysuppressedinthedelay section.) Thegoal Hamiltonian [26] (in units of (cid:126)) of our work is to find optimal time-dependent choices for η(W)(t), η(W)(t) and η(R)(t), such that A(R)(t) and H=ωca†a+ 21(cid:88)N ωjσjz+i(cid:88)N (cid:2)gjσj−a†−gj∗σj+a(cid:3)− A(|1R(cid:105))(|0t(cid:105)) have m|1(cid:105)inimal overlap to be able to dis|0c(cid:105)riminate j j between them. −i(cid:2)η(t)a†e−iωpt−η(t)∗aeiωpt(cid:3) . (1) In order to develop an efficient framework we employ a description in terms of Volterra equations that relate Hereσ±, σz arethePaulioperatorsassociatedwitheach cavity amplitudes A(t) and pump profiles η(t) [1] j j individualspinoffrequencyω anda†, aarecreationand j (cid:90) t annihilationoperatorsofthesinglecavitymodewithfre- A(t)=Ω2 dτ K(t−τ)A(τ)+D(t) , (2) quency ω . An incoming signal is characterized by the c 0 carrierfrequencyω andbytheenvelopeη(t)whosetime p variation is slow as compared to 1/ω . The interaction where D(t) depends on the time integral of the driving p part of H is written in the dipole and rotating-wave ap- signal(seeSupplementaryNote1)andontheinitialcon- proximation (terms ∝ aσ−, a†σ+ are neglected), where ditions for the cavity amplitude and for the state of the j j (cid:82) g isthecouplingstrengthofthej-thspin. Thedistance spin ensemble. K(t−τ) = dωρ(ω)S(ω,t,τ) is a mem- j between spins is assumed to be large enough such that ory kernel functionthatcharacterizesthespin-cavitydy- the direct dipole-dipole interactions between spins can namics (see Supplementary Note 1). This description is be neglected. Furthermore, the large number of spins valid in the regime of weak driving powers and for a suf- allows us to enter the strong-coupling regime of cavity ficientlylargespinensemblethatcanbecharacterizedin QED with the collective coupling strength, Ω2 =(cid:80)Ng2 terms of a continuous spectral density ρ(ω). Due to the j j [27], whichleadstotheenhancementofasinglecoupling linearity of the Volterra equations, rescaling their solu- √ strength g , by a factor of N. tions by a global prefactor leaves them perfectly valid. j We are aiming at the transfer of states from the cav- Here we take the amplitude of the write pulses, η(W)(t), |0/1(cid:105) ity to the spin ensemble, its storage over a well-defined such that the net power injected into the cavity corre- period of time, and its transfer back to the cavity. Our sponds to the power of a coherent driving signal with an control scheme thus consists of a write and readout sec- amplitudeequaltothecavitydecayrate,η =κ. Thelat- tion, with a variable delay section in between. Specifi- ter pumping prepares, on average, a single photon in the cally, we will construct (i) write pulses that encode the empty cavity for stationary transmission experiments. logicalstates|0(cid:105)and|1(cid:105)inthespinensemble,(ii)apulse Another consequence of the linearity of the govern- that acts during the delay between writing and readout, ing equations is that, given two pump profiles η (t), 1/2 and(iii)areadoutpulsethatmapsthetwologicalstates resulting in the two corresponding cavity amplitudes of the spin ensemble on two mutually orthogonal states A (t), any coherent superposition of pump profiles 1/2 of the cavity field. Note that the write pulses (i) are c η (t)+c η (t) will result in the corresponding cavity 1 1 2 2 specific for input states |0(cid:105) and |1(cid:105), but pulses (ii) and amplitudes c A (t)+c A (t). As was demonstrated in 1 1 2 2 (iii) are generic as they must be designed without prior [22,25]thecavityspindynamicsundertheactionofclas- knowledge of the information stored in the ensemble. sical pulses can immediately be transferred to the quan- 3 2 2 ) ) t 0 t 0 ( ( η η 2 2 2 2 ) ) t 0 t 0 ( ( η η 2 2 1 ×10−3 5 ×10−3 250 42.0 ×10−4 2 2 ω) t)| | | ( ( ) ) ρ A (t (t | A A 0 0.0 | | 0.992 1.0 1.008 1.07 t[µs] 1.15 ω/ω s 0 0 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 t[ns] t[µs] FIG. S2. (Color online.) Preparation of the spin ensemble configurations, |0(cid:105) and |1(cid:105), for the case of a spin density ρ(ω) following a q-Gaussian distribution centered around the cavity frequency ω [1]. In the right column two additional holes were c burnt into ρ(w) at frequencies ω ±Ω (see green arrows in the inset) to suppress decoherence [22, 23] and to make room for a s delaysection[whitearea]betweenthewrite[green(lightgray)area]andreadout[grayarea]sections. Notethevastlyextended time scale on the x-axis in the right column. First row: real [blue (dark gray)] and imaginary parts [cyan (light gray)] of the optimizedwritepulseη(W)(t)forstate|0(cid:105)andofthereadoutpulseη(R)(t)(blackandgray). Theendofthewritesection[green |0(cid:105) (light gray) area] is indicated by the vertical dashed lines. Second row: real [orange (light gray)] and imaginary parts [brown (dark gray)] of the optimized write pulse η(W)(t) for state |1(cid:105). The read pulse η(R)(t) is the same as for state |0(cid:105). Third row: |1(cid:105) cavityamplitude|A(t)|2 fortheresultingnon-overlappingtime-binnedstates|0(cid:105)[blue(darkgray)]and|1(cid:105)[orange(lightgray)]. The carrier frequency of all pulses matches the resonance condition, ω = ω = 2π·2.6915GHz, and the coupling strength p c Ω=2π·12.5MHz. The ratio of the powers between the readout and write pulses is 0.068 (0.013) for the case without (with) hole burning. The amplitudes of the write and read pulses are presented in units of κ=2π·0.4MHz (see also the main text and Supplementary Note 4). tum case because of the linearity of our model. astime-binnedstatessincethisscenarioisreminiscentof As a first step we need to find optimal write and read- the concept of time-binned qubits where information is out pulses for the logical states. We do this through the stored in the occupation amplitudes of two well distin- optimization of the following functional guishable time bins [17, 29]. (cid:10)(cid:12) (cid:12)(cid:11) Thebottleneckforextendedinformationstoragetimes η(W),mη(Win),η(R) (cid:12)A|0(cid:105)(t)−A|1(cid:105)(t)(cid:12) t , (3) in the ensemble is its inhomogeneous broadening, as de- |0(cid:105) |1(cid:105) terminedbythecontinuousspectraldensityρ(ω)appear- involving both the write pulses and the readout pulse ing in our theoretical description. Specifically, the total to minimize the overlap between the cavity amplitudes decoherencerate,Γ,inthelimitofstrongcoupling(when of the logical states |0(cid:105) and |1(cid:105) in the readout section. Ω > Γ) can be estimated as Γ ≈ κ + πΩ2ρ(ω ± Ω) s In practice we expand all involved driving pulses in a [1, 30], indicating that the dominant contribution to Γ basis of trial functions with unknown coefficients and stems from the spectral density ρ(ω) at frequencies close construct the functional Eq. (3) with all these coeffi- to the maxima of the two polaritonic peaks, ω =ω ±Ω. s cients contained as unknown variables. We then search To suppress this decoherence rate Γ it is thus advisable the functional’s minima under several constraints using to work with spin ensembles having a spectral density the standard method of Lagrange multipliers (see Sup- thatfallsofffasterthan1/Ω2 initstailssuchthatΓ→κ plementary Notes 1,2). for large Ω. The corresponding “cavity protection ef- A typical result of this optimization (first without a fect” [1, 30, 31] has meanwhile been demonstrated also delay section) is depicted in Fig. S2 (left column), where experimentally [25], but has the drawback of requiring the amplitudes of all optimized pulses as well as that of prohibitively large coupling strengths to take full effect. the resulting cavity responses are depicted. One can in- Alternatively, one can burn two narrow spectral holes at deed see that the two different configurations stored in frequencies close to ω ±Ω. This technique [22–24] was s the spin ensemble, |0(cid:105) and |1(cid:105), are retrieved by the same recentlyshowntobebotheasilyimplementableandvery readout pulse in the form of two well-separated cavity efficient in suppressing the decoherence rate Γ even be- responses. In the following we refer to these responses low the bare cavity decay rate κ [23]. Incorporating this 4 ideally, the corresponding superposition of time-binned cavityresponseswouldbeobservedundertheapplication of the readout pulse η(R)(t). Since the cavity response is oftheformA(R)(t;α,β)=α·A˜(R)(t)+β·A˜(R)(t)+A˜(R)(t) |0(cid:105) |1(cid:105) where the two cavity responses A˜(R) (t) only depend on |0/1(cid:105) the stored spin configurations |0(cid:105) and |1(cid:105), and A˜(R)(t) is the response induced by the readout pulse, the desired superpositionofcavityoutputsisobtainedifthereadout pulsesatisfies(α+β)η(R)(t)=η(R)(t). Togetherwiththe normalization|α|2+|β|2 =1thisimpliesthatfortheam- (cid:112) (cid:112) plitudesα =1−x±i x(1−x)andβ =x∓i x(1−x) x x with x ∈ [0,1] the desired cavity response will be ob- tained. As a result, the proposed storage sequence does not only work for the two basis states |0/1(cid:105), but, indeed foraone-dimensionalsetofcoherentsuperpositions,such as for a rebit [29, 32]. Note that when being only interested in reading out the parameters α and β (and not in further processing the resulting cavity response) one is not restricted by FIG.S3. (Coloronline.) (a)Retrievedcavityamplitudeinthe the above rebit parametrization, but has the full qubit readoutsection,resultingfromasuperpositionofwritepulses, parameter space at one’s disposal. As we show in Sup- αx·A(|0R(cid:105))(t)+βx·A(|1R(cid:105))(t)(normalizedtoamaximumvalueof1) plementaryNote2,αandβ canbeunambiguouslydeter- ianndthβe a(bsseeentceexot)ffnrooimse.xW=e0utsoed1tinhestreepbsitopfa0r.2a5m(eptrairzatmioenteαrxs mined through the overlaps O0/1 = (cid:82)ττacdtA(R)(t;α,β)· as inxFig. S2, left column). (b,c) Retrieved average values A|(0R/)1∗(cid:105)(t), where A(|0R/)1(cid:105)(t)=A˜(|0R/)1(cid:105)(t)+A˜(R)(t). (cid:104)αR(cid:105)and(cid:104)βR(cid:105)(onlyrealpartsareshown)fromtheresulting In principle, this information retrieval is exact, but solution A(R)(t;α,β) in the presence of noise. The averag- noise (which is not included in the previous theoretical ing is performed with respect to 200 noise realizations for a modelling) affects the readout if it reaches values com- noise amplitude δη/η = 0.05, where η is the write ampli- 0 0 parable to the cavity amplitudes. Therefore in the next tude without noise. The input configurations are parameter- izedasα=cos(ϑ/2)andβ =sin(ϑ/2)eiϕ withϑ∈[0,π]and line of our study we examine the robustness of our op- ϕ∈[0,2π]. (d)ReconstructedBlochspherewithspatialcom- timal control scheme against possible noise. For that ponents r =((cid:104)α (cid:105)∗,(cid:104)β (cid:105)∗)·σ ·((cid:104)α (cid:105),(cid:104)β (cid:105)) for i=x, y, z, purpose,wesubjectthepreviouslyestablishedoptimized i R R i R R evaluatedfromtheretrievedaveragedparameterstakenfrom pulses, η(W)(t), η(W)(t) and η(R)(t), to a small pertur- (b,c), where σ is the ith Pauli matrix. The colored symbols |0(cid:105) |1(cid:105) i bation by adding Gaussian white noise as an additional in (b-d) emphasize the rebit encoding from (a). The refer- driving term in our Volterra equations (see Supplemen- ence states, |0/1(cid:105), are taken from the left column of Fig. S2. (e) Maximum of the absolute errors, (cid:15) = |α−(cid:104)α (cid:105)| and tary Note 3). We treat the problem numerically using α R (cid:15) = |β −(cid:104)β (cid:105)|, in retrieval of the input configurations for well-established methods for integrating stochastic dif- β R differentnoiseamplitudesδη0. Verticaldashedlineshowsthe ferential equations (see, e.g., [3]) and accumulate statis- levelofnoiseatwhichthecalculationsfor(b-d)areperformed. tics by integrating many trajectories for different noise realizations. Wethenaveragetheresultingretrievedval- ues with respect to noise realizations and calculate the hole burning protocol in the present analysis allows us absolute retrieval errors as the deviation from the input to increase the dephasing time from 1/Γ ∼ 75ns [see configuration, (cid:15) =|α−(cid:104)α (cid:105)| and (cid:15) =|β−(cid:104)β (cid:105)|. The α R β R Fig. S2 (left column)] to microsecond time scales [see typicalresultsofourcalculationsaredisplayedinFig.S3. Fig. S2 (right column)] for which we can now meaning- It turns out that (cid:15) and (cid:15) scale approximately linearly α β fully introduce a delay section in between the write and withthenoiseamplitudeand,e.g.,themaximalabsolute the readout section. In Fig. S2 (right column) we show error of retrieval shown in Fig. S3 is at most 0.02 for 200 thatwithparameterstakenfromrecentexperiments[23] noise realizations when taking the noise amplitude to be we extend the storage time and thereby our method’s 5% of the incoming amplitude of the write pulse. These temporal range of control beyond one micro-second. results confirm the robustness of our approach with re- With these long coherence times we can now proceed spect to possible noise in a real physical system. to the main goal of storing coherent superpositions of Beyond the prospect of substantially advancing the the two spin configurations, |0(cid:105) and |1(cid:105). Those can be storagecapabilitiesforquantuminformation,ourpresent created by the corresponding superposition η(W)(t) = workalsohasfarreachingconsequencesforquantumcon- α·η(W)(t)+β·η(W)(t)oftherespectivewritepulses,and, trol. Given the extremely unfavorable scaling properties |0(cid:105) |1(cid:105) 5 ofcompositequantumsystemswithparticlenumber,any [10] M. Afzelius, C. Simon, H. de Riedmatten, and N. Gisin, theoretical description of a quantum many-body system Phys. Rev. A. 79, 052329 (2009). is an extremely challenging task. Since the identifica- [11] J.Verdu´, H.Zoubi,Ch.Koller, J.Majer,H.Ritsch, and J. Schmiedmayer, Phys. Rev. Lett. 103, 043603 (2009). tionofoptimalcontrolstrategiesistypicallymuchharder [12] W. Tittel, M. Afzelius, T. Chaneliere, R. L. Cone, S. than the mere description of a system’s dynamics (the Kro¨ll, S. A. Moiseev, and M. Sellars, Laser & Photon. latter is naturally required for the former), optimal con- Rev. 4, 244 (2010). trol is typically a viable option for rather small systems [13] B. Kraus, W. Tittel, N. Gisin, M. Nilsson, S. Kro¨ll, and only. Our efficient method to control a large quantum J. I. Cirac, Phys. Rev. A 73, 020302 (R)(2006). system even in the presence of inhomogeneous broaden- [14] B. Julsgaard, and K. Mo¨lmer, Phys. Rev. A 88, 062324 ing thus provides not only a vital step towards working (2013). quantummemories,butitalsoopensupanewavenuefor [15] C. Grezes, B. Julsgaard, Y. Kubo, M. Stern, T. Umeda, J. Isoya, H. Sumiya, H. Abe, S. Onoda, T. Ohshima, V. thedesignofarchitecturesfortheprocessingofquantum Jacques, J. Esteve, D. Vion, D. Esteve, K. Mo¨lmer, and information in general. P. Bertet, Phys. Rev. X, 4, 021049 (2014). We would like to thank Himadri Dhar and Matthias [16] N. Gisin, S. A. Moiseev, and C. Simon, Phys. Rev. A, Zens for helpful discussions and acknowledge support 76, 014302 (2007). by the Austrian Science Fund (FWF) through Project [17] M. U. Staudt, S. R. Hastings-Simon, M. Nilsson, M. No. F49-P10 (SFB NextLite). F.M. acknowledges sup- Afzelius, V. Scarani, R. Ricken, H. Suche, W. Sohler, W. port by the European Research Council (ERC) grant Tittel, and N. Gisin, Phys. Rev. Lett. 98, 113601 (2007). [18] H.deRiedmatten,M.Afzelius,M.U.Staudt,C.Simon, “odycquent”. and N. Gisin, Nature 456, 773 (2008). [19] M. Gu¨ndogan, P. M. Ledingham, K. Kutluer, M. Mazzera, and H. de Riedmatten, Phys. Rev. Lett. 114, 230501 (2015). [1] T.P.Harty,D.T.C.Allcock,C.J.Ballance,L.Guidoni, [20] P. Jobez, I. Usmani, N. Timoney, C. Laplane, N. Gisin, H. A. Janacek, N. M. Linke, D. N. Stacey, and D. M. and M. Afzelius, New J. of Phys., 16 083005 (2014). Lucas, Phys. Rev. Lett., 113, 220501 (2014). [21] I. Novikova, N. B. Phillips, and A. V. Gorshkov, Phys. [2] J.L.O’Brien,G.J.Pryde,A.G.White,T.C.Ralphand Rev. A. 78, 021802(R) (2008). D. Branning, Nature 426, 264 (2003). [22] D. O. Krimer, B. Hartl, and S. Rotter, Phys. Rev. Lett. [3] J. H. Plantenberg, P. C. de Groot, C. J. P. M. Harmans 115, 033601 (2015). and J. E. Mooij, Nature 447, 836 (2007). [23] S. Putz, A. Angerer, D. O. Krimer, R. Glattauer, W. J. [4] P.J.J.OMalley,R.Babbush,I.D.Kivlichan,J.Romero, Munro,S.Rotter,J.Schmiedmayer,andJ.Majer,Nature J.R.McClean,R.Barends,J.Kelly,P.Roushan,A.Tran- Phot. 11, 36 (2017). ter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, [24] D. O. Krimer, M. Zens, S. Putz, and S. Rotter, Laser & A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Photonics Review 10, 1023 (2016). Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, [25] S. Putz, D. O. Krimer, R. Amsu¨ss, A. Valookaran, T. D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. No¨bauer, J. Schmiedmayer, S. Rotter, andJ. Majer, Nat. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. Phys. 10, 720 (2014). M. Martinis, Phys. Rev. X 6, 031007 (2016). [26] M. Tavis and F. W. Cummings, Phys. Rev. 170, 379 [5] P.L.Stanwix,L.M.Pham,J.R.Maze,D.LeSage,T.K. (1968). Yeung, P. Cappellaro, P. R. Hemmer, A. Yacoby, M. D. [27] K.Sandner,H.Ritsch,R.Amsu¨ss,C.Koller,T.No¨bauer, Lukin,andR.L.Walsworth,Phys.Rev.B82,201201(R) S.Putz,J.Schmiedmayer,andJ.Majer,Phys.Rev.A85, (2010). 053806 (2012). [6] R.Amsu¨ss,Ch.Koller,T.No¨bauer,S.Putz,S.Rotter,K. [28] D.O. Krimer, S. Putz, J. Majer, and S. Rotter, Phys. Sandner,S.Schneider,M.Schrambo¨ck,G.Steinhauser,H. Rev. A 90, 043852 (2014). Ritsch, J. Schmiedmayer, and J. Majer, Phys. Rev. Lett. [29] J. Brendel, N. Gisin, W. Tittel, and H. Zbinden, Phys. 107, 060502 (2011). Rev. Lett., 82, 2594 (1999). [7] Y. Kubo, F. R. Ong, P. Bertet, D. Vion, V. Jacques, D. [30] I.Diniz,S.Portolan,R.Ferreira,J.M.G´erard,P.Bertet, Zheng, A. Dr´eau, J.-F. Roch, A. Auffeves, F. Jelezko, J. and A. Auff´eves, Phys. Rev. A 84, 063810 (2011). Wrachtrup, M. F. Barthe, P. Bergonzo, and D. Esteve, [31] Z.Kurucz,J.H.Wesenberg,andK.Mølmer,Phys.Rev. Phys. Rev. Lett. 105, 140502 (2010). A 83, 053852 (2011). [8] S. Probst, H. Rotzinger, S. Wu¨nsch, P. Jung, M. Jerger, [32] J.Batle,A.R.Plastino,M.Casas,andA.Plastino,Op- M. Siegel, A. V. Ustinov, and P. A. Bushev, Phys. Rev. tics and Spectroscopy 94, 700 (2003). Lett. 110, 157001 (2013). [33] R. Toral, and P. Colet, Stochastic Numerical Methods: [9] P. Rabl, D. DeMille, J. M. Doyle, M. D. Lukin, R. J. An Introduction for Students and Scientists, Wiley-VCH Schoelkopf, and P. Zoller, Phys. Rev. Lett. 97, 033003 (2014). (2006). 6 Supplementary Note 1. Volterra equation for the cavity amplitude OurstartingpointistheHamiltonian(1)ofthemainarticlefromwhichwederivetheHeisenbergequationsforthe cavity and spin operators, a˙(t) = i[H,a(t)]−κ·a(t), σ˙−(t) = i[H,σ˙−(t)]−γ·σ−(t), respectively. Here a stands for k k k the cavity annihilation operator and σ− are standard downward Pauli operators associated with the k-th spin. κ and k γ are the dissipative cavity and individual spin losses, respectively. (All notations are in tact with those introduced in the main article.) During the derivations we use the following simplifications and approximations valid for various experimental realizations: (i) k T (cid:28) (cid:126)ω (the energy of photons of the external bath, k T, is substantially smaller B c B than that of cavity photons, (cid:126)ω ); (ii) the number of microwave photons in the cavity remains small as compared c to the total number of spins participating in the coupling (limit of low input powers of an incoming signal), so that the Holstein-Primakoff-approximation, (cid:104)σ(z)(cid:105)≈−1, always holds; (iii) the effective collective coupling strength of the k spin ensemble, Ω2 = (cid:80)N g2 (g stands for the coupling strength of the kth spin), satisfies the inequality Ω (cid:28) ω , k=1 k k c justifying the rotating-wave approximation; (iv) the spatial size of the spin ensembles is sufficiently smaller than the wavelength of a cavity mode. Having introduced all these assumptions, we derive the following system of coupled first-order linear ODEs for the cavity and spin amplitudes in the ω -rotating frame p N A˙(t)=−[κ+i∆ ]A(t)+(cid:88)g B (t)−η(t), (4) c k k k=1 B˙ (t)=−[γ+i∆ ]B (t)−g A(t), (5) k k k k where A(t) ≡ (cid:104)a(t)(cid:105) and B (t) ≡ (cid:104)σ−(t)(cid:105). ∆ = ω −ω and ∆ = ω −ω are the detunings with respect to the k k c c p k k p probe frequency ω . p By formally integrating Eqs. (5) with respect to time for the spin operators and inserting them into Eq. (4) for the cavity operator, we get ∞ t N (cid:90) (cid:90) A˙(t)=−[κ+i∆c]A(t)+(cid:88)gkBk(T1) e−[γ+i∆k](t−T1)−Ω2 dωρ(ω) dτA(τ)e−[γ+i∆ω](t−τ)−η(t), (6) k=1 0 T1 where ∆ =ω−ω , B (T ) are the initial spin amplitudes at t=T and ρ(ω)=(cid:80)N g2δ(ω−ω )/Ω2 stands for the ω p k 1 1 k=1 k k continuous spectral spin distribution (see [1] for more details about this transformation). Next we formally integrate Eq. (6) in time and simplify the resulting double integral on the right-hand side by partial integration. We also consider the case when the cavity is initially empty, A(T ) = 0, and all spins are in the 1 ground state, B (T ) = 0. To speed up our numerical calculations and to separate different time sections from each k 1 other (see the main text and Supplementary Note 2 for details), we divide the whole time integration into successive subintervals, T ≤ t ≤ T , with n = 1,2,... (see Fig. 1). This allows us to derive the recurrence relation for the n n+1 cavity amplitude for the n-th time interval, A(n)(t), which depends on all previous events at t<T . Finally, we end n up with the following expression for A(n)(t) t (cid:90) A(n)(t)= dτK(t−τ)A(n)(τ)+D(n)(t)+F(n)(t), (7) Tn where the non-Markovian feedback within the n-th time interval is provided by the kernel function K(t−τ) ∞ (cid:90) e−[γ+i∆ω](t−τ)−e[κ+i∆c](t−τ) K(t−τ)=Ω2 dωρ(ω) . (8) [γ+i∆ ]−[κ+i∆ ] ω c 0 The driving term D(n)(t) in Eq. (7), t (cid:90) D(n)(t)= − dτη(n)(τ)e−[κ+i∆c](t−τ), (9) Tn includes an arbitrarily shaped, weak incoming-pulse η(n)(t), defined in the time interval [T ,T ]. The memory n n+1 contributions from all previous time intervals for t < T are given both through the amplitude A(n−1)(T ) and n n 7 S 1. Schematics of the time-divisions of the cavity-amplitude A(n)(t). The input field η(n)(t) is applied to the system in the time interval [T ,T ] and drives the corresponding cavity-amplitude A(n)(t) (indicated by vertical blue arrows). The n n+1 non-Markovian contributions from previous time intervals [T ,T ] are indicated by horizontal green arrows. n−1 n through the memory integral I(n)(ω), which are contained in the function ∞ (cid:90) e−[γ+i∆(ω)](t−Tn)−e−[κ+i∆c](t−Tn) F(n)(t)= A(n−1)(T )e−[κ+i∆c](t−Tn)+ Ω2 dωρ(ω) ·I(n)(ω) , n [γ+i∆ ]−[κ+i∆ ] ω c 0 (10) where (cid:90)Tn I(n)(ω)=I(n−1)(ω)e−[γ+i∆ω](Tn−Tn−1) + dτA(n−1)(τ)e−[γ+i∆ω](Tn−τ), (11) Tn−1 In accordance with the initial conditions introduced above at t = T , A(0)(T ) = 0 and I(1)(ω) = 0, so that F(1)(t) 1 1 vanishes in the first time interval, F(1)(t)=0 (T ≤t≤T ). 1 2 Supplementary Note 2. Optimal control based on the Volterra equation In the main text of the manuscript, we split our time interval into two parts, a write and readout section, with a variable delay section in between. In the write section, two independent optimized write pulses η(W)(t) prepare two different configurations of the spin ensemble, which are referred to as logical states |0(cid:105) and |1(cid:105) of the spin ensemble. It is followed by the delay section characterized by almost completely suppressed cavity responses, and finally by the readout section where two logical states of the spin ensemble are retrieved and mapped on two mutually orthogonal states of the cavity field by means of the readout pulse η(R)(t) (see Fig. 2). Note that the optimized readout pulse is generic being the same for both |0(cid:105) and |1(cid:105) states. For the sake of simplicity we do not explicitly specify the delay pulsebutimposeonthereadoutpulseη(R)(t)aconstraintsuchthatthecavityresponsesaremaximallysuppressedin the delay section [T ,τ ], see Fig. 2. Thus, the write and readout pulses are defined within the time intervals, [T ,T ] 2 a 1 2 and [T ,T ], respectively, in terms of the notations introduced in the Supplementary Note 1 and the delay section is 2 3 formally absorbed into the readout section. We then expand η(W)(t) and η(R)(t) in terms of sine functions (cid:88)N1 η(W)(t)= ξ ·sin(kω (t−T )), (12) k f 1 k=1 (cid:88)N2 η(R)(t)= ζ ·sin(lω (t−T )), (13) l f 2 l=1 where ξ and ζ are the expansion coefficients and ω is the fundamental frequency. (Without loss of generality, we k l f choose ω to be the same in both sections.) The linear property of the Volterra equation (7) allows us to expand f the cavity amplitude in the write section, A(W)(t), in a series of time-dependent functions with the same expansion 8 coefficients ξ as in Eq. (12), k A(W)(t)=(cid:88)N1 ξ ·a(W)(t). (14) k k k=1 Here a(W)(t) are solutions of the following Volterra equations k t t (cid:90) (cid:90) a(W)(t)= dτK(t−τ)a(W)(τ)− dτ sin(kω (τ −T ))e−[κ+i∆c](t−τ), (15) k k f 1 T1 T1 where the kernel function K(t−τ) is given by Eq. (8). The solution in the readout section, A(R)(t), in turn, consists of two contributions A(R)(t)=(cid:88)N2 ζ a(R)(t)+(cid:88)N1 ξ ψ(R)(t). (16) l l k k l=1 k=1 Similar to the ansatz for the write section, the first term in Eq. (16) also contains the same expansion coefficients ζ l as the corresponding driving signal in the readout section (see Eq. (13)) with the time-dependent functions, a(R)(t), l obeying the following Volterra equations (T ≤t≤T ) 2 3 t t (cid:90) (cid:90) a(R)(t)= dτK(t−τ)a(R)(τ)− dτ sin(lω (τ −T ))e−[κ+i∆c](t−τ). (17) l l f 2 T2 T2 Additionally, the second term in Eq. (16) describes the non-Markovian memory and appears in the readout section due to the energy stored both in the cavity and spin ensemble during the time interval T ≤ t ≤ T (write section). 1 2 Therefore, it depends only on the coefficients ξ of the write pulse (12) and the time-dependent functions ψ(R)(t), k k which can be found by substituting the expressions (14,16) into Eqs. (7-11) for n = 2. It can be shown that these functions satisfy the following Volterra equations t (cid:90) ψ(R)(t)= dτK(t−τ)ψ(R)(τ)+f(R)(t), (18) k k k T2 with the feedback from the previous write section defined by f(R)(t)=a(W)(T )e−[κ+i∆c](t−T2)+Ω2(cid:90)∞dωρ(ω)e−[γ+i∆ω](t−T2)−e−[κ+i∆c](t−T2)·(cid:90)T2dτa(W)(τ)e−[γ+i∆ω](T2−τ).(19) k k 2 [γ+i∆ ]−[κ+i∆ ] k ω c 0 T1 Note that a(W)(t) in Eq. (19) are defined in the write section only and are known solutions of Eq. (15). k Inthemaintextofourmanuscriptweusetwodifferentpulsesη(W)(t)=(cid:80)N1 ξ|0(cid:105)·sin(kω (t−T ))andη(W)(t)= |0(cid:105) k=1 k f 1 |1(cid:105) (cid:80)N1 ξ|1(cid:105)·sin(kω (t−T )) in the write section, which are characterized by two sets of expansion coefficients from k=1 k f 1 Eq. (12). As a result, the cavity amplitudes in the write section are also represented by these sets of expansion coefficients and are given by Eq. (14), namely A(W)(t)=(cid:88)N1 ξ|0(cid:105)·a(W)(t), A(W)(t)=(cid:88)N1 ξ|1(cid:105)·a(W)(t). (20) |0(cid:105) k k |1(cid:105) k k k=1 k=1 Note that by injecting these pulses into the cavity, we create two independent configurations (denoted as |0(cid:105) and |1(cid:105)) of the spin-cavity system at the beginning of the readout interval, t=T . 2 Next, we perform a readout by applying a single optimized readout pulse (13), which is the same for the states |0(cid:105) and |1(cid:105). The cavity amplitudes in the readout section, in turn, are governed by Eq. (16) as A(R)(t)=(cid:88)N2 ζ a(R)(t)+(cid:88)N1 ξ|0(cid:105)ψ(R)(t), A(R)(t)=(cid:88)N2 ζ a(R)(t)+(cid:88)N1 ξ|1(cid:105)ψ(R)(t), (21) |0(cid:105) l l k k |1(cid:105) l l k k l=1 k=1 l=1 k=1 (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) A˜(R)(t) A˜(R)(t) A˜(R)(t) A˜(R)(t) |0(cid:105) |1(cid:105) 9 S 2. Time-divisions for the optimization scheme of the time-binned cavity-responses |0(cid:105) and |1(cid:105). The write section, [T ,T ], is 1 2 followed by the variable delay section, [T ,τ ], and the readout section [τ ,τ ]. The states |0(cid:105) and |1(cid:105) reside the first, [τ ,τ ], 2 a a c a b and the second half, [τ ,τ ], of the readout section, respectively. b c where A˜(R)(t) describes the contribution from the readout pulse only which is the same for both cavity responses and the two other terms, A˜(R)(t) (i=0,1), explicitly depend on the states |0(cid:105) and |1(cid:105) created in the write section. |i(cid:105) Thus, the cavity amplitude is determined at every moment of time by Eqs. (14-19) (and, as a consequence, all spin configurations), if all expansion coefficients, ξ|0(cid:105), ξ|1(cid:105), and ζ are provided. k k l As a next step, we develop an optimization scheme aiming at achieving two well-resolved cavity responses in the readout section, A(R)(t) and A(R)(t), as is sketched in Fig. 2. (The results of numerical calculations are presented |0(cid:105) |1(cid:105) in Fig. 2 of the main paper.) For this purpose we use the standard method of Lagrange multipliers by introducing the functional F(ξ|0(cid:105),ξ|1(cid:105),ζ ) subject to several constraints listed below, and search for its minima with respect to the k k l expansion coefficients of all three pulses. Namely, we write the following expression for the functional F(ξk|0(cid:105),ξk|1(cid:105),ζl)=(cid:90) τcdt|A(|0R(cid:105))(t)|2+(cid:90) τbdt|A(|1R(cid:105))(t)|2+(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) τcdtA|(0R(cid:105))(cid:63)(t)A|(1R(cid:105))(t)(cid:12)(cid:12)(cid:12)(cid:12)−λd|0e(cid:105)lay·(cid:90) τadt|A|(0R(cid:105))(t)|2− τb τa τa T2 (cid:90) τa (cid:18)(cid:90) τb (cid:19) λ|1(cid:105) · dt|A(R)(t)|2−λ|0(cid:105)·|A(R)(τ )|2−λ|1(cid:105)·|A(R)(τ )|2−λ|0(cid:105) · dt|A(R)(t)|2−S − delay |1(cid:105) T |0(cid:105) a T |1(cid:105) a ∆T |0(cid:105) T2 τa λ|1(cid:105) ·(cid:18)(cid:90) τcdt|A(R)(t)|2−S(cid:19)−λ|0(cid:105)·(cid:32)(cid:88)|ξ|0(cid:105)|2−P(cid:33)−λ|1(cid:105)·(cid:32)(cid:88)|ξ|1(cid:105)|2−P(cid:33), (22) ∆T |1(cid:105) P k P k τb k k where λ-s are the Lagrange multipliers. The first three terms in Eq. (22) are the functions to be minimized which ensure that the overlap between the time-binned states in the readout section is negligibly small. The rest of the terms are constraints which additionally guarantee the following conditions to be simultaneously fulfilled: (i) the cavity responses within the delay section are maximally suppressed; (ii) the cavity at the beginning of the readout section is almost empty for both states; (iii) the integral taken with respect to the time-binned cavity amplitudes squared within the readout section has the same value S; (iv) a net power P of the write pulses per fundamental period 2π/ω is the same. f In our numerical calculations we used the sequential Least Squares Programming (SLSQP) minimization method [2] embedded in the internal python library scipy.optimize to find the minima of the functional F(ξ|0(cid:105),ξ|1(cid:105),ζ ). k k l In the main text we create an arbitrary superposition of write pulses (each of which separately prepares the logical state |0(cid:105) or |1(cid:105)) by applying the superimposed write pulse η(W)(t)=α·η(W)(t)+β·η(W)(t) (23) |0(cid:105) |1(cid:105) aiming to extract the encoded information (given by complex numbers α and β) from the solution for the cavity amplitudeinthereadout sectiondesignatedbygrayareainFig.2. (Notethatthereadingpulseη(R)(t)isalwayskept the same.) The solution in the readout section can be written as A(R)(t;α,β)=α·A˜(R)(t)+β·A˜(R)(t)+A˜(R)(t), (24) |0(cid:105) |1(cid:105) whereallthreepreviouslyestablishedwell-knownamplitudesA˜(R)(t), A˜(R)(t)andA˜(R)(t)areintroducedinEq.(21). |0(cid:105) |1(cid:105) 10 We then project our resulting solution (24) onto the functions A(R)(t) and A(R)(t) from Eq. (21), namely, we write |0(cid:105) |1(cid:105) (cid:90)τc O = dtA(R)(t;α,β)·A(R)∗(t)=α·F +β·F +F , (25) i |i(cid:105) i,0 i,1 i,R τa where the overlap integrals F = (cid:82)τcdtA˜(R)(t)·A(R)∗(t) with i,q =0,1 and F = (cid:82)τcdtA˜(R)(t)·A(R)∗(t). Since F i,q |q(cid:105) |i(cid:105) i,R |i(cid:105) i,q τa τa and F are known we finally end up with the following set of two algebraic equations i,R O =α·F +β·F +F , (26) 0 0,0 0,1 0,R O =α·F +β·F +F , (27) 1 1,0 1,1 1,R from which the retrieved values α and β can be evaluated. R R Supplementary Note 3. Retrieval of encoded parameters in the presence of noise Here we study the influence of noise on the quality of our optimization scheme presented in the main article and introduced in the Supplementary Note 2. For that purpose, we subject the previously established optimal driving amplitudes, η(W)(t), η(W)(t) and η(R)(t) (see Supplementary Note 2), to a small perturbation represented by the |0(cid:105) |1(cid:105) driving term, δη (t) = δη·υ(t), where δη is the amplitude of perturbation and υ(t) stands for a Gaussian white noise noiseofmeanandcorrelationsgivenby,respectively,(cid:104)υ(t)(cid:105)=0and(cid:104)υ(t(cid:48))υ(t)(cid:105)=δ(t−t(cid:48)). Wethennumericallyintegrate the Volterra equation (6) from Supplementary Note 2 with respect to time by adding the perturbation δη (t) to noise the corresponding deterministic optimal driving amplitudes η(t), which in our specific case are represented by the known writing and readout amplitudes, η(W)(t), η(W)(t) and η(R)(t). We treat the problem numerically by using the |0(cid:105) |1(cid:105) well-established numerical method for integrating stochastic differential equations (see e.g. [3]). In a nutshell, the stochasticcontributiontothecavityamplitudeistakenintoaccountaftereachtimestepofnumericalintegrationinthe √ followingway: A(t )→A(t )+ dt·δη (t ),whereA(t )afterthearrowcorrespondstothedeterministic m+1 m+1 noise m m+1 part of the cavity amplitude at t = t obtained using the standard Runge-Kutta method and δη (t ) is the m+1 noise m stochastic drive taken from the previous time step. We then accumulate statistics by integrating many trajectories fordifferentnoiserealizations. Next,weextracttheencodedparametersα andβ inthepresenceofnoisereplacing R R the overlap integrals in Eqs. (26, 27) for the case without noise by the corresponding overlap integrals evaluated for different noise realizations. The result of calculations for the average retrieval values of (cid:104)α (cid:105) and (cid:104)β (cid:105) and their R R absolute errors, (cid:15) =|α−(cid:104)α (cid:105)| and (cid:15) =|β−(cid:104)β (cid:105)|, with respect to the encoded values are depicted in Fig. 3 of the α R β R main paper. Supplementary Note 4. Numerical values for the optimized readout pulse coefficients Here we present numerical values of the coefficients ξ|0(cid:105), ξ|1(cid:105) and ζ of the optimal readout pulses η(W)(t), η(W)(t) k k l |0(cid:105) |1(cid:105) and η(R)(t) defined by Eqs. (12-13), which are presented in the main text. We take the amplitude of the write pulses such that the net power injected into the cavity, P(W) = 1 T(cid:82)fdt|η(W)(t)|2 = κ2, with i = 0,1, such that it |i(cid:105) Tf |i(cid:105) 0 corresponds to the power provided by a coherent driving signal with the amplitude equal to the cavity decay rate, η = κ. Specifically, using the expansion (12) for the write pulses η(W)(t), we obtain the following expression for the power of the write pulses per fundamental period T : f P(W) =η(W)2· 1(cid:88)N1 (cid:12)(cid:12)ξ|i(cid:105)/η(W)(cid:12)(cid:12)2 =κ2, (28) |i(cid:105) 2 (cid:12) k (cid:12) k=1 whereη(W) =κand1/2·(cid:80)N1 |ξ|i(cid:105)/κ|2 =1duetotheconstraintimposedontheexpansioncoefficients. Ontheother k=1 k hand the power of the readout pulse is substantially smaller than that of the write pulses and for the case without hole burning (see left column of Fig. 2 in the main text) we obtain 1(cid:88)N2 (cid:12) (cid:12)2 P(R) =η(R)2· (cid:12)ζ /η(R)(cid:12) =0.068·κ2, (29) 2 (cid:12) l (cid:12) k=1 where η(R) =0.26·κ and again we use as the constraint 1/2·(cid:80)N2 |ζ /η(R)|2 =1. l=1 l