Fractional Brownian Gyrator

When a physical system evolves in a thermal bath at a constant temperature, it arrives eventually to an equilibrium state whose properties are independent of the kinetic parameters and of the precise evolution scenario. This is generically not the case for a system driven out of equilibrium which, on the contrary, reaches a steady-state with properties that depend on the full details of the dynamics such as the driving noise and the energy dissipation. How the steady state depends on such parameters is in general a non-trivial question. Here, we approach this broad problem using a minimal model of a two-dimensional nano-machine, the Brownian gyrator, that consists of a trapped particle driven by fractional Gaussian noises -- a family of noises with long-ranged correlations in time and characterized by an anomalous diffusion exponent $\alpha$. When the noise is different in the different spatial directions, our fractional Brownian gyrator persistently rotates. Even if the noise is non-trivial, with long-ranged time correlations, thanks to its Gaussian nature we are able to characterize analytically the resulting nonequilibrium steady state by computing the probability density function, the probability current, its curl and the angular velocity and complement our study by numerical results.

Abstract. When a physical system evolves in a thermal bath at a constant temperature, it arrives eventually to an equilibrium state whose properties are independent of the kinetic parameters and of the precise evolution scenario. This is generically not the case for a system driven out of equilibrium which, on the contrary, reaches a steady-state with properties that depend on the full details of the dynamics such as the driving noise and the energy dissipation. How the steady state depends on such parameters is in general a non-trivial question. Here, we approach this broad problem using a minimal model of a two-dimensional nano-machine, the Brownian gyrator, that consists of a trapped particle driven by fractional Gaussian noises -a family of noises with long-ranged correlations in time and characterized by an anomalous diffusion exponent α. When the noise is different in the different spatial directions, our fractional Brownian gyrator persistently rotates. Even if the noise is non-trivial, with long-ranged time correlations, thanks to its Gaussian nature we are able to characterize analytically the resulting nonequilibrium steady state by computing the probability density function, the probability current, its curl and the angular velocity and complement our study by numerical results.

