text:
# Cosmic ray acceleration at supergalactic accretion shocks: a new upper energy limit due to a finite shock extension

## 1 Introduction

Formation of large-scale extragalactic structures (galaxy clusters and super-clusters) involving accretion of plasma to the gravitational potential wells is accompanied by formation of strong shock waves (Ryu et al. 1993). These shocks were considered as viable sources of ultra high energy cosmic rays (UHECRs), up to energies $`E_{max}10^{1819}`$ eV (Kang et al. 1997, Hillas 1984). Evaluation of this maximum energy is not reliable due to inadequate knowledge of physical conditions at the shocks considered. The first estimates show that maximal particle energies accelerated in such conditions must be limited because of the extensive time scales involved. The magnetic fields in the considered flows are not expected to be stronger than $`1`$ $`\mu `$G for the accretion in rich galaxy clusters and $`0.1`$ $`\mu `$G for supergalaxy scale accretion (Kronberg 2001). At the same time the former flows can involve plasma velocities $`1000`$ km/s, while the latter a few times lower ones (Juszkiewicz et al. 2000). The ordinary diffusive Fermi acceleration considered (Drury 1983) is fastest for highly turbulent conditions with the diffusion coefficient comparable to the Bohm diffusion coefficient, $`\kappa _B`$. In the considered situations, the involved acceleration time scale $`\tau \kappa _B/V_{shock}^2`$ becomes longer than the age of the universe for particle energies greater than $`10^{18}`$ eV.

However, as pointed out by Kang et al. (1997; cf. Jokipii 1987 and Ostrowski 1988), even faster acceleration can take place at perpendicular shocks occurring when large-scale magnetic field structures are present in the intergalactic medium and the accretion proceeds perpendicular to the field direction. Then, the cosmic ray cross-field diffusion coefficient governs the first order Fermi acceleration process at the shock. With the involved spatial scales along the shock normal of respectively 1 Mpc for a cluster, or a few Mpc for a super-cluster flow, there is enough space and time to accelerate particles up to $`10^{21}`$ eV, three orders of magnitude greater than in the Bohm diffusion limit (cf. Siemieniec-Oziȩbło & Ostrowski 2000).

Below, in section 2, we present evaluations of the acceleration time scales in perpendicular shocks. The derived time scale is used in section 3 to discuss the particle diffusion and drifts along the shock front, perpendicular to the plasma flow direction. We show that particle escape to the sides of the shock structure may substantially limit $`E_{max}`$ derived for perpendicular shocks. In section 4 we briefly discuss the importance of this constraint for large-scale extragalactic shocks.

## 2 Particle acceleration at perpendicular shock waves

The acceleration time scale derived within the diffusion approximation for oblique shocks (Ostrowski 1988) becomes formally invalid for perpendicular shocks which are always superluminal. Nevertheless, Jokipii (1987) used the diffusion equation to derive the energetic particle acceleration time scale for such shocks:

$$\tau =\frac{3}{U_1U_2}\left(\frac{\kappa _1}{U_1}+\frac{\kappa _2}{U_2}\right),$$
$`(2.1)`$
where $`\kappa _1`$ and $`U_1`$ ($`\kappa _2`$ and $`U_2`$) are the normal diffusion coefficient and shock velocity as measured upstream (downstream) of the shock. With the scaling assumed by Jokipii (1987)

$$\frac{\kappa _{}}{\kappa _{}}=\frac{1}{1+\left(\frac{\lambda _{}}{r_g}\right)^2},$$
$`(2.2)`$
where symbols $``$ and $``$ indicate quantities perpendicular and parallel to the mean magnetic field, respectively, and $`\lambda _{}`$ is a mean free path for particle scattering along the magnetic field ($`\kappa _{}=\frac{1}{3}\lambda _{}w`$, where $`w`$ is a particle velocity), and $`r_g=pc/eB`$ is the particle gyroradius. The equation (2.2) describes the cross field diffusion, where, at the mean free time, the particle moves on average $`\lambda _{}`$ along the magnetic field and is scattered at a distance $`r_g`$ across it. Continuing the acceleration process at perpendicular shocks requires sufficiently effective particle cross-field diffusion. With $`\eta \lambda _{}/r_g`$ this requirement can be expressed (Jokipii 1987) as

$$\eta <\frac{w}{U_1}.$$
$`(2.3)`$
## 3 Maximum energy of accelerated particles

### 3.1 $`\tau `$: a diffusive approach

In order to evaluate the maximum energy attainable in the acceleration process acting at the perpendicular shock wave of a finite extension $`\mathrm{\Delta }X_{max}`$ along the magnetic field one has to consider particle escape from the shock region, represented by free escape boundaries of figure 1. Let us assume Eq. (2.1) to be a valid acceleration time estimate. Additionally, we consider conditions with the constant $`\eta `$ upstream and downstream of the shock i.e. $`\eta =\lambda _{,1}/r_{g,1}=\lambda _{,2}/r_{g,2}`$ . With the shock compression ratio $`rU_1/U_2`$ and a weakly perturbed upstream magnetic field $`B_1`$, the downstream magnetic field $`B_2=rB_1`$ and the particle gyroradius becomes $`r_{g,2}=r_{g,1}/r`$ . Fast acceleration at perpendicular shocks appears only for weakly turbulent conditions with $`\lambda _{}r_g`$. From Eq. (2.2), in this limit $`\kappa _{}=\kappa _{}\left(\frac{r_g}{\lambda _{}}\right)^2`$. For simplicity let us consider a purely perpendicular shock with the mean magnetic field parallel to the shock surface and, respectively, $`\kappa _1=\kappa _{,1}`$ and $`\kappa _2=\kappa _{,2}`$, and relativistic cosmic ray particles with $`w=c`$. Then, with the Bohm diffusion coefficient $`\kappa _B\frac{1}{3}r_gc`$, the diffusion coefficients

$$\kappa _i=\frac{1}{\eta }\kappa _{B,i}(i=\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}),$$
$`(3.1)`$
and the equation (2.1) transforms to

$$\tau =\frac{2r}{r1}\frac{1}{\eta }\frac{r_{g,1}}{c}\left(\frac{c}{U_1}\right)^2.$$
$`(3.2)`$
With the constraint (2.3) and the strong shock with $`r=4`$ the acceleration time scale must be longer than $`c/U_1`$ times a particle gyroperiod.

### 3.2 $`\tau `$ : a lower limit for the perpendicular shock

An independent lower limit for the acceleration time scale can be derived by considering physical conditions approaching a weak scattering limit. Then, in the absence of scattering, let us consider the shock-drift acceleration process of upstream particles at the perpendicular shock propagating in a uniform magnetic field. A particle with the perpendicular momentum component $`p_{}`$ compressed in the shock will increase it to

$$p_{,2}=p_{,1}\sqrt{r}$$
$`(3.3)`$
due to the approximate magnetic momentum conservation. For strong shocks, on average particles roughly double their momenta at individual transmissions through the shock, which take approximately

$$\mathrm{\Delta }t=\frac{2r_{g,1}}{U_1}=\frac{2r_{g,2}}{U_2}.$$
$`(3.4)`$
Thus the characteristic acceleration time cannot be shorter than $`\tau \mathrm{\Delta }t`$. The same result can be obtained by considering particle drift along the shock, perpendicular to the weakly perturbed magnetic field. Addition of the diffusive process that is required to allow downstream particles to return to the shock to continue acceleration can only make $`\tau `$ longer. One may note that the minimum acceleration time in the Jokipii (1987) derivation approximately coincides with the above one. Thus, the time (3.4) can be considered to be an absolute lower limit for the acceleration time in the perpendicular shock and thus it defines the ultimate limit for $`\kappa _{}`$ reduction (or the upper limit for the considered $`\eta `$) in derivation of the expression (3.2).

### 3.3 New constraints for $`E_{max}`$

Cosmic ray removal from the acceleration region near the shock proceeds due to particle advection with the plasma flow downstream of the shock and due to diffusion or drifts towards the escape boundaries around the shock region. In the present paper we discuss the role of such escape in limiting the maximum energy of accelerated particles.

#### 3.3.1 A limit due to escape along the magnetic field

First let us discuss the diffusive escape along the magnetic field (cf. Fig. 1) limiting the maximum particle energy at the shock of a finite extension. To evaluate this limit we consider a solution of the 1-D spatial diffusion equation for particles moving along the magnetic field (along the $`x`$-axis of our coordinate frame), with the diffusion coefficient $`\kappa _{}`$. For this evaluation we neglect particle transport perpendicular to this field, both along the shock normal and to the sides. For $`N_o`$ particles injected at $`x=0`$ and $`t=0`$, their distribution along the $`x`$-axis is given by

$$n(x,t)=\frac{N_o}{\sqrt{4\pi \kappa _{}t}}\mathrm{exp}\left[\frac{x^2}{4\kappa _{}t}\right].$$
$`(3.5)`$
If a dispersion of this distribution $`\sqrt{2\kappa _{}t}`$ becomes larger than the shock sideway extension $`\mathrm{\Delta }X_{max}`$, the accelerated particles tend to efficiently escape to the sides, leading to the spectrum steepening and, finally, to formation of the energy cut-off. Such modification of the particle spectrum is possible when the time required for escape becomes comparable to the acceleration time scale $`\tau `$, i.e. when the particle dispersion along the $`x`$-axis over time $`\tau `$ becomes comparable to the shock extension $`\mathrm{\Delta }X_{max}`$:

$$\mathrm{\Delta }X_{max}2\sqrt{2\kappa _{}\tau }.$$
$`(3.6)`$
Applying the above Eq. (3.2) and the expression for $`\kappa _{}`$, we obtain the requirement for diffusively accelerated particles at such finite-extent shocks as

$$r_g\frac{\mathrm{\Delta }X_{max}}{2}\frac{U}{c}.$$
$`(3.7)`$
As the particle gyroradius $`r_g=pc/eB`$, this expression provides the upper energy limit for accelerated particles due to particle escape along the mean magnetic field. For nonrelativistic shocks with the factor $`U/c1`$, Eq. (3.7) provides a severe constraint for the maximum accelerated particle energy.

#### 3.3.2 A limit due to drift across the magnetic field

As discussed in section 3.2, particle acceleration at the perpendicular shock with the weakly perturbed magnetic field involves particle drift along the shock surface, across the magnetic field. Thus, the shock extension along this direction (along the $`y`$-axis of our coordinate frame) also provides a limit for the maximum energy of accelerated particles. The energy gain of a particle drifting along the shock is

$$\mathrm{\Delta }E=|qL|,$$
$`(3.8)`$
where the ”V$`\times `$B” electric field $`=(𝒰_\mathcal{1}\times _\mathcal{1})/𝒸=(𝒰_\mathcal{2}\times _\mathcal{2})/𝒸`$, $`q`$ is the particle electric charge and $`L`$ is the distance the particle drifted along the shock ($`\mathrm{\Delta }E>0`$). With the above expressions for the perpendicular shock the equation (3.8) can be rewritten as

$$\frac{\mathrm{\Delta }E}{E}\frac{U}{c}\frac{L}{r_g}.$$
$`(3.9)`$
In order to double the original energy in this process, $`\mathrm{\Delta }E=E`$, an average particle has to move a distance $`L=(c/U)r_g`$ along the shock surface. Thus, if the shock extension available for drifting particles $`L=\mathrm{\Delta }Y_{max}/2`$, the maximum energy of particles accelerated by the shock drift mechanism is provided by the condition

$$r_g\frac{\mathrm{\Delta }Y_{max}}{2}\frac{U}{c},$$
$`(3.10)`$
analogous to the above condition (3.7).

## 4 Conclusions

