Stochastic resonance in a suspension of magnetic dipoles under shear flow

We show that a magnetic dipole in a shear flow under the action of an oscillating magnetic field displays stochastic resonance in the linear response regime. To this end, we compute the classical quantifiers of stochastic resonance, i.e. the signal to noise ratio, the escape time distribution, and the mean first passage time. We also discuss limitations and role of the linear response theory in its applications to the theory of stochastic resonance.


I. INTRODUCTION
The dynamics of periodically-driven stochastic systems has been an active field of research over the last years [1]. This kind of systems arise commonly in fields of physics, chemistry and biology. Examples are found in problems involving transport at the cellular level [2,3], optical and electronic devices [4] and signal transduction in neuronal tissue [5,6], to cite just a few.
A particularly interesting phenomenon, occurring in periodically-driven non-linear noisy systems, is stochastic resonance (SR) [7]. This term refers to the enhancement of the response of the system to a coherent signal when the intensity of the noise grows, instead of the degradation that one naively expects. The mechanism leading to this phenomenon is quite simple. Imagine a system which exhibits an energetic activation barrier. In the presence of noise, the system can be assumed to surmount this barrier with a rate proportional to e −∆E/D , where ∆E is the height of the barrier and D is the intensity of the noise acting on the system. The inverse of this rate defines the average waiting time, T (D), between two noise-induced transitions. In the presence of a periodic forcing, the height of the barrier is periodically raised an lowered. When the period of the external force synchronizes with 2T (D), the barrier surmounting will be enhanced by the cooperative effect of the noise and the periodic forcing.
Although originally proposed for systems in a double-well potential, this original scheme has been overcome. In fact, it is known that SR is exhibited by several classes of monostable systems, among which one could mention excitable and threshold systems [8][9][10][11] or systems which do not follow an activated dynamics but a relaxational dynamics [12,13].
In this paper we will show that a magnetic dipole immerse in a shear flow exhibits stochastic resonance when a weak oscillating magnetic field is acting on it. The presence of this flow takes the system out of equilibrium impressing certain peculiarities in the behaviour of this system. In order to treat this problem we will analyze the response of the system in the linear regime. Previous study of the dynamics of a dipole under an oscillating magnetic field has revealed that linear response theory (LRT) predicts a monotonically decreasing behaviour for the ratio between the output signal and the output noise or signal to noise ratio (SNR), i.e. for very weak applied fields the dipole does not exhibit SR [13]. In the present case there is a new ingredient, absent in [13]: the presence of the shear flow, which is determinant for many interesting aspects of the dynamics of this system. Additionally, although we show the existence of SR in the linear regime, we discuss limitations and role of the LRT in its application to the theory of SR, mainly related to questions concerning the fluctuation-dissipation theorem.
The paper is organized as follows: in section II we analyze the dynamics of a dipole in a shear flow and find the fixed points. Section III is devoted to study the response of the system to an oscillating magnetic field by computing the susceptibility. In section IV we calculate the power spectrum and the signal-to-noise ratio. In section V we compute the escape time distribution and from it the mean first passage time. Finally, in section VI we discuss our results. withˆ z being the unit vector along the z-axis, and of an oscillating field H = He −iωtˆ x, withˆ x being the unit vector along the x-axis. The dynamics of these dipoles is governed by the following equation of motion where I is the moment of inertia of the particles, ξ r = 8πη 0 a 3 is the rotational friction coefficient, with η 0 the solvent viscosity, and a the radius of the particle. For t ≫ τ r ,with τ r = I/ξ r being the inertial time scale, the motion of the particle enters in the overdamped regime. This time scale defines a cut-off frequency ω r = τ −1 r , such that the condition for overdamped motion is equivalent to ω ≪ ω r . In this case eq. (1) yields the balance condition between the magnetic and hydrodynamic torques acting on each particle which, together with the rigid rotor evolution equation Here λ(t) ≡ (m s H/ξ r ω 0 )e −iωt , with Ω p being the angular velocity of the particle. The computation of the fixed points of eq. (4) when the magnetic field is held constant, i.e. λ(t) = λ 0 = cnt., and their linear stability analysis are given in detail in Appendix A. After some algebra eq. (4) becomes For λ 0 ≥ 1, this equation has only a linearly stable stationary state. The orientation of the suspended particles is fixed toˆ This means that in this regime the hydrodynamic torque, which tends to cause the rotation of the particles is insufficient to overcome the magnetic torque which maintains constant their orientations. For λ 0 < 1, which is the case we are interested in, the particles undergo a rotation around a fixed axis lying in the y-z plane, being the director of this axis given bŷ In this case the hydrodynamic torque is strong enough to make the dipole to precess around the orientationˆ R s eq. (7) (see Appendix A).

