Resilience of $\mathcal{PT}$ symmetry against stochasticity in a gain-loss balanced oscillator

We investigate the effects of dichotomous noise added to a classical harmonic oscillator in the form of stochastic time-dependent gain and loss states, whose durations are sampled from two distinct exponential waiting time distributions. Despite the stochasticity, stability criteria can be formulated when averaging over many realizations in the asymptotic time limit and serve to determine the boundary line in parameter space that separates regions of growing amplitudes from those of decaying ones. Furthermore, the concept of $\mathcal{PT}$ symmetry remains applicable for such a stochastic oscillator and we use it to distinguish between an underdamped symmetric phase and an overdamped asymmetric phase. In the former case, the limit of stability is marked by the same average duration for the gain and loss states, whilst in the the latter case, a higher duration of the loss state is necessary to keep the system stable. The overdamped phase has an ordered structure imposing a position-velocity ratio locking and is viewed as a phase transition from the underdamped phase, which instead displays a broad and more disordered, but nevertheless, $\mathcal{PT}$ symmetric structure. We also address the short time limit and the dynamics of the moments of the position and the velocity with the aim of revealing the extremely rich dynamics offered by this apparently quite simple mechanical system. The notions established so far may be extended and applied in the stabilization of light propagation in metamaterials and optical fibres with randomly distributed regions of asymmetric active and passive media.


