Protecting operations on qudits from noise by continuous dynamical decoupling

Reginaldo de Jesus Napolitano,1 Felipe Fernandes Fanchini,2, ∗ Adonai Hilario da Silva,1 and Bruno Bellomo3 São Carlos Institute of Physics, University of São Paulo, PO Box 369, 13560-970, São Carlos, SP, Brazil Faculdade de Ciências, UNESP Universidade Estadual Paulista, 17033-360 Bauru, São Paulo, Brazil Institut UTINAM, CNRS UMR 6213, Université Bourgogne Franche-Comté, Observatoire des Sciences de l’Univers THETA, 41 bis avenue de l’Observatoire, F-25010 Besançon, France


I. INTRODUCTION
The development and implementation of effective quantum computers is of great interest to the scientific community as well as to the world economy as a whole [1]. Indeed, quantum computers promise to revolutionize many important tasks, even with a reduced number of algorithms known to be more efficient than their classical analogues [2]. In this sense, reducing errors while keeping quantum algorithms simple is an important aspect to be addressed. Several strategies have been proposed to contrast decoherence effects, such as reservoir engineering methods [3][4][5][6][7][8], optimal quantum control protocols [9,10], measurement-based control [11][12][13], and pulsed dynamical decoupling of qubits [14][15][16][17]. Some of us worked at continuous dynamical decoupling strategies, applied to single-and two-qubit systems [18][19][20][21].
Continuous dynamical decoupling techniques have been theoretically investigated and experimentally implemented in several contexts. For example, these kinds of techniques have been applied in the case of nitrogenvacancy centres to separate a single nuclear spin signal from the bath noise [22], to extend the coherence time for the electron spin [23,24] while protecting quantum gates [24], to provide single-molecule magnetic resonance spectroscopy [25], and in the context of sensing high frequency fields [26]. Furthermore, they can be used to create a dephasing-insensitive quantum computation scheme in an all-to-all connected superconducting circuit [27], and for engineering an optical clock transition in trapped ions, robust against external field fluctuations [28].
In this paper, we present a complete theoretical prescription for a generalized continuous dynamical decoupling (GCDD) of an ensemble of arbitrary qudits from environmental noise in the presence of an arbitrary quantum operation. Our prescription consists of a continuously-varying control Hamiltonian and a modification of the intended quantum operation. We first develop our procedure for the case of an arbitrary qudit and apply it to the case of a qutrit implemented with a specific model where we explicitly show how to realize all the possible SU(3) group operations which are, in general, needed for our scheme. We finally show that our scheme can be extended to the case of an ensemble of qudits, identical or not. Since our procedure works for an arbitrary number of levels and can be generalized to the case of many qudits, it extends the range of applicability of continuous dynamical decoupling strategies to more complex systems with respect to previous studies.
The paper is organized as follows. In Sec. II we present the GCDD procedure for the specific case of a qudit. In Sec. III we give our prescription to build up the control fields necessary for the GCDD procedure. In Sec. IV we show how to apply our GCDD procedure in the case of an atomic qutrit realized with some states of 87 Rb. In Sec. V we present our numerical simulations showing the protection realized by our GCDD procedure in the case of the atomic qutrit previously presented, in the presence of a noise taking into account both damping and dephasing effects. In Sec. VI we comment on the possibility to extend our GCDD to the case of an ensemble of many qudits, identical or not, by referring to Appendix D where the case of two qudits is explicitly treated. Some parts of our analysis are reported in various Appendices.

II. THE GCDD PROCEDURE
Here, we present our GCDD procedure for the specific case of an arbitrary qudit of dimension d.
Let H G be the Hamiltonian generating the intended evolution of an arbitrary qudit state, in the ideal, noisefree case. Notice that the free Hamiltonian of the qudit can be included in H G or assumed to be eliminated before the GCDD procedure. The procedure is then directly applied to a system which, in the absence of H G , is degenerate since all qudit levels have the same energy. After a gate operation time τ , the desired evolution operator acting on an initial qudit state is then given by U G = exp(−iH G τ / ). Instead of H G , we aim to use external fields whose interaction with the qudit is described by a non-autonomous Hamiltonian H lab (t), acting continuously during the time interval τ and only on the qudit Hilbert space, and generating, despite the presence of noise, an effective evolution of the qudit that, at least up to a high enough fidelity, is at time τ the same as the one provided by U G .
To describe H lab (t), we first split this control field Hamiltonian into two terms, H lab (t) ≡ H c (t) + H gate (t), where H c (t) is the control Hamiltonian that will continuously decouple the qudit evolution from the interference of the environment, while H gate (t) will provide the modified gate Hamiltonian that in the end will effectively reproduce the action generated by H G . Associated with H c (t) there is a unitary operator U c (t) that we require to be periodic with a period t 0 and that should satisfy (but non necessarily) the dynamical-decoupling condition [54] t0 where I E is the identity operator of the environmental Hilbert space and H int is the Hamiltonian interaction term coupling the qudit with its environment. It is important to emphasize that a weaker, but sufficient, condition would also work to obtain dynamical decoupling. This weaker condition can be written aŝ where I d is the identity operator of the qudit Hilbert space of finite dimension d, c is a constant, and B is some operator that acts on the environment. Indeed, if the integral of Eq. (1) results to be proportional to a tensor product between the qudit identity and an environment operator, the system will also be protected. To explicitly illustrate this aspect, we consider in detail this situation in Appendix A, where we show that also in this case the qudit dynamics is protected from the effects of the coupling with the environment. In the following we show that our protection scheme satisfies this weaker condition. We finally observe that we choose t 0 such that τ results to be an integer multiple of t 0 itself, for reasons that we explain below. The total Hamiltonian of the system and the environment is then written as where H E is the free Hamiltonian of the environment. In the picture obtained by unitarily transforming Eq. (3) using U c (t), we obtain the Hamiltonian in what we henceforth call the control picture: where, as we explain shortly, we have chosen H gate (t) as Note that at time τ , the qudit state in the control picture coincides with the one in the original picture at the same time, explaining why we have chosen Eq. (5) and t 0 such that τ is an integer multiple of it. Since H gate (t) is used in the presence of continuous dynamical decoupling, the evolution proceeds as if only U G governed the qudit evolution in the control picture (see Appendix A for details). At time τ , even in the original picture the qudit state is then the one that the ideal noise-free evolution would produce, up to a high-enough fidelity. The dissipative dynamics is assumed to be resulting from a perturbing interaction between the qudit and its environment described by the very general Hamiltonian: where B r,s , for r, s = 0, 1, 2, . . . , d − 1 are operators that act on the environmental states and |k , for k = 0, 1, 2, . . . , d − 1, are d normalized state vectors forming a basis set for a Hilbert space of dimension d, henceforth called the qudit space. This is also to be considered the logical basis. Notice that if the free Hamiltonian has not been eliminated before the application of the GCDD procedure, it can be thought as included in the above interaction Hamiltonian in such a way that the protection procedure in the end eliminates also the action of this free Hamiltonian.