III. RESPONSE TO AN OSCILLATING MAGNETIC FIELD
The analysis of Section II has been carried out for the deterministic dynamics of a magnetic dipole in a shear flow. Fluctuations are introduced by means of a Brownian torque. The corresponding Langevin equation is where F B (t) is a Gaussian white noise of zero mean and correlation function The Fokker-Planck equation corresponding to eq. (8) is given by where L 0 and L 1 are operators defined by with D r = k B T /ξ r being the rotational diffusion coefficient, and R =ˆ R × ∂/∂ˆ R the rotational operator. Notice that the first and second terms on the right hand side of eq. (11) 1 , correspond to convective and diffusive term, respectively. Moreover eq. (10) which, according to eq. (8), rules the Brownian dynamics in the case of overdamped motion, is valid in the diffusion regime. This regime is also characterized by the condition t ≫ τ r , or equivalently ω ≪ ω r , which implicitly involves the white noise assumption.
To solve the Fokker-Planck equation (10) we will assume that λ 0 ≡ |λ(t)| constitutes a small parameter such that this equation can be solved perturbatively. Thus up to first order in λ, the solution of the Fokker-Planck equation (10) is Here Ψ 0 (t ′ ) = e (t ′ −t0)L0 Ψ 0 (t = t 0 ) is the zero order solution at time t ′ , and withˆ R 0 being an arbitrary initial orientation. As follows from eq. (11) 1 , the unperturbed operator L 0 is composed of the operators R z and R 2 , which are proportional to the orbital angular momentum operators of quantum mechanics L z and L 2 , respectively, and, therefore, their eigenfunctions are the spherical harmonics [14] R z Y l m (ˆ R) = imY l m (ˆ R), Given that we know how R acts on the spherical harmonics, it is convenient to expand the initial condition in series of these functions, since the spherical harmonics constitute a complete set of functions which are a basis in the Hilbert space of the integrable functions over the unit sphere [15] Using this expansion in eq. (12), for the first order correction to the probability density, ∆Ψ ≡ Ψ − Ψ 0 , we obtain Notice that the integral of ∆Ψ(ˆ R, t) over the entire solid angle is zero, in agreement with the fact that the unperturbed solution Ψ 0 (ˆ R, t) is normalized. Since, we are interested in the asymptotic behavior we will set t 0 → −∞. In this limit, eq. (12) becomes where now and Ψ 0 (ˆ R, t) = 1 4π (19) is the uniform distribution function in the unit sphere.
From eq. (18) the contribution of the AC field to the mean value of the orientation vectorˆ R can be obtained aŝ This equation can be written in the more compact form where the response function [16] has been defined as for τ > 0. By causality, we can write t → ∞ in the upper limit of the integral in eq. (21); hence, this equation becomeŝ where χ i (ω) is the generalized susceptibility, which is the Fourier transform of χ i (τ ) From this equation we obtain the components of the susceptibility The quantities χ x , and χ y , have poles at ω = ±ω 0 ± 2D r i. The inverse of the imaginary part of these poles, (2D r ) −1 , defines the Brownian relaxation time.