Introduction
The Brownian gyrator model, that consists of a pair of coupled, standardly defined Ornstein-Uhlenbeck processes each evolving at its own temperature, is one of the simplest models exhibiting a non-trivial out-of-equilibrium dynamics and, for this reason, has received much interest in the last two decades. It was originally introduced as a solvable model of dynamics with two temperatures [1], and it was realized only years later that a systematic torque is generated on the object [2], so that it can be seen as a minimal nano-machine that, on average, undergoes steady gyration around the origin. Following this remarkable observation, various facets of the model have been studied in different contexts. In addition to the theoretical calculation of the evolution of some relevant properties (see e.g. Refs. [3][4][5][6][7]), several more general questions have also been addressed including the role of cross-correlations between different degrees of freedom coupled to different thermostats on the production of entropy [8], the directionality of interactions between cellular processes [9], statistical properties of the energy exchanged between two heat baths in electric circuits [10,11] and their spectral properties [12], dynamics in suspensions of two types of particles exposed to different heat baths [13] and in laser-cooled atomic gases [14], as well as the synchronization of out-of-equilibrium Kuramoto oscillators [15]. Furthermore, exerting external forces on a Brownian gyrator has lead to derive a special fluctuation theorem and to identify associated effective temperatures [12,16,17]. The latter setting with external forces was found to be mathematically equivalent to the bead-spring model studied via Brownian dynamics simulations in Ref. [18] and analytically in Refs. [19,20]. Finally, different versions of a Brownian gyrator were studied experimentally, especially from the point of view of fluctuation relations [10,11,21].
All the above-mentioned analyses concentrated exclusively on situations in which the noise acting on both components of the gyrator is thermal and delta-correlated in time. However, nano-machines can be encountered or engineered in a variety of environments including active suspensions [22,23] or complex biological mediums [24] where the noise is athermal and correlated in time. It is thus important to understand the effects of correlated noise on the Brownian gyrator. A first step was made recently in Ref. [25] which considered persistent noise that is exponentially correlated in time and the associated memory kernel for the damping so that, individually, the noise on each component remains thermal. Even for such short-ranged correlations, essential departures from the standard white noise model emerge which raise the question of the effect of other types of noise, athermal and/or with long-ranged time correlations. From a broader perspective, because of its simplicity, the Brownian gyrator model is an ideal playground to understand the effects of different types of noise correlations on the properties of a nonequilibrium steady state.
In the present paper, we tackle the question formulated above by generalizing the Brownian gyrator model to a fractional Brownian gyrator (FBG) in which the noise acting on each component is the so-called fractional Gaussian noise (FGN) [26]. The latter, in fact, represent a family of colored noises parametrized by an anomalous diffusion exponent α such that an overdamped particle in free space subject to a FGN performs the standard fractional Brownian motion [26,27]. When α < 1 the increments are anti-correlated and the dynamics is sub-diffusive whereas for α > 1 the increments are positively correlated which entails a super-diffusive motion. The standard Brownian gyrator model is recovered in the borderline case when α = 1. Our analysis here covers all three situations within a unified framework, so that we can inquire for each type of correlations if they are beneficial or, on the contrary, detrimental to rotation. In any case, rotation can appear only if the symmetry between the two spatial components is broken. Here, to assess independently the different features of the noise, we discuss especially two situations: In the first case, both noises have the same α but different amplitudes while in the second case the amplitudes are equal while the values of α differ. We calculate exactly the probability density function P (x, y; t) and present its explicit form in the limit t → ∞. We further compute the probability current, its curl, and eventually calculate the ensemble-averaged angular velocity. For each quantity, we compare our results to numerical simulations.
The paper is outlined as follows: We first introduce the dynamics and detail the two variants of the model that breaks the symmetry using different parameters in Sec. 2. In Sec. 3 we compute analytically the first moments of the position before giving the full probability density function and considering numerical applications. We further compute some important quantities characterizing the nonequilibrium steady state in Sec. 4: the probability current, its curl and the mean angular velocity. Finally, we summarize our findings in Sec. 5. Details of intermediate calculations are relegated to Appendices.

Definition
Let us first present the standard Brownian gyrator (BG) as introduced in Ref. [1]. We denote by x(t) and y(t) the Cartesian coordinates of an overdamped particle evolving in a potential U (x, y) according to the Langevin equation where γ is the damping coefficient, U (x, y) = (1/2)(k x x 2 +k y y 2 )−uxy is a generic harmonic potential and ξ x (t) and ξ y (t) are additive Gaussian white noises with zero mean. Adopting appropriate space and time units, we can choose without loss of generality k x = k y = 1 and γ = 1 so that the dynamics rewritė with noise variance We are interested in the case when the particle is trapped, which requires |u| < 1. When T x = T y , the system is out of equilibrium and exhibits a persistent rotation in steady state [2].
In this paper, we generalize the BG model by considering the situation in which ξ x (t) and ξ y (t) are two fractional Gaussian noises (FGNs) with different parameters so that their correlation function reads The kernel N α j (t 1 , t 2 ) is defined such that in free space,ẋ free = ξ x (t) yields the so-called fractional Brownian for x free (see, e.g., [31]): In particular, the mean squared displacement exhibits an anomalous diffusion where α ∈ (0, 2) is the anomalous diffusion exponent. α = 1 corresponds to the standard Brownian motion while α < 1 and α > 1 generate sub-diffusion and super-diffusion respectively. The explicit expression of the kernel N α j that satisfies Eq. (5) for arbitrary α ∈ (0, 2) is not needed for the calculations presented in this paper. The interested reader can find detailed discussions on this topic in Refs. [28,29,31,32].
Let us note that, although in the following we will call abusively T x and T y "temperatures", one should keep in mind that they are not physical temperatures but rather only the noise amplitudes. Indeed, even if the components are independent (u = 0), each component is out of equilibrium since the noise and damping do not satisfy the fluctuation-dissipation theorem, except when α = 1.
Numerical simulations are performed by solving Eq. (2) using an Euler-Maruyama algorithm. The FGN is generated by using the Davies and Harte method [33] and is provided by the Stochastic python module.