III. OUR PRESCRIPTION
Now, we prescribe how to construct the required U c (t) and we address an arbitrary quantum gate operating on the qudit state, showing how to protect its action against general noise. Also, we illustrate how to generate the fields in the laboratory whose action allows one to implement the control Hamiltonian H c (t) and the timedependent gate Hamiltonian H gate (t), the sum of which generates the evolution of the qudit state driven by such external fields.

A. Constructing Uc(t)
We begin defining H L as the Hermitian operator whose action on the logical basis states gives being ω 0 the control frequency corresponding to the dynamical-decoupling period t 0 . The quantum-Fourier transform of the logical basis is given by [2] for n = 0, 1, . . . , d − 1. We define the Hermitian operator H F by its action on the quantum-Fourier transformed basis, which is for n = 0, 1, . . . , d−1. The required control unitary transformation is then given by where we have defined a real constant ω r as Using Eqs. (7)(8)(9)(10)(11), it is now straightforward to obtain where the time integration is not zero if and only if (n − m) + (p − q)d = 0. Given the constraints on the various parameters, this happens only if m = n and p = q. It follows from Eqs. (6) and (13) that a weaker condition than Eq. (1) is obtained since the integral there appearing results to be equal to t0 This weaker condition has the same form of Eq. (2). As said above, this condition is enough to obtain dynamical decoupling.

B. The laboratory Hamiltonian
To implement the GCDD, we need a prescription for H c (t) and H gate (t). Using Eq. (11), we can calculate the control Hamiltonian, where, for simplicity, we have defined U L (t) ≡ exp(−iH L t/ ). Equation (5) gives H gate (t) in terms of U c (t) and H G . Hence, in the laboratory we need to generate external fields such that they interact with the qudit according to the following Hamiltonian: The term proportional to the unit matrix is immaterial to the dynamics, since it only gives rise to a global phase factor multiplying the evolved state vector.
The action of the gate Hamiltonian H G is what we want to effectively perform, after the time interval τ , through H lab (t). The Hamiltonian H G can be expanded in the computational basis by where g * s,r = g r,s because H G is Hermitian. If some of the eigenvalues of H G are positive, let's choose, of these, the one with the highest absolute value, say, g 0 , with g 0 0, where the equal sign is chosen if there are no positive eigenvalues of H G . Then, let's define the Hermitian operator G such that Therefore, −G does not have positive eigenvalues and, thus, it's a non-positive operator. Notice that G, however, is a non-negative operator [we have introduced the minus sign appearing in Eq. (17) just for convenience)]. Now, let us take a look at H c (t) of Eq. (14). It is easy to see, from Eqs. (7) and (10), that if we define Hermitian operators H L and H F by and respectively, then H L and H F are both non-negative operators. Now, using Eqs. (17), (18), and (19), we obtain from Eq. (15) that where we have also used Eq. (8). Since a unitary transformation of a non-negative operator is still a non-negative operator and H L , H F , and G, as we have defined them, are all non-negative, it follows that the last term within square brackets on the right-hand side of Eq. (20) is a non-negative operator. Then, we can define where and rewrite Eq. (20) as where Equation (23) gives the explicit laboratory Hamiltonian that can be used to implement the protective scheme. Also, it is important to emphasize that once the target state of the gate is achieved, thanks to our GCDD procedure, we can preserve the final state in the absence of the gate but under the noise. This is done by simply turning the protection on, without the control fields that generated the gate operation (this is obtained by using H G = 0). In this way the memory of the final state is protected.
In the next section, we apply the GCDD method to the case of a qutrit based on the three magnetic hyperfine states of the ground energy level of 87 Rb. However, in principle, we could use the same kind of twophoton atomic transitions, employed to implement this qutrit, and involve any number of Zeeman hyperfine states. Thus, although the model implementation we use is limited to the simple case of the ground state hyperfine states of the 87 Rb atom, other atomic systems could be used to obtain control over systems of qudits of dimension d = 2 or d > 3. Let us, then, proceed with the illustration of the application of the GCDD method.