IV. POWER SPECTRUM
In order to discern whether or not SR is present in the relaxation process of the system under consideration we compute the power spectrum, which, following the Wiener-Khinchine theorem, is given by the Fourier transform of the correlation function [17,1]. Since we will take as output signal the projection ofˆ R parallel to the magnetic field, i.e.R x , we only compute the correlation function of this quantity.
The correlation function ofR x is defined by where the initial condition is taken as Ψ(ˆ R, t 0 ) = δ(ˆ R −ˆ R 0 ). The above quantity can be calculated from the solution Fokker-Planck equation simply by recalling the following properties of a Markov process [17].
where t 1 > t 2 > ... > t n . By combination of these two properties eq. (28) becomes To proceed further, we compute in first place the integral overˆ u in eq. (30). From eq.
thus we have The result of these integrals are (for the detailed derivation, see Appendix B) After introducing these expressions into eq. (30) we will obtain three terms corresponding to an expansion of the correlation function in powers of λ(t), of zeroth, first and second order respectively. The presence of this driving yields an explicit dependence of the correlation function on the time t, instead of depending only on the time difference, as occurs in the stationary case. The method for removing this dependence on the initial time is to average the correlation function over a period of the driving [1]. After doing this the first order term vanishes, and consequently we do not worry about it and compute only those whose average gives a non-zero contribution, i.e. the zeroth and second order terms. Taking this into account, and by applying eq. (12) where the sign ∼ indicates that the terms which vanish after averaging over the period of the driving has been neglected (although the average has not been performed yet). After introducing the corresponding expressions for Ψ(ˆ v, t|ˆ R 0 , t 0 ) and by using eq. (33) we obtain (the details of this computation are given in Appendix B) where we have defined At this stage, and before applying the Fourier transform to the correlation function to obtain the power spectrum of the processR x (t), we average eq. (35) to obtain This computation has been carried out in the assumption that τ is a positive quantity. To extend our computation to τ < 0 we have to use the backwards Fokker-Planck equation. The operator which generates the backwards evolution of the probability distribution −L † [18] being L the Fokker-Planck operator and L † its adjoint operator. Consequently the formal solution of the backwards Fokker-Planck, equivalent to eq. (12) is given by (t < t 0 ) As in the case of the operator L 0 , the spherical harmonics are eigenfunctions of −L † 0 with eigenvalues given by Thus the process to follow in the calculation of the correlation function for τ < 0 is identical to the corresponding computation for τ > 0 changing the eigenvalues of the operator L 0 by the ones of −L † 0 , which yields for τ < 0. We now apply to this averaged correlation function (now defined for −∞ < τ < ∞) the Wiener-Khinchine theorem, which states that the power spectrum and the correlation function are related through a Fourier transform. Thus, Since our purpose is to discern whether or not SR is present in the relaxation process of the quantityR x (t), we proceed to compute the signal to noise ratio, R, i.e. the ratio between the weight of the delta function in (41) 1 and the noisy part of Q(Ω) computed at the frequency of the driving. From eqs. (25) and (41) we achieve This quantity has been plotted in Fig. 1 as a function of the inverse of the Péclet number, P e −1 = D r /ω 0 , which measures the ratio between the time scales associated to diffusion (thermal noise) and flow. The presence of a maximum in R for a non-zero value of this parameter shows the existence of stochastic resonance in the relaxation process of a dipole in a shear flow. In addition to the slow relaxation to the single attractor of the dynamics, our model includes another effect, which hides, to some extent, the SR profile. To understand this, note that even though the signal is too weak, it nevertheless causes the position of the attractor of the dynamics to vary, and so the output will always have a nonzero component at the signal frequency. This fact causes the SNR to go to infinity in the zero noise limit [19].