The maximum energy a particle can gain in the drift - diffusion acceleration expressed by the famous Hillas (1984) formula needs more specification concerning the characteristic spatial scales of the shock. Apart from the shock normal extension (Bhattacharjee & Sigl 2000) one also needs information about its sideway extent ($`\mathrm{\Delta }X_{max}`$ and $`\mathrm{\Delta }Y_{max}`$) influencing the maximum energy and, at highest energies, the spectral shape of the accelerated particles. As an illustration let us consider an accretion flow with velocity $`U=300`$ km/s, the size of $`h3`$ Mpc along the flow and of $`L10`$ Mpc to the sides, harbouring a shock wave inside with, say, $`\mathrm{\Delta }X_{max}`$ and/or $`\mathrm{\Delta }Y_{max}L`$ . In the perpendicular shock configuration the plasma normal extent provides the upper limit for particle energy expressed as $`r_g<h`$. On the other hand particle escape to the sides provide the limits (3.7) and (3.10). If the magnetic field in the plasma equals $`0.1`$ $`\mu `$G, the former constraint requires approximately that $`E_{max}<10^{20}`$ eV, while the others provide the limit $`E_{max}10^{17}`$ eV. When varying the above accretion flow parameters within the reasonable ranges one finds that the upper energy limit for particles accelerated at such shocks can hardly exceed a scale of EeV, and in most cases is close to $`10^{17}`$ eV. This low $`E_{max}`$ results from the fact that for the weakly perturbed magnetic field required for fast acceleration at the considered perpendicular shock, the accelerated particles can diffuse and drift far along the finite shock to escape to its sides. For highly turbulent conditions and/or the parallel shock configuration a similar value for $`E_{max}`$ arises due to the slow acceleration constrained by radiation losses (cf. Kang et al. 1997).

One should note that obtaining the same expressions for constraints (3.7) and (3.10) is not a coincidence. The upper limit (3.7) can be rewritten as $`\lambda _{}\mathrm{\Delta }X_{max}/2`$, while the limit (3.10) means $`L\mathrm{\Delta }Y_{max}/2`$ . So, for the highest energy particles their dispersion along the shock can be at most comparable to their mean free path or their maximum drift length. Such particles stream along the magnetic field nearly free, advancing their positions a distance $`r_g`$ in each gyroperiod, and they drift across the magnetic field changing their $`y`$ coordinate $`r_g`$ at each gyroperiod, as well, provided $`B_2B_1`$. Thus for a given minimum acceleration time (3.2) or (3.4) such highest energy particles can move comparable distances along and across the magnetic field at the shock.

###### Acknowledgements.

We are grateful to the anonymous referee for useful suggestions. The work was supported by the Komitet Badań Naukowych through the grant PB 258/P03/99/17 (MO) and within the project 2 P03B 112 17 (GSO).


format: text

text:
# Aging in ferromagnetic systems at criticality near four dimensions

## I Introduction

The time evolution of a system relaxing towards equilibrium is characterized by two different regimes: a transient behavior with off-equilibrium evolution, for $`t<t_{\mathrm{eq}}`$, and a stationary equilibrium evolution for $`t>t_{\mathrm{eq}}`$. In the former a dependence of the behavior of the system on initial conditions is expected, while in the latter homogeneity of time and time reversal symmetry (at least in absence of external fields) are recovered; dynamics of fluctuations are thus described in terms of “equilibrium” dynamics, with a characteristic time scale diverging at the critical point (critical slowing down) .

Consider a ferromagnetic model in a disordered state for the initial time $`t=0`$, and quench it at its critical temperature. During the relaxation a small external field $`h`$ is applied at $`𝐱=0`$ after a waiting time $`s`$. At time $`t`$ the order parameter response to $`h`$ is given by the response function $`R_𝐱(t,s)=\delta \varphi _𝐱(t)/\delta h(s)`$, where $`\varphi `$ is the order parameter and $``$ stands for the mean over stochastic dynamics. Since the system does not reach the equilibrium this function will depend both on $`s`$ (the “age” of the system) and $`t`$. This behavior is usually referred to as aging and was firstly noted in spin glass systems .

To characterize the distance from equilibrium of an aging system, evolving at a fixed temperature $`T`$, the fluctuation-dissipation ratio (FDR) is usually introduced :

$$X_𝐱(t,s)=\frac{TR_𝐱(t,s)}{_sC_𝐱(t,s)},$$
(1)
where $`C_𝐱(t,s)=\varphi _𝐱(t)\varphi _\mathrm{𝟎}(s)`$, the two-time correlation function. When the waiting time $`s`$ is greater than $`t_{\mathrm{eq}}`$ the dynamics is homogeneous in time and the fluctuation-dissipation theorem leads to $`X_𝐱(t,s)=1`$. This is no longer true in the aging regime .

In recent years, several works  have been devoted to the study of the FDR for systems exhibiting domain growth , or for aging systems such as glasses and spin glasses, showing that in the low-temperature phase $`X(t,s)`$ turns out to be a non-trivial function of its two arguments. In particular analytical and numerical studies indicate that the limit

$$X_{𝐱=0}^{\mathrm{}}=\underset{s\mathrm{}}{lim}\underset{t\mathrm{}}{lim}X_{𝐱=0}(t,s),$$
(2)
vanishes throughout the low-temperature phase both for spin glass and simple ferromagnetic system .

Only recently  attention has been paid to the FDR, for non-equilibrium, non-disordered, and unfrustrated systems at criticality. It has been argued that the FDR (2) is a novel universal quantity of non-equilibrium critical dynamics. Correlation and response functions were exactly computed in the simple cases of a random walk, a free Gaussian field, and a two-dimensional XY model at zero temperature and the value $`X_{𝐱=0}^{\mathrm{}}=1/2`$ was found . The same problem has been addressed for the $`d`$-dimensional spherical model , for the one dimensional Ising-Glauber chain  and Monte Carlo simulation was done for the two dimensional Ising model . In all cases $`X_{𝐱=0}^{\mathrm{}}`$ has values ranging between $`0`$ and $`\frac{1}{2}`$ while for some urn models a different range has been found . Also the scaling form for $`R_{𝐱=\mathrm{𝟎}}(t,s)`$ was rigorously established using conformal invariance .

In this work we investigate the non-equilibrium correlation and response functions and the associated FDR for the $`O(N)`$ ferromagnetic model with purely dissipative relaxation dynamics (Model A of Ref. ) both at the critical point and in the high-temperature phase, using a field-theoretical approach (never applied so far — see concluding remarks of Ref. ), at first order in an $`ϵ`$-expansion.

The paper is organized as follows. In Section II we briefly introduce the model, the scaling forms and the Gaussian (mean-field) result. In Section III we derive the $`ϵ`$-expansion for the FDR for all values of $`s`$ and $`t`$. Finally in Section IV we summarize our results and discuss some points needed further investigation. In Appendix we report some useful details on the one loop calculation.

## II The model

Let us consider the purely dissipative relaxation dynamics of a $`N`$-component field $`\phi (𝐱,t)`$ described by the stochastic Langevin equation (Model A of Ref. )

$$_t\phi (𝐱,t)=\mathrm{\Omega }\frac{\delta [\phi ]}{\delta \phi (𝐱,t)}+\xi (𝐱,t),$$
(3)
where $`[\phi ]`$ is the Landau-Ginzburg Hamiltonian

$$[\phi ]=\mathrm{d}^dx\left[\frac{1}{2}(\phi )^2+\frac{1}{2}r_0\phi ^2+\frac{1}{4!}g_0\phi ^4\right],$$
(4)
$`\mathrm{\Omega }`$ the kinetic coefficient, and $`\xi (𝐱,t)`$ a zero-mean stochastic Gaussian noise with

$$\xi _i(𝐱,t)\xi _j(𝐱^{},t^{})=2\mathrm{\Omega }\delta (𝐱𝐱^{})\delta (tt^{})\delta _{ij}.$$
(5)
The equilibrium correlation functions, generated by the Langevin equation (3) and averaged over the noise $`\xi `$, can be obtained by means of the field-theoretical action

$$S[\phi ,\stackrel{~}{\phi }]=dt\mathrm{d}^dx\left[\stackrel{~}{\phi }\frac{\phi }{t}+\mathrm{\Omega }\stackrel{~}{\phi }\frac{\delta [\phi ]}{\delta \phi }\stackrel{~}{\phi }\mathrm{\Omega }\stackrel{~}{\phi }\right].$$
(6)
where $`\stackrel{~}{\phi }(𝐱,t)`$ is the response field.

In Ref.  this formalism was extended in order to incorporate a macroscopic initial condition into Eq. (6): one has also to average over the initial configuration $`\phi _0(𝐱)=\phi (t=0,𝐱)`$ with a weight $`e^{H_0[\phi _0]}`$ given by

$$H_0[\phi _0]=\mathrm{d}^dx\frac{\tau _0}{2}(\phi _0(𝐱)a(𝐱))^2.$$
(7)
This specifies an initial state $`a(𝐱)`$ with correlations proportional to $`\tau _0^1`$. In this way all response and correlation functions may be obtained, following standard methods , by a perturbative expansion of the functional weight $`e^{(S[\phi ,\stackrel{~}{\phi }]+H_0[\phi _0])}`$.

The propagators (Gaussian two point correlation and response functions) of the resulting theory are

$`\stackrel{~}{\phi _i}(𝐪,s)\phi _j(𝐪,t)_0=`$ $`\delta _{ij}R_q^0(t,s)=`$ $`\delta _{ij}\theta (ts)G(ts),`$ (8)
$`\phi _i(𝐪,s)\phi _j(𝐪,t)_0=`$ $`\delta _{ij}C_q^0(t,s)=`$ $`{\displaystyle \frac{\delta _{ij}}{q^2+r_0}}\left[G(|ts|)+\left({\displaystyle \frac{r_0+q^2}{\tau _0}}1\right)G(t+s)\right],`$ (9)
where

$$G(t)=e^{\mathrm{\Omega }(q^2+r_0)t}.$$
(10)
It has also been shown that $`\tau _0^1`$ is irrelevant (in the renormalization group sense) for large times behavior .

### A Scaling forms

When a ferromagnetic system is quenched from a disordered initial state to its critical point, the correlation length grows as $`t^{1/z}`$, where $`z`$ is the dynamical critical exponents . So in momentum space, applying standard scaling arguments, all the universal functions depend only on the two products $`q^zt`$ and $`q^zs`$.

In particular we expect the scaling forms

$`R_q(t,s)=q^{2+\eta +z}\left({\displaystyle \frac{t}{s}}\right)^\theta F_R(\mathrm{\Omega }q^z(ts),t/s),`$ (11)
$`C_q(t,s)=q^{2+\eta }\left({\displaystyle \frac{t}{s}}\right)^\theta F_C(\mathrm{\Omega }q^z(ts),t/s),`$ (12)
where $`\theta `$ is the initial-slip exponent of response function related to the initial-slip exponent of the magnetization $`\theta ^{}`$ and to the autocorrelation exponent $`\lambda _c`$  by the relation

$$\theta ^{}=\theta +z^1(2z\eta )=z^1(d\lambda _c).$$
(13)
The functions $`F_R(y,x)`$ and $`F_C(y,x)`$ are universal apart the normalizations for small arguments. These functions are regular functions of both arguments, and for large $`x`$ they behave as:

$`F_R(y,x)=F_R^{\mathrm{}}(y)+O(x^1),`$ (14)
$`F_C(y,x)={\displaystyle \frac{F_C^{\mathrm{}}(y)}{x}}+O(x^2),`$ (15)
so that, for $`s0`$, these scaling forms reduce to ones of Ref. . We would also mention that transforming Eqs. (11,12) in the real $`𝐱`$ space (to this end we have to assume that $`F_R`$ and $`F_C`$ are rapidly decreasing functions of $`y`$ for $`y\mathrm{}`$) and setting $`𝐱=0`$ one easily obtain the scaling forms reported in Ref. . They also reduce to the equilibrium ones  when $`ts1`$.

