Geometric Phase in Quantum Synchronization

We consider a quantum limit-cycle oscillator implemented in a spin system whose quantization axis is slowly rotated. Using a kinematic approach to define geometric phases in nonunitary evolution, we show that the quantum limit-cycle oscillator attains a geometric phase when the rotation is sufficiently slow. In the presence of an external signal, the geometric phase as a function of the signal strength and the detuning between the signal and the natural frequency of oscillation shows a structure that is strikingly similar to the Arnold tongue of synchronization. Surprisingly, this structure vanishes together with the Arnold tongue when the system is in a parameter regime of synchronization blockade. We derive an analytic expression for the geometric phase of this system, valid in the limit of slow rotation of the quantization axis and weak external signal strength, and we provide an intuitive interpretation for this surprising effect.


I. INTRODUCTION
In a seminal paper, Berry showed that a quantum system initialized in an eigenstate of its parameterdependent Hamiltonian acquires a so-called geometric phase factor upon adiabatic transport around a closed path in parameter space [1].Unlike the familiar dynamical phase, acquired by the system due to its time evolution, the geometric phase (GP) depends solely on the curvature of the parameter space and the path taken through it, and is thereby a purely geometric quantity.Pancharatnam had already discovered a similar phase in classical optics earlier [2] and Hannay generalized the concept to classical mechanics afterwards [3].A prominent example of a GP in classical mechanics is provided by the Foucault pendulum, whose plane of oscillation rotates by an angle that depends only on the latitude of the pendulum if the period of oscillation is much shorter than a day [4].GPs appear in diverse settings including light propagation in an optical fiber [5], the Aharonov-Bohm effect [6][7][8], and the quantum Hall effect [9].GPs have also been proposed [10,11] as a way to implement quantum gates that are robust against certain pulse imperfections and parameter uncertainties [12][13][14], and such gates have been experimentally demonstrated in a number of systems [15][16][17][18][19].
Building on the work of Pancharatnam, the concept of a quantum GP has been extended to nonadiabatic evolution [8], noncyclic evolution [20], and mixed states [21][22][23] (including periodic fermionic systems in mixed states [24]).A so-called kinematic approach has been formulated [25], which is based on the (time dependent) density matrix along a path in state space and enables the concept of a GP to be generalized to nonunitary evolution [26].In classical nonlinear dynamics, Kepler et al. [27,28] extended the concept of a GP to classical dissipative systems with self-sustained oscillations (i.e., limit cycles).They showed that, even though these systems are nonconservative, cyclic adiabatic deformations of the limit cycle lead to GP shifts which could potentially be observed in certain chemical reactions.
In this paper, we analyze geometric phases in quantum limit-cycle oscillations using the definition of a GP in nonunitary quantum evolution proposed in Ref. 26.We implement a numerically stable algorithm to calculate the GP of a quantum system undergoing dissipative time evolution.We use this algorithm to demonstrate that a quantum limit-cycle oscillator implemented in a spin-1 system acquires a GP that depends only on the trajectory through parameter space if the direction of its quantization axis changes slowly compared to the timescales of its coherent and dissipative dynamics, similar to the classical case [27,28].
We then consider the more general case of a quantum limit-cycle oscillator subject to an external signal and show that the GP has a tongue-like structure very similar to the well-known Arnold tongue in synchronization [29].Surprisingly, this tongue-like structure of the GP vanishes together with the conventional Arnold tongue of synchronization if the system is in a regime of interference-based quantum synchronization blockade [30].To understand this surprising effect, we derive an analytical formula for the GP, which is valid if the quantization axis is rotated sufficiently slowly.This paper is structured as follows: In Sec.II, we summarize the kinematic approach to the GP in an open quantum system before we introduce the numerical algorithm to compute the GP in a generic open quantum system in Sec.III.In Sec.IV, we focus on the specific example of a van der Pol (vdP) oscillator subject to a weak external signal.We demonstrate the surprising similarities between the GP and the Arnold tongue of synchronization in this system.To gain better insight into this phenomenon, we derive an analytical expression for the GP in an arbitrary quantum limit-cycle oscillator with nondegenerate populations.Finally, we conclude in Sec.V.