V. MEAN FIRST PASSAGE TIME
In this section we study the behaviour of the escape time distribution and the mean first passage time of the magnetic dipole immersed in a shear flow. To this end, we have to account for the fixed point orientations of eq. (4) in the case λ 0 < 1. In this situation there is a single fixed point corresponding to an orientation contained in the plane x = 0 or, equivalently, φ = π/2. However, when λ(t) > 0 this stationary orientation is in the subspace z > 0 (cos θ > 0) and in the subspace z < 0 (cos θ < 0) if λ(t) < 0. Therefore, we are going to study the escape from the region z > 0 (cos θ > 0) assuming that the initial orientation of the dipole is contained in this region. Consequently, we have to solve the Fokker-Planck equation (10) with absorbing boundary conditions in the plane cos θ = 0 [20,21], i.e.
Since this escape problem will be treated perturbatively, the first step is to analyze the eigenvalue problem of the operator L 0 eq. (11) under the boundary condition (43). It is easy to check that the eigenfunctions and the eigenvalues are the same with the restriction that only those spherical harmonics which vanish at cos θ = 0 are solution of this eigenvalue problem. From the parity properties of the associated Legendre functions [15] one can see that (43) selects only the spherical harmonics such that l + m = 2n + 1 with n = 0, 1, .... Thus, we have where l, m denotes that the sum is carried out over 0 ≤ l < ∞ and −l ≤ m ≤ l restricted by l + m = 2n + 1.
In order to evaluate the mean first passage time (MFPT), we have to compute previously the survival probability, S(ˆ R 0 , t), and the escape time distribution (ETD) which are related through where S(ˆ R 0 , t) is defined by with R, the region from which we are studying the escape problem (in the present case cos θ > 0), andˆ R 0 ∈ R, the initial orientation of the dipole. The probability distribution Ψ(ˆ R, t|ˆ R 0 ) is obtained from eq. (12) with the boundary conditions (43) and the initial condition Before proceeding to obtain the survival probability, there are some facts to consider which will facilitate further computation. Looking at eq. (46) one can realize that, due to the integration over the azimuthal angle, only terms with m = 0 contribute to S(ˆ R 0 , t). Consequently, the selection rule l + m = 2n + 1 reduces to keep only the odd values of l. In addition, we are interested only in the modes with greater relaxation times. Therefore, from the whole series eq. (44) we are only interested in the term l = 1, m = 0. Thus, our purpose is to obtain the coefficient a 10 (t) up to first order in λ(t) from eqs. (10) and (43). Up to zeroth order, we have a (0) The obtention of the first order contribution, a 10 (t), requires a little bit more elaborated calculation. To proceed further with this computation, the operator L 1 acting on Y lm (ˆ R) yields where the action of R z on Y lm (ˆ R) is given by eq. (14) 1 and From eqs. (49) and (50) together with the rules for the addition of angular momenta familiar from quantum mechanics [14] and the selection rule l + m = 2n + 1 imposed by the boundary condition (43), one can deduce that only the terms Y 2±1 (ˆ R) contributes to a 10 (t). The mentioned rules of addition of angular momenta implies that the product of two spherical harmonics, Y lm Y pq has projection onto a third spherical harmonic, Y rs , only when m + q = s. On the other hand these same rules imposes the restriction that the product Y lm Y pq only projects onto the subspaces such that |l − p| ≤ r ≤ l + p. By using these restrictions one can see that when in eq. (49) one takes l = 1 and m = 0 one obtains a vanishing contribution and only when l = 2 and m = ±1 the contribution to a Taking these considerations into account and by using the following results .
we obtain the first order correction to the coefficient a 10 (t).
a (1) Eqs. (48) and (52) together with (44) and (46) allows us to obtain the survival probability S(ˆ R 0 , t), which is given by This quantity is directly related to the MFPT, since where we have used eq. (45). Consequently, the MFPT is given by where we have takenˆ R 0 =ˆ R s . In Fig. 2 we have plotted the quantity ∆T /T 0 . This figure shows that this quantity exhibits a minimum, as required for the appearance of SR.
The knowledge of the survival probability allows us to obtain the ETD, ρ(ˆ R 0 , t). From eqs. (45) and (53) the ETD is given by By takingˆ R 0 =ˆ R s , we finally obtain where ∆ρ ≡ ρ − ρ 0 , being ρ 0 the corresponding ETD when the amplitude of the oscillating field is set to zero. The succession of maxima in the ETD (see Fig. 3) is indicating that the dynamic of the magnetic dipole suspended in a shear flow under a periodic field exhibits SR.