Introduction
Given that physical systems are in general not conservative but rather tend to dissipate energy, some external forces are always necessary in order to reactivate their dynamics. There exist many systems, such as the simple pendulum, a motor engine or electric circuit that undergo transitions between states in which energy is gained and dissipated. If the gain is tunable enough in order to compensate the loss, then the resulting device simulates perfectly a conservative system and thus preserves the time reversal symmetry T but for many reasons, mostly technical (regulation or automatism), the compensation is not always perfect so that the loss and gain have to be treated separately.
However, even an imperfect control of the energy balance can result in extraordinary new properties. Indeed, let us consider a system that switches between two possible states, one in which energy is gained and the other where energy is lost. It is clear that such a system breaks the time reversal invariance T and looses its nice property of energy conservation. What is more obvious is that it is not invariant under the operation P which consists in a swap between the loss and gain states. Yet, in order to preserve some symmetry, one idea is to combine the two operations and impose on the non-conservative system a weaker requirement of invariance under the PT transformation.
Works involving PT symmetry have been initiated in the context of quantum mechanics using non Hermitian Hamiltonians [1] in which P refers to the parity operator, which reverses the position. The investigation of a generic class of PT symmetric Hamiltonians has shown that their energy spectrum remains real below a critical point but becomes imaginary beyond it, manifesting a transition to new 'exotic' quantum states. Subsequent works [1][2][3][4] have led to a reformulation of the use of these concepts in the framework of a non-Hermitian model involving only two quantum states for which the parity operator corresponds to the swap operation between the two states.
Since these seminal works [1,2,[5][6][7], this new field has emerged in other contexts such as classical optics and electrical circuits in order to better understand the interplay between active and passive transmission, but also in tight binding systems [8][9][10]. In optical fibers, the simultaneous use of active and passive components displays very interesting properties such as transient wave amplification in an array of coupled waveguides with an arbitrary space distribution of gain and loss [11]. Furthermore, there are experiments which demonstrated that PT symmetric materials can exhibit power oscillations, non-reciprocal light propagation and tailored energy flow [6,7,12]. In addition, the existence of giant amplifications is predicted, meaning that a passive medium may be helpful to enhance the gain effect of an active medium [13]. Similar problems were studied in another experiment with a pair of coupled oscillators in the form of an LRC circuit [12]. Instead of considering a single oscillator that switches between gain and loss states, the authors of [12] examined an electronic dimer made of two coupled oscillators, one with gain and the other with loss. The experiment succeeds in displaying all the phenomena encountered in systems with generalized PTsymmetries.
In order to understand the basics of a PT -symmetric gain and loss process, a very simple one-dimensional harmonic oscillator was considered. The prototype model consisted of two separate states of frictional and gain forces linearly proportional to the velocity that alternate periodically in time [14]. This model contains only the oscillator frequency, the damping coefficient and the alternating period as parameters. Quite remarkably, it provides a complete analysis with a phase diagram that distinguishes the stable from the unstable regimes according to the parameter values.
However, as the dynamics might be even less controllable in the presence of randomness, a natural question arises on how it affects, or rather breaks the PTsymmetries. In this paper, we investigate how an effective PT symmetry persists in the presence of dichotomous noise introduced by replacing the fixed time periods with random intervals in the simple generic model developed in [14]. More precisely, the oscillator switches randomly in time between a damping state in which energy is dissipated or lost and an anti-damping state in which energy is accumulated or gained. This oscillator can represent one electromagnetic mode in a cavity that is amplified randomly in order to compensate the losses.
There exist earlier studies of the effects of random damping on the stability of harmonic oscillators [15,16]. The stabilities of the first two moments of the oscillator position and velocity have been analyzed, but only for uncorrelated Gaussian and colored noise and not for dichotomous noise. In this context, we also mention the work on the PT symmetric coupler in [17] with Gaussian white noise, where amplification occurs despite the perfect balance of gain and loss. In contrast to these previous works, besides determining the moments, we are also able to characterize in the asymptotic limit the exact nature of the probability distribution generated by the random noise and thus predict the oscillator energy distribution. Furthermore, we also introduce an alternative notion of stability based on the energy logarithm of the system which we motivate through the properties of the probability density function of the state of the system. In addition to the PT symmetric states, we also found regimes in which this symmetry is broken even though the average durations for the loss and gain states were equal. This observation confirms the known statement that energy amplification occurs even when the system is predominantly dissipative over time [11,17] and can be formally established using a mathematical framework based on the master equation. Finally, we succeed in pointing out the analogy with phase transitions in thermodynamics, in which beyond a certain critical value, the stochastic oscillator breaks its PT symmetry towards an ordered phase. This paper is organized as follows. In section 2, we formulate the stochastic oscillator problem in terms of the master equation and define an asymptotic stability criterion. In section 3, we present the results for both simulations and analytics and show how they can be related to a phase transition. Section 4 concerns a more restricted stability criterion involving the position and velocity averages. We discuss the short time behaviour and stability involving higher order moments of the velocity in the limit of zero frequency oscillation in section 5 before ending with the conclusion in section 6.

General consideration
We consider a simple harmonic oscillator that randomly switches between a damping and anti-damping phase. The equation of motion of such an oscillator with natural angular frequency ω has the form with a time-dependent damping coefficient θ(t) that can take only two constant values, either +γ or −γ, i.e. θ(t) is a piece-wise constant function in the form of dichotomous noise (see figure 1). In the former case, the oscillator undergoes damping and therefore loses energy (loss state) whilst in the latter case it gains energy (gain state). We introduce stochasticity through the damping coefficient so that the system fluctuates between the gain (g) and loss (l) phases with residence times τ . For the gain or loss phase the residence time is sampled from a distinct exponential distribution of the form τ −1 g/l exp(−τ /τ g/l ). The two phases are therefore characterized by well defined average residence times, τ g in the case of gain and τ l in the case of loss.

Master equation
In order to deal with the stochastic system, we begin by examining the time evolution of p(x, v, t), the probability to find the oscillator in position x with velocity v =ẋ at time t. To this end we write down the master equation of such a system, keeping in mind that the oscillator can also be in any of the two states g or l. Therefore, we have the following system of coupled partial differential equations: is the probability for the oscillator to be in the gain/loss state at time t. Furthermore, the probability is conserved so that P g (t) + P l (t) = 1. The first term on the left-hand-side of (2) is the deterministic Liouville term, whilst the second one is the gain/loss term that arises from the nonconservative nature of equation (1). The term on the right-hand-side describes the stochastic switching rate of the oscillator between the gain and loss phases.
In as such, one would have to solve the coupled pair of equations in (2) for p g (x, v, t) and p l (x, v, t) in order to completely characterize the stochastic system. However, such a task is extremely heavy and unnecessary for our purpose. As mentioned in the introduction, we are interested in the stability of the stochastic oscillator and a reliable criterion for it. A first simplification arises by noticing that the action-angle variables are more suitable for handling the master equation. These polar-type coordinates lead to a separation of variables in the master equation in (2) (see Appendix A.5 for details).
In the optics terminology, J and ϕ correspond respectively to the amplitude and phase while x and v correspond to the quadrature components. Subsequent Laplace and Mellin transforms allow us to eliminate the time and J derivatives in the resulting master equation. Indeed, these transformations defined respectively asf simplify the master equation into: d dϕ (ω ∓ γ sin 2ϕ) + s ∓ 4γk cos 2 ϕ p g/l (k, ϕ, s) = ∓ p g/l (k, ϕ, s) τ g −p g/l (k, ϕ, s) τ l +p g/l (k, ϕ, t 0 ), (5) where p g/l (J, ϕ, t 0 ) is the starting distribution (initial condition), which we assume to be a delta function. Still however, the last form cannot be solved analytically exactly. Nevertheless, essential information can be derived about the asymptotic time limit of the solution (see Appendix A.1 and Appendix A.2). For large times, we deduce indeed that the variable ln J follows a normal distribution by showing that any moment of the cumulant expansion of ln J scales linearly with time (see Appendix A.3). As a consequence, the average value ln J depends linearly on time and the relative square root variance has the scaling δ 2 ln J / ln J → 1/ √ t so that asymptotically ln J becomes a deterministic variable (see Appendix A.4). On the other hand, the angle ϕ remains generally distributed over a broad value range. The numerical simulations of the evolution of an ensemble of stochastic oscillators confirm these theoretical results: Figure 2 shows the linear time dependence for the average around which the square root variance remains small in comparison; figure 3 shows the histogram of the distribution of ln J, which converges to a Gaussian in the asymptotic time limit. In this particular example ω = 1, γ = 0.9, τ g = 5 and τ l = 1. The average is taken over an ensemble of 10 4 oscillators, each with the initial condition (x 0 , v 0 ) = (1, 1) at t = 0. The positive slope of ln J(t) for those particular parameter values implies that the system is unstable in the asymptotic limit. The dashed black line, whose slope was obtained from (6), corresponds to the theoretical result.

The stability criterion
In order to assess the stability of the dynamics of the stochastic oscillator, we can use the following results established in Appendix A.4. If the asymptotic marginal distribution P g/l (ϕ) = lim t→∞ ∞ 0 p g/l (J, ϕ, t) dJ exists, then the asymptotic constant η associated to the linear evolution of the first moment is given by: This real-valued constant is the basis of the stability criterion that we shall employ in the next section. A positive η corresponds to a diverging first moment of ln J implying that the system is unstable while a negative η corresponds to a stable system. In order to apply the stability criterion defined in (6), we need to determine P g (ϕ) and P l (ϕ) asymptotically.
3. Stability results of the stochastic oscillator in the asymptotic limit

Simulated results compared with the theory
We used numerical simulations of the evolution of an ensemble of stochastic oscillators from which we extract the asymptotic first moment of ln J(t) and thus determine η.  The results are plotted in figure 4. The green area corresponds to negative values of the slope of the ensemble average of ln J whilst the red corresponds to positive ones. For γ < ω, the two regions are separated by the line of symmetry τ g = τ l . This result corresponds to what one might expect -if the amount of time spent in the gain state is on average longer than in the loss state, then the average value of the energy diverges over time. On the other hand, if the system on average spends more time in a loss state, then its average energy decays to zero over time. What comes as a surprise is that when γ > ω the system energy can diverge even when τ l > τ g . In order to see this better it is worth theoretically studying the dependence of η in terms of γ.
In the limit of large τ g and τ l , we calculate explicitly the formula (6) using the expressions for P g/l (ϕ) derived in Appendix A.5 and Appendix A.6 to obtain the simple analytic forms: From (7), a necessary condition for stability is that τ g < τ l , independent of the values of ω and γ. Furthermore for γ < ω, the asymptotic expression is in good agreement with the simulation results, whereby the oscillator is at the edge of stability for τ g = τ l .  (6). The green color indicates the stable region, where the average of ln J is negative and the red color indicates the unstable region, where it is positive. Panels A and B: In this case γ < ω. It can be seen that the line of symmetry τ g = τ l separates the stable region from the unstable. Panel C: This is the critical case where ω = γ. Since the system is not in the asymptotic regime, the symmetry line τ g = τ l does not separate the two regions perfectly. Panel D: An example of the case where γ > ω. It shows that the symmetry line τ g = τ l is well below the line that separates the two regions. Consequently there exist cases where τ g is well below τ l and yet the system is still unstable. The inset shows that for small values of τ g and τ l symmetry is regained, as discussed at the end of Appendix A.6. In all four cases an ensemble of 1000 oscillators that evolved up to t = 300 were used, each oscillator having the initial condition (x 0 , v 0 ) = (1, 1) at t = 0.  (7). The INSET shows the dependence of η on γ in the intermediate τ g , τ l limit. The blue curve is the result of a simulation with τ l = τ g = 5. In this limit, symmetry is broken even for γ < ω. In all simulations an ensemble of 10 5 oscillators evolving up to t = 300 were used, each oscillator having the initial condition Asymptotically, the energy logarithm ln J can be considered as a deterministic quantity, which does not decay nor diverge when there is a perfect balance between gain and loss, when τ g = τ l . On the other hand, when γ > ω a gain-loss balance (τ g = τ l ) does not induce stability in the system. On the contrary, the system can remain active with a growing or constant energy even when the loss states dominate over the gain states, i.e. when τ l > τ g . This is illustrated by the green curve in figure 5; there is a value of γ above which the system's energy diverges no matter how large τ l becomes compared to τ g .

The effective transition to a broken PT symmetry phase
We can more precisely formalize what was said above by investigating further the system by analyzing the symmetry properties of the probability density functions p g and p l . We define the space reflection, or parity operator P as the exchange between the gain and loss probability densities of ϕ. In other words, P has the effect of swapping g and l so that we have the exchange P g ↔ P l . A similar definition is encountered in [2,3] where for simplicity the real space is represented by a two-valued position (let us say ±1) and the parity operation represent the exchange between +1 and −1. This redefined operation has been used as a swap operation between two quantum states and subsequently for the swap between the gain and loss states in [14]. We can also include the time-reversal operation T , where t → −t, v → −v and x → x so that ϕ → −ϕ. If we apply both operations at the same time, the distributions P l (ϕ) and P g (ϕ) remain invariant so that the PT symmetry is fulfilled. The ensemble of stochastic oscillators is effectively PT -symmetric in the so-called underdamped regime when γ < ω and under the condition that τ g = τ l , although neither the stochastic equation (1) nor the master equation (2) obey such a symmetry. Indeed, the phase probability densities P g (ϕ) and P l (ϕ) obtained by simulation and represented in the first graph of figure 6 have a mirror symmetry with respect to the origin (ϕ → −ϕ). For comparison, the analytic expression for P g/l (ϕ) is obtained by solving the master equation in (2) in the large time limit (numerical integration and asymptotic expression in the large τ g/l limit in Appendix A.6).
However, this symmetry is only effective if we compare values of ln J up to its square root variance given that ln J keeps diffusing normally. But if we compare the square root variance relatively to any non trivial average of ln J, it shrinks to zero in the large time limit. Therefore these considerations have only a strict sense in the asymptotic limit viewed here as the analog of the thermodynamic limit where the concept of a large particle number of a thermodynamic system is replaced by one of large time, and where the so-called normal quantities are the average and the variance of ln J that both scale linearly with time (see Appendix A.3 and Appendix A.4).
The mirror symmetry of P g/l (ϕ) is maintained only for γ < ω. Once this condition is no longer satisfied, the distribution of ϕ initially broadly distributed in the symmetric case condenses by forming two delta-like peaks. It is in this sense that the system becomes deterministic once γ is greater than ω. At the same time, however, the mirror symmetry of the probability densities is broken as can be seen in the second graph of figure 6 leading to a PT -symmetry violation. In the limit of large τ g and τ l we are able to calculate a simple expression that determines the values ϕ g/l at which the two peaks in P g/l (ϕ) occur (see Appendix A.6 for details): The new "phase" obtained is ordered in the sense that it corresponds to a ratio locking of the velocity over the position with different fixed values for the gain state and the loss state. This result may also been obtained more intuitively by noticing that for γ ≥ ω the oscillator is damped with no oscillations. It corresponds to the overdamped regime as opposed to the underdamped regime where oscillations persist. We can indeed solve (1) using the ansatz x(t) = e λt x(0) and find that λ = ±γ ± γ 2 − ω 2 is real only in the overdamped regime. Hence, in the large time limit only one eigenvalue is dominant and therefore using ωx(t)/v(t) = ω/λ for the dominant eigenvalue, we recover (8) accordingly. Such a relation could not have been used in the underdamped regime since the phase of oscillations would have randomized the trajectories. We interpret this observation as a phase transition from a disordered state to an ordered state with symmetry breaking in analogy to what happens in phase transition phenomena in thermodynamics. It can therefore be concluded that the PT -symmetry breaking occurs at the point of critical damping (γ = ω). In analogy to the Ising model [18], we start from a symmetric state with no ordering above a critical point, the broad angle distribution in our case (or the spin distribution in the magnet), and go towards a broken symmetry state with a well defined order with two possible opposite angle values e.g. ratio locking (or a well defined value of spin).
We end this section by adding that the symmetry breaking established in the limit of large τ l,g is essentially valid also in the intermediate regime, despite a little bias (τ l > τ g ) for γ around ω shown in panel C of figure 4 and in the inset of figure 5. The symmetry is totally restored, however, in the limit of small τ l,g whatever the value of γ as can been seen from the inset in figure 4.

Dynamic equations for the position and velocity average
The stability criterion obtained in the previous section does not mean that all physical quantities of interest are stable. This statement can be illustrated by considering the evolution of the first moments of p g and p l . The fact that ln J is stable does not necessarily mean that averages involving position and speed average are stable. Indeed, if f (ϕ) is a function to average, from the Feymann-Gibbs inequality, we deduce: On the contrary, the stability of position and velocity averages implies stability of ln J. Therefore, there exist additional requirements that enhance the stability of the stochastic oscillators, adding to the richness of their dynamics. By multiplying both equations in (2) by x and v and integrating over the entire space and all the velocities, we obtain a system of coupled first order differential equations: -π -π/2 0 π/2 π ϕ P g (ϕ) P l (ϕ) Figure 6. The distribution of the phases P g (ϕ) and P l (ϕ). TOP: The system exhibits mirror symmetry and dispersed phases in the case where γ < ω and τ g = τ L . In this example γ = 0.5, ω = 1 and τ g/l = 20. An ensemble of 10 6 oscillators was allowed to evolve up to t = 1500. The black curve corresponds to the numerical solution of the analytical result in (A.31). BOTTOM: The symmetry is broken when γ > ω and the phases become localized. When γ becomes larger than ω a transition occurs from a disordered to an ordered phase with a velocity-position ratio locking. In this example γ = 0.5, ω = 1.5 and τ g/l = 10. An ensemble of 10 6 oscillators was allowed to evolve up to t = 300, each oscillator having the initial condition (x 0 , v 0 ) = (1, 1) at t = 0. where (11) The linear dynamical system described here is of the formż = M · z and therefore an asymptotically stable condition is verified when all of the real parts of the roots of the characteristic polynomial associated with M are negative. It is straightforward to determine the four eigenvalues which we shall denote as λ = (λ 1 , λ 2 , λ 3 , λ 4 ) (see Appendix B for details). Consequently the stability is also asymptotic in time.
In a way similar to what was presented in the previous sections, the dynamical system governing the first moments can be separated into two regimes, namely, the case where γ < ω and γ > ω. We proceed by examining the stability conditions for such a system and without loss of generality we assume that ω = 1. We show that in this case also, depending on the value of γ, it is possible to have a situation in which the average quantities we studied diverge even though the system, on average, spends more time in the loss state than in the gain state. In contrast to the results obtained regarding the stability of ln J in the previous section, the average values in (10) diverge when τ l = τ g since there exists at least one eigenvalue that is positive for any combination of the other parameters (see Appendix B for details).

Underdamped case: γ < ω
In this case, all four eigenvalues λ i are complex and they come in conjugate pairs (see Appendix B). The real parts of two of them are equal, Reλ 1 = Reλ 2 , and negative for all combinations of (τ g ,τ l ) whilst the real parts of the other two, which are also equal (Reλ 3 = Reλ 4 ), can be either positive or negative. Since the imaginary parts of all four eigenvalues are always non-zero, all the solutions of the system (10) are oscillatory, whether they decay or grow, for any value of τ l and τ g . For a fixed γ we numerically determine the regions in the τ l − τ g parameter space for which every eigenvalue of the system has a negative real part. This corresponds to the stable region and is indicated in green in panels A and B of figure 7. Moreover, the relation lim τ l →∞ λ i = γ − 1/τ g for i = 3, 4 derived using B.2 imposes the critical value τ * g = 1/γ beyond which stability is never attained no matter how large the value of τ l is.

Overdamped case: γ ≥ ω
Once the stochastic oscillator is overdamped, more of the eigenvalues can attain positive real parts, rendering the dynamics more elaborate. In particular, in contrast to the underdamped case, Reλ 2 can also be positive, depending on the values of τ g and τ l . Nevertheless, the stability does not differ qualitatively from the underdamped case. The boundary that separates the stable from the unstable region continues to shift towards lower values of τ g when γ increases further. Another interesting property deduced from the eigenvalues is that in the overdamped case there exist oscillatory solutions to (10) in contrast to the simple damped harmonic oscillator that is monotonically damped under those conditions. These oscillations occur when at least one of the eigenvalues has a non-zero imaginary part. This is the case for values of (τ l , τ g ) for which |Imλ 1 | + |Imλ 2 | + |Imλ 3 | + |Imλ 4 | = 0. These results are presented in panels C and D of figure 7 where, although the monotonic motion is predominant in the overdamped case, there persists a region in the τ g − τ l plane for which the solutions are oscillatory.

Full time analysis at zero frequency
Until now, we have analyzed the asymptotic behavior of the oscillator in the large time limit and also the first moment of the position and the velocity. For the sake of completeness, it remains in principle to discuss about the short and intermediate time regimes and the stability of the higher order moments. A full analysis is beyond the scope of this paper but we can address the particular case of ω = 0, so that the problem reduces to the simpler dynamical equation in which the position coordinate is eliminated: This dynamical equation describes the evolution of a wave function that is amplified or undergoes loss randomly. A similar study where the frequency takes randomly two possible values has been shown useful in the context of superconducting qubits [19]. This type of model has already been used to determine the first passage time at which, for instance, the speed exceeds a critical value [20][21][22][23]. The interest here is to illustrate how the exact solution complements the results obtained so far. We start from the initial condition: p g/l (v, t = 0) = p g/l δ(v − v 0 ) and solve this equation using the variable change: r = ln(v/v 0 )/(2γ) (see Appendix C). The short time analysis differs from the asymptotic analysis because the distribution is not normal anymore. If we start with a well defined speed v 0 and study its subsequent evolution for short time (t τ g/l ), we obtain the spread of the velocity distribution but confined within a cone: where 1 + (x) is the Heaviside function. The short time behavior is characterized by a propagation of a delta distribution along this cone inside which the probability densities develop. In the opposite case for large time (t τ g/l ), we find a normal distribution for the variable r with average and variance: Thus, the stability of the average is satisfied for τ l ≥ τ g in contrast to the previous statements in (7), which always predict instability. The apparent contradiction is resolved by remembering that according to (9) the stability of r does not imply the stability of any moment of the velocity and/or any moment of the position. Indeed for the n th moment, we determine the following more restrictive stability criterion: Therefore, there always exists an order n above which the stability criterion is not satisfied in accordance with the variance in (14), which always increases. However, if we restrict to the v variable only without the position x then the simplified system becomes effectively always PT symmetric when τ l = τ g .

Conclusions and perspectives
We have studied the dynamic evolution of stochastic oscillators subject to dichotomous noise made of alternating gain and loss states random in time and we have unveiled an intimate connection of this non conservative system with PT symmetry. We established a useful criterion that fixes the boundary line between a stable regime with a likely decaying amplitude and an unstable regime with a likely growing one. Although the oscillator evolution becomes more stochastic with time, it is nevertheless possible to effectively define the useful concept of PT symmetry in the asymptotic time limit. In other words, despite the breaking of time reversal invariance due to noise, the oscillator can still remain resilient so as to preserve at least the PT symmetry. Application of this invariance property allows to distinguish between different regimes or phases: a) an underdamped regime (or weakly damping-amplifying oscillator ) for which the boundary lines between stable and unstable regions satisfy this symmetry; b) an overdamped regime (or strongly damping-amplifying oscillator) for which this boundary line becomes asymmetric. We interpret these results in analogy to thermodynamics as a phase transition from a symmetric disordered state consisting of a broad distribution to an ordered state with a restricted distribution imposing a ratio locking of the position over the velocity separately for both the gain and loss states.
To complete the panorama, we also examined the time evolution of the position and velocity averages of the oscillator. We showed that the stability of the oscillator does not necessarily imply bounded dynamics of these averages. It appears indeed that the stability diagrams are more elaborate, illustrating the much richer structure of this apparently very simple system. For instance, despite the absence of oscillations in the deterministic case, in the overdamped regime the presence of a random gain can restimulate them. Higher order moment analysis together with a study of the short and intermediate time limits confirm this broad range of different regimes with different physics such as the cone-like propagation of the velocity distribution.
The formalism developed here in the particular case of a stochastic oscillator is quite general and may be applied to other situations where dichotomous noise is present such as the stabilization of light propagation in metamaterials and optical fibres with random regions of asymmetric active and passive media [7]. After reducing the master equation (2) to a pair of coupled ordinary differential equations by means of integral transforms, we write the master equation (5) under a matrix form with p g/l (J, ϕ, t 0 ) as the initial distribution: Note that because of the inversion symmetry x → −x and v → −v, the stochastic oscillator is invariant under the transformation ϕ → ϕ + π. Therefore, we can restrict the angle to the interval ϕ ∈] − π/2, π/2]. Writing the solutionp(k, ϕ, s) as a linear combination of the eigenfunctionsp i (k, ϕ) of M with eigenvalue −s i we obtain In that case, the master equation becomes Multiplying from the left with the eigenfunctionp † j (k, ϕ) of the adjoint problem with eigenvalue s j and using the bi-orthogonality property we obtain Taking the inverse Laplace transform, we arrive at Resilience of PT symmetry against stochasticity in a gain-loss balanced oscillator 19 Therefore, the inverse Laplace transform of (A.5) giveŝ In the asymptotic time limit where t → ∞, assuming that s 0 is the dominant eigenvalue in the sense that the real component Re s 0 has the lowest value among all eigenvalues, we simplify the dynamics into:

. Perturbation theory
We progress further using the perturbation expansion applied to the eigenvalue problem: In the leading order, s (0) i is equal to the eigenvalue of the unperturbed operator M (0) while, in the first order, we find using (A.7) the first correction: i,l ) T is the solution to the unperturbed eigenvalue problem (i (th) eigenvector) that involves only M (0) . Since the operator is not self-adjoint, we must usẽ p (0) i , which is the solution to the zeroth order adjoint problem defined by where the adjoint operator has the form It is easily verified that in the unperturbed case, the adjoint problem has the trivial solutionp (0) 0 = (1 1) T with the eigenvalue s (0) 0 = 0. This trivial eigenvalue is the most dominant because an eigenvalue with a real positive value would lead to a probability that is not conserved in time. In fact, this eigenvalue is precisely associated to the probability conservation. Therefore, up to first order in k, the dominant solution of the eigenvalue problem (A.14) is rewritten using (A.2) more simply as: −π/2 4γ cos 2 ϕ (p g (0, ϕ, 0) −p l (0, ϕ, 0)) dϕ π/2 −π/2 (p g (0, ϕ, 0) +p l (0, ϕ, 0)) dϕ . (A.17) Appendix A.3. Evolution of the central moments of ln J We start with the Mellin transform shown in (4) and write it down as a characteristic function of the total probability for the variable ln J: If we write down the exponential as a power series and integrate over ϕ, we obtain an expression with the moment generating function of ln J on the right-hand-side: Therefore, in order to obtain the evolution in time of the n th central moment of ln J, we take the n th derivative of the logarithm of the moment-generating function above (i.e. derive the cumulant generating function) with respect to k According to perturbation theory, we can expand the dominant eigenvalue s 0 in powers of k. By substituting this expansion into the Mellin-transformed probability distribution obtained in the asymptotic time limit in equation (A.11) we obtain: Finally, for the n th central moment with respect to the quantity ln J, we identify simply in the asymptotic time limit: is the n th expansion coefficient in k of the eigenvalue s 0 = ∞ n=0 s (n) 0 k n /n!. Note that the asymptotic result does not depend anymore on the initial condition through A 0 (k).