II. GEOMETRIC PHASE IN OPEN QUANTUM SYSTEMS
The GP of a quantum system can be defined as the difference between the global phase acquired during the time evolution and the local phase changes accrued along the way [31].This subtraction is equivalent to enforcing a parallel-transport condition [31].To apply this definition to a quantum system in a mixed state undergoing nonunitary evolution, one has to consider the GP of a purification of the system, which is measurable in principle [26] by an interferometric measurement of the purified state.Note, however, that the interferometric measurement requires unitary evolution of an enlarged system comprising the system of interest and ancillary degrees of freedom because the value of the GP depends on the chosen purification [32] (see App.A for more details).Given an open quantum system undergoing evolution along a path P in the space of density matrices, where ρ(t) is the time-dependent density matrix of the system, p k ≥ 0 are its populations (which we assume to be nondegenerate functions for t ∈ [0, τ ]), |φ k (t) are the corresponding eigenvectors, and N is the Hilbert-space dimension, Tong et al. [26] proposed the following definition of the GP γ[P]: Intuitively, Eq. ( 2) is the sum over the Pancharatnam phases of each eigenstate |φ k (t) of the density matrix ρ(t), weighted by the corresponding population p k (t) at the start and the end of the path P, where the exponential factors subtract the local phase changes accrued along P.

III. ALGORITHM TO CALCULATE THE GEOMETRIC PHASE NUMERICALLY
Solutions to dissipative quantum systems in closed analytical form are rare: beyond the simplest models, it is impossible to solve the differential equations arising from a quantum master equation (QME) analytically.We therefore use a numerical approach to evaluate Eq. ( 2) and calculate the time evolution of a given quantum system by solving its QME numerically exactly using the QuantumOptics package [33] in Julia [34].This provides us with the density matrix ρj = ρ(t j ) of the system at equidistant discrete time steps t j = j∆t, j ∈ {0, . . ., n step }, where ∆t = τ /n step .Using this set Algorithm 1: Numerical calculation of the geometric phase using Eq. ( 2).
For each time step t j , we diagonalize ρj numerically to find the populations p k (t j ) and the associated eigenvectors |φ k (t j ) defined in Eq. ( 1).Note that for the specific problem considered later, the populations p k (t) are constant and distinct, such that we can order them ascendingly and there is no ambiguity in the labeling of the eigenvectors in different time steps.Since the eigenvectors are only defined up to a global phase factor, we use the convention that the kth entry of the kth eigenvector is real and positive.This is equivalent to choosing a particular gauge and does not affect the GP since Eq. ( 2) is gauge-invariant.
In step 2 of Alg. 1, the exponential phase factors in Eq. ( 2) are calculated by numerically differentiating {|φ k (t j ) | j = 1, . . ., n step } with respect to time using a symmetric difference quotient that is of fourth order in ∆t [35].The overlaps { φ k (t j ) φk (t j ) | j = 1, . . ., n step } are then numerically integrated over time using an extended Simpson rule in step 3, which is also of fourth order in ∆t [36].The specific choice of the gauge of the eigenvectors |φ k (t j ) ensures that the overlaps are smooth functions of time and that the numerical integration is stable.
We benchmarked this algorithm using the exactly solvable case of a qubit subject to dephasing [26], see App.B for more details.