Let us introduce

$$𝒳_𝐪=\frac{\mathrm{\Omega }R_𝐪(t,s)}{_sC_𝐪(t,s)}$$
(16)
then, from the scaling forms above is simple to show that, assuming $`F_R(y,x)=O(y)`$, $`x`$, the FDR may be written as

$$𝒳_{𝐪=0}^{\mathrm{}}=\underset{s\mathrm{}}{lim}\underset{t\mathrm{}}{lim}𝒳_{𝐪=0}(t,s)=(1\theta )^1\frac{F_R^{\mathrm{}}(0)}{F_C^{\mathrm{}}(0)}.$$
(17)
### B Gaussian FDR

For the Gaussian model we know exactly the response and correlations functions, so we can evaluate the FDR (in  the related quantity $`X_𝐱`$ has been considered, see Sec. IV). From Eqs. (8, 9) and definition (16) we have

$$𝒳_q^0(t,s)=\left(\frac{_sC_q^0}{\mathrm{\Omega }G_q^0}\right)^1=\left(1+e^{2\mathrm{\Omega }(q^2+r_0)s}+\mathrm{\Omega }q^2\tau _0^1e^{2\mathrm{\Omega }(q^2+r_0)s}\right)^1.$$
(18)
If the theory is off-critical ($`r_00`$) the limit of this ratio for $`s\mathrm{}`$ is $`1`$ for all values of $`q`$, in agreement with the idea that in the high-temperature phase all modes have a finite equilibration time, so that equilibrium is recovered and as a consequence the fluctuation-dissipation theorem holds. For the critical theory, i.e. $`r_0TT_c=0`$, if $`q0`$ the limit ratio is again equal to one, whereas for $`q=0`$ we have $`𝒳_{q=0}^0(t,s)=1/2`$. This analysis clearly shows that the only mode characterized by aging, i.e. that “does not relax” to the equilibrium, is the zero mode in the critical limit.

## III One-loop FDR

The aim of this section is the computation of the non-equilibrium response and correlation functions for the purely dissipative dynamics of the $`N`$-vector model at one-loop order. We use here the method of renormalized field theory in the minimal subtraction scheme. The breaking of homogeneity in time gives rise to some technical problems in the renormalization procedure in terms of one-particle irreducible correlation functions (see  and references therein) so our computation is done in terms of connected functions.

At one-loop order we have to evaluate, taking also into account causality , the three Feynman diagrams in Figure 1, one for the response function and two for the correlation one. In terms of these diagrams we have

$`R_q(t,s)=R_q^0(t,s){\displaystyle \frac{N+2}{6}}g_0(a)+O(g_0^2),`$ (19)
$`C_q(t,s)=C_q^0(t,s){\displaystyle \frac{N+2}{6}}g_0((b)+(c))+O(g_0^2).`$ (20)
In order to evaluate the FDR at criticality we have to set in this perturbative expansion $`r_0=0`$ (massless theory). We also set $`\tau _0^1=0`$, since it is an irrelevant variable , and $`\mathrm{\Omega }=1`$ to lighten the notations. The first step in the calculation of the diagrams is the evaluations of the critical “bubble” $`B_c(t)`$, i.e. their common one-particle irreducible part. We have, in generic dimension $`d`$,

$$B_c(t)=\frac{\mathrm{d}^dq}{(2\pi )^d}C_q(t,t)=\frac{1}{d/21}\frac{(2t)^{1d/2}}{(4\pi )^{d/2}}=N_d\frac{\mathrm{\Gamma }(d/21)}{2^{d/2}}t^{1d/2},$$
(21)
where $`N_d={\displaystyle \frac{2}{(4\pi )^{d/2}\mathrm{\Gamma }(d/2)}}`$. Note that the equilibrium contribution to $`B_c(t)`$ is zero for $`d>2`$.

Let us consider $`t>s`$ in the following. We may write

$`(a)={\displaystyle _0^{\mathrm{}}}dt^{}R_q^0(t,t^{})B(t^{})R_q^0(t^{},s),`$ (22)
$`(b)={\displaystyle _0^{\mathrm{}}}dt^{}R_q^0(t,t^{})B(t^{})C_q^0(t^{},s),`$ (23)
$`(c)={\displaystyle _0^{\mathrm{}}}dt^{}R_q^0(s,t^{})B(t^{})C_q^0(t^{},t),`$ (24)
where we set $`r_0=0`$ in $`R_q^0`$ and $`C_q^0`$.

Performing the integration and expanding in powers of $`ϵ`$ we find for the response function

$$R_q(t,s)=G(ts)\left(1+\stackrel{~}{g}_0\frac{N+2}{24}\mathrm{ln}\frac{t}{s}\right)+O(ϵ^2,\stackrel{~}{g}_0^2),$$
(25)
and for the correlation function

$$C_q(t,s)=\frac{G(ts)G(t+s)}{q^2}\left(1+\stackrel{~}{g}_0\frac{N+2}{24}\mathrm{ln}\frac{t}{s}\right)\stackrel{~}{g}_0\frac{N+2}{24}\frac{G(t+s)}{q^2}f(2q^2s)+O(ϵ^2,\stackrel{~}{g}_0^2),$$
(26)
where

$$f(v)=2\left[_0^vd\xi \mathrm{ln}\xi e^\xi +(1e^v)\mathrm{ln}v\right],$$
(27)
and $`\stackrel{~}{g_0}=N_dg_0`$. Note that $`f(0)=0`$, $`f^{}(0)=2`$ and $`f(v)`$ has the following asymptotic expansion, for $`v1`$

$$f(v)=2\frac{e^v}{v}\left(1+\frac{1}{v}+\frac{2}{v^2}+\mathrm{}+\frac{k!}{v^k}+\mathrm{}\right).$$
(28)
In order to obtain the critical functions we have to set the renormalized coupling equal to its fixed point value. At first order in $`ϵ`$

$$\stackrel{~}{g}_0=\stackrel{~}{g}^{}=\frac{6}{N+8}ϵ+O(ϵ^2).$$
(29)
Finally we get (called $`P_N=\frac{N+2}{N+8}`$)

$$R_q(t,s)=G(ts)\left(1+ϵ\frac{P_N}{4}\mathrm{ln}\frac{t}{s}\right)+O(ϵ^2),$$
(30)
$$C_q(t,s)=\frac{G(ts)G(t+s)}{q^2}\left(1+ϵ\frac{P_N}{4}\mathrm{ln}\frac{t}{s}\right)ϵ\frac{P_N}{4}\frac{G(t+s)}{q^2}f(2q^2s)+O(ϵ^2),$$
(31)
that are fully compatible with the scaling form given in the previous section, with

$$F_R(y,x)=e^y+O(ϵ^2),$$
(32)
and

$$F_C(y,x)=e^y\left[1+ϵ\frac{P_N}{4}f\left(\frac{2y}{x1}\right)\right]e^{y\frac{x+1}{x1}}+O(ϵ^2).$$
(33)
In particular we recognize the exponent $`\theta =P_Nϵ/4+O(ϵ^2)`$ in agreement with Ref. , $`z=2+O(ϵ^2)`$, $`\eta =O(ϵ^2)`$ as expected, and that $`F_R(y,x)`$ is not affected by $`O(ϵ)`$ corrections. It is also easy to find that

$$F_C^{\mathrm{}}(y)=2y\left(1+ϵ\frac{P_N}{2}\right)e^y+O(ϵ^2).$$
(34)
Computing the derivative with respect to $`s`$ of the two-time correlation function and taking its ratio with the response function we have

$$𝒳_q^1(s)=1+e^{2q^2s}\frac{P_Nϵ}{4}e^{2q^2s}\left[\frac{e^{2q^2s}1}{q^2s}f(2q^2s)+2f^{}(2q^2s)\right]+O(ϵ^2).$$
(35)
Note that, at least at this order, the result is independent of the observation time $`t`$. Using the large $`v`$ behavior of $`f(v)`$, cf. Eq. (28), we find that the limit of the FDR for $`s\mathrm{}`$ is equal to 1 for all $`q0`$. Instead for $`q=0`$ we have (using (27))

$$𝒳_{q=0}^{\mathrm{}}=\frac{1}{2}\left(1ϵ\frac{P_N}{4}\right)+O(ϵ^2),$$
(36)
in agreement with Eq. (17) and with the scaling forms (32) and (34).

Taking into account the effect of the mass $`r_0`$ (deviation from critical temperature) in the previous computations, one obtains for the non-critical bubble (contributing to the mass renormalization)

$$B(t)=N_d\left[\frac{\pi }{2\mathrm{sin}d\pi /2}\frac{1}{2}\mathrm{\Gamma }(d/2)\mathrm{\Gamma }(1d/2,2r_0t)\right]r_0^{d/21},$$
(37)
where $`\mathrm{\Gamma }(x,y)`$ is the incomplete $`\mathrm{\Gamma }`$ function . Using this expression it is possible to determine, as previously done, correlation and response functions. We report the basic formulas in the Appendix. The final result is obtained computing the ratio $`𝒳_q`$ in terms of the renormalized parameters of the theory. It is then trivial, but algebraically cumbersome, to show that $`𝒳_q^{\mathrm{}}`$ is equal to $`1`$ for all $`q`$ in the high temperature phase.

## IV Discussion

In this work we considered the off-equilibrium properties of the purely dissipative relaxational dynamics of an $`N`$-vector model in the framework of field theoretical $`ϵ`$-expansion. We computed at first order in $`ϵ`$ the FDR, as defined in (16) as a function of the waiting time $`s`$ and of the observation time $`t`$ both at criticality and in the high-temperature phase. The main result is that the ratio $`𝒳_q^{\mathrm{}}`$ is always $`1`$ unless at criticality for $`q=0`$, when it takes the value

$$𝒳_{𝐪=0}^{\mathrm{}}=\frac{1}{2}\left(1\frac{ϵ}{4}\frac{N+2}{N+8}\right)+O(ϵ^2).$$
(38)
To compare our result with some particular limit considered in the literature  we have to relate this quantity to the analog in the real $`𝐱`$ space. The following heuristic argument may be useful to realize that the two ratios are exactly equal, i.e.

$$X_{𝐱=0}^{\mathrm{}}=𝒳_{𝐪=0}^{\mathrm{}}.$$
(39)
We may rewrite the FDR in real $`𝐱`$ space as a mean value of that in momentum space with a weight given by $`R_𝐪`$:

$$X_{𝐱=0}^1\frac{\mathrm{d}^dq_sC_𝐪(t,s)}{T\mathrm{d}^dqR_𝐪(t,s)}=\frac{\mathrm{d}^dqR_𝐪(t,s)\frac{_sC_𝐪(t,s)}{TR_𝐪(t,s)}}{\mathrm{d}^dqR_𝐪(t,s)}=𝒳_𝐪^1_{R_𝐪}.$$
(40)
Now, since we expect $`R_𝐪e^{q^2(ts)}`$, in the limit $`s,t\mathrm{}`$ (in the right order) $`X_{𝐱=0}^1`$ will take contributions only for the $`q=0`$ mode, i.e. apart a normalization, the weight function $`R_𝐪`$ is a $`\delta (𝐪)`$. However we note that at the first order in $`ϵ`$ the equality (39) is identically satisfied by our result since $`R_𝐪e^{q^2(ts)}`$ , cf. Eq. (25).

In the limit $`N\mathrm{}`$ Eq. (38) reduces to $`X^{\mathrm{}}=1/2ϵ/8+O(ϵ^2)`$ that is the same as the expansion of the result for the spherical model near four dimension .