Two variants
Although all our analytical results can be derived for the general FBG introduced above, it will help to discuss specifically two variants of the model in which we vary different parameters. In the first one (Model 1) the anomalous diffusion exponent is unique, α x = α y ≡ α, but the temperatures T x and T y are different, so that the noise correlations read In Model 2, the temperature T x = T y ≡ T is the same for both components but the anomalous diffusion exponents are different, so that

Spatial distribution
We first solve the equations of motion for a given realization of the noise which will allow us to obtain any moment of the position distribution. Thanks to their linear structure, the equations of motion Eq. (2) are readily integrated to give with the kernels Q c (t) = e −t cosh(ut) , Q s (t) = e −t sinh(ut) .
For simplicity, and since we are interested only in the steady-state properties, we chose the initial condition x(0) = y(0) = 0. Since the noises have zero mean it follows that x(t) = y(t) = 0 at all time.
In the following, we first compute the second moment of the position in Sec. 3.1 and the full probability density function in Sec. 3.2 and analyze quantitatively the effect of the FGN on the spatial distribution in Sec. 3.3.

Mean squared displacement and cross correlations
The second moments of the position distribution can be computed at all time from Eq. (9) and the expression of the noise correlation Eq. (5). They take the form where A j (t), B j (t), and C j (t) are given by the following integrals where j = x, y and to compute C j (t), we have used the symmetry N α j (τ 1 , τ 2 ) = N α j (τ 2 , τ 1 ).
Since we are interested in the nonequilibrium steady state (NESS), we evaluate the quantities in Eqs. (14)-(16) in the limit t → ∞. Details of the derivation can be found in Appendix A. Eventually, we find where A j , B j , C j denote the limiting values of A j (t), B j (t), C j (t) in the limit t → ∞.
One can notice that the asymptotic values of A j , B j and C j are finite (exist) only for |u| < 1. This is the only case of interest to us. When |u| ≥ 1, the physics is very different since the particle is not trapped, a fact which is reflected in the diverging mean-squared displacements along both directions as t → ∞.

Probability density function and characteristic function
Since the noise is Gaussian and the equations of motion are linear, the probability density function P (x, y; t) is fully determined by the second moments. This is most concisely expressed using the characteristic function where x 2 (t) , y 2 (t) , and x(t)y(t) are given by Eqs. (11)- (13). The calculation of the characteristic function Eq. (20) is performed in Appendix B.
The probability density function is then obtained as the inverse Fourier transform of Eq. (20), which gives 2 is satisfied as long as the moments remain finite.
Integrating Eq. (21) over one variable, the marginal distributions P (x; t) and P (y; t) take the expected form

Quantitative analysis of the effect of the FGN
We discuss here the effect of the FGN on the position moments in the special cases of Model 1 and 2.
In Model 1, which has α x = α y and T x = T y , the second moments Eq. (11)-(13) in the steady state simplify to and y 2 is obtained by interchanging T x and T y in the expression for x 2 . Putting α = 1, one recovers the results of Ref [3,14]. Fig. 1 shows the prediction of Eq. (23) for varying α. Comparing with numerical simulations confirms that the predictions are exact. We observe that, varying α, the shape of the curves is not affected. However, quantitatively, α has a strong effect. In particular, the dependence on u becomes much steeper when α increases, with the amplitude of all quantities increasing by more than an order of magnitude for all |u| > 0.5 when increasing α from 0.2 to 1.6. When u = 0, we find, as expected, that the two components decouple, i.e. xy = 0. The variance becomes x 2 = Γ(1 + α)T x so that we can interpret Γ(1 + α)T x ≡ T eff,x as an effective temperature and the marginal distribution P (x) ∝ e −x 2 /(2T eff,x ) as an effective Boltzmann distribution although the system remains out of equilibrium as discussed in Sec. 2.
For Model 2, which has T x = T y and α x = α y , Figure 2 shows the moments x 2 , y 2 and xy as functions of the anisotropy parameter u when increasing the difference in exponent while keeping α x = 2 − α y . We see that increasing the difference in exponent from 0 to 1.2, the amplitude of all quantities increases, more significantly as |u| gets larger, although the effect remains moderate.