IV. APPLICATION OF THE GCDD
To illustrate the application of the GCDD method, we describe a possible implementation of a qutrit, exploiting Rb (not to scale). We are using two-photon transitions for three different detunings ∆s, for s = 1, 2, 3, and, for each of these detunings, we use σ ± -and π-polarized laser light. The wavelength 780.241 nm corresponds to a frequency of the order of 384.230 THz. The qutrit space comprises the subspace spanned by the three magnetic states of the F = 1 ground level of 87 Rb, with magnetic quantum numbers mF = −1, 0, 1. We represent these degenerate states by the kets |m , for m = −1, 0, 1, respectively. the three magnetic hyperfine states of the ground energy level of 87 Rb. In the following, we use atomic data from Ref. [55]. Figure 1 shows the relevant D 2 -line hyperfine states of 87 Rb. In the following, we show that this model allows one to generate continuously all the possible SU(3) group operations which are, in general, needed to apply the GCDD procedure.
The 87 Rb atom, in the absence of external magnetic fields, has a ground-state manifold of three degenerate magnetic states. This is so because 87 Rb has a nuclear spin equal to 3/2 and a fundamental electronic manifold of states with symmetry 5 2 S 1/2 . This amounts to a hyperfine ground state with total angular momentum F = 1, so that there are three magnetic states whose projections along the quantization axis have quantum numbers m F = −1, 0, 1. These three ground states are degenerate in the absence of magnetic fields and we denote them by |m , for m = −1, 0, 1 (the notations 1 and +1 are both used in the following). The 5 2 S 1/2 ground manifold of states (including also the five magnetic states of the F = 2 ground level, besides the already-mentioned three F = 1 states) can be excited to states of the 5 2 P 3/2 excited manifold by absorbing photons with wavelengths of about 780 nm (called the D 2 spectral line of 87 Rb). The lowest-energy magnetic hyperfine state of the 5 2 P 3/2 manifold is not degenerate and has a totalangular-momentum quantum number F = 0, whose pro-jection is m F = 0. We denote this state by |e and we name ω g the energy of the ground states |m , for m = −1, 0, 1, and ω e the energy of the excited state |e .
If we use only frequencies corresponding to virtual transitions with wavelengths greater than the optical 780 nm, that is, if we use only photons with frequencies ω s that are red-detuned from the F = 1 ↔ F = 0 transition (this means that the detunings ∆ s ≡ ω s − ω e + ω g are negative), then we can approximate the relevant set of atomic states to be the one involving only the ground states |m , for m = −1, 0, 1, and the excited state |e . For the restricted Hilbert space of these four atomic hyperfine states, denoted by H 4 , we have the identity operator If the photons are detuned far enough to the red of the transitions |m ↔ |e , for m = −1, 0, 1, then the excited state is not going to be effectively populated, avoiding spurious transitions to the F = 2 ground states through spontaneous emission from |e . The effective qutrit, therefore, consists of the states |m , with m = −1, 0, 1, whose Hilbert space we denote by H 3 .
As explained in Appendix B, the control over the states in H 3 using the GCDD method is accomplished through two-photon transitions. These transitions are used to couple these states among themselves in a controlled way. In particular, in order to be able to generate continuously all the possible SU(3) group operations we need three independent detunings ∆ s , with s = 1, 2, 3, one for each transition |m ↔ |e , for m = −1, 0, 1. In particular, each laser beam is red-detuned from |e , i.e., ∆ s < 0 for s = 1, 2, 3. For each of the three different laser colors we can use the linear polarization and both circular polarizations, thus obtaining a total of nine independent Rabi frequencies. These independent control parameters are enough to effectively emulate the action of any 3 × 3 Hermitian matrix used to represent a generic modified single-qutrit quantum gate, H gate (t) [cf. Eq. (5)], together with the control fields described by H c (t), that are required for the GCDD, as explained in detail in Appendix B. There, we present all the details and values of the parameters that we can use, in principle, to implement the above atomic qutrit and the needed effective control Hamiltonian. Here, it suffices to say that using the rotating-wave approximation and adiabatic elimination of |e [56], we get an effective Hamiltonian based on two-photon interactions given by [cf. Eq. (B76)] where Ω s,−m (t), for s = 1, 2, 3 and m = −1, 0, 1, are adiabatically time-varying Rabi frequencies allowing one to emulate the time dependent control Hamiltonian of Eq. (15) up to an immaterial term proportional to I 3 [cf. Eq. (23)] (this term, which is immaterial for the dynamics, is not even needed if the atomic energies are shifted in such a way that ω g = ω l ). As quantum gate to implement for the above qutrit we choose the Hadamard one. Starting from the ideal definition of the Hadamard unitary quantum gate [57] in the ground-state subspace basis {| − 1 , |0 , |1 }, namely, we can invert the equation to obtain the gate Hamiltonian: where τ is the characteristic gate time. For example, starting from the state |0 , after the action of the Hadamard operation, the output state after a time τ becomes where ϕ = 2π/3. In the next section, in our numerical simulations, we consider two paradigmatic noises due to bosonic thermal baths, chosen to disturb our intended gate operation. We consider the amplitude damping noise which simulates dissipation involving, respectively, | − 1 and |0 , and |1 and |0 , and the dephasing noise which destroys the relative coherences among these states.