IV. GEOMETRIC PHASE OF A QUANTUM LIMIT-CYCLE OSCILLATOR
A. Quantum van der Pol oscillator The numerical algorithm introduced in Sec.III can be applied to any quantum system whose density matrix ρ(t) is known as a function of time and has distinct populations p k (t).In the following, we focus on a specific example, namely, a quantum limit-cycle oscillator implemented in a spin-1 system whose quantization axis changes slowly.
A classical limit-cycle oscillator is a nonlinear dynami- FIG. 1.(a) We consider a quantum limit-cycle oscillator, implemented in a spin-1 system, whose quantization axis is rotated with angular frequency ω on the surface of a cone with opening angle α and symmetry axis n(α), as defined in Eq. ( 9).The quantization axis rotates counter-clockwise about n(α) and points along the z-direction at time t = 0. (b) Trajectory of the instantaneous eigenstate |φ+1(t) of the density matrix ρ(t) on the spin-1 Bloch sphere in the laboratory frame, obtained by numerical integration of Eq. ( 18).In the limit-cycle state, i.e., without any external signal, the eigenstate follows the rotation of the quantization axis (red line).
If an external signal Ĥsig(t) = 0 is applied to synchronize the quantum limit-cycle oscillator, |φ+1(t) precesses about the instantaneous quantization axis (blue line).The eigenstate |φ−1(t) follows a similar trajectory on the opposite part of the Bloch sphere (not shown).For presentation purposes, ω and the signal strength T have been chosen much larger than in the numerical examples.
Here, we consider a quantum limit-cycle oscillator implemented in a spin-1 system, which is convenient since it allows us to work with a finite-dimensional Hilbert space, N = 3.We follow the framework introduced in Ref. 30, which defines synchronization based on the phase-space dynamics of the spin system.A quantum limit-cycle oscillator can be modeled by a QME of the form ( = 1) where pator.The Hamiltonian Ĥ0 determines the natural frequency of oscillation ω 0 of the limit-cycle oscillator.We choose the quantization axis to be the z-direction and set The spin operators obey the commutation relation [ Ŝj , Ŝk ] = iε jkl Ŝl , where j, k, l ∈ {x, y, z}, and they are the generators of the rotation group SO(3).A basis of the Hilbert space is given by the joint eigenstates |S, m of Ŝ2 = Ŝ2 x + Ŝ2 y + Ŝ2 z and Ŝz , where S = 1 and m ∈ {+1, 0, −1}.The set of jump operators { Γ1 , . . ., ΓM } determines how the amplitude of the limit-cycle oscillator is stabilized in phase space and should not introduce any phase preference of the oscillation.For our numerical examples, we consider a spin-1 implementation of a quantum vdP oscillator [42,43], such that M = 2 and Here, γ g and γ d denote the gain and damping rates, respectively, and Ŝ± = Ŝx ± i Ŝy are the raising and lowering operators.The specific form of the jump operators ( 5) and ( 6) can be motivated as follows.In the quantum regime, γ g γ d , the bosonic quantum vdP oscillator populates only the lowest three Fock states [42,43].Thus, the bosonic system can be approximated by a spin-1 system whose jump operators have the same matrix representation as the corresponding bosonic jump operators restricted to the lowest three Fock states [30].

B. Demonstration of a geometric phase in a quantum van der Pol oscillator
Kepler et al. [27] demonstrated that a classical limitcycle oscillator acquires a geometric phase if its limit cycle is adiabatically deformed along a closed path in parameter space.To generate a similar effect in the quantum limit-cycle oscillator defined in Eq. ( 3), we choose to rotate the direction of the quantization axis along a path P lab in parameter space, i.e., where the rotation operator is Here, Ŝ = ( Ŝx , Ŝy , Ŝz ) is the vector of spin operators, and the unit vector n(α) = (sin α, 0, cos α) defines the symmetry axis of a cone with opening angle α.The quantization axis rotates on the surface of this cone, as shown in Fig. 1(a), which constitutes the path P lab in parameter space.The time evolution of the quantum system along the path P lab induces the path P in the space of density