A Monte Carlo simulation of the two dimensional Ising model gave for the FDR the value $`X^{\mathrm{}}=0.26(1)`$ , qualitatively in agreement with our result for $`ϵ=2`$, $`X^{\mathrm{}}=5/12<1/2`$. To have a reliable quantitative prediction the knowledge of higher loop contributions is required.

Setting $`ϵ=1`$ for $`N=1`$, one obtains $`11/24`$. This number, that is a rough estimate of the actual three-dimensional value, may be measured in numerical or experimental works.

This work may be easily extended to more realistic models than those previously considered in literature, contributing to the understanding of out-of-equilibrium dynamic phenomena, currently under intensive investigation, by means of the powerful tools of perturbative field theory.

## Acknowledgments

We would thank S. Caracciolo, A. Pelissetto and E. Vicari for a critical reading of the manuscript and useful comments. We acknowledge in a special way A. Pelissetto for having introduced us to these problems.

## Details of computations

We summarize here the main analytical results useful for the computation of correlation and response functions for $`r_00`$ at one loop. Again one has to perform all the needed integrations over the times, as in Eq. (23) with the free field correlator and response function (89). At variance with the critical theory a renormalization of the parameter $`r_0`$ is now required to cancel dimensional poles both in $`R_q`$ and $`C_q`$.

Let us introduce the function

$$Y(t)_0^td\tau B(\tau ),$$
(41)
and

$$W(t)_0^td\tau G(2\tau )B(\tau ),$$
(42)
where $`G(t)`$ and $`B(t)`$ are given in Eq. (10) and (37), respectively. In terms of $`Y`$ and $`W`$ we obtain (for $`t>s`$ and $`\tau _0^1=0`$, $`\mathrm{\Omega }=1`$)

$$(a)=G(ts)[Y(t)Y(s)],$$
(43)
and

$`(b)+(c)={\displaystyle \frac{1}{q^2+r_0}}\left\{G(ts)[Y(t)Y(s)]G(t+s)[Y(t)+Y(s)]+2G(t+s)W(s)\right\}.`$ (44)
(45)
In the following with $`Y`$ and $`W`$ we mean also their analytic continuation in $`d`$.

An explicit computation leads to

$$Y(t)=\frac{r_0^{d/22}}{2(4\pi )^{d/2}}\left\{(2r_0t+d/21)[\mathrm{\Gamma }(1d/2)\mathrm{\Gamma }(1d/2,2r_0t)]+(2r_0t)^{1d/2}e^{2r_0t}\right\},$$
(46)
and

$`W(t)=`$ (48)
$`={\displaystyle \frac{1}{2(4\pi )^{d/2}}}{\displaystyle \frac{r_0^{d/21}}{q^2+r_0}}\left\{G(2t)[\mathrm{\Gamma }(1d/2)\mathrm{\Gamma }(1d/2,2r_0t)](q^2/r_0)^{d/21}\mathrm{\Delta }(1d/2,2q^2t)\right\},`$
where we introduced

$$\mathrm{\Delta }(v,w)_0^wd\tau \tau ^{v1}e^\tau $$
(50)
(for $`v0`$ its analytic continuation has to be considered).

Expanding (46,LABEL:expW) in $`ϵ=4d`$, we obtain

$$2(4\pi )^{d/2}Y(t)=\frac{2}{ϵ}(2r_0t+1)(2r_0t+1)[\gamma (2r_0t)\mathrm{ln}r_0]+1+\frac{e^{2r_0t}}{2r_0t}+O(ϵ),$$
(51)
and

$`2(4\pi )^{d/2}(q^2+r_0)W(t)=`$ (53)
$`={\displaystyle \frac{2}{ϵ}}[r_0G(2t)+q^2]+q^2[\mathrm{ln}q^2\delta (2q^2t)]r_0G(2t)[\gamma (2r_0t)\mathrm{ln}r_0]+O(ϵ),`$
where

$$\gamma (v)1+e^v\left(\mathrm{ln}v+\frac{1}{v}\right)+_0^vd\xi \mathrm{ln}\xi e^\xi ,$$
(54)
and

$$\delta (v)1+e^v\left(\mathrm{ln}v\frac{1}{v}\right)_0^vd\xi \mathrm{ln}\xi e^\xi .$$
(55)
It is easy to find that $`f(v)`$ in Eq. (27) is related to $`\delta (v)`$ by

$$f(v)=2\left[1+\mathrm{ln}v\delta (v)\frac{e^v}{v}\right].$$
(56)
Plugging Eqs. (51) and (53) into Eq. (43) and (45) and then into Eqs. (20) it is easy to realize that to cancel dimensional poles both in $`R_q(t,s)`$ and $`C_q(t,s)`$ a renormalization of the bare mass $`r_0`$ is sufficient (at least in the case $`\tau _0^1=0`$ we are considering)

$$r_0=Z_rr\text{with}Z_r=1+\frac{N+2}{3}\frac{g_0}{(4\pi )^{d/2}}\frac{1}{ϵ}+O(g_0^2),$$
(57)
in agreement with what one would expect from the corresponding static field-theory (see, for instance, Ref. ). All the previously stated results easily follow from explicit expressions given above.


text:
# 1 Introduction

## 1 Introduction

A particle in a viscous medium, in the presence of a spatially periodic potential, can be driven by a time-periodic force at a nonzero average velocity $`v_0`$ determined by the resonance condition,

$$2\pi m/\omega _0=l/v_0,$$
(1)
where $`\omega _0`$ is the driving frequency, an integer $`m`$ is the order of the resonance, and $`l`$ is the period of the potential. This phenomenon was experimentally observed in the form of phase-rotating states (Shapiro steps) in small Josephson Junctions (JJ) driven by ac bias current . Later, it was shown that the same effect can also be realized for a fluxon (magnetic flux quantum) moving in a periodically modulated long JJ under the action of the ac bias current . The latter effect has been recently experimentally observed, in the form corresponding to both the fundamental resonance \[$`m=1`$ in Eq. (1)\] and higher-order ones, in a long circular JJ with an effective spatial modulation induced by uniform dc magnetic field .

In terms of the phase difference $`\theta `$ of the superconducting wave function across the junction, an ac-driven small JJ is described by a pendulum equation (for a detailed derivation see the book )

$$\frac{d^2\theta }{dt^2}+\mathrm{sin}\theta =\alpha \frac{d\theta }{dt}+ϵ\mathrm{cos}\left(\omega _0t+\psi (t)\right)+j(t).$$
(2)
Here, time is measured in units of the inverse Josephson plasma frequency, $`\alpha `$ is the normalized JJ conductance, $`ϵ`$ is the ac-current amplitude normalized to the Josephson critical current, $`\omega _0`$ is the normalized driving frequency, $`\psi `$ is an arbitrary phase, and $`j(t)`$ represents intrinsic thermal noise in JJ. The equation of motion for the fluxon in the above-mentioned long circular JJ in the magnetic field takes exactly the same form in the “nonrelativistic” limit, i.e., if the fluxon’s velocity is much smaller than the Swihart velocity of the junction. The resonant relation (1), where $`l=2\pi `$, implies that, neglecting the noise, the ac drive in Eq. (2) may support phase rotation of the pendulum at the average velocity $`d\theta /dt=v_0`$ in the presence of the friction. Below, we assume that the spatial modulation period is always normalized so that $`l2\pi `$, i.e., $`v_0=\omega _0/m`$. We note that $`d\theta /dt`$ is proportional to the voltage across JJ, hence the phase rotation in the ac-driven JJ gives rise to a nonzero dc voltage, a feature that is used in Josephson voltage standards (see, e.g., Ref.  and references therein).

In real applications, the driving ac signal is always slightly nonmonochromatic, having a finite width $`\delta \omega `$ in the spectral domain; in other words, $`\psi `$ in Eq. (2) is not a constant phase, but rather a slowly varying function of time, representing random phase fluctuations of the driving signal. This introduces a finite lifetime $`\mathrm{}`$ of the phase-locked state, which is of direct relevance to applications, affecting the stability of the Josephson voltage standards.

Intrinsic thermal fluctuations, represented by the term $`j(t)`$ in Eq. (2), also contribute to dephasing of the ac-driven motion. In terms of the small ac-driven JJ, the effect of thermal fluctuations was considered in earlier works , where Eq. (2) with the monochromatic drive and thermal noise was reduced to a Langevin equation for a particle driven by a random force in a periodic potential. The phase-locked state is then represented by the particle trapped at a minimum of the potential, and the dephasing implies that the particle is extracted by the random force from the trapped state. The corresponding dephasing time was taken as the inverse Kramers’ escape rate , i.e.,

$$\mathrm{}\mathrm{exp}\left(\mathrm{\Delta }U/T\right),$$
(3)
where $`\mathrm{\Delta }U`$ is the difference between maximum and minimum values of the effective potential, and $`T`$ is the temperature (see exact definitions below).

In this work, we focus on effects of the frequency fluctuations in the ac drive. In section 2, we demonstrate that Eq. (2) with a finite linewidth of the ac signal can be transformed, at the first order of the perturbation theory, into an equation driven by a strictly monochromatic signal and an additive random force. However, in contrast to the random force representing the intrinsic thermal noise \[$`j(t)`$ in Eq.(2)\], the correlator of the effective additive noise generated by the ac-drive’s nonmonochromaticity does not contain the friction coefficient $`\alpha `$, as this correlator, which has a non-thermal origin, does not obey the fluctuation-dissipation theorem. Next, using the energy-balance technique , in section 3 we reduce the monochromatically driven equation with the random force to a Langevin equation in a periodic potential $`U(\theta )`$, which makes it possible to predict the dephasing time by means of the Kramers’ expression (3).

In section 4, we specially consider the situation when the small amplitude $`ϵ`$ of the driving signal is close to the phase-locking threshold $`ϵ_{\mathrm{t}hr}`$ . By means of the Fokker-Planck equation corresponding to the above-mentioned effective Langevin equation , we demonstrate that a critical value $`\delta \omega `$ of the driving-signal’s linewidth, at which a sharp transition from the locking to unlocking (in the form of random $`2\pi `$ phase slips) takes place, may be represented as a universal function of the drive’s amplitude $`ϵ`$.

In section 5, we present results of direct numerical simulations of Eq. (2) with the nonmonochromatic drive, which are reported in the form showing the logarithm of the dephasing time $`\mathrm{}`$ as a function of $`ϵ`$ and of the linewidth $`\delta \omega `$. If $`ϵ`$ is above the threshold value $`ϵ_{\mathrm{t}hr}`$, and is not too large, the numerically found lifetime $`\mathrm{}`$ is found to be quite close to that predicted by the perturbation theory. However, at large values of $`ϵ`$ the simulations reveal a new effect, which cannot be predicted by the perturbative analysis: $`\mathrm{}(ϵ)`$ reaches a maximum value and then decreases. As the nonmonotonic character of the dependence $`\mathrm{}(ϵ)`$ and the existence of the maximum in it are quite important features, in section 6 we report results of direct simulations of the mode with the strictly monochromatic drive and intrinsic thermal (additive) noise. We conclude that the dependence $`\mathrm{}(ϵ)`$ in this case has the same nonmonotonic character. Although the latter model was studied earlier , this feature was not reported.

It is relevant to mention that, besides small and long JJs, essentially the same dynamical model as the one considered in this work applies to ensembles of oscillators coupled via a mean field, which may be laser arrays or biological oscillators . As is known, the global coupling may synchronize the phase-rotation states of the oscillators, each of them being driven by the mean field. On the other hand, various perturbations affecting the mean field make it a slightly nonmonochromatic drive . Thus, desynchronization of the globally coupled rotating oscillators is another manifestation of the problem considered in this work.

