Quantum control operations with fuzzy evolution trajectories based on polyharmonic magnetic fields

We explore a class of quantum control operations based on a wide family of harmonic magnetic fields that vary softly in time. Depending on the magnetic field amplitudes taking part, these control operations can produce either squeezing or loop (orbit) effects, and even parametric resonances, on the canonical variables. For these purposes we focus our attention on the evolution of observables whose dynamical picture is ascribed to a quadratic Hamiltonian that depends explicitly on time. In the first part of this work we survey such operations in terms of biharmonic magnetic fields. The dynamical analysis is simplified using a stability diagram in the amplitude space, where the points of each region will characterise a specific control operation. We discuss how the evolution loop effects are formed by fuzzy (non-commutative) trajectories that can be closed or open, in the latter case, even hiding some features that can be used to manipulate the operational time. In the second part, we generalise the case of biharmonic fields and translate the discussion to the case of polyharmonic fields. Using elementary properties of the Toeplitz matrices, we can derive exact solutions of the problem in a symmetric evolution interval, leading to the temporal profile of those magnetic fields suitable to achieve specific control operations. Some of the resulting fuzzy orbits can be destroyed by the influence of external forces, while others simply remain stable.

A typical quantum control programme consists of a setup of oscillating external fields irradiating a micro-object (qubit) confined in a very small cavity either to induce a class of driven motion or to manoeuvre its degrees of freedom toward a special configuration. One shall be particularly cautious with some other details that cannot be disregarded, in that the settled presence of environmental noise, the radiative pollution or the emergence of electric shocks, amongst other issues, could result in an inoperative control protocol. Then, how to achieve the desired effects, at least averagely, in a simple but realisable way?
Customarily, some of these drawbacks can be successfully circumvented if a reasonable approximation is considered. There are control techniques based on magnetic nuclear resonance 1,2 , that accomplish acceptable results even if one ignores the intrinsic inhomogeneities as well as any electrical component conveyed by the sequence of variable magnetic fields, although beforehand the approximation relies on the quality of the generated magnetic pulses. Other schemes addressed to systems with discrete spectrum 3 can suppress unwanted effects if the pulses of electrical forces depend only on the temporal domain. It is also the case for the effective description of an atom that moves across an optical device and interacts with a photon, which can be represented as an oscillating electric field 4,5 as long as the dimensions of the cavity are small enough such that the approximation is not drastically far from the accurate description. Nonetheless, whenever the particle is not strictly localised, another class of control techniques shall be considered. Think about the motion of an electron, in an unbounded state, that is controlled with a protocol of simple Rabi rotations. Hence, the implementation of ample cavities would be necessary due to the electron's free propagation. These operations are usually implemented via ion traps [6][7][8][9] or optical devices 10,11 , and receive special attention in a variety of applications where the unbounded motion shall be guided in a particular way, e.g. particle accelerators.
The interest in this kind of problems has motivated a number of works not limited to ion traps, but also in the realm of squeezed states 12 -useful to produce efficient non-demolition measurements 13,14 . Although the method is, once again, based on an approximation. It basically relies on the application of harmonic oscillator potentials Scientific Reports | (2020) 10:22256 | https://doi.org/10.1038/s41598-020-79309-8 www.nature.com/scientificreports/ with time-dependent frequencies, producing an effective oscillator field inside the walls of a hyperbolic structure, which is in turn connected to a pulsating electric potential 6,15 . Our discussion below is addressed to these types of control problems, focusing our attention on time-dependent magnetic fields that are modulated in intensity by way of control currents. We shall design magnetic field pulses whose temporal profile is such that a charged particle under their influence evolves in time with either a stable (ion traps) or unstable (squeezing) motion, even neither of them, which yields a particular operation that produces parametric resonance.
Whenever the magnetic pulse is smooth enough, we will be allowed to look for exact solutions of the evolution problem. However, the application of soft operations will not be absent from very tiny discontinuities in the evolution process 16 due to the concatenation of pulses. Indeed the presence of control currents produces a circular electric field enclosing the magnetic field, which could yield a sudden electrical jump. In our proposal, such discontinuities will be partially avoided by observing that the unitary matrices of evolution satisfy a simple Toeplitz algebra.
Even more, these unitary matrices of evolution truly reproduces the evolution loop phenomenon within the structure of fuzzy orbits 17,18 , which makes the particle to return to its initial values after a certain period of time (i.e. the interval of operation), albeit according to the type of motion taking part the fuzzy orbits can be open or broken. From our viewpoint, these orbits have a semiclassical character in that we are to consider only periodic, quadratic Hamiltonians, yet the loop phenomenon in itself also arises in problems with non-quadratic Hamiltonians 19,20 . This interesting aspect entices a variety of quantum control applications 16,21,22 as well as some theoretical aspects in quantum tomography [23][24][25] .
In agreement with causality, the change in the field intensity due to time-varying currents will not be registered instantaneously in every part of the laboratory, instead the effects will be sensed depending on the dimensions of the laboratory and on the velocity of the signals. In that regard, an accurate treatment of the problem should include the retarded field effects, nevertheless, for simplicity we have restricted our analysis to the nonrelativistic regime neglecting the amount of time in which the fields propagate. Our simplification deserves special attention for a practical implementation, and we shall survey the consequences as well.
For a self-contained discussion, let us start summarising the known facts 26 . Consider a spinless particle of charge e and fixed mass m moving in a uniform magnetic field B(x) ∈ R 3 , generated by a class of vector potential A(x) through the relation B(x) = ∇ × A(x) . In particular, if the magnetic field is homogeneous a possible choice of A(x) is where B(x) will be aligned in the direction Oz, for the sake of convenience. As well, it is worthy of note how this gauge condition does not define uniquely A(x) when B(x) is given.
Let us start with a time-independent, quadratic Hamiltonian of the form: where V = 1 m P − e c A(X) is the operator that accounts for the velocity of the particle and P is the operator of momentum. At this level A(X) is also an operator and depends on the observable of position X . Besides, we see that the observables X and P obey the usual commutation relations [X j , P k ] = i δ jk , [X j , X k ] = 0 and [P j , P k ] = 0 . Whereas the components of the operator V satisfy [V j , V k ] = ie m 2 c ǫ jkl B l (X) , where ǫ jkl is 1 (-1) if j, k, l is and even (odd) permutation of 1,2,3, otherwise it is zero. Likewise, the commutation relations between the components of X and V are effortlessly obtained given that A(X) commutes with X j , it follows that [X j , V k ] = i m δ jk . There are special features that shall not escape attention. Unlike the classical problem where the trajectory described by the particle is fully localisable, in the quantum analogue a type of abstract space localisations will enter into the picture. The quantum particle will describe a trajectory whose centre (X,Ȳ ) is a fuzzy point, i.e. it is not fully localisable given that the operators associated with its position do not commute [X,Ȳ ] � = 0 -even either commuting or not there is nothing at that point, despite its geometric origin 27,28 . For instance, take the gauge (1), then the Hamiltonian will separate into two components H = H ⊥ + H � , the transversal motion (the one of our interest) will occur in the xOy plane whereas the parallel motion will be aligned to Oz. The components of the operator V are: where ω c = e mc B is the cyclotron frequency and B = |B(X)| . The respective rotation centre results in a fuzzy point with coordinates: H can be interpreted as a fuzzy surface spanned by the transversal orbits. One can conclude that the www.nature.com/scientificreports/ surface πρ 2 is a constant of the motion, yet it cannot change continuously 29 since it is quantised. This fact has a number of applications in the framework of loop quantum gravity 30,31 , for instance. Indeed the fuzzy centres are not limited to appear in quantum Hall-effect models 32 under the influence of a fixed magnetic field, but still occur in time-dependent problems. In our case, we will discuss how the timevarying magnetic operations become responsible for a number of noncommutative aspects supported on closed or open semiclassical trajectories.
In turn, we are to define the time-dependent version of the Hamiltonian (2), actually it can be written in a similar fashion but considering a magnitude-controlled current I(t) that modulates A(X) and B(X) , in this way their corresponding time-dependent versions are written as A(t, X) = I(t)A(X) and B(t, X) = I(t)B(X) . Such definitions imply that the time-dependent, quadratic Hamiltonian describing the motion of the non-relativistic, spinless, charged particle reads as 33 : Given that H leads to the well known free particle eigenfunction 1 √ 2π e ip z z/ , we are to separate this part from our analysis, and only pay attention to H ⊥ (t) = H osc (t) + H rot (t) , where the oscillatory term H osc (t) is truly equivalent to a two dimensional harmonic oscillator, whereas the easily integrable term H rot (t) represents the rotations caused by L z -the component of the angular momentum projected onto Oz-which is a conserved quantity. The former picture resembles, for example, a cylinder whose symmetry axis is parallel to Oz and carries a homogeneous current density on its surface.
To proceed it is convenient to express the perpendicular terms in (Eq. 5) in dimensionless variables, in order, we symbolise with T = 2π ω the time scale or period of operation. Our new variables become: in consequence the intensity of the oscillatory field is written as: The function β(t) is, in some extent, the cornerstone of our study: We shall implement a continuous, real function that solves the evolution problem. One is much in need to be specially careful with this task, for if β(t) is bounded and piecewise continuous, the evolution of X and P will be continuous as well, although the same is not true for the components of the kinetic momentum provided there are sudden jumps that can be interpreted as electric shocks 34,35 E = − e c ∂ ∂t A(t, X). The substitution of Eqs. (6) and (7) into the perpendicular terms of the Hamiltonian (Eq. 5) yields the simplified, dimensionless expression: in the special case that B(t) oscillates periodically (a regular Floquet problem) with frequency ω , the substitution H ⊥ (tω) → H ⊥ (t) ω corresponds also to a dimensionless representation. One must be aware that in such case the stability regions of the motion described by the older Hamiltonian will no be valid for the dimensionless Hamiltonian any more.
Even though our Report is wholly concerned with magnetic operations, we would like to remark that despite the elementary and minimalistic structure of the Hamiltonian (Eq. 8), it can also describe the evolution of charged particles in hyperbolically shaped Paul's traps 6 of radius r 0 . Under such circumstances, the time-dependent electric potentials of the type �(t, X) = eφ(t) . One immediately notes that the evolution problem can be handled in a similar fashion to the magnetic problem-in both cases yielding evolution matrices u(t, t 0 ) identical for classical and quantum dynamics.
Below we are to solve the evolution problem described by Eq. (8). The already studied biharmonic amplitudes 9,34 β(t) occupy the first part of our survey, but we shall reformulate the stability map where the magnetic control operations are comprised. In the second part, we shall generalise the biharmonic approach in terms of polyharmonic amplitudes β(t) , that at the same time play the role of exact solutions. For this case we cannot generate a stability map such as the Strutt diagram that governs the harmonic and biharmonic cases. The systematic classification of fuzzy orbits in the polyharmonic case remains open.