Geometric phase ±γ[P]
Opening angle Demonstration of a geometric phase (GP) in a quantum van der Pol limit-cycle oscillator whose quantization axis is slowly rotated on the surface of a cone with opening angle α along the path P lab ,as shown in Fig. 1.If the rotation frequency ω is much smaller than the natural frequency of oscillation ω0 and the dissipation rates γ g,d , the system acquires a purely geometric phase whose sign depends only on the direction of the path P lab in parameter space (red lines and markers).For larger ω, the evolution is no longer adiabatic (blue lines and markers) and the GP for the respective path P in the space of density matrices depends also on the velocity at which P lab is traversed.The black dashed line shows an analytical formula for the GP in the limit ω/ω0 → 0, given by Eq. ( 35).The parameters are ω0/γ d = 10, γg/γ d = 0.1, and nstep = 10 8 .
matrices defined in Eq. ( 1), for which we can calculate the geometric phase γ[P] using Eq. ( 2).As discussed by Aharonov and Anandan [8] in the context of unitary evolution, γ[P] can be viewed equally well as a geometric phase of the path P lab in parameter space in the limit of adiabatic evolution.In our scenario, adiabatic evolution corresponds to a rotation frequency ω of the quantization axis that is much smaller than the relaxation rates and the natural frequency of the system, ω γ g , γ d , ω 0 [55].To demonstrate the existence of a GP γ[P], which (in addition to being a geometric quantity of the path P in the space of density matrices) depends only the geometry of the path P lab through parameter space, we simulate a full rotation of the quantization axis for different values of the rotation frequency ω, and compute the phase acquired along this path P lab with Alg. 1.For slow rotation of the quantization axis, i.e., adiabatic evolution, we expect the resultant phase to depend only on the geometry of the path traced out by the rotation.Therefore, reversing the rotation direction, P lab → −P lab , should only flip the sign of the phase.Figure 2 demonstrates that this is indeed the case if the rotation frequency ω is small enough, thus, the limit-cycle oscillator indeed acquires a GP in the regime ω γ g , γ d , ω 0 that depends only on the rotation of the quantization axis in parameter space.Note that, in the nonadiabatic case (i.e., for fast rotations of the quantization axis compared to the relaxation rates and natural frequency of the system), the geometric phase along the path P can no longer be viewed as a geometric quantity of the path P lab in parameter space, such that the geometric phases obtained for P lab and −P lab differ.These numerical results in Fig. 2 can be explained by the following heuristic argument.In the limit of infinitely slow rotation, i.e., ω → 0, we expect the system always to remain in the steady state along the current direction of the quantization axis, i.e., ρ(t) is well approximated by the steady-state solution of Eq. ( 3) rotated by R(α, t).Sjöqvist et al. [22] showed that the GP of a mixed state ρ(t) undergoing unitary evolution is the weighted sum of the GPs acquired by each eigenvector |φ k (t) , i.e., we find (10) where the populations are The phase factors given by the exponential functions measure the solid angle traced out by each eigenvector |φ k (t) and are a generalization of Berry's result to a spin-1 system [1].As shown in Fig. 2, this formula for the GP in the limit ω → 0 matches perfectly with the numerical results.Note that we will provide a more rigorous derivation of the GP in Sec.IV D, where we show that the heuristically motivated result given by Eqs. ( 10) and ( 11) is a limiting case of a more general calculation.
C. Arnold tongue of the geometric phase in the presence of an external signal So far, we have shown that a single isolated quantum limit-cycle oscillator acquires a geometric phase upon adiabatic rotation of its quantization axis.Quantum limit-cycle oscillators are of particular interest because they can be synchronized to an external signal at frequency ω.In this so-called entrainment phenomenon, the external signal causes the limit-cycle oscillator to deviate from its natural frequency ω 0 .The degree to which the frequency of oscillation is modified depends on the detuning ∆ = ω − ω 0 and the strength of the signal.In general, the entrainment is strongest on resonance, ∆ = 0, and the range of detuning where synchronization occurs grows with increasing signal strength T .This gives rise to the so-called Arnold tongue of synchronization, a roughly triangular-shaped region in the ∆-T parameter space [29].To describe the presence of an external signal, we replace Ĥ0 → Ĥ0 + Ĥsig (t) in Eq. ( 3) with the signal Hamiltonian which aims to rotate the state of the limit-cycle oscillator about an axis in the equatorial plane at an angle φ with respect to the positive x-axis.In the remaining parts of this paper, we focus on the GP of a quantum limit-cycle oscillator subject to the external signal given by Eq. ( 12), and we discover striking similarities between the Arnold tongue of synchronization and a corresponding plot of the GP as a function of ∆ and T .
As a preparation, we first calculate the Arnold tongue of synchronization of the system for a fixed quantization axis, i.e., for ω = 0.In a frame rotating at the signal frequency ω and using a rotating-wave approximation, Eq. ( 3) becomes time-independent and, to first order in the small signal strength T , its steady state has the form For the quantum vdP oscillator considered here, the populations p m are given in Eq. ( 11) and the coherences T c m,m are Following Refs.46 and 30, we define a phasespace quasiprobability distribution of the limitcycle oscillator by calculating its Husimi-Q function Q(θ, φ|ρ) = θ, φ| ρ |θ, φ , where |θ, φ = e −iφ Ŝz e −iθ Ŝy |S = 1, m = +1 are coherent spin states [56].Since we are working in a frame rotating at ω, the variable φ determines the relative phase between the limit-cycle oscillator and the applied signal.From the Q function, we obtain the shifted phase distribution of ρ, which is zero if the relative phase is uniformly distributed (no synchronization) and nonzero if there is a preferred relative phase.A single-number measure of synchronization can be obtained by considering the maximum of S(φ|ρ), S is positive (zero) if there is (no) synchronization and the particular value of φ max = − arg(c 1,0 + c 0,−1 ) maximizing S(φ|ρ) determines the relative phase lag between the limit-cycle oscillator and the signal.
The Arnold tongue of quantum synchronization can now be visualized by plotting S(ρ ss ) as a function of ∆ and T , which is shown in Fig. 3(a).In contrast to the classical case, quantum noise smears out the synchronization transition and leads to a smooth crossover from no synchronization at large detunings and weak signal strength (dark colors) to a roughly triangularshaped region of synchronization for sufficiently large signal strength around resonance (bright colors).
We now analyze the GP of a limit-cycle oscillator with an applied signal whose quantization axis is slowly rotated, ω = 0, as shown in Fig. 1(a).Numerically integrating the QME and using the algorithm described in Sec.III to calculate the GP γ[P], we find the GP shown in Fig. 3(b).Comparing Figs.3(a) and (b), we find a striking similarity between the GP and the Arnold tongue of synchronization: Both quantities take a constant value at large detuning and small signal strength, and vary strongly in a triangular region around resonance whose width grows with increasing signal strength.
One may argue that this coincidence is not surprising since both S(ρ) and γ[P] depend on the density matrix of the system.For T → 0 or |∆| → ∞, the external signal cannot significantly affect the limit-cycle oscillator and its density matrix is essentially independent of T and ∆ and equivalent to that of an unperturbed vdP oscillator.Close to resonance and for large enough T , however, the signal will significantly affect the oscillation dynamics and changes in both the synchronization measure and the GP are to be expected.
However, the similarities between the GP and the synchronization measure do not end here.In Fig. 3(c), we plot S(ρ ss ) for parameters in a so-called interferencebased quantum-synchronization-blockade regime [30].In this regime, the gain and damping rates, γ g and γ d , respectively, are chosen such that the coherences entering the definition of S(ρ) in Eq. ( 17) have the same magnitude but opposite signs, c 1,0 = −c 0,−1 .On resonance, this relation takes the following form In this regime, each coherence is nonzero, i.e., the signal does modify the dynamics of the limit-cycle oscillator appreciably, but an interference effect prevents phase localization such that S(ρ ss ) = 0. Therefore, the Arnold tongue of synchronization vanishes, as shown in Fig. 3(c).Surprisingly, the Arnold-tongue-like structure in the plot of the GP vanishes in the synchronization-blockade regime, too, even though the S(ρ ss ) and γ[P] measure very different properties of the density matrix.This result suggests a deeper connection between synchronization and the GP.