## 2 Transformation of frequency fluctuations into an additive noise

We begin our analysis from Eq. (2) without intrinsic thermal noise, i.e., with $`j=0`$. White-noise fluctuations of the ac-drive’s frequency, $`\omega (t)d\psi /dt+\omega _0`$, are assumed to be subject to the Gaussian correlations:

$$\frac{d\psi (t)}{dt}\frac{d\psi (t^{})}{dt^{}}=2\mathrm{\Omega }\delta (tt^{}),$$
(4)
$`\mathrm{\Omega }`$ being intensity of the fluctuations. The relative linewidth of the ac drive, $`\delta \omega \left(\omega \omega _0\right)/\omega _0`$ (calculated at $`3`$ dB level), can be then estimated as $`\delta \omega 2.4\mathrm{\Omega }/\omega _0`$. We note that even low-quality sources of radiofrequency radiation, that may be used as the ac drive for JJ, have $`\delta \omega _{}^<\mathrm{\hspace{0.17em}\hspace{0.17em}10}^3`$, while the dissipative constant $`\alpha `$ in JJ, although being small, is normally in the range $`\alpha _{}^>\mathrm{\hspace{0.17em}\hspace{0.17em}10}^2`$, therefore we hereafter assume $`\delta \omega \alpha `$. In other words, we may assume that a characteristic time of the variation of the random ac-drive’s phase shift $`\psi (t)`$ is much larger than the relaxation time $`1/\alpha `$.

To convert the frequency fluctuations into an effective additive noise, which is more convenient for the subsequent analysis, we define a new time variable which includes a slowly varying stochastic term,

$$\tau t+\chi (t),with\chi (t)\omega _0^1\psi (t).$$
(5)
Then, transforming the time derivatives $`d/dt`$ into $`d/d\tau `$ according to this definition, making use of the above relations $`\delta \omega \alpha 1`$, and keeping newly appearing small perturbations at two lowest orders, we cast the underlying Eq. (2) into a form

$$\frac{d^2\theta }{d\tau ^2}+\mathrm{sin}\theta =\alpha \frac{d\theta }{d\tau }+ϵ\mathrm{cos}\left(\omega _0\tau \right)2\frac{d\chi }{d\tau }\frac{d^2\theta }{d\tau ^2}\alpha \frac{d\chi }{d\tau }\frac{d\theta }{d\tau }.$$
(6)
In fact, the last term in Eq. (6) is much smaller than the previous one, as $`\alpha `$ is small, and, in the first approximation, one may substitute, in the latter term, $`d^2\theta /d\tau ^2\mathrm{sin}\theta `$. Thus, the final form of the perturbed equation, in which the frequency fluctuations were converted into the effective additive random force, is

$$\frac{d^2\theta }{d\tau ^2}+\mathrm{sin}\theta =\alpha \frac{d\theta }{d\tau }+ϵ\mathrm{cos}\left(\omega _0\tau \right)+2\frac{d\chi }{d\tau }\mathrm{sin}\theta .$$
(7)
The time transformation (5) affects the Gaussian correlator (4). It is easy to find that, in terms of the renormalized time and renormalized random phase $`\chi `$ (see Eq. (5)), an exact form of the correlator is

$$\frac{d\chi (\tau )}{d\tau }\left(1\frac{d\chi (\tau )}{d\tau }\right)^1\frac{d\chi (\tau ^{})}{d\tau ^{}}\left(1\frac{d\chi (t^{})}{d\tau ^{}}\right)^1=\frac{2\mathrm{\Omega }}{\omega _0^2}\frac{\delta (\tau \tau ^{})}{\left|1\frac{d\chi }{d\tau }\right|}.$$

However, in view of the smallness of the frequency fluctuations, in the lowest approximation we may adopt a simple form of the correlator,

$$\frac{d\chi (\tau )}{d\tau }\frac{d\chi (\tau ^{})}{d\tau ^{}}=\frac{2\mathrm{\Omega }}{\omega _0^2}\delta (\tau \tau ^{}).$$
(8)
Equations (7) and (8) will be a basis for further analysis, while numerical simulations will be run for the full underlying equation (2).

## 3 An effective Langevin equation and estimate for the lifetime of the phase-locked state

In the zero-order approximation, which implies $`ϵ=\alpha =d\chi /d\tau =0`$, Eq. (7) has a known solution,

$$\theta _0(t)=2\mathrm{a}m((\tau \tau _0)/k;k),$$
(9)
corresponding to the phase rotation at a nonzero average frequency (phase velocity)

$$\omega _0=\pi /kK(k).$$
(10)
Here, am is the Jacobi’s elliptic amplitude with the modulus $`k`$ ($`0<k<1`$ ), $`K(k)`$ and $`E(k)`$ are complete elliptic integral of the first and second kinds, and $`\tau _0`$ is an arbitrary constant. In this approximation, Eq. (2) conserves the energy, the value of which for the law of motion (9) is determined by the value of $`k`$,

$$\frac{1}{2}\left(\frac{d\theta _0}{d\tau }\right)^2\mathrm{cos}\theta _0=\frac{2}{k^2}1.$$
(11)
A possibility to support persistent phase rotation at a velocity $`\omega _0`$ in the presence of the friction by the monochromatic ac drive is predicted by the energy-balance equation . To this end, we calculate the net rate of the change of energy due to the action of the friction and drive, averaged over the rotation period $`2\pi /\omega _0`$, under the resonance condition given by Eq. (1):

$$\overline{\frac{d}{d\tau }}=\alpha \overline{\left(\frac{d\theta _0}{d\tau }\right)^2}+ϵ\omega _0\omega _1\mathrm{cos}(\omega _0\tau _0),$$
(12)
where the overbar stands for the time average, $`\tau _0`$ being the same constants as in Eq. (9), and $`\omega _1`$ is the amplitude of the resonant harmonic in the Fourier decomposition of the time-dependent velocity (instantaneous frequency) $`d\theta _0/d\tau `$, taken as per the unperturbed law of motion (9). An elementary calculation yields

$$\overline{\left(\frac{d\theta _0}{d\tau }\right)^2}=\frac{4E(k)}{k^2K(k)},\omega _1=\frac{2Q}{1+Q^2},Q\mathrm{exp}\left[\pi K\left(\sqrt{1k^2}\right)/K\right]$$
(13)
($`Q`$ is called the Jacobi parameter). Then, the phase-locked ac-driven regime, corresponding to $`\stackrel{\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}}{d/d\tau }`$ $`=0`$, is possible at two constant values (one stable and one unstable) of the phase difference between the ac drive and the rotating pendulum,

$$\omega _0\tau _0=\pm \mathrm{cos}^1\left[\frac{\alpha }{ϵ\omega _0\omega _1}\overline{\left(\frac{d\theta _0}{d\tau }\right)^2}\right]$$
(14)
, provided that the amplitude $`ϵ`$ exceeds the threshold value

$$ϵ_{\mathrm{t}hr}=\frac{\alpha }{\omega _0\omega _1}\overline{\left(\frac{d\theta _0}{d\tau }\right)^2}.$$
(15)
In the presence of the additive random force in Eq. (7), a perturbed equation of motion for the rotating pendulum can again be obtained from the averaged energy-balance equation, which has the same form of Eq. (12) with a difference that now $`\tau _0`$ is a slowly varying function of the time $`\tau `$ \[roughly speaking, varying as slowly as the random phase $`\psi (t)`$ in the underlying equation (2)\]. Notice that the time dependence of $`\tau _0`$ determines a change $`\delta \omega _0=\omega _0\left(d\tau _0/d\tau \right)`$ of the average phase-rotation velocity, and the latter may be related to the change of the energy (11), through its kinetic part, as $`\delta =^{}\delta \omega _0`$, where $`^{}`$ stands for $`d/d\omega _0`$, calculated for the unperturbed law of motion given by Eq. (9), $`^{}=4K^2/[\pi k(K+kdK/dk)]`$). Thus, $`\stackrel{\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}}{d/d\tau }=^{}\omega _0(d^2\tau _0/d\tau ^2)`$, and the balance equation (12) takes the form

$$\omega _0^{}\frac{d^2\tau _0}{d\tau ^2}=\alpha \overline{\left(\frac{d\theta _0}{d\tau }\right)^2}(1+2\frac{d\tau _0}{d\tau })+ϵ\omega _0\omega _1\mathrm{cos}\left(\omega _0\tau _0\right)+2\omega _0\frac{d\chi }{d\tau }\mathrm{sin}\left(\theta (\tau )\right),$$
(16)
the second term in $`(1+2d\tau _0/d\tau )`$ being a contribution to the energy dissipation rate due to the small change of the average velocity, while a similar correction to the last term in Eq. (16) is negligible. Eq. (16) can be transformed into a more convenient form by defining $`\omega _0\tau _0(\tau )\zeta (\tau )`$:

$$\omega _0^{}\frac{d^2\zeta }{d\tau ^2}+\mathrm{\hspace{0.17em}2}\alpha \overline{\left(\frac{d\theta _0}{d\tau }\right)^2}\frac{d\zeta }{d\tau }=[\omega _0\alpha \overline{\left(\frac{d\theta _0}{d\tau }\right)^2}+ϵ\omega _0^2\omega _1\mathrm{cos}\zeta ]+2\omega _0^2\mathrm{sin}\left(\theta (\tau )\right)\frac{d\chi }{d\tau }.$$
(17)
Thus, we have arrived at an effective Langevin equation  for a particle driven in a viscous medium by the sum of a regular force, represented by the terms in the square brackets, and a stochastic force, represented by the last term in the equation. The subsequent step is, following a well-known procedure , to introduce the Fokker-Planck (FP) equation corresponding to this Langevin equation, taking into regard the correlator (8). An essential feature of the thus derived FP equation is the presence of the extra multiplier $`\mathrm{sin}^2\left[\theta (\tau )\right]`$ in front of the diffusion (second-derivative) term in it, due to the multiplier $`\mathrm{sin}\left[\theta (\tau )\right]`$ in the stochastic-force term in Eq. (17). The coefficient $`\mathrm{sin}^2\left[\theta (\tau )\right]`$ may be averaged in time as per the unperturbed law of motion (9). It is easy to calculate the average value:

$$\overline{\mathrm{sin}^2\theta _0(\tau )}=\frac{4}{3k^4}\left[\left(2k^2\right)\frac{E(k)}{K(k)}2\left(1k^2\right)\right]$$
(18)
\[note that, in the limit $`k0`$, the expression (18) approaches an obvious value $`1/2`$\].

The FP equation takes essentially the same form as it would take in the known problem  of the dephasing of the ac-driven JJ phase rotation under the action of the additive thermal noise represented by the term $`j(t)`$ in Eq. (2) \[the most important difference of the effective Langevin equation (17) from its counterpart in the thermal-noise problem is the presence of the multiplier $`\mathrm{sin}\left[\theta (\tau )\right]`$ in the last term of the equation\]. With the FP equation taking the usual form, one can directly use the Kramers’ expression (3) to determine the lifetime of the ac-driven state, as it was done for the case of the thermal noise in Refs. . In particular, an effective potential corresponding to the potential force, i.e., the expression in square brackets in Eq. (17), is

$$U(\zeta )=ϵ\omega _0^2\omega _1\mathrm{sin}\zeta +\omega _0\alpha \overline{\left(\frac{d\theta _0}{d\tau }\right)^2}\zeta ,$$

hence the potential-barrier height $`\mathrm{\Delta }U`$, which should be substituted into Eq. (3), can be easily found as a difference of values of the potential taken between two points where the above-mentioned potential force vanishes. As a result, we obtain