Currents and rotation
A NESS is characterized by the existence of a probability current. We first compute this steady-state current in Sec. 4.1. We then use it to characterize the rotational motion of our gyrator through the angular velocity in Sec. 4.2 and the curl of the current in Sec. 4.3.

Steady-state current
Computing the current for our FBG involves subtleties due to the fact that for the noises with long-ranged temporal correlations, we cannot write a closed-form differential equation for the probability density function which would generalize the Fokker-Planck equation that can be written when α = 1 and give us an expression of the current. We circumvent here this fundamental difficulty by starting directly from the definition of the current J(x, y, t) = (J x (x, y, t), J y (x, y, t)) as an ensemble-average over noise realizations † It is easy to check that this is a good definition of the current, in the sense that it satisfies the continuity equation ∂ t P (x, y, t) + ∇ · J(x, y, t) = 0.
Using the integral representation of the Dirac delta function, J x can be expressed as which becomes, using the equation of motion Eq. (2), The average containing the first two terms on the r.h.s. can be computed by using the characteristic function with Φ(t, λ x , λ y ) (we have omitted the dependence on the other parameters for brevity) given in Appendix B.1. The last term requires different characteristic functions which are given in Appendix B.2.
Putting everything together (see details in Appendix C) we arrive at the following expression for the current in steady state We recognize the effective temperatures T eff,x/y = Γ(1 + α x/y )T x/y introduced in Sec. 3.3 which control the probability distribution at u = 0 when the two components are decoupled. In addition, we see from Eq. (30) that for any u the relation between current and probability is the same as the one that would be obtained for a standard Brownian gyrator driven by Gaussian white noises at the effective temperatures T eff,x/y .
The current field is shown in Fig. 3 for the special case of Model 1 with T x = 1, T y = 2 and u = 0.3 and two values of α. Interestingly, one observes that when α becomes larger, the symmetry axis of the current rotates.

Angular velocity
One can use the expression Eq. (30) for the current to compute the steady-state angular velocity Since the current vector belongs to the xy plane, we write ω = ω ẑ and, after straightforward algebra, we obtain the explicit, although lengthy, expression: where As for the standard BG, one needs a double symmetry breaking (geometric and in the noise) to observe a rotation. Indeed, one sees from Eq. (33) that the angular velocity vanishes if either u = 0 or if α x = α y and T x = T y .
On the contrary, as soon as both symmetry breakings are present, one observes a rotation in general. This is most clearly seen by performing a series expansion ω in powers of u to simplify Eq. (33) to At this order, the second symmetry breaking occurs when  Let us now discuss the predictions of Eq. (33) and compare with numerical measurements. In practice, we do not use Eq. (32) to compute the angular velocity in simulations but rather compute the time derivative of the polar angle. The mean of ω is displayed in Fig. 4 for Model 1 and Model 2. For Model 1, we observe that ω becomes larger as the asymmetry |u| increases and as α decreases. For Model 2, ω shows two extrema with |u| and increases in amplitude when the difference between the exponents increases.

Curl of the current
The rotational motion of the particle in the NESS can also be characterized by the curl of the current ∇ ∧ J which can be computed from Eq. (30). In particular, let us give here its expression at the origin in the case of Model 1 Again one can see the double symmetry breaking that is necessary to obtain a rotational motion since ∇ ∧ J(0, 0) in Eq. (39) vanishes when either T x = T y or u = 0. The known result for the Brownian gyrator [3], is easily retrieved by observing that G + (2u, u, 1) = u and G − (1 + u 2 , u, 1) = u 2 . Fig. 5 shows the curl at the origin for various temperatures and two exponents. One observes that the extrema depend little on the temperature but strongly on α, increasing by an order of magnitude when going from α = 1.75 to α = 0.25. The latter trend is also observed on the current field and its curl at arbitrary positions, displayed in Fig. 6. In the case of Model 2, unsurprisingly, the curl at the origin increases when the difference between the exponents increase, as illustrated in Fig. 7.