D. Approximate analytic expression for the geometric phase
To get a better understanding for the numerical results shown in Fig. 3, we now derive approximate analytic expressions for the GP of a quantum limit-cycle oscillator whose quantization axis is slowly rotated as described by the QME (18).In a first step, we use the fact that the rotation R(α, t) is adiabatic, i.e., the timescale 2π/ω on which the direction of the quantization axis changes is much longer than any other timescale of the system.An approximate solution for ρ(t) can thus be obtained by calculating the steady state of the system for a fixed orientation of the quantization axis, and then rotating this steady state according to R(α, t).This motivates us to introduce the frame which co-rotates with the quantization axis, Naively, one may now attempt to simplify Eq. ( 21) by switching to a rotating frame with respect to the signal, and by performing a rotating-wave approximation.However, this approach eliminates the Ŝx correction to the quantization axis and leads to incorrect results.To preserve this term, we first diagonalize the modified Hamiltonian Ĥ0 + Ĥaxis to leading order in ω/ω 0 using a Schrieffer-Wolff transformation χSW = e Ŵ χe − Ŵ (22) with the generator The QME for the density matrix χSW in this new frame is The Schrieffer-Wolff transformation rotates the spin basis states such that Ĥ0 + Ĥaxis becomes diagonal up to quadratic corrections, whereas the signal acts now along a combination of the x and z directions, We can now finally switch to a frame rotating at the signal frequency ω, and perform a rotating-wave approximation.The resulting QME is which differs from the naive approach to simplify Eq. ( 21) outlined above by the fact that χrot is given in a basis which is rotated by e Ŵ compared to the basis of χ in Eq. ( 21).In this way, we retain information on the Ŝx term in Ĥaxis .
Assuming the populations {p +1 , p 0 , p −1 } to be nondegenerate, we can now diagonalize χrot,ss perturbatively in the signal strength T and find the eigenvectors Their corresponding eigenvalues p +1 , p 0 , and p −1 remain unchanged up to corrections of O(T 2 ).To obtain the lab-frame density matrix ρ, we now undo the transformations ( 27), (22), and (20).The frequency ω appears twice in these transformations, namely, as the small expansion parameter ω/ω 0 in the Schrieffer-Wolff transformation, and as the potentially large rotation angle ωt of the quantization axis.To separate these different roles of ω clearly, we introduce two new parameters for the backtransformation, ω → ω SW in Eq. ( 22) and ω → ω R in Eq. ( 20), such that we can work perturbatively in ω SW /ω 0 while keeping all orders of ω R /ω 0 , given in App. C. For cyclic evolution, τ → 2π/ω, they reduce to γ[P] = arg(z) with If no signal is applied, T → 0, this result reduces to our guess for the GP of an unperturbed quantum limit-cycle oscillator given in Eq. (10).These analytical formulas are in excellent agreement with the numerical results shown in Fig. 3(b) and (d) if ω and T are small, as demonstrated in Fig. 4. We stress that the derivation presented in this section and the result (10) are not specific to the quantum van der Pol oscillator considered in our numerical examples.Any spin-1 limit-cycle oscillator subject to a semiclassical signal of the form (12) will have a density matrix of the form (13) [30].Only the specific formulas for the populations p m and the coherences T c m,m will differ depending on the dissipators Γk stabilizing the limit cycle.The results can thus be directly applied to other limitcycle oscillators (as long as all nonzero populations p k are distinct), and they can be easily generalized to other spin numbers S.