$$\mathrm{\Delta }U=2ϵ\omega _0^2\omega _1\left[\sqrt{1\left(ϵ_{\mathrm{t}hr}/ϵ\right)^2}\left(ϵ_{\mathrm{t}hr}/ϵ\right)\mathrm{cos}^1\left(ϵ_{\mathrm{t}hr}/ϵ\right)\right],$$
(19)
where the definition (15) for the threshold value of the amplitude has been used to simplify the expression.

However, an additional difference of the present case from the thermal-noise problem, that must be taken into regard before applying the expression (3), is that the frequency fluctuation intensity appearing in the correlator (8) does not obey the fluctuation-dissipation theorem, and therefore it does not include the intrinsic dissipative constant $`\alpha `$ of the pendulum (JJ). By properly defining the effective temperature $`T_{\mathrm{e}ff}=4\omega _0^2\overline{\mathrm{sin}^2\theta _0(\tau )}\mathrm{\Omega }/\alpha _{\mathrm{e}ff}`$, where $`\alpha _{\mathrm{e}ff}2\alpha \overline{\left(d\theta _0/d\tau \right)^2}`$ is the effective friction constant from Eq. (17), and using the potential-barrier height $`\mathrm{\Delta }U`$ (19) and the expression (18), we can rewrite the Kramers’ expression (3) as:

$$\mathrm{}\mathrm{exp}\left(\frac{6ϵ\alpha k^2EQ}{\left[\left(2k^2\right)E2\left(1k^2\right)K\right]\left(1+Q^2\right)\mathrm{\Omega }}\left[\sqrt{1\left(ϵ_{\mathrm{t}hr}/ϵ\right)^2}\left(ϵ_{\mathrm{t}hr}/ϵ\right)\mathrm{cos}^1\left(ϵ_{\mathrm{t}hr}/ϵ\right)\right]\right).$$
(20)
This is the main prediction of the analytical consideration, which will be compared to results of direct simulations of Eq. (2) in section 5.

## 4 Dephasing the phase-locked state near the locking threshold

In this section we investigate the system described by the Langevin equation (17), by explicitly solving the associated FP equation. However, we can first simplify Eq. (17), recalling the fundamental physical condition according to which the frequency fluctuations are much smaller than $`\alpha `$, or, in other words, the random force varies on a time scale $`1/\alpha `$ . Consequently, the acceleration term on the left-hand side of Eq. (17) may be neglected as compared to the velocity term, which yields a simplified Langevin equation

$$\dot{\zeta }=F_0+F_1\mathrm{cos}\zeta +f(t),$$
(21)
where $`F_0\omega _0/2`$, $`F_1ϵ\omega _0^2\omega _1[2\alpha \overline{\left(d\theta _0/dt\right)^2}]^1`$, and

$$f(t)\omega _0^2\left[\alpha \overline{\left(\frac{d\theta _0}{d\tau }\right)^2}\right]^1\frac{d\chi }{d\tau }\mathrm{sin}\left[\theta (\tau )\right].$$
(22)
The FP equation (in this case, it is, in fact, the Smoluchowski equation ) for a probability distribution function $`P(\zeta ,t)`$, corresponding to Eq. (21) with the Gaussian correlator (8), is

$$P_t=F_1\mathrm{sin}\zeta P+\left(F_0F_1\mathrm{cos}\zeta \right)P_\zeta +\mathrm{\Xi }P_{\zeta \zeta },$$
(23)
where the subscripts stand for the partial derivatives, and

$$\mathrm{\Xi }\overline{\mathrm{sin}^2\theta _0(\tau )}\left[\alpha \overline{\left(\frac{d\theta _0}{d\tau }\right)^2}\right]^2\omega _0^2\mathrm{\Omega },$$
(24)
where $`\mathrm{\Omega }`$ is the same as in Eqs. (4) and (8), and the average value $`\overline{\mathrm{sin}^2\theta _0(\tau )}`$ is given by Eq. (18).

Information about the distribution of the phase $`\zeta `$ can be obtained from the stationary version of Eq. (23),

$$D\frac{d^2P}{d\zeta ^2}=(1b\mathrm{cos}\zeta )\frac{dP}{d\zeta }b\mathrm{sin}\zeta P,$$
(25)
where the final set of notation is $`bF_1/F_0`$ and $`D\mathrm{\Xi }/F_0`$. These two parameters are interpreted , respectively, as the ratio of the drive’s amplitude to the friction coefficient, and an FP diffusion coefficient, which is proportional to an effective drive’s linewidth.

In the absence of the diffusion ($`D=0`$), a solution to Eq. (25) is

$$P(\zeta )=\left(2\pi \right)^1\sqrt{1b^2}/(1b\mathrm{cos}\zeta ),$$
(26)
where the normalization $`_\pi ^{+\pi }P(\zeta )𝑑\zeta =1`$ is imposed. The solution (26) is regular at $`b<1`$, while the singularity at $`b=1`$ exactly corresponds to the drive’s amplitude attaining the threshold value (15), i.e., to the onset of the phase-locking regime.

Collecting results produced by the numerical solution of Eq. (25), we have concluded that it is possible to define a critical value $`D_{\mathrm{c}r}`$ of $`D`$, at which a sharp transition from the phase-locked state at $`D<D_{\mathrm{c}r}`$ to an unlocked one at $`D>D_{\mathrm{c}r}`$ takes place. The transition is still better illustrated by consideration of the probability flux, $`J`$ $`\left[DP^{}+(1b\mathrm{cos}\zeta )P\right]`$, in terms of which the time-dependent FP (Smoluchowski) equation (23) is written as $`F_0^1P_t+`$ $`J/\zeta =0`$. At the points of a minimum of the stationary distribution function, $`|J|`$ gives the phase-slippage rate, i.e., a rate of the transition from a vicinity of a phase-locked point to a point differing by a phase shift $`\pm 2\pi `$. In Fig. 1, we display the phase-slippage rate vs. $`D`$ at different constant values of $`b`$, as found from the numerical solution of Eq. (25). The existence of critical values $`D_{\mathrm{c}r}`$, such that virtually no phase slippage takes place at $`D<D_{\mathrm{c}r}`$, is evident. Note that for $`b<1`$, when the locking is impossible even to the monochromatic drive in the absence of the noise ($`D=0`$), $`D_{\mathrm{c}r}`$ does not exist. For $`b=1`$, i.e., exactly at the locking threshold (15 ), $`D_{\mathrm{c}r}=0`$, which cannot be seen on the logarithmic scale used in Fig. 1. It is necessary to mention that a picture which may be interpreted as showing the flux $`J`$ as a function of $`b`$ at different fixed values of $`D`$ is available in the book ; nevertheless, we find it relevant to present Fig. 1 here, as we need to display $`J(D)`$ at different fixed values of $`b`$ (otherwise the existence of the critical values $`D_{\mathrm{c}r}`$ is not obvious).

The sharp unlocking transition can also be seen in terms of the ratio of the aforementioned minimum value of $`P(\zeta )`$ to its maximum value, which is attained fairly close to the unperturbed (i.e., pertaining to the monochromatic drive) locking point. These data (not displayed here) show that the ratio is virtually equal to zero at $`D<D_{\mathrm{c}r}`$, and abruptly begins to increase exactly at $`D=D_{\mathrm{c}r}`$. In order to quantify $`D_{\mathrm{c}r}`$ we define it as the value of $`D`$ for which $`J=10^3`$.

Figure 2 shows the most important characteristic of the unlocking transition, viz., $`D_{\mathrm{c}r}`$ vs. the effective drive’s amplitude $`b`$. The dependence is nearly linear, except for the region $`0<b11`$, i.e., just above the locking threshold (15). It is difficult to show the dependence in this region directly, therefore, in Fig. 2 we instead display $`r(b)`$ for $`b11`$, where $`r`$ is defined so that $`D_{\mathrm{c}r}`$ is approximated by an expression $`C_1(b1)^r`$ with a suitable constant $`C_1`$.

An asymptotic value of $`r`$ for $`b10`$ can be found analytically, as one can use the exact solution (26) to describe an approximate form of the distribution function, except in a sensitive region of small $`\zeta `$: $`P(\zeta )1/\mathrm{sin}^2(\zeta /2)`$. At small $`\zeta `$, one can expand Eq. (25), taking into regard that $`b1`$ and $`D`$ are now small too. This yields an equation,

$$\frac{d^2P}{d\eta ^2}=\left(\frac{1}{2}\eta ^2\frac{b1}{D^{2/3}}\right)\frac{dP}{d\eta }\eta P,$$
(27)
where $`\eta D^{1/3}\zeta `$. Thus, the solution in the sensitive region (which is $`\eta 1`$, or $`\zeta \sqrt{b1}D^{1/3}`$) depends on the single parameter $`(1b)/D^{2/3}`$. Although a solution to Eq. (27) can be matched to the aforementioned approximation $`P(\zeta )1/\mathrm{sin}^2(\zeta /2)`$, valid at larger $`\zeta `$, only numerically, it is obvious that, as $`b10`$, the dependence $`D_{\mathrm{c}r}(b)`$ must be $`D_{\mathrm{c}r}=C_1(b1)^{3/2}`$, with $`C_10.17`$ found from numerical data. The value $`r=3/2`$, obtained for $`b10`$, is in good agreement with the numerical results displayed in Fig. 2.

The linearity of the dependence $`D_{\mathrm{c}r}(b)`$ at large $`b`$ can be explained in a very simple way: neglecting in this case $`1`$ in the expression $`(1b\mathrm{cos}\zeta )`$ in Eq. (25), we immediately conclude that the asymptotic solution depends on the single parameter $`b/D`$, hence the dependence must take the form $`D_{\mathrm{c}r}=C_2b`$, with a constant $`C_20.14`$ found numerically. The latter result can also be interpreted in another way: the minimum (threshold) value of the ac-drive’s amplitude necessary to support the rotation of the pendulum grows nearly linearly with the linewidth $`\delta \omega `$, so that it may be approximated by $`ϵ_{\mathrm{t}hr}\left(\delta \omega \right)ϵ_{\mathrm{t}hr}^{(0)}(1+\mathrm{c}onst\delta \omega `$), where $`ϵ_{\mathrm{t}hr}^{(0)}`$ is given by Eq. (15) and the constant is roughly $`1/C_2`$.

These analytical results, which comply well with the numerical findings, justify the introduction of the very concept of the critical value $`D_{\mathrm{c}r}`$ of the phase-diffusion constant in the Smoluchowski equation, which was originally defined above in a phenomenological way, just by looking at Fig. 1.

## 5 Numerical results

To check the limits of validity of the analytical results obtained above, we have performed numerical simulations of the full stochastic equation (2), using a simple Euler scheme (the use of this scheme is considered in the book ). As usual, we halved the time step until the results would converge to a stable value within few-per-cent accuracy. In this section, we will first focus on the case of the parametric noise, so we now set $`j(t)=0`$, keeping the random term $`\psi (t)`$ in Eq. (2). The basic phenomenon sought in the previous analysis was the escape from the state synchronized with the external drive as per Eq. (1). However, since the JJ is driven by the ac-term alone, once the system is no longer phase-locked, it cannot sustain progressive motion and will therefore quickly decay to the zero-voltage state. So the prediction of Eq.(20 ) actually refers to the lifetime of the phase-locked state; an abrupt transition from this state to the zero-voltage one (after about $`300`$ time unit) is evident in Fig. 3, which displays the time dependence of the phase velocity, found in a typical run of the simulations of the stochastic equation (2).

To estimate the lifetime of the phase-locked state at a given “temperature” (spectral linewidth of the drive), we have run the simulations for many different realizations of the random phase $`\psi (t)`$, and averaged the results for the lifetime. The number of the realizations was determined by the condition that the computed average has to converge to an established value. In Fig. 4, the logarithm of the thus computed average lifetime is plotted versus the inverse linewidth, so that the potential barrier (19) could be estimated from the slope of the linear part of this dependence.