Results
We start with the operator equations satisfied by the unitary evolution operators U(t, t 0 ): www.nature.com/scientificreports/ nonetheless, for convenience and without loss of generality, we chose deliberately t 0 = 0 in order to introduce the abbreviations u(t) = u(t, 0) and U(t) = U(t, 0). As well let us represent with the column-vector Q of dimension 2N the set of ordered pairs of observables X and P as Q = (Q 1 , . . . , Q 2N ) T = ((X, P) 1 , . . . , (X, P) 2N ) T (the notation A T indicates the transpose of the matrix A). In our particular case Q = (X, P x , Y , P y ) T , recall that we are exclusively interested in the motion projected into the plane xOy, described by the Hamiltonian H ⊥ (t) . Then a set of time-dependent observables X(t) and P(t) in the Heisenberg picture will evolve as a linear combination of their initial values X and P according to the rule: where u(t) is a 2N × 2N time evolution matrix that determines uniquely the evolution operator U(t), and vice versa, iff Q spans a complete set of observables in a Hilbert space H.
We can write down the j-th coordinate of the fuzzy centre (Eq. 4): The direct application of the matrix u(t) to the initial vector Q will generate the whole trajectory of motion. Furthermore, since H ⊥ (t) is quadratic in its variables, the evolution operator U(t) will be the same in the classical and quantum regimes. Thus, the classical motion trajectories U(t) will be analogously interpreted in the quantum case, where, the mean position X will evolve according to i d dt �X� = �[X, H ⊥ (t)]� (Ehrenfest's theorem). Note that if two unitary operators U 1 (t) and U 2 (t) lead to the transformation and U 2 (t)U † 1 (t) commutes simultaneously with any function depending on X and P . Moreover, the observables X and P generate an irreducible algebra in L 2 (R) , therefore U 1 (t)U † 2 (t) must be a phase factor provided these are unitary operators. As a result we have where ϕ is a real number. Such detail becomes relevant in quantum control problems, for if any two unitary operators differ by a phase factor, they generate the same transformation of quantum states, so both are equivalent in this sense.
The condition of evolution loop will be satisfied inasmuch as the whole set of observables Q return to their initial conditions after a time interval T. In consequence the analysis is considerably simplified if we assume that Eq. (8) possesses a periodic temporal dependence, such that H ⊥ (T + t) = H ⊥ (t) , that is a Floquet Hamiltonian. As well, the periodicity of Eq. (8) implies that U(t) = U(T + t) , concluding that any observable Q in the Heisenberg picture will evolve periodically even if it does not depend explicitly on time. It follows from Eq. (10) that Q(T + t) = U † (T + t)QU(T + t) = Q(t) , or for a fixed period of time evolution Q(T) = U † (T)QU(T) = Q , which means that the loop condition indicates U(T) = e iϕ 1.

Biharmonic fields.
We shall now analyse the evolution loop trajectories (flowing inside of a solenoid with a time-varying, homogeneous current density on its surface) as well as the squeezing effects achieved by biharmonic oscillations of the type 34 : where the selection of either β 2 = 0 or ω 2 = 0 enables the harmonic case.
The temporal profile of β(t) is not restricted to a particular form, in fact it can be quite arbitrary. Nevertheless any reasonable physical implementation of β(t) shall preferably avoid any sudden jump such as the ones furnished by kicked protocols 16,21,34,35 , since the resulting electric delta shocks could compromise an even evolution. Accordingly, an alternative is to choose a sufficiently smooth pulse in terms of harmonic functions.
An interesting consequence of our cylindrical model is that the commutation relation [H osc (t), H rot (t ′ )] = 0 is invariably satisfied, thus the evolution operator U(t) can be split into two distinct components, , describing the evolution around Oz and the evolution of the magnetic oscillator U osc (t) , respectively, each one satisfying Eq. (9). The resulting evolution matrix u(t), therefore, will consist of two steps of evolution starting from the initial condition Q to Q osc (t) to Q(t) unfolded by the consecutive application of the oscillatory and rotational evolutions, in that order.
The matrix of rotations u rot (t) will be generated by the operator U rot (t) , which is straightforwardly integrated as whereas the matrix u osc (t) of dimension 4 × 4 will be constructed from U osc (t) . Fortunately, the integration of the oscillatory part can be further simplified due to u osc (t) can be reduced to a single matrix h(t) that will make each pair of observables Q x = (X, P x ) T and Q y = (Y , P y ) T to evolve simultaneously, i.e. www.nature.com/scientificreports/ will act at once in each of the subspaces spanned by Q x and Q y . In that regard, for Q j with j = x, y , we have: where concluding that h(t) simply obeys the differential equation: while the integration of this equation becomes standard for stationary fields β(t) = β 0 , yielding solutions of the form e t , in general we want a computer to do the job.
Observe that the determinant of the symplectic matrix h(t) is an integral of the evolution in that det(h(t)) = 1 . This attribute turns essential to classify the particle's motion as well as to find exact solutions, as we shall show later.
The matrix h(t) is at the core of any evolution process generated by the class of periodic fields β(T + t) = β(t) , and we must discuss the ascribed dynamical features to such matrix. First, h(t) is symplectic, it has two eigenvalues + , − such that + = 1 − , hence the algebraic structure of h(t) is wholly determined by the scalar � = tr(h(t)) , as suggested by the characteristic polynomial D( ) = 2 − � + 1 and its roots ± = 1 2 (� ± √ �) , where = 2 − 4 . We are granted, therefore, to classify the outperformed motion through | | into three categories 34 : . The motion is completely stable and restricted to an oscillatory evolution. The annihilation and creation operators, a − and a + , can be obtained from the eigenvectors of h(t). Control operations of this type are useful in ion traps. | | = 2. The eigenvalues of h(t) are merely ±1 , defining a separatrix between the stable and unstable regions.
The motion described by field amplitudes in the threshold region might originate an effective parametric resonance, which supports a variety of interesting applications, e.g. a charged particle could be attracted by repulsive forces. |�| > 2. The matrix h(t) has a pair of real eigenvalues, + = e +σ and − = e −σ , σ ∈ [0, π 2 ] . This class of motion is unstable, producing a squeezing effect over the annihilation and creation operators a − and a + , i.e. expanding a − at the cost of a + , or conversely as well.
Despite the classification given above is entirely based on the Heisenberg's evolution of canonical observables, regarding material particles, the parametric resonance region, | | = 2 , can be understood as an analogue of the parametric amplification studied by Mollow and Glauber 36 in the context of coherent photon states. Even though the trajectory picture commonly attracts the concern of many researches as for the design of ion traps.
As an example, a well known ion trap is the radio-frequency, quadrupole device engineered by Paul 6 in which the oscillating elastic forces β(t) are written in terms of Mathieu functions. Indeed, the stability regions for such trap are simply identified from a Strutt map, useful in a variety of time-dependent problems with elastic potentials, e.g. tomography 23 .
We shall now illustrate how this diagram works, not only for ion traps, but for another quantum control operations. To this aim, we have devised a computer routine to integrate numerically Eq. (17) in terms of a biharmonic field (Eq. 12) with fixed β 0 = 0 . A scanning process shall be meticulously undertaken, to track the different values of | | in the amplitude domain. We have chosen β 1 , β 2 ∈ [−15, 15] with angular frequencies ω 1 = 2π , ω 2 = 4π and t ∈ [0, 1] . As a result, we have drawn the biharmonic map in Fig. 1, where the three regions according to | | , can be located in the amplitude space spanned by {β 1 , β 2 } , in other words, depending on the values of such amplitudes, a type of quantum control operation will be furnished.
First we are to survey the stable motion assured by the set of amplitudes {β 1 , β 2 } lying in the clear areas of the map in Fig. 1. To explain how the evolution loops are originated, we have integrated Eq. (17) up to T = 6 field periods, finding the closed trajectory shown in Fig. 2-A, whose corresponding dimensionless amplitudes are β 1 = π 4 and β 2 = −10 . In this example, the particle starts its evolution with velocities (p x , p y ) = (−5, 20) at the point (x, y) = (10, −20) and finishes at the same point with velocities (p x , p y ) = (6.25, −2.41) . Any evolution loop with symmetry under parity reflection, such as this one, will have a vanishing fuzzy point (X,Ȳ ) indicating immunity to the effect of external driven forces. In fact, whenever the elastic field β(t) is biharmonic, the potential A(t, X) will evanesce at the beginning t = 0 and at the end t = T of the process, then the canonical momenta and the kinetic momenta coincide, simplifying the interpretation of the wave packet at these points.
The oscillatory motion is not restricted to circuits of evolution, it truly can be broken into open trajectories under some circumstances. These cases are particularly interesting since they involve a kind of free evolution embedded in the particle's history of motion. We recall that in the Schrödinger picture the free evolution along a time interval τ corresponds to the operator e − iτ P 2 2 , which transforms the canonical observables into X j → X j + τ P j and P j → P j via a matrix of evolution of the form (Eq. 14), although in this case h(t) → h(τ ) has the simple structure  www.nature.com/scientificreports/ therefore, should h(τ ) represents a segment of evolution loop the rest of the trajectory merely corresponds to the inverse operation e + iτ P 2 2 , that is h(−τ ) , such that h(τ )h(−τ ) = 1 . In other words, the action of h(−τ ) reverts the effects produced by h(τ ) over the set of observables, concentrating the wave packet instead of contributing to its spread across the interval of operation.
Such effects are produced by allowing the field β(t) to take amplitudes {β 1 , β 2 } living in the threshold region, | | = 2 , allowing three classes of temporal manipulations in accordance with the effective time of operation τ : (1) Accelerated free evolution, τ > T , the system will get older faster than the actual time of operation; (2) Retarded free evolution, 0 < τ < T , the system will experience a retarded progress in its temporal motion; (3) Inverted free evolution, τ < 0 , the effect is generated by the successive application of two biharmonic periods T = 2 as long as = −2 (i.e. h 11 = h 22 = −1 ) hence τ = −2h 12 with h 12 > 0.
Finally, we are to review the squeezing effects granted by those points pertaining to the coloured areas of the Strutt map, Fig. 1. Our comments ahead, however, are not addressed to the squeezed-states philosophy, rather we focus our attention on the source of the squeezing operations regarding the periodic-evolution model (Eq. 8).
Once again the symplectic structure of the matrix h(T), permits to simplify the analysis, and formulate the squeezing condition in terms of its eigenvalues as: fulfilled inside the region |�| > 2 . However, the time evolution outlined by Eq. (17) has to be calculated within a time interval where β 2 (t) is antisymmetric around the interval centre, or the squeezing effects would not occur 37 . Moreover, depending on the actual value of , there will be two possible squeezing transformations: the purely positive transformation for � > 2 and the parity transformation for � < 2 . One can freely choose between these two options via the points (β 1 , β 2 ) , based on the control purposes.
In the light of the condition (Eq. 19), it follows that the matrix elements h 12 and h 21 are equal to zero, and therefore any canonical observable Q will be transformed as: which read as the amplification of Q ′ at the cost of compressing Q ′′ , or vice versa. Such transformations are, in fact, embedded in the eigenvectors of h(T).
It might also be noted that if a fuzzy point (Eq. 4) does not vanish, the squeezing operations will sense the influence of any external force but in an orthogonal direction to its exertion, otherwise the control operations become resistant to external perturbations, such as radiation pollution 38 .
Nonetheless, in general the kind of transformations (Eq. 20) regarding eigenvalues x j → x j , p j → 1 p j in intervals [nT, (n + 1)T] , can only be produced at those separatrix points where h 12 = h 21 = 0 is satisfied, namely at the intersection points of the matrix trajectories in the separatrix belt. The squeezing transformations (Eq. 20) are realisable if the amplitudes β 1 , β 2 take values inside the unstable regions specified by the Strutt map. We have traced some of these points at which an intense squeezing = h 11 is achieved, particularly for = 2 and = 4 . For instance, a biharmonic field with the amplitudes (β 1 , β 2 ) = (−10.3, −6.9) would attain a squeezing/ amplification factor = 4.

Polyharmonic fields: exact operations.
There exists a miscellany of methods to generate quantum control operations of the form (Eq. 20), amongst them the programmes of kicked (or discrete) pulses that disrupt a continuous evolution process 16,21,34 . Unfortunately, as we have already mentioned, this kind of operations result technically impractical leading to the necessity of a more method based on soft evolution operations. If the discrete pulses are replaced by smooth pulses, e.g., biharmonic fields β(t) , the imperfections are partially circumvented. Indeed, the design of such operations could be devised considering two different smooth pulses 35 , both in the stable region of the Strutt diagram, but demanding that the product between them belongs to the instability region, yet the successive application of both pulses makes the particle to absorb an amount of energy 39 , that shall be less than the difference between the neighbouring energy levels of the particle. How can we really avoid any sudden jump in the composition of quantum control operations? We shall not answer that question with utter certainty, rather we would like to sketch an attempt in the following discussion.
From our viewpoint, the composition approach consists of taking portions of the time evolution course described by the Hamiltonian (Eq. 8), with the exception that this time the elastic fields (Eq. 7) shall be transformed as β 2 (t) → β(t) to facilitate the quest of exact solutions. Hence, in the subspace of oscillations, the matrix h(t) becomes: noting that it adopts the form of a symplectic matrix of rotations.
In particular, for any ωt = nπ 2 (n is an integer) we obtain the squeezing transformations: www.nature.com/scientificreports/ now, the successive application of two matrices of this type leads to: leading again to the set of operations in Eq. (20) with the requirement of two distinct frequencies ω 1 and ω 2 at different times, say t 1 and t 2 , such that ω 1 t 1 = ω 2 t 2 = π 2 . The sudden jumps have not been removed already, for example, should the protocol be applied in the void background there would be at least three jumps involved: 0 → ω 1 → ω 2 → 0 . What if we translate this procedure to the continuum?
We have shown that the set of operations (Eq. 20) are actually generated by symplectic matrices h(t) with the property h 11 = h 22 = 1 2 , i.e. h(t) has the form of a Toeplitz matrix. This agreeable property grants that if η and ξ are Toeplitz matrices, their anti-commutator algebra ηξ + ξη , as well as their symmetric products ηξ η and ξ ηξ , furnish matrices of the same class. It follows that the squeezing transformations (Eq. 20) can be built from the contribution of a big number of symmetric products between symplectic matrices η(t) of the form (Eq. 21), each one at different time subintervals t j with a definite elastic field amplitude β(t j ) , where j = 0, 1, 2, . . . , namely: once again, conveying the property h 11 = h 22 = 1 2 . To obtain the continuous analogue we shall assume that the infinitesimal jumps dh(t) are assembled by the infinitesimal contributions dη(t) = �(t)dt , from the right to the left sides as in Eq. (24), each one depending on a symmetric field β(t) = β(−t) around t = 0 , arriving at the matrix differential equation in the expanded interval [−t, t]: which is explicitly written as nonetheless, the actual trajectory determined by the whole evolution process, requires the integration of Eq. (17) along an asymmetric interval, since the character of Eq. (25) is rather auxiliary as will be clear below. On top of that the anti-commutative structure of Eq. (25) defines a β(t) in terms of a smooth enough, real function θ(t) . What we are to discuss now is precisely how to obtain an exact solution of the inverse evolution problem.
Even without any specialised consideration on adiabatic invariants 40,41 , β(t) as expressed in Eq. (29) is an exact solution of the inverse evolution problem for h(t) given a function θ(t) , which is, in principle, arbitrary though constricted to satisfy non-trivial conditions at singular points. Nevertheless, for any asymmetric interval [t 0 , t] the dependence of h(t) on β(t) must be determined by the integration of Eq. (17), as we have already shown for the case of biharmonic fields.
Some simple relations between β(t) and θ(t) deserve observation. Any field amplitude given by Eq. (29) in a symmetric interval [−t, t] , is accomplished by a sufficiently smooth function θ(t) to assure continuity and differentiability. Specifically, at any point t where θ(t) = 0 , there must be dθ(t) dt = ±2 . As well, if θ(t) = 0 but dθ(t) dt = 0 then Eq. (25) becomes a matrix of squeezing transformations such as the one in Eq. (23). Moreover, if dt 3 = 0 it follows necessarily that dβ(t) dt = 0. The only missing piece is the construction of θ(t) . As we said, this function is truly arbitrary, although from an empirical viewpoint a natural choice would be a θ(t) in terms of harmonic functions to induce soft control www.nature.com/scientificreports/ operations and prevent the evolution process from any leap. Amongst the most elementary cases, let us consider the polyharmonic function since this function is antisymmetric around t = 0 the corresponding β(t) as for Eq. (29) will be symmetric around that point. Particularly, if we fix the frequencies as ω 1 = 1, ω 3 = 3, ω 5 = 5 and ω 7 = 7 , at the endpoints of the symmetric interval − π 2 , π 2 a matrix of squeezing (Eq. 23) will emerge with h 12 = ±ω = b . Accordingly, the following initial conditions must be fulfilled: where c = 0 is a real parameter that we shall carefully fix to achieve the desired control operation. This set of rules gives additional information regarding the nontrivial relation between β(t) and θ(t) . The condition θ π 2 = b determines the magnitude of the squeezing transformation depending on the whole trajectory, whereas θ ′ (0) = 2 endows β(t) with non-singularity at t = 0 . In turn, the coefficients of θ(t) are: see the numerical examples of the generated pulses β(t) in Fig. 3, in which we have selected β − π 2 = β π 2 = 0. Unlike the method based on biharmonic fields, where we could generate a map such as the one presented in Fig. 1, in the current case we are unable to pursue an analogue survey to scan the type of control operations produced by a set of amplitudes {β 1 , β 2 } . Even though, in the light of polyharmonic pulses the functions (Eq. 29) and (Eq. 30) truly offer additional information regarding the effects achieved at the end of the symmetric interval of operation: If β(t) > 0 the magnetic field would induce a squeezing effect over the canonical observables in the concatenated intervals − π 2 , π 2 and π 2 , 3π 2 , otherwise, the magnetic field would reproduce an evolution loop effect, either closed or broken.
We shall remark that the trajectory smeared in the interval − π 2 , π 2 ∪ π 2 , 3π 2 is not yet determined at this stage, which calls for a separate computer routine to integrate (Eq. 17) in the asymmetric interval where the initial time will be t 0 = − π 2 . That is, the integration of Eq. (17) will be split into two subintervals, each one with a definite pulse β(t) . The evolution process will start at the time − π 2 with the first pulse β(t) , then continuing the integration in the following subinterval at π 2 with the second pulse β(t) and finishing at 3 2π , in order to draw a congruent trajectory.
Please note, as for the pulses in Fig. 3, the resulting control operation would be absent of any abrupt discontinuity even when it was regulated by the concatenation of two different fields, for example, the amplitudes β(t) with solid and dashed curves can achieve a squeezing transformation in the fashion of (Eq. 20), whose actual trajectory in the xOy plane is plotted in Fig. 4(A) for a charged particle with initial conditions (x, y) = 1 2 , 10 (30) θ(t) = a 1 sin(ω 1 t) + a 3 sin(ω 3 t) + a 5 sin(ω 5 t) + a 7 sin(ω 7 t), θ π 2 = a 1 − a 3 + a 5 − a 7 = b, θ ′ (0) = a 1 + 3a 3 + 5a 5 + 7a 7 = 2, . The canonical observables of position are transformed as X → x X and Y → y Y in agreement with (Eq. 20), our numerical computation indicates that such particle would be subject to a squeezing transformation with x = −0.403 and y = −1.125 , meaning that the momenta have been amplified and compressed, respectively. These effects, furthermore, can become significantly boosted by external driven forces either time-dependent or not. In such scenario the evolution process must account for the unperturbed Hamiltonian (Eq. 8) plus the external field. For instance, we want to consider a force aligned in the Ox direction, F = (F, 0, 0) , leading to the perturbed Hamiltonian H(t) = H ⊥ (t) + FX , such that the corresponding evolution operator is constituted as U(t) = U(t)W(t) , where U(t) is the unperturbed evolution operator, whereas W(t) represents the contribution of the perturbative force that obeys the differential equation where X(t) is the time-dependent observable fostered by the Heisenberg evolution (Eq. 10).
Even more, since the evolved observable X(t) is linear in the canonical variables, the commutators [X(t ′ ), X(t)] are reduced to numbers. The continuous Baker-Campbell-Hausdorff formula 42 then implies that Eq. (33) has a solution of the form e iF t 0 dt ′ X(t ′ ) up to a phase factor e iφ , where φ is real. It follows that if t = T is the loop period in the unperturbed evolution operator, then any initial coordinate X j of the fuzzy centre (Eq. 4) will be reconstructed from the integral t 0 dt ′ X(t ′ ) , therefore U(T) = U(T)W(T) = e iφ e iTFX j . Suggesting that X j will be unaffected by the transit of W(t), namely U † (T)X j U(T) = e −iTFX j X j e iTFX j =X j . However any other coordinate of the fuzzy centre could suffer a drift U † (T)X k U(T) =X k − iTF[X j ,X k ] , of course j = k.
We conclude that any fuzzy point (X j , . . . ,X k ) with noncommutative coordinates will result in a broken loop every time it is manoeuvred by an external force and the coordinates of its centre will drift in a transversal direction with respect to the force. On the contrary, if the fuzzy centre is such that its centre [X j ,X k ] = 0 , then at the end of the period T the point (X j , . . . ,X k ) will return to its initial value.
To illustrate our discussion, we have examined the inverted squeezing operation referred in Fig. 4(A) in the presence of a time-dependent perturbation sin(t) aligned in the Ox direction to distort the operational trajectory, see Fig. 4(B). As we mentioned before, the evolution process (17) first sweeps in the interval t ∈ − π 2 , π 2 with a given polyharmonic pulse β 1 (t) , and continues with a different β 2 (t) in the interval t ∈ π 2 , 3π 2 . To invert the operation it is required, therefore, a second application of β 2 (t) with t ∈ 3π 2 , 5π 2 and then a second application of β 1 (t) with t ∈ 5π 2 , 7π 2 . In our example, note that the perturbed evolution returns to its initial position after completing the period of application. Also the effect of the driven external force sin(t) in this case produces a set (33)  The evolution process starts at t 0 = − π 2 in terms of a polyharmonic amplitude (Eq. 29) defined by the parameters b = 2, c = −3 finalising its application at t = π 2 -where the amplitude vanishes. Immediately, at the very same time, the evolution process continues its development now in terms of the pulse defined by b = 9 5 , c = − 7 2 that will end its application at t = 3π 2 . The resulting operation attains a squeezing transformation with x = −0.403 and y = −1.125 . (B) To invert the operation and recover the initial observables's values, the second pulse has to be applied once again from t = 3π 2 to t = 5π 2 , where the evolution process is now conducted by the first pulse during from t = 5π 2 to t = 7π 2 . The dotted line represents the distorted evolution generated by the presence of an external harmonic force sin(t) in the direction Ox, even though note the original configuration is entirely recovered by reversing the operation. www.nature.com/scientificreports/ of transformations (Eq. 20) at t = 3π 2 with x = 3.92 and y = −0.96-unlike the unperturbed case, this time we have amplified the observable X and compressed Y at the cost of compressing P x and amplifying P y .
We have commented before that those polyharmonic pulses (Eq. 29) greater than zero convey the adequate oscillatory motion to produce squeezing operations, given that the corresponding matrix h(t) fulfils |�| > 2 . However, it has not been surveyed yet a more general type of pulses that can be either greater or less than zero, as the amplitude β(t) with dotted line in Fig. 3, and whose associated matrix h(t) furnishes a stable motion in the region |�| < 2 . Accordingly, this class of pulses β(t) does bring the alternative of evolution loop operations by merely implementing the programme outlined for biharmonic fields.
In this regard, please refer to the loops shown in Fig. 5. In panel A, we present a symmetric polyharmonic loop that is resistant to the application of an external potential. By resistant we mean that the loop is absolutely stable-as the examples shown in Fig. 2-and cannot be broken under the influence of forces coming from the outside, in fact, under these circumstances any solution of Eq. (33) will account as a phase factor. On the other hand, the loop in panel B has been deformed by an external field of constant magnitude (see the figure caption for details). It truly has a fuzzy point (Eq. 4) with noncommutative coordinates, then its fuzzy centre will suffer a drift in the transversal direction and eventually the loop could be broken.

Discussion
Our discussion above is fundamentally focused on the time dependent family of matrices h(t) responsible for a number of quantum control operations, although our scheme does merely constitute an approximation. Usually this type of imperfections are not entirely disappointing: In Paul's ion traps 6 , for instance, the propagation of the electromagnetic signals through the device's interiors is usually disregarded due to the pretty small size of the trap. However, adopting a rigorous attitude, even those softly changing potentials sensed on the trap surfaces must produce field corrections starting from 1 c (post Newtonian) terms in the Einstein-Infeld-Hoffmann (EIH) approximation 43 . Below we are to briefly discuss a similar problem for softly changing magnetic fields β(t) such as Eq. (12) or (29).
Throughout this report, we have been considering a time dependent, homogeneous magnetic field B(t, x) in a cylindric solenoid. It truly does not fulfil the Maxwell's equations, yet it obeys a sequence of EIH approximations. To evaluate the conveyed errors, let us look for the exact time dependent vector potentials of a cylindrical solenoid in the form (Eq. 1): where the magnetic field B instead of depending only on t, could also depend on the radius r of the solenoid's cross section, and n is a unit vector in the field direction. To assure the relativistic sense of Eq. (34) we must assume A(t, x) = 4π c J(t, x) = 0 , with the notation = 1 c 2 ∂ 2 ∂t 2 − ∇ 2 (the d' Alembert operator). As easily seen, www.nature.com/scientificreports/ the application of the Laplacian to the right hand side of Eq. (34) is equivalent to the acts of the operator D = ∂ 2 ∂r 2 + 3 r ∂ ∂r to B(t, r) alone. Hence, A(t, x) means: The magnetic field B(t, r) must be analytic around Oz, then we preferably look for a solution in the form: where B 0 (t) = B(t) is the homogeneous quasi-static approximation. Since Dr 2n = 4n(n + 1)r 2(n−1) , the substitution of Eq. (36) into Eq. (35) yields: We now introduce a dimensionless time τ = t T , where T is some conventional time unit corresponding to the laboratory observation, and by writing Eq. (37) in terms of the time derivatives ∂ ∂τ one can reduce it to: We kept here the B(t) depending on the actual time t, in order to assure that all derivatives ∂ n ∂τ n B(t) will be expressed in magnetic field units. The curious property of this formula is the absence of terms proportional to 1 c -in fact, the field propagation law (Eq. 35) is solved exclusively in terms of extremely small contributions proportional to 1 c 2 . It suggests that the superpositions of delicate wave fronts running towards the solenoid centre create a good approximation of the quasi-static theory.
Even though, this does not explain how to create (or at least approximate) the first magnetic step B(t), to wake up the whole iterative series (Eq. 38). (That is, how to induce the homogeneous surface currents which do not depend on z, but depend on time in any desired way.) In the static case, the magnetic field inside the solenoid is generated by a stationary current I circulating around the surface, i.e., B = 4π c �I �z . The question is, how to produce the circulating currents depending on time, but homogeneous (independent from z) on every section of the surface? Should the solenoid was constructed as a single spiral wire around the cylindrical surface, connected at its both extremes to the potential difference �(t) , then even a subtle change δ�(t) would be propagated along the solenoid as a current pulse, creating a softly changing but z-dependent field, instead of the desired quasi-static B(t).
There is an alternative. One rather could consider a cylindrical surface of non-conducting material (e.g., glass) of radius a, charged uniformly with a surface density σ , so that each circular belt of 1cm of height contains a charge aσ . The experimental challenge is not extraordinary. Suppose the cylinder has a radius a = 20cm, and it is rotating with an angular velocity ω = 1s −1 around its symmetrical axis aligned to Oz and has 1C of charge at every horizontal belt of 1cm, then inside the cylinder will be generated a homogeneous magnetic field nB of intensity: at least in the post-post-Newtonian approximation. Hence, by employing the softly changing angular velocity ω → ω(t) one truly can generate the practically homogeneous magnetic field nB(t) of the quasi-static environment described by Eq. (38). Will such technique work?
As a remark, the actual operation time t T → t given in Eq. (6) can be arbitrarily large (or short) and the fields (Eq. 7) arbitrarily strong (or weak). To form an idea of the real orders of magnitude required by our control operations, we shall present the Table 1 below.
To this aim, we have also decided to check the magnitude of the Abraham-Lorentz radiative force. While the ordinary force of the variable oscillator trajectory is simply F osc (t) = mẍ(t) , the hypothetical radiative force is expressed as F rad (t) = mγ ...
x (t) , where γ = 2 3 e 2 mc 3 is the particle dependent characteristic time. Using now the definitions of dimensional quantities (Eq. 6) we can compare the magnitudes of the conventional and radiative forces for the squeezing operations (see the last row in Table 1) finding the radiative ones extremely small albeit slowly increasing as T decreases (higher frequencies).
Our preliminary evaluations become valid whenever the physical size of the solenoid is realisable. In many laboratories the techniques to cool and trap ions are adequate to hardly study the atomic structure but insufficient for wider purposes, since the time dependent oscillator potential is created only in a strictly local scale, such is the case of a quadrupole trap whose immediate vicinity from the central axis is conformed, in some cases, just by four metal bars 15 .
In our calculations we have assumed pure states, evolving under the influence of slowly changing external fields, without taking corrections for traces of matter across the ion traps or solenoids. Moreover, we have neglected the possible packet reflection or absorption by the laboratory walls. Suggesting an apparatus whose dimensions make any of these considerations to become insignificant. Actually, in case that particle would be absorbed by the (surface) walls, the problem leads to the fundamental question of the time operator which persists without a truly convincing solution even in case of flat surfaces. We thus hope that for wide traps, our proposal brings something of interest to the quantum control theory.   www.nature.com/scientificreports/ As well, we have noticed that in works 41,44 using the Ermakov-Milne invariants the operations are supposed to be faster that the times given here. In that regard, our contribution is slightly different: Its main goal is the use of a simple Toeplitz algebra to obtain the exact solutions (Eq. 29) whereas other approaches use the suggestive idea of frictionless driving, which transports the states without changing the eigenvalues of certain invariants. But could the more general variational methods be applied at the level of a θ(t) function? The question remains open.
Received: 6 November 2020; Accepted: 7 December 2020 Table 1. Physical conditions (in cgs) to achieve the squeezing transformations x = −0.43 and y = −1.125 on a proton moving inside a cylindrical solenoid. The reported quantities have been calculated regarding the example in Fig. 4. The physical magnitudes of q, p and v, with q = x = y and p = p x = p y correspond to the dimensionless values x ′ = y ′ = p ′ x = p ′ y = 1 according to Eq. (6). Note how the incredibly modest strength of the magnetic field grows as the operation time T becomes shorter, to form an idea, a control operation of T = 10 −2 s. would require half the strength of Earth's magnetic field on the equator. In the last row it is reported the average ratio of the Abraham-Lorentz radiative force to the time-dependent oscillator forces, at different operation intervals T.