V. NUMERICAL SIMULATIONS
In this section we present our numerical simulations in the case of the qutrit depicted in Fig. 1 when both damping and dephasing are simultaneously present. In particular, we perform our simulations by following the prescription of Sec. III and we assume the presence of two identical baths with Ohmic spectral density, characterized by an exponential cutoff function with an angular cutoff frequency ω c = 2π/τ c , where τ c is the bath correlation time. The effectiveness of our protective scheme is studied by means of the fidelity measuring how close are the states obtained, respectively, by the dissipative dynamics induced by Eq. (3) (treated by means of a Redfield master equation as explained below) and by the noisefree dynamics governed by the Hamiltonian H lab (t) of Eq. (15). We recall that this latter dynamics gives after a time τ the same output state of the ideal noise-free , with |0 as initial state. Here, the bath correlation time is τc = τ /4, ωc/(kBT ) = 1, and the involved coupling constants are all equal leading to an effective coupling parameterλ = 0.1 (see Appendix C for details). The solid line represents the fidelity with no protection, while the dotted, dot-dashed, and dashed ones refer to the protective scheme with ω0/(2π) equal to 2/τc, 4/τc, and 16/τc, respectively. In the insets, we represent the gate fidelity (fidelity at time τ ) as a function of n = ω0τc/(2π) (the interpolated curve just guides the reading) and we report its numerical values.
dynamics governed by the Hadamard gate Hamiltonian H G of Eq. (29) (i.e., the case without protective scheme). The dissipative dynamics in the absence of the protective scheme is obtained by turning off the control Hamiltonian H c (t). The fidelity is defined for two arbitrary states ρ and σ as [Tr{ Here, the open quantum system dynamics is obtained by means of the Redfield master equation which describes the dissipative dynamics in the general case. In the Appendix C we show in details how to calculate the dynamics governed by this master equation. Figure 2 shows that if we do not use the GCDD method during the time τ in which the Hadamard gate operates and let the noise affect the dynamics, starting from the state |0 , the fidelity rapidly decreases. When the GCDD protection is turned on, we obtain better results by increasing the control frequency ω 0 , with the final gate fidelity moving towards one. In the insets, we show the gate fidelity, i.e., the fidelity at time τ , as a function of n = ω 0 τ c /(2π) [this implies t 0 = τ c /n, according to Eq. (8)] and we report its numerical values. In particular, the smaller is t 0 with respect to τ c , the more effective is the decoupling procedure. We stress out that, by construction, the time at which to look for a state close to the original target is exactly the time τ at which the original gate would have produced that state in the absence of the environment and of the control fields. The actual value of the gate time τ is not specified in these simulations, the other quantities being given in units of it. Its value must just be such that the derivation of the effective Hamilto-nian of Eq. (26) can be coherently performed. We state this condition for a question of coherence in our analysis, even if in our simulations we do not make use of Eq. (26), but we directly follow the prescriptions given in Sec. III. We also observe that the results shown in Fig. 2 have been obtained when the various coupling constants involved in the interaction Hamiltonian with the environment are all equal leading to an effective coupling parameterλ = 0.1 (see Appendix C for details). We have also tested some configurations with the various coupling constants not all equal, finding similar results.
The illustration of the GCDD method shown in Fig. 2 for our qutrit model using 87 Rb and a modulated set of laser beams, can, in principle, be realistically implemented in the laboratory. Even if only the qutrit case has been considered, our results show that quantum computation could be implemented using laser light and atomic systems, which are available in setups with trapped ions, for example. This kind of implementation is attractive because it already presents long coherence times, implying high efficiency of our procedure [56,58,59]. We remark that the implementation of our procedure can, in principle, be extended to the case of a qudit with more levels.

VI. EXTENSION TO THE CASE OF A MULTI-QUDIT SCENARIO
The domain of applicability of our GCDD procedure can be extended to the case one wants to protect the action of a multi-qudit gate on an ensemble of qudits, identical or not, subject to general noise, which can act locally or non-locally on them. This can be obtained by extending the procedure presented in detail in Appendix D for the case of two qudits. There, it is explicitly shown how to build up the control fields necessary to protect the qudits from the action of a general noise. The extension to the case of more than two qudits, identical or not, is straightforward, as indicated in Appendix D.
A simple but important direct application of this procedure is the possibility to preserve the memory of a multi-qudit entangled state against noise.

VII. CONCLUSIONS
In conclusion, here we have presented a generalized continuous dynamical decoupling procedure to decouple an ensemble of qudits from any possible noise and still apply an arbitrary many-qudit quantum gate on them. Our study extends the domain of applicability of dynamical decoupling strategies. Indeed, we provide a general procedure of this kind to protect the action of a general quantum gate on one or many qudits, identical or not, against general noise.
Importantly, we have explicitly provided a detailed analysis of a specific example, employing a Rubidium atom, which, in principle, could be experimentally implemented, where it is explicitly shown how to realize the operations which could be in general needed to apply our GCDD procedure.
We think then that our approach represents an important step towards the protection of quantum information, especially when many-level quantum systems are employed. We stress that our procedure appears to be directly implementable in experiments with atomic systems, nitrogen-vacancy centers, or other setups where current technology permits to generate the control fields required for the protection scheme. Here, we show that dynamical decoupling is still obtained if a weaker condition than the one of Eq. (1) is satisfied.
Let us start observing that Eq. (1) is stronger than necessary to achieve dynamical decoupling. A weaker, but sufficient, condition would be that, instead of being equal to zero, the integral in Eq. (1) results in a tensor product between the qudit identity and an environment operator. To see this, we recall that the idea of dynamical decoupling comes from the usual Magnus expansion [54] of the total propagator for the whole system, the qudit and the environment, that is, where we have defined the time average of H(t) as and H(t) is the total Hamiltonian given in Eq. (4).
Hence, we obtain: Equation (A1) is not exact because, as usual [54], we have neglected terms of order equal and superior to t 0 in the Hamiltonian part appearing in the argument of the exponential. However, it becomes exact if the number of periods within the time interval τ tends to infinity. We now assume that a weaker condition than the one of Eq. (1) is satisfied, namely [the same condition given in Eq.
where c is a constant and B is some operator that acts on the environment. It follows that Eq. (A3) becomes: Therefore, with this average Hamiltonian, Eq. (A1) gives which shows that the qudit is decoupled from the environment, since the interactions get effectively eliminated, even in the case in which c = 0 in Eq. (A4).

Appendix B: GCDD for an atomic qutrit
In this Appendix we show how to effectively implement an atomic qutrit with the system described in Sec. IV and how to realize within this system an effective control Hamiltonian like the one of Eq. (15) (in general up to an immaterial term proportional to I 3 ), needed to apply our GCDD procedure.

The interaction Hamiltonian between the atom and the laser beams
We introduce nine laser beams, whose electric-field vectors, each being the resultant with a different polarization, can be written as and (B2) where the polarization versors are chosen, in terms of a space-fixed system of Cartesian coordinates, aŝ representing, respectively, the σ ± polarizations, and representing the π polarization. Here, the z-axis of this system is chosen to represent the quantization axis. It is noteworthy that in Eqs. (B1) and (B2), for each polarization, there are three different superposed amplitudes, E s,± (t) and E s,0 (t), each corresponding to a different polarization-independent frequency, ω s , for s = 1, 2, 3. Figure 1 shows the scheme we are describing. The amplitudes E s,± (t) and E s,0 (t), as we discuss below, must follow a prescribed relatively slow time-dependent modulation. It is worth mentioning that we treat the driving electric fields of Eqs. (B1) and (B2) as classical, intense laser fields. We are justified to use such a semiclassical approach because of the relatively high intensities and detuning magnitudes used, so that quantum fluctuations of the number of photons is completely negligible in the regime we consider here. The laser beams of Eqs. (B1) and (B2) interact with the atom according to the Hamiltonian since we have the three resultant laser fields continuous and simultaneously present, each one with a different polarization, where d is the atomic electric-dipole operator, which is Hermitian. In Cartesian coordinates, we write and, using Eq. (25) in d = I 4 dI 4 , we obtain where we have used the fact that the electronic excited state has a parity that is opposite to the parity of the ground states, that is, for m, m = −1, 0, 1. Now, we can write the operator d in terms of its spherical-tensor components: that is, where we have used Eqs. (B3) and (B4) and defined its spherical components as usual: Because |e has zero angular momentum, from Eq. (B10) it follows that since total angular momentum is conserved. From the Wigner-Eckart theorem [60], we have where D is a reduced matrix element of the dipole operator and is, thus, independent of m or q. Hence, we rewrite Eq. (B13), using Eq. (B14), as where we have used the rotating-wave approximation [61], which is justified since we will choose detunings ∆ s = ω s − ω e + ω g , for s = 1, 2, 3, and Rabi frequencies Ω s,q (t), defined as for s = 1, 2, 3 and q = −1, 0, +1, much smaller in modulus than the transition frequency ω e − ω g and the laser frequencies ω s . We have also used Eqs. (B3) and (B4) to calculate the scalar products between polarization vectors. Values of |Ω s,q (t)|/2π of the order of a few MHz, let us say, roughly are routinely obtained in the context of optical manipulation of rubidium [58,59]. These independent control parameters are enough to emulate the action of any 3 × 3 Hermitian matrix used to represent a generic modified single-qutrit quantum gate [cf. Eq. (5)] together with the control fields required for the continuous dynamical decoupling, as we explain below. Given Eq. (B18), we can rewrite Eq. (B17) as This is the interaction Hamiltonian whose effective version, for large detunings to the red of the D 2 line, allows us to realize in the laboratory the GCDD Hamiltonian of Eq. (23), in general up to an immaterial term proportional to the identity operator in the qutrit Hilbert space. In the following we show how to do this through adiabatic elimination of the excited state |e .

Effective implementation of the GCDD Hamiltonian for the atomic qutrit
The unperturbed atomic Hamiltonian is written as where the energy ω g of the ground states |m , for m = −1, 0, 1, and the energy ω e of the excited state |e are such that (ω e − ω g ) is equal to the energy corresponding to the D 2 line, with wavelength given by 780.241 nm, which corresponds to a frequency of the order of 384.230 THz: Using Eq. (B21) as the unperturbed Hamiltonian in the usual interaction picture with the interaction Hamilto-nian of Eq. (B20), we have where we recall that the detunings are defined by The coherence times involved in superpositions of atomic quantum states are typically of the order of a second or longer [56,58,59]. Thus, because the quantum-gate operation τ is an integer multiple of t 0 and should be shorter than these typical coherence times, we can take, roughly, Hence, because we take Eq. (B25) as valid, we see that is the corresponding rough estimate we can take for ω 0 [cf. Eq. (8)]. As we show below, about ten Hz for ω 0 /(2π) are enough for the GCDD method to work. In other words, about ten Hz corresponds to the order of magnitude of the Hamiltonians we need to emulate the H L /(2π ) and H F /(2π ) operators of Eqs. (18) and (19) [see also Eqs. (7), (8), and (10)].
To implement the GCDD method in the context of laser control of an atomic qutrit, here we show how to use two-photon transitions. We need detunings that are much greater in magnitude than the typical few MHz of the Rabi frequencies in modulus [cf. Eq. (B19)], so that we can make the adiabatic elimination of the state |e [56]. One can easily implement this, given the relatively large difference in energies in the transitions indicated in Fig. 1. As we can see in this figure, the detunings can be as large as a few GHz, and still the states of the ground F = 2 level do not get involved in the transitions (they are at the energy corresponding to 6.8 GHz above the F = 1 states we use). Since the detunings we use are negative, meaning that the photons excite a virtual level well below the |e excited state, the higher excited states are not going to interfere with our transition scheme. Moreover, the D 2 -line natural line-width for 87 Rb is of the order of 6 MHz [55], so that laser photons detuned to the red of the D 2 transition frequency by a few GHz will not practically populate the excited state |e . As we show in detail below, the magnitudes involved in the effective two-photon Hamiltonian are proportional to the square of Rabi frequencies divided by the detuning, which can be substantially higher than the few tens of MHz (at least) required for an efficient GCDD implementation. Typically, if we use, roughly, using Eq. (B19) we find that Hence, using large detunings as in Eq. (B27), we end up with an effective Hamiltonian (explained below) that can have its magnitude as in Eq. (B28),flexibly above the minimal requirement of Eq. (B26) for the GCDD method to work, as we have discussed above. We notice that the values for |∆ s | associated to Eq. (B27) are consistent with the rotating wave-approximation used to obtain Eq. (B17) since they are much smaller than the transition frequency ω e − ω g and the laser frequencies ω s . Now, let us write the interaction-picture state as since |ψ I (t) ∈ H 4 . We can introduce the following projection operators: and P e ≡ |e e| .
We immediately see that |ψ I (t) = (P g + P e ) |ψ I (t) = P g |ψ I (t) + P e |ψ I (t) .  (B39) Our intention is to start with the atom in the groundstate subspace, that is, the population of the excited state is initially zero. Thus, using this fact, that is, in Eq. (B39), we obtain where we have used the fact that P e is a projector operator and, therefore, P 2 e = P e . From Eqs. (B23), (B29), (B30), and (B31), we see that  for m = −1, 0, 1, where we have defined the kernel function: We can also arrange Eqs. (B45) and (B46) in matrix for-mat: where we have defined and Iteration of Eq. (B47) yields: (B50) Let us calculate a generic element of the first kernel integral in Eq. (B50): where we have used Eq. (B46). If we initially focus on the integral over t 2 in Eq. (B51), we must make an assumption about the time dependence of Ω s ,q (t 2 ). Our aim here is to use Eq. (B23) to effectively emulate the GCDD Hamiltonian of Eq. (23). As we show in the following, although we are going to assume our Rabi frequencies Ω s ,q (t 2 )/2π with magnitudes of a few MHz [see Eq. (B19)], its time dependence is to be modulated in such a way that the corresponding spectral density results to be centered at about ω 0 /2π, whose value is here assumed to be of the order of 10 Hz [see Eq. (B26)], with a width much smaller than |∆ s |. Let G s ,q (ω) be the Fourier transform of Ω s ,q (t ), namely, whose inverse is given by Using Eq. (B53), we obtain where we have changed the order of the integrals. Whatever form G s ,q (ω) might have, it is assumed to be characterized by a central value ω 0 > 0 and a width both much smaller than |∆ s |, so that we can writê where we have used Eq. (B53). With the result of Eq. (B56) we can now tackle Eq. (B51): We already see that all the terms in the second double sum on the right-hand side of Eq. (B57) are of second order in the quotient between the order of magnitude of the Rabi frequencies [see Eq. (B19)] and the order of magnitude of the detunings. Let us define this order of magnitude more rigorously by assuming that, for all t ∈ [0, t 0 ], s = 1, 2, 3, and q = −1, 0, 1, we define η as the maximum absolute value of Ω s,q (t)/∆ s , that is, Using the rough estimates of Eqs. (B19) and (B27) we see that η can be even less than 10 −3 . Using Eq. (B53), we can now calculate the following integral: Here we have two situations: s = s and s = s . Hence, taking these two cases into account, we obtain We now substitute Eq. (B60) back into Eq. (B59). We have assumed that the functions G * s,−m (ω 2 ) and G s ,q (ω 1 ) only contribute in a frequency region around ω 0 much smaller than |∆ s |. Thus, if we choose the detunings such that, for s = s , the absolute difference |∆ s − ∆ s | is of the same order of magnitude of the max{|∆ s |} s∈{1,2,3} , that is, which we estimate as about a few GHz [see Eq. (B27)], then we can assume ∆ s − ∆ s + ω 2 − ω 1 ≈ ∆ s − ∆ s in the denominators of Eq. (B60) when this equation is substituted back into Eq. (B59), obtaininĝ After substituting Eq. (B62) into Eq. (B57) we end up witĥ t t dt 1ˆt We conclude, therefore, that defining the kernel matrix including terms with s = s , as in Eq. (B49), amounts to producing contributions of second order in η that are negligible when compared with the terms with s = s , which are of first order in η, as we can directly verify in Eq. (B63) [cf. Eqs. (B58) and (B61)]. Therefore, in the present problem, we keep only the s = s terms in Eq. (B63) and neglect any other terms of second order in η:ˆt Now we can differentiate Eq. (B64) with respect to t and (B65) We see from the above discussion and Eqs. (B47) and (B50) that a time-local approximation of Eq. (B45) is of first order in η [see Eq. (B58)], and, using Eq. (B65), we can write it as Equation (B66) can be arranged in a matrix representation: where, for s = 1, 2, 3, we define We see, in Eq. (B67), that, effectively, we have found a Hamiltonian in the interaction picture given by 3 s=1 H I,s (t). Now, where |ψ S (t) is in the Schrödinger picture [see Eq. (B23) for the definition of the interaction picture in our problem]. Using Eqs. (B21) and (B30), we see that since in this case P g = I 3 , the identity operator acting on H 3 [cf. Eq. (B30)]. The above equation simply means that once computed the evolution of the state C(t), the corresponding state in Schrödinger picture state is obtained by multiplying it for the immaterial global phase factor exp(−iω g t). It follows that the dynamics in the qutrit subspace H 3 is governed by the effective Hamiltonian To make this scheme work, we know that ∆ s ∈ R and ∆ s < 0. Therefore, we adopt the convention that so that where √ −∆ s ∈ R and √ −∆ s > 0. Thus, Eq. (B68) becomes Ωs,0(t) √ ∆s Ωs,−1(t) √ ∆s , (B74) since, from Eq. (B73), it follows that Based on Eqs. (B71) and (B74) we can now express H eff (t) in a way that is analogous to Eq. (23), allowing us to connect the Rabi frequencies of this Appendix with the elements of the operator Υ of Sec. III B: Ωs,0(t) √ ∆s where we have defined Then, if we choose the Rabi frequencies and detunings appearing in Eq. (B77) so that Θ(t) is Hermitian, we can identify it with Υ(t) of Eq. (23) and this is how we can implement, up to an immaterial term proportional to I 3 , the GCDD method for an atomic qutrit manipulated using two-photon transitions. We notice that the immaterial terms proportional to I 3 in Eqs. (23) and (B76) can be made equal if the atomic energy scale can be modified in such a way that ω g = ω l in the new scale [see Eq. (24) for the definition of ω l ]. Accordingly, thus, we choose, along the diagonal of Eq. (B77): Using Eqs. (B72) and (B73), we then obtain: For the off-diagonal elements of Eq. (B77), we choose: Now, we are able to identify the remaining independent elements of Θ(t) with those of Υ(t). By imposing Eqs. (B79) and (B80), and that Θ(t) = Υ(t), we obtain: , , , , , We observe that the above derivation leading to the effective Hamiltonian of Eq. (B76) could be coherently obtained also for values of t 0 different from the one choosen in Eq. (B25) but satisfying the conditions required for performing the various approximations involved in the derivation. In this sense, we do not fix a specific value for t 0 in the numerical simulations based on Appendix C and depicted in Fig. 2. Consequently, the value of the gate time τ is not specified in these simulations and the other quantities are given in units of it.
Appendix C: Environmental noise due to bosonic thermal baths Here, we explain how we use two baths of thermal bosons to simulate the perturbations caused by the noisy environment considered in Sec. V (see also some previous works of some of us on continuous dynamical decoupling of qubit systems [18][19][20][21]).
In the picture obtained by unitarily transforming the total Hamiltonian of the system and the environment of Eq. (3) using U c (t), we obtain the total Hamiltonian in the control picture given in Eq. (4): As explained in Sec. V, in our example we consider a qutrit subject to independent amplitude damping and dephasing noises. We divide the interaction Hamiltonian and the environment free Hamiltonian in two parts, i.e., int and H E = H E , where the superscripts 1 and 2 refer, respectively, to the amplitude damping bosonic bath and to the dephasing one. Here, we suppose that the two baths are identical, besides being independent. The above Hamiltonian terms are given in terms of the usual lowering and raising operators a (i) k and a (i) † k for each mode k of the i-th bath, with i = 1, 2. In particular, the first interaction-Hamiltonian term (which introduces the damping noise) is given by where B (1) = k g k a k . In a similar way, the second interaction-Hamiltonian term (which introduces the dephasing noise) is given by where B (2) = k g k a (2) k and Λ (2) = λ , while the bath Hamiltonian associated to this class of error is given by H To obtain a master equation governing the three-level system dynamics, we transform the total Hamiltonian to the interaction picture. It is written as wherẽ H (s) . With this transformation, the Redfield master equation is written as [62] and ρ E is the environment density matrix, here given by a thermal state, that is , k B is the Boltzmann constant, and T is the temperature of the baths. We remark that the above Redfield master equation is derived under the Born approximation, linked to the assumption of weak coupling between the system and the environment, and to a part of the global Markovian approximation, so that the master equation is not Markovian [62]. Now, substitutingH I (t) into the master equation, we finally obtain where the correlation functions G 1 (t, t ) and G 2 (t, t ) are given by It is important to emphasize that, since we suppose identical baths, the correlation functions are equivalent for B (1) andB (2) . Thus, the expressions for G 1 (t, t ) and G 2 (t, t ) are given by is the average number of photons in a mode with frequency ω k . Finally, in the continuum limit, the sums become integrals and we obtain, using s = t − t , where we have exploited the fact that the correlation functions are homogeneous in time, n(ω) is the continuous frequency version of n k , namely, n(ω) = 1/[exp(β ω) − 1], and J(ω) is the spectral density.
In the numerical simulations of Fig. 2, we apply the GCDD to the qutrit considered in Sec. IV and we choose for the spectral density J(ω) = α 2 ω exp(−ω/ω c ), where α is a dimensionless constant prefactor and ω c is the angular cut-off frequency. We have also set equal all the λ coupling constants, introducing as effective coupling parameterλ = αλ is, they have dimensions d 1 and d 2 , both integers, which may be different. Without loss of generality, we assume d 2 d 1 > 1. As for the one-qudit case, the free Hamiltonians of the qudits can be included in H G or assumed to be eliminated before the GCDD procedure.
The starting point is then a generalization of Eq. (3), which reads where, now, H gate (t) corresponds to the action of a twoqudit gate H G , I d1d2 ≡ I d1 ⊗ I d2 , and the control Hamiltonian comprises only local one-qudit operations, that is, H c (t) = H and H (2) c (t) in general are different, as described below. As for the noise, we consider a generalized interaction Hamiltonian as: where B p,q,r,s are operators that act on the environmental states. Let us notice that this Hamiltonian may contain, in general, both one-qudit perturbations (local noise) as well as two-qudit perturbations (non-local noise). Then, we have B p,q,r,s = C p,q δ r,s + D r,s δ p,q + E p,q,r,s , so that where the terms involving the operators C p,q and D r,s describe local noises, while the terms involving the operators E p,q,r,s describe non-local noise. As for the case of one qudit treated in Sec. II, the condition there expressed in Eq. (1) is not going to be satisfied in our procedure, since the integral there appearing is going to be proportional to the identity of the two qudits, tensorially multiplied by an operator that only acts on the environment. Analogously to what has been done in Appendix A, it is possible to show that this result is sufficient to dynamically decouple the qudits from the noisy perturbations. In the end, the final state is the same that, up to first order of perturbation, one would obtain in the absence of interaction with the environment. In other words, Eq. (1), applied to the two-qudit case treated here, is sufficient, but not necessary, since if the integral is proportional to the qudits-system identity, instead of zero, this is sufficient for guaranteeing dynamical decoupling. This is due to the fact that the integral appearing in Eq. (1) is present in the exponential in the Magnus expansion of the total propagator at first order, so that, if it is proportional to the system's identity, dynamical decoupling is still obtained (see Appendix A).
It is then enough to obtain for the integral of Eq. (1) something proportional to the two-qudit identity. To this aim, we use control Hamiltonians H (D4) where j = 1, 2, runs over the two qudits, for k j = 0, 1, . . . , d j − 1, for n j = 0, 1, . . . , d j − 1, and The main point is that, differently from the one-qudit case where we have used ω 0 , here we make use of two different frequencies, ω (1) 0 and ω (2) 0 . Then, we define the control propagator for the two qudits as In order to show that the integral of Eq. (1) in the case considered here of two qudits is proportional to the twoqudit identity, we consider the integral for a generic term in the decomposition of H int in Eq. (D2) and we exploit the fact that the frequencies we have chosen amount to different time scales. Proceeding analogously to what has been done in Sec. III A, we obtain in the two-qudit version of Eq. (13)  1) is here proportional to the two-qudit identity: We emphasize again that this result is enough to eliminate any noise involved in Eq. (D2), even when the environment involves non-local terms acting on the qudits, and that this result is only possible, in general, if we can really separate all the time scales of each control propagator, so that the final results of Eq. (D10) can be given in terms of factored Kronecker deltas. The choice made in Eq. (D6) is motivated by this aim. In the case of a single qudit, treated in Secs. II and III, we need only two sets of time scales, determined, respectively, by {kω 0 , k = 0, 1, . . . , d − 1} and {kdω 0 , k = 0, 1, . . . , d − 1}. Now, having two qudits, we need four different sets of time scales in order to have in Eq. (D10) the time integral giving rise to several factored Kronecker deltas. According to our choice, we have {kω 0 , k = 0, 1, . . . , d 1 − 1} and {kd 1 ω 0 , k = 0, 1, . . . , d 1 − 1} for one of the qudits, while we use {kd 2 2 ω 0 , k = 0, 1, . . . , d 2 − 1} and {kd 3 2 ω 0 , k = 0, 1, . . . , d 2 − 1} for the other one. We finally stress that when d 1 = d 2 all that has been said above is still valid and we can, therefore, protect two identical qudits using the same procedure.
The problem simplifies in the case of only local noises, that is, when in Eq. (D3) the terms with E p,q,r,s are missing and then H int reduces to In this case, we can define each U It follows that the integral in Eq. (1) is still proportional to the two-qudit identity: We notice that in the above simplified case with only local noises, U c (t) in the case of identical qudits, since in this case d 1 = d 2 .
Now, we address the problem of the two-qudit gate, irrespective whether the qudits are identical or not. Let us start noticing that since the control Hamiltonians are local operators acting on each qudit, they do not change any possible entanglement between the two qudits. We proceed in an analogous way as we did for the case of one qudit, that is, we use where now, of course, H G is the desired two-qudit gate and U c (t) is the two-qudit control propagator defined in Eq. (D9). It might be very difficult to be able to implement such an orchestrated two-qudit Hamiltonian, but, once it is done, we have a two-qudit gate protected against a very general kind of perturbing environment by using H lab (t) = H c (t) + H gate (t), where H c (t) is given by where U (j) A simple and direct application of the procedure explained in this Appendix is represented by the case when one does not want to protect the action of a gate on the qudits but only preserve their state since, for instance, it is an entangled state. In this case, the state of the two qudits can be preserved from the action of the noise, as a protected memory state, by simply considering the case H G = 0, that is by using H lab (t) = H c (t), where H c (t) is given in Eq. (D19).