E. Interpretation of the Arnold tongue of the geometric phase
The derivation of the approximate expression (35) for the GP of a quantum limit-cycle oscillator provides an intuitive understanding of the results shown in Fig. 3.As shown in Fig. 1(a), the quantization axis slowly rotates on the surface of a cone.In the absence of an external signal, T = 0, each eigenstate |S, m traces out a path with a solid angle 2π(1 − cos α)|m| subtended from the origin of the Bloch sphere, as shown by the red curve in Fig. 1(b).The time-dependent signal Hamiltonian Ĥsig (t) in Eq. ( 18) tries to tilt the states |S, m away from the instantaneous quantization axis and causes them to precess, as shown by the blue curve in Fig. 1(b).The action of the drive is counteracted by the dissipative limitcycle stabilization mechanism, which attempts to relax the system to a state without any coherences in the basis defined by the instantaneous direction of the quantization axis.In a reference frame that corotates with Ĥsig (t) about the instantaneous quantization axis [see Eq. ( 27)], the density matrix thus acquires constant coherences, as shown in Eq. ( 13).
Since their magnitude increases with the signal strength T and decreases with increasing detuning |∆|, the GP changes in an Arnold-tongue-like region around resonance.
In the case of cyclic evolution shown in Eq. ( 35), only the imaginary part of the coherences modifies the phase factor of each eigenstate |φ m .For noncyclic evolution, the coherences also change the overlap φ m (0)|φ m (t) between eigenstates, see App.C for details.In both cases, however, the functional dependence of γ[P] on the coherences T c m,m is more complex than the simple sum of coherences encountered in the synchronization measure S(ρ).Therefore, the suppression of the Arnold tongue of the GP for parameters in the synchronization-blockade regime must have a different origin than the destructive interference of coherences that causes the Arnold tongue of synchronization to vanish.
The surprising disappearance of the Arnold tongue of the GP in Fig. 3(d) can be traced back to a more general suppression of coherences in the synchronization blockade regime.As shown in Eqs. ( 14) and ( 15), the coherences are in general different functions of the gain and dissipation rates γ g and γ d .Focusing on the resonant case ∆ = 0, we find where we introduced the ratio γ = γ g /γ d and considered the limit of large γ.For generic values of γ g and γ d , the coherences T c +1,0 and T c 0,−1 will have different magnitudes and will not interfere destructively.However, the expressions in the first lines of Eqs. ( 36) and (37) show that the coherences tend to zero with increasing dissipation rates γ d or γ g because the limit-cycle stabilization scheme dominates over the influence of the signal.Therefore, one can find specific combinations of dissipation rates for which the destructive interference occurs, e.g., by increasing the ratio γ g /γ d .This fine-tuning of the ratio γ g /γ d , however, comes at the cost of an overall reduction of the magnitude of the coherences (about an order of magnitude for our parameters), which leads to a slower change of the GP when T or ∆ are changed.This effect causes the disappearance of the Arnold tongue of the GP in Fig. 3(d).