VI. DISCUSSION
We have shown that the relaxation process of a dipole immersed in a shear flow exhibits SR upon application of a weak periodic field. To this end we have computed three quantities typically used to characterize SR, namely the signal-to-noise ratio, the escape-time distribution and the mean first passage time. All of them behave as expected for a process in which SR occurs.
Previous works devoted to analyze whether or not SR is present in the relaxation process of an overdamped dipole in a fluid at rest have revealed that this phenomenon does not occur in the linear regime [13]. Effectively, linear response theory (LRT) predicts a maximum in the signal, i.e. in the susceptibility, as a function of the noise level. Unlikely, the SNR decreases monotonically with the noise level. This behaviour can be easily understood. In the limit of zero noise the output of the system has a small component (proportional to the applied field) at the frequency of the signal whereas the background noise vanishes at zero noise level, being this behaviour the responsible for the monotonic dependence of the SNR on the noise intensity.
In our case, the situation is completely different. When the fluid in which the dipole is suspended is submitted to a pure rotation (vortex flow), both output signal and output background noise exhibit a peak at the same value of D r (see Fig. 4). Consequently, although the background noise vanishes when D r goes to zero, the characteristic SR profile of the SNR can not be completely hidden, as shown in Fig. 1. This feature arises as a consequence of the presence of shear acting on the suspension, thus, the appearance of SR in the system studied in this paper is a non-equilibrium feature.
In a sense, the mechanism yielding SR in this system is similar to the one operating in SR in threshold devices [9,10]. Due to the presence of noise, the dipole can eventually acquire enough energy to get out from its stable orientation by crossing the absorbing barrier (the threshold) cos θ = 0. After this, the system is driven to its stable position. This process produces a short spike in the magnetization. Of course, the time that the system takes to return to the fixed point has to be smaller than the semiperiod of the oscillating magnetic field. Thus, SR in this system can be understood in the same way as, for example, the SR in level crossing detectors [9].
LRT has been one of the widest used tools in the study of stochastic resonance [22]. When the system, in the absence of the external periodic force, is in thermal equilibrium a very adequate form of describing stochastic resonance is in terms of the susceptibility. This is because the noisy part of the power spectrum is given directly by the susceptibility through the fluctuation-dissipation theorem, This result is correct when the fluctuations whose spectral density is given by N (Ω) have the thermal equilibrium state as reference state [16]. However, in the present case we are dealing with a system which is maintained in an out-of-equilibrium steady state due to the presence of a shear flow. It is evident from Fig. 5, where we have plotted the imaginary part of the susceptibility corresponding toR x and the noisy part of the power spectrum, that these two quantities are clearly different. Note, that if ω 0 = 0, i.e. the system in the absence of the periodic field is in equilibrium, the relation (58) is fulfilled. Thus, we have shown that although we can define a susceptibility which describes the response of our system to a small perturbation we can not describe SR by means of LRT. The reason can be found in the fact that due to the non-equilibrium nature of the attractor of the dynamics the fluctuation-dissipation theorem fails to be valid.