Conclusion
We have considered the two-dimensional motion of a particle in a harmonic potential. The particle undergoes Brownian motion in the overdamped limit under the action of different fractional Gaussian noises along the two spatial directions. As such it is out of equilibrium and display a non-vanishing current in steady state. Most strikingly, the particle rotates around the origin. This model can be seen as a minimal model of a nonequilibrium steady state driven by long-range correlated noise that can be tuned to be sub-diffusive or superdiffusive. We have computed analytically the full position probability distribution and the current in steady state, from which we deduced the angular velocity and the curl of the current characterizing the rotational motion of the particle. Numerical simulations confirmed these derivations in two special cases where the rotation is driven either by different amplitudes of the noise or by different exponents. The effect of the fractional noise, compared to the white noise case α = 1 is rather paradoxical: When decreasing α the position distribution becomes more and more insensitive to the asymmetry u but at the same time the current and angular velocity increase. For large values of α, the particle probability becomes much more spread and sensitive to the asymmetry but rotates slower around the origin.
Of course, the properties of the nonequilibrium steady-state that we have determined here depend on the details of the dynamics. For example, we expect that using a generalized Langevin equation as in Ref. [25] but including power-law correlations instead of exponential ones will lead to a different steady state compared to our result even though the free dynamics of each component is exactly the same in both models (see e.g. [28][29][30]). Long-ranged correlated non-Gaussian noise could also yield different results. However, there is hope that some universality is at play and that, qualitatively, our conclusions capture the effect of long-ranged time correlations beyond a specific model. It would be interesting to question this universality by studying other models in the future or even by comparing to experimental measurements. acknowledgments A.S. acknowledges the warm hospitality of LPTMC, Sorbonne Université during his research stay in October-December 2021.

Appendix A. Auxiliary functions
The calculation of A j (t), B j (t), and C j (t) involves integrals of the following form the choice of the kernel ω j depending on the specific function one is considering. We use the fact that is the autocorrelation function in free space. Integrating by parts twice we find an expression for g α (t) where ω (u) ≡ dω(u)/du.
For the function A j (t) the kernel is given by ω 1 (t) = ω 2 (t) = Q c (t). The first term in the right hand side of Eq. (A.4) yields the contribution 2t α because Q c (0) = 1. The terms involving a single integration are identical and their contribution is denoted A (I) j (t). Then, the contribution stemming from the term involving a double integral is denoted A Recalling that Q c (t) = [u sinh(ut) − cosh(ut)]e −t , the change of variables t − t 1 = zt gives It is convenient to express the integral in terms of the function Since we are interested in lim t→∞ A j (t), it is sufficient to consider the infinite-q limit of the function F (q, α) for positive q. In such a limit we can use the asymptotic expansion Taking 0 < u < 1, it follows that Let us consider now the contribution from the double integral We further split A (II) j (t) into a contribution due to the terms t α 1 + t α 2 in G α (t 1 , t 1 ) and the remainder due to |t 1 − t 2 | α . This splitting is denoted as follows and The integrals in Eq. (A.14) are factorized, therefore The integration yields When z → +∞ we have Focusing on those terms which in the large-t expansion are not exponentially suppressed, we find A We can now consider the integral Eq. (A.15), which can be written in the form which follows by applying the change of variables t 1 → t−t 1 and t 2 → t−t 2 in Eq. (A.15).
Since the correlation function G α (t 1 , t 2 ) is invariant under the interchange of t 1 and t 2 , we have Again, since we are interested in the large-t behavior, we can simply perform a Laplace transform with respect to t and apply the final value theorem. The Laplace transform is written as follows where s is the Laplace parameter. The integrations in Eq. (A.23) can be computed analytically and the corresponding result admits the following Taylor expansion for small s where a n (u, α) are functions that can be determined analytically. The final value theorem states that lim Let us mention that it is also possible to invert the Laplace transform in Eq. (A.22) and obtain the formal expansion for large t, including also the subleading corrections for large times. Collecting the results obtained so far, we find The calculation of this function, which has kernels ω 1 (t) = ω 2 (t) = Q s (t), turns out to be easier than the calculation of A j (t) because Q s (t) vanishes at t = 0. As a result, We proceed by following the very same guidelines which we followed for the calculation of A j (t). Using the same notation as before, we split Eq. (A.28) by writing The first integral gives When t → ∞ the function φ α reduces to a constant while the overall t-dependence The second integral in Eq. (A.30) is first written in the form We apply the Laplace transform of both sides and we get (A.36) The final value theorem then allows us to write Now using kernels ω 1 (t) = Q c (t) and ω 2 (t) = Q s (t), we write as before where C (I) The first integral gives and asymptotically it reduces to We split the second integral as already done in the previous sections writing Up to exponentially small terms, We can now consider C (II,2) j (t), which is given by In contrast with the functions A (II,2) j (t) and C (II,2) j (t), the integrand is not symmetric under the exchange of t 1 and t 2 . This time, we approach the calculation by performing directly the Laplace transform and then integrating with respect to t 1 and t 2 . The application of the final value theorem then gives Assembling all the pieces together, we finally obtain and therefore the characteristic function Eq. (B.2) is Ψ(t, λ x , λ y ; T j ; α j ) = exp 2 −1 G(t, λ x , λ y ; T j ; α j ) .