V. CONCLUSION
Based on the kinematic approach to the geometric phase (GP) proposed by Tong et al. [26], we have developed a numerically stable algorithm to calculate the GP in an open quantum system.We used it to demonstrate the existence of a GP in a spin-1 implementation of a quantum vdP oscillator whose quantization axis is slowly rotated.We have shown that if the quantum vdP oscillator is synchronized to an external signal, the GP plotted as a function of the detuning and signal strength exhibits a structure similar to the Arnold tongue of synchronization: the GP changes strongly in a roughly triangular region around resonance whose width increases with increasing signal strength.Surprisingly, this Arnold tongue of the GP vanishes if the system is in a parameter regime where an interference-based quantum synchronization blockade occurs.
These striking similarities between the Arnold tongue of the synchronization measure S(ρ) and the structure of the GP naturally lead to the question if there is a deeper connection between GPs and quantum synchronization.For instance, could quantum synchronization be an indicator of a nonzero GP or vice versa?In general, a nonzero GP does not imply that a quantum system is synchronized because, even in the absence of an external signal, the quantum vdP limit-cycle oscillator shows a GP (see Sec. IV B) but it is clearly not synchronized.Moreover, GPs occur even in completely unitary evolution [1], i.e., in quantum systems that are no limit-cycle oscillators at all.
Conversely, the numerical data presented in Sec.IV C suggests that a nonzero synchronization measure S(ρ) could be an indicator of changes in the GP relative to its value in an unperturbed limit-cycle oscillator.
Using perturbation theory in the small frequency ω of the rotation of the quantization axis and in the small signal strength T , we have derived an approximate analytical expression for the GP, which reveals that the mechanism leading to the suppression of the GP (namely, a suppression of the coherences compared to a regime of regular synchronization) is different from the mechanism leading to a suppression of the synchronization measure (namely, destructive interference of the coherences).Despite these differences, the synchronization measure S(ρ) and the GP show qualitatively the same behavior in a quantum vdP oscillator.It is an exciting direction for further research to understand if this is a generic feature that holds for arbitrary limit-cycle oscillators and external signals.Since the assumptions in our derivation of a perturbative analytic formula for the GP in Sec.IV D are very general, the same technique could be applied to other quantum limit-cycle oscillators to address this open question.
Kepler et al. [27,28] developed a general approach to the GP in classical limit-cycle systems, but the deformations of the classical limit cycle they analyzed differ from the rotation of the quantization axis we considered here.It would therefore be interesting to connect and compare these results by analyzing the classical equivalent of a quantum vdP oscillator whose quantization axis slowly rotates.
Appendix A: Interferometric measurement of the geometric phase In this appendix, we comment on the possibility to measure the GP of a mixed state undergoing nonunitary evolution in an interferometric measurement.
In the case of a pure initial state and unitary evolution, the GP is uniquely defined and can be measured in an interferometric measurement, as shown in Fig. 5.The first beam splitter transforms the initial pure state into a superposition of pure states propagating along the two arms of the Mach-Zehnder interferometer (MZI).The state in the upper arm undergoes adiabatic unitary evolution along the desired path P and acquires a GP γ[P], whereas the state in the lower arm acquires only a controllable U (1) reference phase χ.Having passed through the two arms, the states interfere at the second beam splitter and the probabilities p 0 and p 1 of finding the system in the two output ports are measured.These probabilities depend on the relative phase γ[P] − χ between the two arms and show an interference profile of the form The dynamical phase can be eliminated from the measurement by enforcing the parallel-transport condition in the upper arm or by choosing a specific path along which the dynamical phase vanishes.Other possibilities include cancelling the dynamical phase using spin-echo techniques, or comparing different paths where the relative signs between the dynamical and GPs differ.One can then use the controllable phase χ to map out Eq. (A1) and determine the GP γ[P].
The GP of a mixed state is defined via a purification of the state in a larger Hilbert space [21][22][23].Since a given mixed state can be purified in many different ways, one may worry that the GP can no longer be uniquely defined.However, Sjöqvist et al. [22] showed that, for mixed states undergoing unitary evolution, the GP can still be measured in the MZI setup shown in Fig. 5 and turns out to be the statistical average of the GPs of the pure eigenstates of the density matrix.Such an interferometric measurement has been experimentally demonstrated in an NMR system [57].
In the following, we calculate the measurement probabilities of a mixed state in a MZI undergoing nonunitary evolution and show that the interferometric measurement of the GP of mixed state cannot be extended to nonunitary evolution.
Analogously to the treatment of the unitary case in [22], we model the MZI in a combined Hilbert space H B ⊗ H S , where H B = {|0 B , |1 B } encodes the two where the Hamiltonian Ĥ0 and the Lindblad operators Γ1 , Γ2 are defined in Eqs. ( 4) to (6).Hence, the visibility decays proportional to e −min(γg,γ d ) τ .For large damping and gain rates or for long times τ , the visibility of the interference pattern tends to zero, ν → 0, and detection of any phase shift becomes impossible.Moreover, γ(τ ) has a very different form than the definition of a GP in nonunitary evolution, Eq. ( 2), such that γ(τ ) = γ[P].
In conclusion, for mixed states and nonunitary time evolution, the interferometric measurement considered here cannot be used to determine the GP.Since the GP for nonunitary evolution is defined via a purification procedure [21,26] and since the choice of the purification matters [32,58], one has to enforce this purification during the interferometric measurement.This implies that one has to perform an interferometric measurement us-ing unitary evolution of the combined system and ancilla [58][59][60][61], which is experimentally very demanding for large quantum systems.
The phase γ(τ ) defined in Eq. ( A17) is related to a definition of geometric phases in open quantum systems using a quantum-jump unraveling of the dissipative dynamics [62].A geometric phase can then be defined by adding up the phase changes at quantum jumps and during nonunitary time evolution between quantum jumps, but the value of the phase still depends on the chosen unraveling [63].In our calculation, the effective timeevolution operator Ũ S eff (τ ) describes the dynamics of the system between two quantum-jump events and γ(τ ) thus corresponds to the special (and rare) trajectory where no jumps occur during the entire duration τ .where θ τ = arctan e −Γτ tan θ 0 + π mod π.For a total evolution time τ = 2π/η and cos θ 0 ≥ 0, this result simplifies to Eq. ( 21) of [26].With this analytical formula for the GP at hand, we benchmarked the accuracy of our numerical algorithm described in Sec.III.We solved Eq. (B1) numerically, calculated the GP using Alg. 1, and compared the result with Eq. (B3).For all tested parameters η and Λ, and for all initial conditions r, we observed quartic convergence in the number of time steps, similar to the result shown in Fig. 6.