ACKNOWLEDGMENTS
The authors thank Miguel Rubí for valuable discussions. This work has been supported by DGICYT of the Spanish Government under grant PB98-1258. One of us (T. Alarcón) wishes to thank to DGICYT of the Spanish Government for financial support. From eqs. (3) and (4) one can see that the time derivative of theˆ R vanishes either when Ω p = 0 or when Ω p ×ˆ R = 0. In the first case we have From this equation, and taking into account that |ˆ R s | = 1, we obtain that the stationary orientation iŝ This solution exists only when λ 0 ≥ 1 and corresponds to a fixed orientation of the dipoles, given that the intensity of the magnetic field is high enough to maintain this fixed direction. The second possibility leads tô Eq. (A.3) provides two equations for three unknowns. If one setsR z = 0 one recovers eq. (A.2). If, by contrast one makesR x = 0 then a new stationary orientation is obtained, which exists only when λ 0 ≤ 1. This orientation gives rise to a rotation of the dipoles with angular velocity since, in this case, the field is not strong enough to inhibit the rotation caused by the shear flow. The linear stability of this fixed points is better analyzed in spherical coordinates. Taking into account that we obtain After expressing the components ofˆ R in spherical coordinates, we obtain the following bidimensional dynamical system where θ and φ are the polar and azimuthal angles, respectively. By linearization of eqs. (A.8) around the λ ≥ 1 fixed points we obtain the following matrix which implies that, if λ is positive, the orientation corresponding to choose the sign "+" in eq. (A.2) is stable, while the other one is unstable. The same linearization procedure carried out around the λ < 1 fixed points leads to The eigenvalues of this matrix are given by In this appendix we work in detail some steps of the computation of the power spectrum corresponding to the relaxation process of a dipole under an oscillating magnetic field in a shear flow, in particular we calculate the integrals which yield eqs. dˆ uû x Ψ(ˆ u, t + τ |ˆ v, t) = I 1 + I 2 , To begin with we focus on the integral I 1 , which can be rewritten as where L † 0 is the adjoint operator of L 0 defined by with A and B two arbitrary observables. Explicitly, L † 0 is given by and therefore eq. (B.2) reads leading to eq. (33) 1 . In eq. (B.5) we have used the relation To compute the integral I 2 , we have to use the following representation of the delta function After introducing this expression into eq. (B.1) 3 we obtain and by using eq. (B.6), eq. (B.8) yields eq. (33), i.e.
Once these expressions have been obtained we can compute the correlation function given by eq. (34), being Ψ 0 (ˆ v, t|ˆ R 0 , t 0 ) and ∆Ψ(ˆ v, t|ˆ R 0 , t 0 ) given by where the initial time t 0 has been fixed to zero and eq. (B.7) has been used.
Eqs. (B.11) provide the evolution of the probability distribution under the condition of the system being initially in the stateˆ R 0 . Since a priori nothing is known about this initial condition, we assume thatˆ R 0 is a random variable uniformly distributed over the orientation space, consequently we average the correlation function over the distribution of initial states (eq. (36)). Taking into account that the average of the correlation function over initial conditions dˆ vv x I 1 , After introducing eqs. (B.5) and (B.6) into (B.13) 2 , we obtain where the orthogonality relation for the spherical harmonics, Let us focus our attention on the integral overˆ v. (B.17) From the rules of addition of angular momentum familiar from quantum mechanics, which imply that the product of two spherical harmonics, Y pq Y rs , only has a non-vanishing projection over a third spherical harmonic, Y lm , when the relations |p − r| ≤ l ≤ r + p and q + s = m are fulfilled, it is easy to see that these integrals will give a non-zero result only when l = 0, 1, 2 [14]. In addition, for the integrals containing the products Y 1±1 (ˆ v)Y 1±1 (ˆ v) the parameter m has to be m = ±2 whereas it must be m = 0 for the integrals with Y 1±1 (ˆ v)Y 1∓1 (ˆ v) to yield a non-zero contribution. However, although these integrals gives in principle a non-vanishing contribution, note that when we perform the integral over the variableˆ u in eq. (B.16) the terms introduced by these contributions finally yield, by the orthogonality property of the spherical harmonics, a vanishing result. Thus, only when l = 0 and m = 0 contributes to I 3 . Taking this into account and using eq. (B.15) finally performing the changes of variables t ′ = t − r and t ′′ = t + τ − s and using eq. (25) we obtain In this integral the upper limit goes to infinity by causality reasons.