Appendix A.4. Stability criterion
As already mentioned in Appendix A.2, the master equation can be treated as a perturbed eigenvalue problem with k as the perturbation parameter and s as the eigenvalue to be determined. Given that the dominant eigenvalue for the unperturbed problem is s  −π/2 4γ cos 2 ϕ (p g (0, ϕ, 0) −p l (0, ϕ, 0)) dϕ The linear growth of the first two central moments means that the relative square root variance or the square root variance-to-mean ratio shrinks to zero asymptotically, i.e., provided ln J = 0. This result leads us to conclude that ln J is in general a well defined statistical variable and can be viewed as a deterministic one in the asymptotic sense (see figure 2). Strictly speaking, only the case ln J = 0 has to be considered non deterministic but distributed within the interval of the square root variance.
Appendix A.5. Solution to the unperturbed eigenvalue problem Setting s = k = 0 in (5), we notice that the functions: and X 2 ≡ −γ sin(2φ)(p g +p l ) + ω(p g −p l ). (A.26) satisfy a much simpler system of linear first order differential equations: where, A ± (ϕ) = τ l (ω + γ sin 2ϕ) ± τ g (ω − γ sin 2ϕ) τ g τ l (ω 2 − γ 2 sin 2 (2ϕ)) . (A.28) The first equation, for X 1 , is trivial whose solution is a constant. The second equation, for X 2 , can be solved exactly with the formal solution: Resilience of PT symmetry against stochasticity in a gain-loss balanced oscillator 22 The arbitrary constant of integration C is determined from the condition of a periodic solution, i.e. X 2 (ϕ) = X 2 (ϕ + π) so that we find: Reversing the relations (A.25) and (A.26) in terms ofp g andp l and inserting the results into (A.23), we obtain a simpler expression for the stability criterion: Appendix A.6. Large τ l , τ g limit We can simplify further and solve (A.23) in some particular but relevant cases. In the large τ l , τ g limit, but keeping the ratio τ l /τ g constant, the asymptotic solution developed in Appendix A.5 can be integrated exactly. We have the two cases: 1) Underdamped case: γ ≤ ω In this limit the exponential terms inside (A.29) and (A.30) reduce to unity and the solution simplifies to a constant: where the last line results from integration over ϕ . Inserting this result into (A.31), we obtain after integration the stability criterion: After some straightforward algebra, we find the angle probability distribution: p g = τ g X 1 ω − sin(2ϕ)γ andp l = τ l X 1 ω + sin(2ϕ)γ . (A.34) We immediately see that these expressions are PT symmetric when τ l = τ g . Normalizing these expressions to unity, we fix the constant to X 1 = ω 2 − γ 2 2π(τ g + τ l ) . (A.35) 2) Overdamped case: γ ≥ ω In this case X 2 is also a constant except at the singularities ϕ 0 , which are the zeroes of ω = ±γ sin(2ϕ 0 ). Some of these singularities are dominant in the sense that their weights are much greater than than those of the others. Let us focus our analysis on the interval ϕ ∈] − π/2, π/2]. In order to determine their importance, we notice that around these singularities, the master equation in (5) decouples so that we can make the approximation: (ω ∓ γ sin 2ϕ) d dϕ ∓ 2γ cos 2ϕ + 1/τ g/l p g/l (ϕ) 0, (A.36) We can expand this equation locally around the singularity ϕ 0 to obtain more simply: ∓2γ cos 2ϕ 0 (ϕ − ϕ 0 ) d dϕ + 1 + 1/τ g/l p g/l (ϕ) 0, (A.37) The solution is thenp g/l (ϕ) ∼ |ϕ − ϕ 0 | −1∓(2γ cos(2ϕ 0 )τ g/l ) −1 (A.38) so that the dominant singularities appear for the highest negative power fixed by the condition: ± cos(2ϕ 0 ) > 0. After a little algebra, we find that these singularities impose the ratio locking condition: Therefore we deduce for the probability distribution: p g = τ g τ g + τ l δ(ϕ − ϕ g ) andp l = τ l τ g + τ l δ(ϕ − ϕ l ). (A.40) The relative weight between the two probabilities is determined by noticing that the total probabilities for the gain and loss states should satisfy P g (t)/P l (t) t→∞ = τ g /τ l . Contrary to the underdamped case, these probability expressions are not PT symmetric, no matter the parameter values chosen. Inserting these last results into (A.23), we find after integration Note that a similar reasoning can be done in the small τ l , τ g limit keeping the ratio τ l /τ g constant. In that case, we would recover (A.33) provided τ g is not too far from τ l .
In the symmetric case, where τ l = τ g = τ , the system is never stable since there exists at least one eigenvalue that is positive for any combinations of γ and τ . Precisely, If we assume that λ 4 = α + iβ and make the substitution γ 2 (γ 2 − 1)τ 4 − τ 2 = α + iβ , then by solving for α, it is straightforward to show that Reλ 4 > 0 for every positive γ and τ . Note that γ and τ are real so that α and β cannot be non-zero at the same time.