FIG. 3 .
FIG. 3. (a)Arnold tongue of a quantum van der Pol oscillator subject to a semiclassical signal given by Eq. (12).The limit-cycle oscillator is synchronized to the external signal in the bright region where the synchronization measure S(ρss), defined in Eq. (17), is nonzero.The plot is symmetric about the ∆ = 0 axis and we only show the left half.The dissipation rates have the ratio γg/γ d = 0.5.(c)The right side of the same plot in a parameter regime where an interference-based quantum synchronization blockade occurs (γg/γ d ≈ 2.84).(b) Plot of the GP for the same parameter values as in (a) and nstep = 10 4 , which shows a strikingly similar Arnold-tonguelike structure.Again, the plot is symmetric about the ∆ = 0 axis.(d) For the interference-based quantum synchronization blockade parameters of (c) and nstep = 10 6 , this structure in the GP disappears.For all four cases, the remaining parameters are ω0/γ d = 1, τ ω0 = 200, φ = 0, α = π/4, ω/γ d = 0.05, and ω = ω0 + ∆.Note that the color scale in (b) and (d) is periodic since the GP is a 2π-periodic quantity (unlike the synchronization measure S(ρss)).

2 )
FIG. 4.Linecuts through the different regions shown in Fig.4(data points), compared with the approximate analytic expressions for the GP from App. C (solid lines).The red dots and lines correspond to the synchronization-blockade parameters of Fig.4(d) and are magnified by a factor of 10 2 .The blue dots and lines correspond to the parameters of Fig.4(b).

1 FIG. 5 .
FIG. 5. Mach-Zehnder interferometer (MZI) setup to measure the GP of a mixed state undergoing unitary evolution.The state in the upper arm undergoes adiabatic unitary evolution along the path P and acquires a GP γ[P] whereas the lower arm acquires only a U (1) reference phase χ.The two detectors (green) measure the probabilities of the state exiting in the two output ports of the second beam splitter.

FIG. 6 .
FIG. 6. Deviation of the numerically calculated GP from the exact solution for a qubit undergoing pure dephasing, modeled by the quantum master equation (B1).The GP is computed using Alg. 1 for a fixed evolution time τ = 2π/η and different numbers nstep of time steps.The dashed line indicates a 1/n 4 step scaling.The parameters are Λ/η = 0.2 and θ0 = π/4.