This method closely follows that of Ref. ; the main difference is, as already mentioned, that we are not looking for mere phase-slippage, but for a jump to the state with zero average velocity. As it is clearly seen in Fig. 3, this occurs at a somewhat later time than the phase slips commence, although the difference is, typically, small (for instance, it is less than $`20`$ time units in the example shown in Fig. 3).

The numerically computed effective energy barrier (represented by lines with symbols) and the one predicted by Eq. (19) (the lines without symbols) are shown in Fig. 5(a). Taking into regard that no fitting parameter was employed, the agreement is very good near the threshold value (15), which is $`ϵ_{\mathrm{t}hr}0.16`$ for the values of the parameters corresponding to Fig. 5(a). However, a drastic deviation from the analytical prediction is evident at larger values of the ac-drive’s amplitude: while the analytical formula (20) predicts an almost linear increase of the effective barrier height with the ac drive amplitude, the numerical results show a maximum followed by a substantial decrease of barrier’s height.

It is relevant to mention that the nearly linear increase of the barrier height with $`ϵ`$ was predicted by the power-balance approach, which was employed above for the analytical consideration (see also Ref. ), a different method, based on an harmonic expansion, would result in a nonmonotonic dependence of the barrier height on $`ϵ`$ . In the case of the additive noise, the latter method produced the energy barrier demonstrating a Bessel-functional behavior, typical of the rf-induced current steps in JJs . In our case, however, such a dependence cannot be analytically justified. Moreover, even if the results of Ref.  are formally applied to our case, yielding

$$\mathrm{\Delta }UJ_1\left(ϵ/\omega _0^2\right),$$
(28)
where $`J_1`$ is the Bessel function, the maximum of $`\mathrm{\Delta }U`$ would occur at a much larger value of $`b`$ \[$`b45`$, instead of $`4`$ in Fig. 5(a) for $`\alpha =0.01`$\]. Thus, the phenomenon reported here is an essentially new one, and still remains to be explained.

It should, moreover, be noticed that the Bessel-function approximation similar to Eq. (28) gives, according to Ref. , a good estimate for the energy barrier only in the limit $`\omega _0^2<<1`$, deviations occurring already for $`\omega _0^20.05`$. Therefore, it is not surprising that, for the parameters considered here ($`\omega _0^2=0.25`$), the agreement with the Bessel-function behavior is very poor.

To check that the newly found dependence of the effective barrier height vs. the normalized drive’s amplitude is not due to some particular feature of the parametric noise, we have also performed extensive simulations of the same stochastic equation (2), but with the additive noise and strictly monochromatic driving signal, for the same values of parameters as those used above in the case of the frequency fluctuations. A typical example is shown in Fig. 5(b), together with the theoretical estimate of the energy barrier according to Ref. :

$$\mathrm{\Delta }U=2J\left(b\frac{ϵ_{\mathrm{t}hr}}{\omega _0^2}\right)\left[\sqrt{1b^2}b^1\mathrm{cos}^1\left(b^1\right)\right]$$
(29)
In this case too, a strong deviation of $`\mathrm{\Delta }U(b)`$ from the linear dependence occurs at relatively low values of the drive amplitude ($`b7`$), although they are higher than those in the frequency-noise case \[which are $`b4`$, see Eq. 5(a)\].

## 6 Conclusions

Results reported in this work are of relevance for applications to systems in which Josephson junctions are phase-locked to an external ac source, such as voltage standards. By substituting reasonable experimental values into Eq. (20), we can estimate that an ac source with the relative linewidth better than $`10^4`$ is needed if the lifetime of the order of $`1`$ second is required for the measurement system.

We also note that our approach could be applied to another problem: a pendulum driven by an ac signal whose frequency is subject to a systematic (rather than random) change, i.e., a zero-linewidth but variable-frequency drive. A system of the latter type was considered in Ref.  for a soliton in a perturbed nonlinear Schrödinger equation.

This work was performed in the framework of the bilateral-cooperation agreement between Consiglio Nazionale di Ricerca (Italy) and the Israeli Ministry of Science and Technology, under the project Nonlinear Dynamics of Josephson Networks. The authors appreciate useful discussions with G. Costabile and access to computer facilities at the University of l’Aquila (Italy).

## Figure Captions

Fig. 1. The phase-slippage rate $`J`$ vs. the diffusion parameter $`D`$ at different fixed values of the normalized drive’s amplitude $`b`$, which are indicated near each curve.

Fig. 2. The normalized noise threshold $`D_{cr}`$ vs. the normalized drive’s amplitude $`b`$ (solid line). The dashed line shows the dependence of the exponent $`r`$ defined in the text.

Fig. 3. A typical example of the evolution of the phase velocity $`\dot{\theta }(t)`$ obtained from the numerical integration of the stochastic equation (2). A loss of phase-locking occurs at $`t300`$. The parameters are $`\alpha =0.01`$, $`ϵ=0.2`$, $`\mathrm{\Omega }=0.01`$, $`\omega _0=2`$.

Fig. 4. Lifetime of the phase-locked state, plotted on a logarithmic scale versus the inverse linewidth. The parameters are $`\alpha =0.01`$, $`\omega _0=2`$. The estimated slope is $`0.00031`$ for $`ϵ=0.2`$, and $`0.0012`$ for $`ϵ=0.5`$.

Fig. 5. Dependence of the normalized energy barrier $`\mathrm{\Delta }U`$ on the normalized drive’s amplitude $`b`$, for the case of frequency fluctuations (a) and additive noise (b). Symbols represent the energy barrier estimated numerically on the basis of Eq. (3) for $`\omega _0=2`$ and $`\alpha =0.1`$ (squares) and $`\alpha =0.01`$ (triangles), connecting lines being a guide for the eye. The curves without symbols represent analytical predictions, viz., Eq. (20) for the frequency-fluctuation case \[in the panel (a), the continuous and dashed curves pertain, respectively, to $`\alpha =0.1`$ and $`\alpha =0.01`$\], and Eq. (29) for the additive-noise case.


text:
# Pair of Heavy-Exotic-Quarks at LHC

## I Introduction

One of the major problems to be studied in particle physics concerns the spectrum of elementary fermions. Many models consider the possible existence of new generations of fermions, such as composite models , grand unified theories , technicolor models , superstring-inspired models , mirror fermions , etc., which predict the existence of new particles with masses in turn of the scale of $`1`$ TeV. In this work, we will study the production mechanism for exotic-quarks at LHC. It is assumed here that the mixing between ordinary and exotic leptons are of the same flavour.

The three models that we consider here include new fermionic degrees of freedom, which introduce naturally a number of unknown mixing angles and fermionic masses . These models are: the vector singlet model (VSM) , which consists in the inclusion of new left- and right-handed fermions in singlets, the vector doublet model (VDM) , that can arise at low energies in the $`\mathrm{𝟐𝟕}`$ representation of $`E_6`$ theories and the fermion-mirror-fermion (FMF) model , where the new particles are introduced to restore the right-left symmetry.

Exotic fermions mixed with the standard fermions interact through the standard weak vector bosons $`W^+,W^{}`$ and $`Z^0`$, according to the Lagrangians

$$_{\mathrm{NC}}=\frac{g}{4cos\theta _W}\left[\overline{F}_i\gamma ^\mu (g_V^{ij}g_A^{ij}\gamma ^5)F_j+\overline{F}_i\gamma ^\mu (g_V^{ij}g_A^{ij}\gamma ^5)f_j\right]Z_\mu $$
(1)
and

$$_{\mathrm{CC}}=\frac{g}{2\sqrt{2}}\overline{Q_i}\gamma ^\mu (C_V^{ij}C_A^{ij}\gamma ^5)Q_jW_\mu ,$$
(2)
where $`F`$ and $`f`$ are the exotic and standard fermions, $`Q`$ are the exotic quarks, $`g_V^{ij}`$ and $`g_A^{ij}`$ are the corresponding neutral vector-axial coupling constants, and $`C_V^{ij}`$ and $`C_A^{ij}`$ are the charged vector-axial coupling constants, which can be obtained from Eq. (2) and are given in Table I, for each of the three models that we study here.

We consider here that all mixing angles have the value $`\theta _i=0.1`$, although phenomenological analysis  give an upper bound of $`sin^2\theta _i0.03`$. This means that the value of $`\theta _i`$ can be scaled up to $`0.173`$.

TABLE I. Coupling constants for a charged heavy fermion interaction: for the vector singlet model (VSM), the vector doublet model (VDM) and the fermion-mirror-fermion (FMF) model :

| Cou. | VSM | VDM | FMF |
| --- | --- | --- | --- |
| $`C_V^{QQ^{^{}}}`$ | $`\mathrm{sin}\theta _{iL}\mathrm{sin}\theta _{jL}`$ | $`\mathrm{cos}(\theta _{iL}\theta _{jL})+\mathrm{cos}\theta _{iR}\mathrm{cos}\theta _{jR}`$ | $`\mathrm{sin}\theta _{iL}\mathrm{sin}\theta _{jL}+\mathrm{cos}\theta _{iR}\mathrm{cos}\theta _{jR}`$ |
| $`C_A^{QQ^{^{}}}`$ | $`\mathrm{sin}\theta _{iL}\mathrm{sin}\theta _{jL}`$ | $`\mathrm{cos}(\theta _{iL}\theta _{jL})\mathrm{cos}\theta _{iR}\mathrm{cos}\theta _{jR}`$ | $`\mathrm{sin}\theta _{iL}\mathrm{sin}\theta _{jL}\mathrm{cos}\theta _{iR}\mathrm{cos}\theta _{jR}`$ |

The other coupling constants for the neutral heavy fermion interaction are given in Ref. .

## II Cross Section Production and Results

The production mechanism can be studied through the analysis of the reaction $`\overline{q}q(gg)\overline{Q}Q`$, provided that there is enough available energy ($`\sqrt{s}2M_Q`$). Such process takes place through the exchange of a photon, a $`Z^0`$ and a $`W`$ in the $`s`$ channel, which concerns electroweak interactions, and a quark and a gluon in the $`s`$ and $`t`$ channel, corresponding to the strong interaction.

For the $`pp`$ colliders, the contribution of the $`Z`$ and $`W`$ exchange, despite the high energy, is expected to be of the same order as the photon exchange, since the couplings of exotic quarks to $`Z`$ and $`W`$ have a weak interaction strength but the mixing angles are not so weak, while for the photon we have a electromagnetic strength only.

Using the interaction Lagrangians, Eqs. ($`1`$) and ($`2`$), we first evaluate the cross section for the electroweak case involving a neutral current and obtain:

$`\left({\displaystyle \frac{d\widehat{\sigma }}{d\mathrm{cos}\theta }}\right)_{\overline{Q}Q}=`$ $`{\displaystyle \frac{\beta _Q\alpha ^2\pi }{s}}[{\displaystyle \frac{c_q^4}{s^2}}[2sM_Q^2+(M_Q^2t)^2+(M_Q^2u)^2]`$ (8)
$`+{\displaystyle \frac{1}{16\mathrm{sin}^4\theta _W\mathrm{cos}^4\theta _W(sM_Z^2)^2}}[2sM_Q^2(g_{V}^{QQ}{}_{}{}^{2}g_{A}^{QQ}{}_{}{}^{2})(g_{V}^{q}{}_{}{}^{2}+g_{A}^{q}{}_{}{}^{2})+`$
$`(g_{V}^{QQ}{}_{}{}^{2}+g_{A}^{QQ}{}_{}{}^{2})(g_{V}^{q}{}_{}{}^{2}+g_{A}^{q}{}_{}{}^{2})((M_Q^2t)^2+(M_Q^2u)^2)`$
$`+4(g_V^{QQ}g_A^{QQ}g_V^qg_A^q((M_Q^2u)^2(M_Q^2t)^2)]`$
$`+{\displaystyle \frac{c_q^2}{2\mathrm{sin}^2\theta _W\mathrm{cos}^2\theta _Ws(sM_Z^2)}}[2sM_Q^2(g_V^qg_V^{QQ}+(g_V^qg_V^{QQ}`$
$`g_A^qg_A^{QQ})(M_Q^2t)^2+(g_V^qg_V^{QQ}+g_A^qg_A^{QQ})(M_Q^2u)^2]],`$
where $`\beta _Q=\sqrt{14M_Q^2/s}`$ is the velocity of the exotic-quark in the c.m. of the process, $`c_q`$ is the charge of the quark, $`Q`$ is the exotic quark and $`\overline{Q}`$ the exotic antiquark, $`M_Z`$ is the mass of the $`Z`$ boson, $`\sqrt{s}`$ is the center of mass energy of the $`q\overline{q}`$ system, $`t=M_L^2\frac{s}{2}(1\beta \mathrm{cos}\theta )`$ and $`u=M_L^2\frac{s}{2}(1+\beta \mathrm{cos}\theta )`$, where $`\theta `$ is the angle between the exotic quark and the incident quark, in the c.m. frame and the couplings $`g_V^{QQ}`$ and $`g_A^{QQ}`$ are given in .

For the charged current we obtain:

$`\left({\displaystyle \frac{d\widehat{\sigma }}{d\mathrm{cos}\theta }}\right)_{\overline{Q}Q}^W=`$ $`{\displaystyle \frac{\beta _Q\alpha ^2\pi }{32\mathrm{sin}^4\theta _Ws(sM_W^2)^2}}[2sM_Q^2(C_{V}^{QQ^{^{}}}{}_{}{}^{2}C_{A}^{QQ^{^{}}}{}_{}{}^{2})+`$ (10)
$`(C_{V}^{QQ^{^{}}}{}_{}{}^{2}C_{A}^{QQ^{^{}}}{}_{}{}^{2})(M_Q^2t)^2+(C_{V}^{QQ^{^{}}}{}_{}{}^{2}+C_{A}^{QQ^{^{}}}{}_{}{}^{2})(M_Q^2u)^2],`$
where $`M_W`$ is the mass of the W boson and $`C_V^{QQ}`$ and $`C_A^{QQ}`$ are given in Table I.

The total cross section for the process $`ppqq\overline{Q}Q`$ is related to the subprocess $`qq\overline{Q}Q`$ total cross section $`\widehat{\sigma }`$, by

$$\sigma =_{\tau _{min}}^1_{\mathrm{ln}\sqrt{\tau _{min}}}^{\mathrm{ln}\sqrt{\tau _{min}}}𝑑\tau 𝑑yq(\sqrt{\tau }e^y,Q^2)q(\sqrt{\tau }e^y,Q^2)\widehat{\sigma }(\tau ,s),$$
(11)
where $`\tau =\frac{\widehat{s}}{s}(\tau _{min}=\frac{4M_Q^2}{s})`$, with $`s`$ being the center-of mass energy of the $`pp`$ system, and $`q(x,Q^2)`$ is the quark structure function.

The strong production of heavy standard quarks is very well studied and those results can be used to the exotic quarks. With respect to the gluon fusion there is an additional enhanced process to be considered. This process is regulated by the axial part of the Z boson and a Higgs. The exchange of a photon is not allowed by C conservation (Furry’s theorem). The contribution of the Z boson is given by

$`\left({\displaystyle \frac{d\widehat{\sigma }}{d\mathrm{cos}\theta }}\right)_{\overline{Q}Q}^Z={\displaystyle \frac{g_{A}^{Q}{}_{}{}^{2}N_c\alpha ^2\alpha _s^2}{256\pi \mathrm{sin}^4\theta _W}}{\displaystyle \frac{M_Q^2}{M_W^4}}\beta _Q\left|{\displaystyle \underset{q=u,d}{}}T_3^q\left(1+2\lambda _qI_q\right)\right|^2`$

and of a Higgs by

$`\left({\displaystyle \frac{d\widehat{\sigma }}{d\mathrm{cos}\theta }}\right)_{\overline{Q}Q}^H={\displaystyle \frac{N_c\alpha ^2\alpha _s^2}{1024\pi \mathrm{sin}^4\theta _W}}{\displaystyle \frac{M_Q^2}{M_W^4}}\widehat{s}^2\beta _Q^3\left|\chi (\widehat{s}){\displaystyle \underset{q=u,d}{}}\left(2\lambda _q+\lambda _q(4\lambda _q1)I_q\right)\right|^2,`$where the summations run over all generations. $`T_3^q`$ is the quark weak isospin \[$`T_3^{u(d)}=+()1/2`$\]. The loop function $`I_iI(\lambda _i=m_i^2/\widehat{s})`$, is defined by

$`I_iI_i(\lambda _i)={\displaystyle _0^1}{\displaystyle \frac{dx}{x}}\mathrm{ln}\left[1{\displaystyle \frac{(1x)x}{\lambda _i}}\right]=\{\begin{array}{c}2\left[\mathrm{sin}^1\left(\frac{1}{2\sqrt{\lambda _i}}\right)\right]^2,\lambda _i>\frac{1}{4}\hfill \\ \frac{1}{2}\mathrm{ln}^2\left(\frac{r_+}{r_{}}\right)\frac{\pi ^2}{2}+i\pi \mathrm{ln}\left(\frac{r_+}{r_{}}\right),\lambda _i<\frac{1}{4},\hfill \end{array}`$
with, $`r_\pm =1\pm (14\lambda _i)^{1/2}`$ and $`\lambda _i=m_i^2/\widehat{s}`$. Here, $`i=q`$ stands for the particle (quark ) running in the loop.

We have also defined

$`\chi (\widehat{s})={\displaystyle \frac{1}{\widehat{s}m_H^2+iM_H\mathrm{\Gamma }_H}}`$with $`\mathrm{\Gamma }_H`$ being the width of the Higgs boson.

The total cross section for the process $`ppgg\overline{Q}Q`$ is related to the subprocess $`gg\overline{Q}Q`$ total cross section $`\widehat{\sigma }`$ by

$$\sigma =_{\tau _{min}}^1_{\mathrm{ln}\sqrt{\tau _{min}}}^{\mathrm{ln}\sqrt{\tau _{min}}}𝑑\tau 𝑑yG(\sqrt{\tau }e^y,Q^2)G(\sqrt{\tau }e^y,Q^2)\widehat{\sigma }(\tau ,s)$$
(13)
where $`\tau =\frac{\widehat{s}}{s}(\tau _{min}=\frac{4M_{\stackrel{~}{Q}}^2}{s})`$, with $`s`$ being the center of mass energy of the $`pp`$ system, and $`G(x,Q^2)`$ is the gluon structure function.

In Fig. $`1`$ we show the cross sections for the production of exotic U-quarks, (for the D-quarks, the results for the cross section are similar), $`\overline{p}p\gamma ,Z\overline{U}U`$, we see that for the three models that we are considering the cross section are similar one to another. In all numerical calculations we take $`\mathrm{sin}^2\theta _W=0.2315`$, $`M_Z=91.118`$ GeV, $`M_W=80.33`$ GeV and used the distribution functions given by M. Glück, E. Reya and A. Vogt . Considering that the expected integrated luminosity for the LHC will be of the order of $`10^5pb^1/yr`$, and taking the mass of the quark-U equal to $`500`$ GeV, we expect to have a total of $`10^3`$ exotic quark pairs produced per year.

Fig. $`2`$ shows the cross section for the production of exotic-U quarks, $`\overline{p}pW\overline{U}U`$, for the VDM and for the FMF model. For the quark mass of 500 GeV we expect a total of $`10^4`$ exotic quarks pairs produced per year for the VDM, while for the FMF we estimate a total of $`3.10^3`$ events per year.

In Fig. $`3`$ we present the cross sections for the gluon fusion mechanism $`\overline{p}pgg\overline{U}U`$. For the same mass of 500 GeV we obtain, for the VDM model $`10^4`$ quark pairs produced per year and $`3.10^4`$ for the FMF model. In this case we take the mass of the Higgs equal to $`100`$ GeV.

In Fig. $`4`$ we compare the cross sections as for the electroweak case and as for the strong case. We observe that for large masses of exotic quarks, (more than $`500`$ GeV), the cross sections for the electroweak case begins to be competitive with the strong case, that is, for every $`100`$ events produced by the strong case we have $`15`$ to $`20`$ events produced by the electroweak case.

We make use of $`p_T`$ distributions to try to distinguish the signatures of the strong and electroweak cases. Fig. $`5`$ shows the distribution of $`p_T`$ for the exotic quarks. We see that the exotic quarks signal for the electroweak case is below for the QCD-case. We take in this case for the mass of the exotic quarks equal to $`600`$ GeV.

Fig. $`6`$ shows the distribution of transverse mass $`M_T(qe,\nu )`$. We see from the figure a Jacobian peak near the $`M_T=M_Q`$. The transverse mass is defined by

$`M_T^2(qe,\nu )=\left[(p_{qeT}^2+M_{qe}^2)^{1/2}+p_{\nu T}\right]^2(\stackrel{}{p}_{qeT}+\stackrel{}{p}_{\nu T})^2`$where from here we know some information about the recoil q-jet. The $`Q\overline{Q}`$ events with real W decays are distinguished from the electroweak single and pair production of W’s by the form of their decays, while for the heavy quark we have an odd number of quarks for the W we have an even number respectively.

It is necessary to point out that the process involving a boson Z to gluon-gluon fusion for VSM and VDM-quarks gives a negligible contribution compared to the FMF-quarks, due to the small mixing angles.

###### Acknowledgements.

I would like to thank to Prof. I. P. Calvalcante for a careful reading of the manuscript and Prof. J. O. P. Éboli for some useful suggestions. One of the authors (P.P.Q.F.) would like to thank the Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro for full financial support (contract No. E-26/150.639/99).

FIGURE CAPTIONS

Figure 1: Total cross section for the process $`\overline{p}p\gamma ,Z\overline{U}U`$, as a function of $`M_U`$ at $`\sqrt{s}=14`$ TeV: (a) vector singlet model (dotted line); (b) vector doublet model (dashed line); (c) fermion mirror fermion (solid line).

Figure 2: Total cross section for the process $`\overline{p}pW\overline{U}U`$ as a function of $`M_U`$ at $`\sqrt{s}=14`$ TeV: (a) vector doublet model (dashed line); (b) fermion mirror fermion (dot dashed line).

Figure 3: Total cross section for the gluon-gluon fusion mechanism, $`\overline{p}pgg\overline{U}U`$, as a function of $`M_U`$ at $`\sqrt{s}=14`$ TeV: (a) vector doublet model (dashed line); (b) fermion mirror fermion (dot dashed line).

Figure 4: Total cross section for the process $`\overline{p}p\overline{U}U`$ as a function of $`M_U`$ at $`\sqrt{s}=14`$ TeV: (a) for the strong case (dot dashed line); for the electroweak case: (b) FMF (doted line) and (c) VDM (dashed line).

Figure 5: $`p_T`$ distributions for $`\overline{Q}Q`$ production: for exotic quarks mass equal to $`600`$ GeV.

Figure 6: Normalized distribution of transverse mass $`M_T(qe,\nu )`$ for $`\overline{Q}Q`$ production: (a) for exotic quark mass equal to $`300`$ GeV (solid line); (b) for exotic quark mass equal to $`500`$ GeV (dashed line).