Appendix B.2. The function Ξ j
In order to compute the characteristic function Ξ j it is instructive to consider a preliminary exercise which can be applied to the calculation of Ξ j . Suppose that we want to compute the following characteristic function where ξ x (t) is the FGN, denotes the convolution defined by and X is some kernel. The calculation of Eq. (B.7) is performed by expanding in power series the exponential and using the fact that only even powers of ξ x contributes in the averaging. Therefore The quantity in angular brackets can be computed by applying Wick's theorem, which yields , Thanks to the identity (2n + 1)!! (2n + 1)! = 1 2 n n! , (B. 12) it follows that Eq. (B.11) is the Taylor expansion of the exponential. The summation of the series gives We can now to carry out the calculation of the characteristic function Ξ j . Let us start with Ξ x , Ξ x = ξ x (t)e iλxx(t)+iλyy(t) . (B.14) The quantity in the exponential can be factorized in two terms: the first contains ξ x and a second ξ y The second factor is precisely Ψ(t, λ y , λ x ; T y ; α y ). The first factor has the form of Eq .(B.7) with the kernel X given by The result, Eq. (B.13), can be applied with where U j (t) and V j (t) are the auxiliary functions To recap, Ξ x is given by and from Eq. (B.1) it follows that An analogous reasoning applies to Ξ y giving Ξ y = iT y λ x V y (t) + λ y U y (t) Φ(t, λ x , λ y ; T x , T y ; α x , α y ) .

(B.22)
In order to compute the functions U j (t) and V j(t) in the simplest possible way, we introduce the function thanks to which it is possible to use the linear decomposition The function W j (t, r) has an interesting significance: it is the equal-time correlation function between the noise of a one-dimensional fractional Brownian motion in a fictitious optical trap. More precisely, this identification is established by considering the onedimensional fractional Brownian motion in the parabolic trap such that the system is described by the Langevin equatioṅ where ζ(t) is a fractional Gaussian noise with zero average and autocorrelation function given by ζ(t 1 )ζ(t 2 ) = N α j (t 1 , t 2 ), the parameter r is the stiffness of the potential. In free space, N α j (t 1 , t 2 ) originates the autocorrelation function of the position X f (t 1 )X f (t 2 ) = (t α 1 + t α 2 − |t 1 − t 2 | α ) . It is now clear that the right hand side of Eq. (B.29) is the function W j (t, r). As a result, we can write W j (t, r) = X(t)ζ(t) , (B.30) as we anticipated. We then plug in Eq. (B.30) the expression for ζ(t) provided by Eq. (B.26) to get W j (t, r) = r X 2 (t) + X(t)Ẋ(t) .