format: text

text:
# Transverse phase-locking in fully frustrated Josephson junction arrays: a new type of fractional giant steps

## I INTRODUCTION

Phase locking phenomena are found in a wide variety of nonlinear driven systems in condensed matter physics.bak  It takes place when an internal frequency of the system locks to a rational multiple of the frequency of an external ac-drive. A simple example of this is the case of an overdamped particle moving in a tilted washboard potential, where the frequency of motion of the particle over the periodic potential can be locked to multiples of the frequency of a superimposed ac force for a finite range of the dc force (tilt of the washboard). Since the internal (washboard) frequency is proportional to the mean velocity of the particle, phase locking results in a constant mean velocity for a certain range of dc-force curve when the interference condition is satisfied. A particularly well known realization of this effect is Shapiro stepsshapiro  in the dc current-voltage ($`IV`$) characteristics of a single small area Josephson junction driven by a time periodic current. Within the washboard analogy outlined above, a simple analysis provides expressions for the appearance of Shapiro steps at specific voltages barone  corresponding to integer multiples of the driving frequency.

Driven systems with many degrees of freedom can also exhibit phase locking. This has attracted broad scientific and technological interest since phase-locking in complex systems can either be induced by collective effects, providing for a low dimensional interpretation of the phenomenon, or itself induce collective (low dimensional) behavior in the complex system. Phase-locking experiments have provided information about such dynamical response of non-equilibrium collective states, where dimensionality, thermal fluctuations, quenched disorder, and the magnitude of external fields can be very relevant. A particular well known example is the large Josephson junction array (JJA), with $`N\times N`$ junctions, driven by an external current $`(I_{dc}+I_{ac}\mathrm{cos}\mathrm{\Omega }t)\widehat{x}`$ with frequency $`\mathrm{\Omega }`$ and with an applied magnetic field density $`f=Ha^2/\mathrm{\Phi }_0`$, where $`H`$ is the magnetic field, $`a`$ the lattice period of the Josephson array and $`\mathrm{\Phi }_0`$ the quantum of flux. Giant Shapiro steps at voltages $`V_n=Nn\mathrm{}\mathrm{\Omega }/2e`$, $`n`$ being an integer, have been observed experimentally in zero magnetic field ($`f=0`$). leeman  Fractional giant Shapiro steps at voltages $`V_{n,q}=Nn\mathrm{}\mathrm{\Omega }/2eq`$, were observed experimentally expfgs ; expshap  for strongly commensurate magnetic fields, $`f=p/q`$, with $`p,q`$ integers, and extensively investigated in numerical simulations.lee ; simff ; fourp ; octav ; ciria  Also, subharmonic giant Shapiro steps at voltages $`V_{n,m}=Nn\mathrm{}\mathrm{\Omega }/2em`$ were observed experimentally expshap  for zero magnetic field ($`f=0`$), and attributed to the nucleation of complex collective dynamical states acvs ; giov  induced by disorder or inductance effects. Shapiro-like phase-locking is also observed in the case of driven vortex lattices in bulk superconductors with two-dimensional periodic pinning arrays, as recently reported both experimentally vanlook  and theoretically.reichhardt  Also superconductors, where vortices are driven over a one-dimensional potential generated by thickness modulations,martinoli  or are confined to move through mesoscopic channels, kokubo  show Shapiro-like phase-locking. Moreover, systems with many degrees of freedom in the presence of quenched disorder can also exhibit phase-locking when there is a dynamically induced periodicity, like charge density waves gruner  and vortex lattices in superconductors with random pinning.fiory ; harris ; kolton

In the phase-locking examples mentioned above, the ac-drive is applied parallel to the dc-drive. However, it was recently shown that a different type of phase-locking, distinct from Shapiro phase-locking, is possible in vortex lattices if the ac-force is applied perpendicular to the dc-force.reichhardt2 ; kolton3  In this case the interference effect is due to an effective parametric ac drive in the longitudinal direction, which is induced by the transverse ac drive. In several systems, like charge density waves or single degree of freedom systems (e.g., the single Josephson junction), the dynamical variables are such that perturbations or displacements can be induced in only one direction (i.e., the displacement field is a scalar). An important characteristic of vortex lattices in superconductors is that the displacement field is two dimensional. In particular, the behavior of displacements in the direction perpendicular to the driving force shows phenomena like a transverse critical current koshelev ; moon ; olson ; marconi  and a transverse freezing transition koshelev ; kolton2  at high velocities. Phase locking in ac-driven vortex lattices, where the ac and dc forces are perpendicular, arises as a direct consequence of the nonlinear coupling between the two directions of motions.

Transverse phase-locking has been reported for vortex lattices moving in rectangular or triangular pinning arrays reichhardt2  as well as in arrays of randomly distributed pinning centers.kolton3  In this paper we investigate the possibility of transverse phase locking in a two dimensional (2D) fully frustrated JJA, where the average of the external magnetic field corresponds to one half flux quantum per plaquette, $`f=1/2`$. This system has several attractive properties. The presence of a magnetic field ($`f0`$) breaks the axial symmetry in the direction of the bias current, and 2D-cooperative behavior may come into play. This leads to the well-known fractional giant Shapiro steps expfgs ; expshap ; lee ; simff ; fourp ; octav ; ciria  induced by a longitudinal ac current. It also, as we will demonstrate in this paper, allows for transverse phase locking when the ac current is perpendicular to the dc current, since the two directions of motion become coupled. Non-equilibrium dynamical phases for fully frustrated JJAs driven by a dc current have previously been studied.niels ; marconi2  Phase locking can be used to characterize temporal order in the different dynamical phases of the JJA at high velocities by their ac-response, as was done in Ref. harris, ; kolton,  for bulk superconductors.

Here, we report transverse phase locking steps in the $`IV`$ characteristics, similar to the longitudinal giant Shapiro steps, but with very different characteristic dependencies on external ac-drive amplitude $`I_{ac}`$ and frequency $`\mathrm{\Omega }`$. The critical depinning current in the system with transverse ac drive is larger than the critical current of the dc driven system. For $`I_{ac}/\mathrm{\Omega }1`$ the depinning critical current $`I_c`$ and the phase-locked step width $`\mathrm{\Delta }S_1`$ for $`V=\mathrm{\Omega }\mathrm{}/2e`$ increase quadratically with $`I_{ac}`$. For $`I_{ac}/\mathrm{\Omega }>1`$ we find windows of $`I_{ac}/\mathrm{\Omega }`$ where depinning is hysteretic and the periodic phase-locked dynamical states become unstable. We characterize the dynamical states around the phase-locking interference condition in the $`IV`$ curve with the voltage noise, Lyapunov exponents and Poincaré sections. We find that zero temperature phase-locking behavior in large fully frustrated arrays is well described by an effective four plaquette model.

The remainder of this paper is outlined as follows. In section II we present the model used for simulating the dynamics of the fully frustrated JJA. In section III we develop an analytical framework for predicting critical current and phase-locking properties of the fully frustrated JJA. Section IV presents simulated transverse phase-locking steps in typical $`IV`$ curves obtained from an effective four plaquette model for the JJA. We calculate numerically the dependence of the critical current and the magnitude of phase-locking steps with the amplitude $`I_{ac}`$ and frequency $`\mathrm{\Omega }`$ of the external ac-drive. The results obtained are analyzed in more detail by studying voltage noise, Lyapunov exponents and Poincaré sections for the dynamical states around the phase-locking interference condition. We also compare numerical simulations results for large arrays with those obtained using the effective four plaquette model. The discussions and conclusions of the investigation are presented in sections IV and V.

## II MODEL

We study a current driven JJA with an ac current perpendicular to the dc current, as shown in figure 1a. A magnetic field, $`H`$, is applied such that half a flux quantum, $`f=Ha^2/\mathrm{\Phi }_0=1/2`$ with $`a^2`$ being the area of a plaquette and $`\mathrm{\Phi }_0=h/2e`$ being the flux quantum, penetrates each plaquette; corresponding to the fully frustrated XY model,xyff ; mon  where the ground state is a “checkerboard” vortex lattice, in which a vortex (flux quantum) penetrates every other square grid (see Fig. 1b). In such ground state, current and phase differences in the junctions are described by a repeated two-junction by two-junction ($`2\times 2`$ plaquette) superlattice unit cell.

Numerical simulations simff ; lee  of large driven arrays suggest that this spatial periodicity is preserved when the dynamics is phase-locked to an external ac-perturbation applied in parallel to the dc force. We will show later, in section IV, that this is also a good approximation when the ac current is perpendicular to the dc current. We will therefore consider the simple system of a $`2\times 2`$ superlattice unit cell of the array and the associated gauge-invariant phase differences in a field of $`f=1/2`$, with the dc current (per plaquette) $`I_{dc}`$ perpendicular to an ac current (per plaquette) $`I_{ac}\mathrm{cos}(\mathrm{\Omega }t)`$, as shown in Fig. 2. Flux quantization, total current conservation at the central node, and the applied total currents in the two directions give the following governing equations,

$`\alpha +\beta \delta \gamma `$ $`=`$ $`\pi (1+2n)`$ (1)
$`\dot{\beta }+\dot{\gamma }+\mathrm{sin}\beta +\mathrm{sin}\gamma `$ $`=`$ $`2I_{ac}\mathrm{cos}\mathrm{\Omega }t`$ (2)
$`\dot{\alpha }+\dot{\delta }+\mathrm{sin}\alpha +\mathrm{sin}\delta `$ $`=`$ $`2I_{dc}`$ (3)
$`\dot{\alpha }+\dot{\gamma }\dot{\beta }\dot{\delta }+\mathrm{sin}\alpha `$ $`+`$ $`\mathrm{sin}\gamma \mathrm{sin}\beta \mathrm{sin}\delta =\mathrm{\hspace{0.33em}0},`$ (4)
where $`n`$ is an integer, $`t`$ is the normalized time in units of $`t_0=2eI_0R_N/\mathrm{}`$, $`R_N`$ being the normal state single junction resistance, $`\mathrm{\Omega }`$ is the normalized frequency in units of $`1/t_0`$, $`I_{ac}`$ and $`I_{dc}`$ are the normalized external currents in units of the single junction critical current $`I_0`$. This model was introduced by Benz et al. ccff  (for $`I_{ac}=0`$) to study the dc current-voltage curve of the $`f=1/2`$ array. They obtained analytically that the critical current per junction of this model is $`I_c=(\sqrt{2}1)I_0`$. ccff  The same model was later used in Ref. fourp,  to study the (longitudinal) Shapiro steps. Since the analysis done in the work of Refs. ccff, ; fourp,  did not include the transverse ac current with accompanying transverse voltage drop, an additional constraint of $`\beta =\delta `$ was implied, reducing the model system to two dynamical degrees of freedom. In contrast, our model system of a four plaquette unit cell consists of three effective dynamical variables. We calculate the instantaneous longitudinal $`V_x`$ and transverse $`V_y`$ (normalized) voltages per junction as,

$`V_x`$ $`=`$ $`(\dot{\alpha }+\dot{\delta })/2`$ (5)
$`V_y`$ $`=`$ $`(\dot{\beta }+\dot{\gamma })/2,`$ (6)
and the $`IV`$ characteristics, $`v_x=V_x`$ as a function of $`I_{dc}`$, where $`\mathrm{}`$ is a time average. The total longitudinal voltage $`v_T`$ for an $`N\times N`$ array, built with this $`2\times 2`$ superlattice unit cell, is $`v_T=Nv_x`$. When vortices move with a mean velocity $`u`$ in such $`2a\times 2a`$ superlattice structure, we can obtain the normalized voltage $`v_x`$ using the relation $`2\pi u/2a=v_x`$, where $`a`$ is the array periodicity (see Fig. 1a). The intrinsic washboard frequency for vortices moving with velocity $`u`$ in the periodic potential of the JJA is $`\omega _0=2\pi u/a`$. Phase-locking in the longitudinal direction is obtained when the frequency $`\mathrm{\Omega }`$ of the ac drive locks to a rational multiple of the intrinsic washboard frequency. For the $`n`$-th harmonic this corresponds to $`\omega _0=n\mathrm{\Omega }`$, i.e. $`2\pi u/a=n\mathrm{\Omega }`$. This leads to phase-locking at voltages $`V_{n,2}=(n/2)\mathrm{\Omega }`$ for fully frustrated JJA. In general, for $`f=p/q`$, the ground state has $`qa\times qa`$ superlattice structure, therefore the voltage for vortices moving with velocity $`u`$ is $`2\pi u/qa=v_x`$, and phase-locking for the $`n`$-th harmonic is obtained at voltages $`V_{n,q}=(n/q)\mathrm{\Omega }`$. This is the condition for the so-called “fractional giant Shapiro steps” expfgs ; expshap ; lee ; simff ; fourp ; octav ; ciria .

## III Phase-Locking and Critical Current Analysis

We will in this section assume that the dynamics of the system is represented by the simple two-plaquette degrees of freedom as shown in figure 2. We will apply the following linear transformation of the phase-variables of Eqs. (1)-(4):

$`\mathrm{\Phi }_x`$ $`=`$ $`{\displaystyle \frac{\alpha +\delta }{2}}`$ (7)
$`\mathrm{\Psi }_x`$ $`=`$ $`{\displaystyle \frac{\alpha \delta }{2}}`$ (8)
$`\mathrm{\Phi }_y`$ $`=`$ $`{\displaystyle \frac{\beta +\gamma }{2}}`$ (9)
$`\mathrm{\Psi }_y`$ $`=`$ $`{\displaystyle \frac{\beta \gamma }{2}}.`$ (10)
With these variables we can represent the constraint ($`f=\frac{1}{2}`$) of Eq. (1) as $`\mathrm{\Psi }_x+\mathrm{\Psi }_y=\frac{\pi }{2}`$, and thereby write the relevant three degrees of freedom in either of the two following forms, eliminating $`\mathrm{\Psi }_x`$:

$`\dot{\mathrm{\Phi }}_x+\mathrm{sin}\mathrm{\Psi }_y\mathrm{sin}\mathrm{\Phi }_x`$ $`=`$ $`I_{dc}`$ (11)
$`\dot{\mathrm{\Phi }}_y+\mathrm{cos}\mathrm{\Psi }_y\mathrm{sin}\mathrm{\Phi }_y`$ $`=`$ $`I_{ac}\mathrm{cos}\mathrm{\Omega }t`$ (12)
$`2\dot{\mathrm{\Psi }}_y\mathrm{cos}\mathrm{\Phi }_x\mathrm{cos}\mathrm{\Psi }_y+\mathrm{cos}\mathrm{\Phi }_y\mathrm{sin}\mathrm{\Psi }_y`$ $`=`$ $`0,`$ (13)
or eliminating $`\mathrm{\Psi }_y`$:

$`\dot{\mathrm{\Phi }}_x+\mathrm{cos}\mathrm{\Psi }_x\mathrm{sin}\mathrm{\Phi }_x`$ $`=`$ $`I_{dc}`$ (14)
$`\dot{\mathrm{\Phi }}_y+\mathrm{sin}\mathrm{\Psi }_x\mathrm{sin}\mathrm{\Phi }_y`$ $`=`$ $`I_{ac}\mathrm{cos}\mathrm{\Omega }t`$ (15)
$`2\dot{\mathrm{\Psi }}_x+\mathrm{cos}\mathrm{\Phi }_x\mathrm{sin}\mathrm{\Psi }_x\mathrm{cos}\mathrm{\Phi }_y\mathrm{cos}\mathrm{\Psi }_x`$ $`=`$ $`0.`$ (16)
The normalized voltages are given by, $`V_x=\dot{\mathrm{\Phi }}_x`$ and $`V_y=\dot{\mathrm{\Phi }}_y`$. We will in the spirit of the usual Shapiro analysis assume that

$`\mathrm{\Phi }_y`$ $`=`$ $`{\displaystyle \frac{I_{ac}}{\mathrm{\Omega }}}\mathrm{sin}\mathrm{\Omega }t,`$ (17)
which is the solution to Eqs. (12) and (15) for large $`I_{ac}`$ and $`\mathrm{\Omega }`$.

### III.1 Critical Current

We will here look at Eqs. (14) and (16). Assuming $`\mathrm{\Phi }_x=\mathrm{\Phi }_x^{(0)}`$ and $`\mathrm{\Psi }_y=\mathrm{\Psi }_y^{(0)}+\epsilon (t)`$, where $`\mathrm{\Phi }_x^{(0)}`$ and $`\mathrm{\Psi }_y^{(0)}`$ are constants and $`\left|\epsilon (t)\right|1`$, the static contribution to equation (16) is:

$`\mathrm{cos}\mathrm{\Phi }_x^{(0)}\mathrm{sin}\mathrm{\Psi }_y^{(0)}`$ $`=`$ $`J_0\left({\displaystyle \frac{I_{ac}}{\mathrm{\Omega }}}\right)\mathrm{cos}\mathrm{\Psi }_y^{(0)}`$ (18)
$`\mathrm{tan}\mathrm{\Psi }_y^{(0)}`$ $`=`$ $`{\displaystyle \frac{J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}{\mathrm{cos}\mathrm{\Phi }_x^{(0)}}},`$ (19)
where $`J_n`$ is the $`n`$’th order Bessel function of the first kind. Inserting this into the static part of equation (14) yields

$`I_{dc}`$ $`=`$ $`\mathrm{cos}\left\{\mathrm{tan}^1\left[{\displaystyle \frac{J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}{\mathrm{cos}\mathrm{\Phi }_x^{(0)}}}\right]\right\}\mathrm{sin}\mathrm{\Phi }_x^{(0)}.`$ (20)
This expression provides a unique relationship between the constant phase, $`\mathrm{\Phi }_x^{(0)}`$, and the dc current, $`I_{dc}`$. However, for increasing $`I_{dc}`$, there exists a critical value, $`I_c^{}`$, for which no real $`\mathrm{\Phi }_x^{(0)}`$ can satisfy Eq. (20). This value is given by

$`I_c^{}`$ $`=`$ $`\mathrm{max}\left[\mathrm{cos}\left\{\mathrm{tan}^1\left[{\displaystyle \frac{J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}{\mathrm{cos}\mathrm{\Phi }_x^{(0)}}}\right]\right\}\mathrm{sin}\mathrm{\Phi }_x^{(0)}\right],`$ (21)
which is the predicted critical dc current for static states ($`v_x=\dot{\mathrm{\Phi }}_x=0`$). We notice that the identical expression for the critical current, Eq. (21), can be obtained from the anisotropic dc driven system,

$`\alpha +\beta \delta \gamma `$ $`=`$ $`\pi (1+2n)`$ (22)
$`\dot{\beta }+\dot{\gamma }+\mathrm{\Gamma }\mathrm{sin}\beta +\mathrm{\Gamma }\mathrm{sin}\gamma `$ $`=`$ $`0`$ (23)
$`\dot{\alpha }+\dot{\delta }+\mathrm{sin}\alpha +\mathrm{sin}\delta `$ $`=`$ $`2I_{dc}`$ (24)
$`\dot{\alpha }+\dot{\gamma }\dot{\beta }\dot{\delta }+\mathrm{sin}\alpha `$ $`+`$ $`\mathrm{\Gamma }\mathrm{sin}\gamma \mathrm{\Gamma }\mathrm{sin}\beta \mathrm{sin}\delta =\mathrm{\hspace{0.33em}0},`$ (25)
where the anisotropy, $`\mathrm{\Gamma }`$ (suppression of transverse critical current), is given by the standard Shapiro barone  critical current, $`\mathrm{\Gamma }=J_0(\frac{I_{ac}}{\mathrm{\Omega }})`$.

### III.2 Phase-Locking

We will here use equations (11)-(13), since $`\mathrm{\Psi }_y=0(\mathrm{mod}\pi )`$ for $`\dot{\mathrm{\Phi }}_x0`$ provides for a simple description of the dynamics. We will assume the ansatz, $`\mathrm{\Phi }_x=\mathrm{\Phi }_x^{(0)}+\mathrm{\Omega }t`$ and $`\mathrm{\Psi }_y=A\mathrm{sin}\mathrm{\Omega }t+B\mathrm{cos}\mathrm{\Omega }t`$. The equation for $`\mathrm{\Psi }_y`$, (13), now reads:

$`2\dot{\mathrm{\Psi }}_y\mathrm{cos}(\mathrm{\Phi }_x^{(0)}+\mathrm{\Omega }t)\mathrm{cos}\mathrm{\Psi }_y+\mathrm{sin}\mathrm{\Psi }_y{\displaystyle \underset{k}{}}J_k\left({\displaystyle \frac{I_{ac}}{\mathrm{\Omega }}}\right)\mathrm{cos}k\mathrm{\Omega }t`$ $`=`$ $`0,`$ (26)
where we will use the approximations: $`\mathrm{cos}\mathrm{\Psi }_yJ_0\left(\sqrt{A^2+B^2}\right)`$ and $`\mathrm{sin}\mathrm{\Psi }_y\mathrm{\Psi }_y`$. With the ansatz for $`\mathrm{\Psi }_y`$ above, this equation has no static component. The time varying component, at frequency $`\mathrm{\Omega }`$, yields the coefficients $`A`$ and $`B`$

$`A_1`$ $`=`$ $`A_0J_0\left(\sqrt{A_0^2+B_0^2}\right)={\displaystyle \frac{2\mathrm{\Omega }\mathrm{cos}\mathrm{\Phi }_x^{(0)}\left[J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)+J_2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)\right]\mathrm{sin}\mathrm{\Phi }_x^{(0)}}{4\mathrm{\Omega }^2+J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}}J_0\left(\sqrt{A_0^2+B_0^2}\right)`$ (27)
$`B_1`$ $`=`$ $`B_0J_0\left(\sqrt{A_0^2+B_0^2}\right)={\displaystyle \frac{2\mathrm{\Omega }\mathrm{sin}\mathrm{\Phi }_x^{(0)}+\left[J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)J_2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)\right]\mathrm{cos}\mathrm{\Phi }_x^{(0)}}{4\mathrm{\Omega }^2+J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}}J_0\left(\sqrt{A_0^2+B_0^2}\right),`$ (28)
where $`A_0`$ and $`B_0`$ are the solutions for $`J_0\left(\sqrt{A_0^2+B_0^2}\right)=1`$. Thus, the solution $`(A,B)=(A_0,B_0)`$ is correct up to order $`\mathrm{\Omega }^1`$ and $`(A,B)=(A_1,B_1)`$ is correct up to $`\mathrm{\Omega }^2`$. Inserting this solution $`(A_1,B_1)`$ for $`\mathrm{\Psi }_y`$ into the $`\mathrm{\Psi }_y`$-linearized equation (11) gives the static properties:

$`J_0\left(\sqrt{A_0^2+B_0^2}\right){\displaystyle \frac{1}{2}}\left(A_0\mathrm{cos}\mathrm{\Phi }_x^{(0)}+B_0\mathrm{sin}\mathrm{\Phi }_x^{(0)}\right)`$ $`=`$ $`I_{dc}\mathrm{\Omega }`$ (29)
$`J_0\left(\sqrt{A_0^2+B_0^2}\right){\displaystyle \frac{1}{2}}{\displaystyle \frac{2\mathrm{\Omega }J_2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)\mathrm{sin}2\mathrm{\Phi }_x^{(0)}}{4\mathrm{\Omega }^2+J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}}`$ $``$ $`\mathrm{\Delta }I_1{\displaystyle \frac{1}{2}}\mathrm{\Delta }S_1\mathrm{sin}2\mathrm{\Phi }_x^{(0)}`$ (30)
$`=`$ $`I_{dc}\mathrm{\Omega }.`$
Thus, the locking range for the this step can be found to second order in $`|\mathrm{\Psi }_y|`$. The dominant part of this expression for the range in phase-locking yields:

$`\mathrm{\Delta }S_1`$ $`=`$ $`{\displaystyle \frac{\left|J_2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)\right|}{4\mathrm{\Omega }^2+J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}}\left(1{\displaystyle \frac{3\mathrm{\Omega }^2+\frac{1}{4}\left(J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)+J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)\right)}{\left(4\mathrm{\Omega }^2+J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)\right)^2}}\right).`$ (31)
The expression displays quadratic growth of the phase-locked step size for small $`I_{ac}`$. This is consistent with the particle (pancake) model results reichhardt2  for vortices in rectangular pinning arrays, and it is different from the known longitudinal (Shapiro) phase-locking of JJAs.leeman ; expfgs ; expshap ; lee

In addition to the range of phase-locking, $`\mathrm{\Delta }S_1`$, equation (30) provides information about the offset, $`\mathrm{\Delta }I_1`$, of the phase-locked step relative to the Ohmic (linear) curve (see inset in figure 3). The offset is given by the part of the equation that does not depend on $`\mathrm{\Phi }_x^{(0)}`$. From Eq. (30) we have:

$`\mathrm{\Delta }I_1`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{4\mathrm{\Omega }^2+J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}}\left(1{\displaystyle \frac{1}{4}}{\displaystyle \frac{4\mathrm{\Omega }^2+J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)+2J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}{\left(4\mathrm{\Omega }^2+J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)J_2^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)\right)^2}}\right).`$ (32)
Expressions (31) and (32) provide a second order (in $`\mathrm{\Omega }^1`$) description of the phase-locking step magnitude and location as a function of the system parameters, $`\mathrm{\Omega }`$ and $`I_{ac}`$, for large $`\mathrm{\Omega }`$.

## IV RESULTS AND DISCUSSION

We will here show the results of numerical simulations of the system analyzed in the previous section. The simulations are conducted with numerical parameters corresponding to the model parameters, using a fourth order Runge-Kutta method such that the normalized time step typically is no larger than 1% of the period of the driving frequency, and often smaller. Since we are mostly concerned with dc $`IV`$ characteristics, we choose to acquire data for averaging over many ac-periods of motion (typically 10<sup>2</sup>-10<sup>3</sup>) after a sufficient initial time of interval allowed for transient behavior. Simulated $`IV`$ characteristics are obtained by performing the necessary averages as described, and then changing the dc current, $`I_{dc}`$, slightly to acquire the next point on the $`IV`$ curve. All simulations are conducted for the fully frustrated case of $`f=\frac{1}{2}`$.

Figure 3 shows a simulated $`IV`$ characteristic, $`v_x`$ as a function of $`I_{dc}`$, simulated for $`\mathrm{\Omega }=1`$ and $`I_{ac}=2.3`$, for the simple four plaquette model showed in figure 2, described by the equations (1)-(4). As is obvious from the figure, we obtain clear signatures of critical current(s) and phase-locked steps. Specifically, we observe the $`\mathrm{\Delta }S_1`$ step at $`v_x=\mathrm{\Omega }`$, and steps $`\mathrm{\Delta }S_{\frac{1}{2}}`$ and $`\mathrm{\Delta }S_2`$ at $`v_x=\frac{1}{2}\mathrm{\Omega }`$ and $`v_x=2\mathrm{\Omega }`$, respectively. We observe a critical current larger than the previously predicted ccff  value of $`I_c=\sqrt{2}10.41`$ for a fully frustrated dc-driven system. The characteristics of this plot are very similar to the behavior observed in JJA with parallel ac $`+`$ dc drives, obtained both by simulations simff ; octav  and experiments,expfgs  as well as analytically for the four plaquette model.fourp

Comparing simulations of the full-RSJ dynamical equations marconi ; marconi2  for different large arrays ($`N\times N`$ junctions) we will later (below) demonstrate that the simple model of a $`2\times 2`$ array gives very good description of the JJA dynamics. However, we will first compare the predictions of the analytical treatment of the previous section to numerical results.

### IV.1 Critical Current and Phase-Locking for $`\mathrm{\Omega }1`$

In order to verify the simple theory for critical current and phase-locking behavior developed above, we have conducted numerical simulations of the four plaquette system described by Eqs. (1)-(4) for $`1\mathrm{\Omega }`$.

The first set of simulations are conducted to investigate the critical current, $`I_c`$, as it is described by Eq. (21). This expression provides an estimate of the critical current, $`I_c^{}`$, for which the JJA switches from a zero-voltage state ($`V_x=0`$ and $`I_{dc}<I_c^{}`$) to a non-zero voltage state. The simulations are conducted accordingly, starting the system at rest for small $`I_{dc}`$ and slowly increasing the dc-bias until non-zero average voltage is detected. The results are shown in figure 4a, where the solid curve represents the expression, Eq. (21), and the markers represent the simulation results for several frequencies $`1\mathrm{\Omega }3`$ as a function of the characteristic ratio, $`\frac{I_{ac}}{\mathrm{\Omega }}`$. The size of the markers are larger than the error on the estimated critical current. It is obvious that the agreement is very good for all simulated data sets, and we conclude that the critical current, as given by Eq. (21), is a relevant estimate for $`\mathrm{\Omega }`$ not smaller than 1.

Figure 4b shows the critical current, $`I_c^{}`$, evaluated from numerical simulations when the dynamical system switches from the non-zero voltage state to the zero-voltage state (see inset in figure 3). We have here shown the results of numerical simulations with markers as for figure 4a, together with the solid curve of Eq. (21). However, it is clear from the figure that the critical current, $`I_c^{}`$, for decreasing dc-bias may be smaller than for increasing dc-bias ($`I_c^{}I_c^{}`$). Since this is a multi-dimensional system, the critical current may be hysteretic, such that decreasing the dc-current, $`I_{dc}`$, for non-zero voltage states ($`v_x=\dot{\mathrm{\Phi }}_x0`$) is subject to different critical characteristics. One simple way of investigating this is to assume a non-phase-locked state of voltage $`v_x0`$, such that $`\mathrm{\Phi }_x=v_xt`$. A primitive analysis can provide a hint to this hysteresis.

The critical current analysis of the previous section, is obviously a critical current for a system operated at the $`\dot{x}=0`$ branch of the $`IV`$-curve. We may instead analyze what may happen for a non-phase-locked $`\dot{\mathrm{\Phi }}_x=v_x0`$ state. We will still assume (17) to be an appropriate description of the transverse current. However, equation (13) becomes (for small $`|\mathrm{\Psi }_y|`$ and with no resonance to $`\mathrm{\Omega }`$):

$`2\dot{\mathrm{\Psi }}_y\mathrm{cos}v_xt+J_0\left({\displaystyle \frac{I_{ac}}{\mathrm{\Omega }}}\right)\mathrm{\Psi }_y`$ $`=`$ $`0`$ (33)
$`\mathrm{\Psi }_y`$ $`=`$ $`{\displaystyle \frac{J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)}{J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)+4v_x^2}}\mathrm{cos}v_xt+{\displaystyle \frac{2v_x}{J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)+4v_x^2}}\mathrm{sin}v_xt.`$ (34)
Inserting the $`\mathrm{\Psi }_y`$ solution into equation (11) yields the static component:

$`v_x+\left|J_1\left({\displaystyle \frac{2v_x}{J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)+4v_x^2}}\right)\right|`$ $`=`$ $`I_{dc},`$ (35)
where the second term on the left hand side is the result of the resonant mixing between the propagation, $`\dot{\mathrm{\Phi }}_x=v_x`$ and the transverse oscillation, $`\mathrm{\Psi }_y`$. However, the overdamped dc driven pendulum equation is also subject to the following simple relationship, barone

$`\sqrt{v_x^2+(I_c^{})^2}`$ $`=`$ $`I_{dc}.`$ (36)
Combining the two expressions provides the relationship

$`I_c^{}`$ $`=`$ $`\sqrt{\left[v_x+\left|J_1\left({\displaystyle \frac{2v_x}{J_0^2\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)+4v_x^2}}\right)\right|\right]^2v_x^2}\mathrm{\hspace{0.33em}\hspace{0.33em}0.826591}.`$ (37)
As we have indicated, the critical current, $`I_c^{}`$ has a maximum value at around 0.82 (for $`J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)=0`$ and $`v_x0.33`$). Thus, we can argue that propagating ($`\dot{\mathrm{\Phi }}_x0`$) solutions may exist for $`I_{dc}>I_c^{}0.82`$, which provides for a hysteretic $`IV`$ characteristic switching between zero and non-zero voltages in the range $`I_c^{}I_{dc}I_c^{}`$, when $`I_c^{}I_c^{}`$. Notice that when $`I_c^{}>I_c^{}`$, the relevant critical current for both zero and non-zero voltage states must be $`I_c^{}`$, since no static states exist for $`I_{dc}>I_c^{}`$. However, when $`I_c^{}<I_c^{}`$, the actual critical current for switching into a zero-voltage state may be anywhere in the interval $`[I_c^{};I_c^{}]`$.

Figure 4 clearly indicates the hysteretic switching in the $`IV`$ characteristics when $`I_c^{}I_{dc}I_c^{}`$, which is the case for small $`|J_0(\frac{I_{ac}}{\mathrm{\Omega }})|`$.

We note that the above rather primitive analysis of the hysteresis provides a fairly good agreement with the results of numerical simulations. The results of the analysis are not completely consistent with its assumptions in that the resulting amplitude of $`\mathrm{\Psi }_y`$ for the optimized $`v_x0.33`$ is about 1.5, which is not a small number. However, a more detailed (nonlinear) analysis of a single frequency representation of $`\mathrm{\Psi }_y`$ yields quantitatively similar and qualitatively identical results ($`I_c^{}0.77`$) as the above, and we therefore conclude that the simple explanation for hysteresis presented here is relevant.

For $`J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)=0`$ we can provide an explicit approximate expression for $`I_c^{}`$ by assuming the critical current is given by the value of $`v_x`$ which optimizes $`J_1\left(\frac{1}{2v_x}\right)`$. This leads to $`v_x1/3.6`$, which inserted into the above equation yields the optimized coordinates: $`(I_c^{},v_x)=(0.81,0.28)`$.

It is noteworthy that we observe, as predicted by the expression for $`I_c^{}`$, $`I_c(\mathrm{\Omega },I_{ac})`$ to be larger than the dc-driven system, $`I_c(\mathrm{\Omega },I_{ac})I_c(\mathrm{\Omega },0)0.41`$. Hence, a transverse ac driving leads to an enhancement of the critical current. This is contrary to the case with the ac-current parallel to the dc-current, where the critical current is reduced, i.e. $`I_cI_c(\mathrm{\Omega },0)0.41`$.simff ; octav ; expfgs ; fourp

The predicted range, $`\mathrm{\Delta }S_1`$, in $`I_{dc}`$ of phase-locking, as given by Eq. (31), is investigated through simulations similar to the above study of the critical current. Comparisons between the predicted expression and results of numerical simulations are shown in figure 5a, which displays the largest magnitude of the range in dc current for which phase-locking is observed as a function of the characteristic ratio, $`\frac{I_{ac}}{\mathrm{\Omega }}`$, for different values of $`\mathrm{\Omega }`$ in the interval $`1\mathrm{\Omega }3`$. Markers represent results of numerical simulations and solid curves represent the corresponding predicted results of Eq. (31). It is obvious that the simulated parameter sets provide very good overall validation of the perturbation analysis, with the larger of the simulated frequencies providing better agreement than the smaller, as expected. However, an observation common to all simulated frequencies is that large deviations from the expected behavior are found for parameter values ($`\mathrm{\Omega }`$ and $`I_{ac}`$) leading to $`J_0(\frac{I_{ac}}{\mathrm{\Omega }})0`$. The reason for this discrepancy is likely due to a dynamical instability, which can be explained by the perturbation analysis above. The average equilibrium position, $`\mathrm{\Psi }_y`$, of the variable $`\mathrm{\Psi }_y`$ can be observed from Eq. (26) if we write $`\mathrm{\Psi }_y=\mathrm{\Psi }_y^{(0)}+\psi _y`$, where $`\mathrm{\Psi }_y^{(0)}`$ is varying slowly in time (much slower than $`\mathrm{\Omega }`$), and $`\psi _y`$ represents all high frequency (including $`\mathrm{\Omega }`$) contributions ($`|\psi _y|1`$). The slow evolution of equation (26) can then be written,

$`2\dot{\mathrm{\Psi }}_y^{(0)}+J_0\left({\displaystyle \frac{I_{ac}}{\mathrm{\Omega }}}\right)\mathrm{sin}\mathrm{\Psi }_y^{(0)}`$ $`=`$ $`0.`$ (38)
Thus, we find that the stable position of $`\mathrm{\Psi }_y`$ is:

$`\mathrm{\Psi }_y`$ $`=`$ $`\{\begin{array}{ccc}0& ,& J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)>0\\ \pi & ,& J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)<0\end{array}.`$ (41)
The consequence of this abrupt transition in $`\mathrm{\Psi }_y`$ is that the locking phase, $`\mathrm{\Phi }_x^{(0)}`$, must experience a similar abrupt transition of $`\pi `$, as can be seen from equation (11). We therefore claim that the apparent discrepancy observed between the numerical simulations and the perturbation theory near the roots of $`J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)`$, is a result of dynamical instabilities arising from switching the average phase, $`\mathrm{\Psi }_y`$ between $`0`$ and $`\pi `$.

We finally show the comparisons of the center of the phase-locked step as a function of the characteristic ratio, $`\frac{I_{ac}}{\mathrm{\Omega }}`$, for different values of $`\mathrm{\Omega }`$ in the interval $`1\mathrm{\Omega }3`$. The predicted behavior, Eq. (32), is subject to the same issues as the predicted range of phase-locking since the origin of both expressions is Eq. (30). Figure 5b shows the offset, $`\mathrm{\Delta }I_1=I_1\mathrm{\Omega }`$, between the center of the step, $`I_1`$, and the Ohmic curve. As is the case for the phase-locking range shown in figure 5a, the comparison between numerical simulations (markers) and the corresponding predicted offsets (solid curves) is very good, except for the instabilities near the roots of $`J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)`$. Notice that the phase-locking analysis leading to the predictions Eqs. (31) and (32) does not depend on the sharp $`\mathrm{\Psi }_y`$ transition between $`0`$ and $`\pi `$. The reason is that this transition provides only a sign change in the effective equations of phase-locking, and the magnitudes of locking range and offset are therefore unaffected as long as $`J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)0`$.

Based on the above presented comparisons between the numerical simulations of critical currents, range of phase-locking and position of the phase-locked step in the $`IV`$ characteristics of a transversely ac-driven JJA and the corresponding results from simple perturbation analysis, we conclude that the high frequency behavior is well described by the presented analytical treatment.

### IV.2 Critical Current and Phase-Locking for intermediate and low $`\mathrm{\Omega }`$

In Fig. 6 we show the critical current behavior for intermediate and low frequencies. For intermediate frequencies (figure 6a) we observe how the critical current, $`I_c^{}`$, increasingly deviates from the high frequency behavior outlined above. Even so, we notice that the overall behavior of the critical current is qualitatively well described by the analysis leading to equation (21) for $`\mathrm{\Omega }\frac{1}{2}`$. We have, for comparison, included an example of the critical current for the longitudinally ac-driven JJA as an inset. Not surprisingly, decreasing the frequency further (see figure 6b) results in rather large discrepancy between the high frequency analysis of section IIIA and the numerical simulations, and no universal behavior of the critical current as a function of the characteristic ratio, $`\frac{I_{ac}}{\mathrm{\Omega }}`$, can be found. However, we do observe that the critical current does seem to increase quadratically for small $`I_{ac}`$.

In Fig. 7 we show the phase-locking range, $`\mathrm{\Delta }S_1`$, at $`v_x=\mathrm{\Omega }`$ as a function of $`I_{ac}/\mathrm{\Omega }`$ for intermediate ($`\mathrm{\Omega }>0.5`$) and low frequencies ($`\mathrm{\Omega }<0.5`$). Again, as for the critical current we observe that the intermediate frequency range provides for reasonably good qualitative comparisons between numerical simulations and the high frequency analysis of section IIIA. We have, for comparison, included an example of the comparable range of phase-locking for the longitudinally ac-driven JJA as an inset. A noticeable feature of figure 7a is that the dynamical instability discussed above around $`J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)=0`$ seems to widen as the frequency is lowered. Figure 7b shows how this instability provides for increasing discrepancy between high frequency analysis and numerical simulations. However, we notice that even the very low frequencies retain the basic feature of quadratic growth of the phase-locking range as a function of $`I_{ac}`$ for small $`I_{ac}`$.

### IV.3 Dynamics of Phase-Locking

Let us now analyze in detail the dynamics of the voltage responses (increasing and decreasing dc-current) to elaborate on our previous results: critical current hysteresis and windows without transverse phase locking. We calculate $`IV`$ curves and Lyapunov exponents as a function of $`I_{ac}`$ and $`\mathrm{\Omega }`$. In order to distinguish between periodic or quasi-periodic dynamics and chaotic dynamics we calculate the maximum Lyapunov exponent, $`\lambda `$, following the standard methods of nonlinear dynamics.strog ; liap  This means that a small perturbation, $`\stackrel{}{ϵ}(0)`$, to the initial condition will displace the new trajectory by an amount $`|\stackrel{}{ϵ}(t)||\stackrel{}{ϵ}(0)|e^{\lambda t}`$. The Lyapunov exponent is then defined as

$$\lambda =\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\frac{|\stackrel{}{ϵ}(t)|}{|\stackrel{}{ϵ}(0)|}=\underset{t\mathrm{}}{lim}\lambda (t).$$

To recognize a chaotic trajectory we evaluate the maximum Lyapunov exponent. If $`\lambda >0`$ the trajectory is locally unstable; i.e., initial points that are arbitrarily close to each other are macroscopically separated by the flow after a sufficiently long time and the attractor is chaotic. Negative Lyapunov exponents are obtained when trajectories that start sufficiently close to a subset are attracted to it. Here we will show two particular cases: $`I_{ac}/\mathrm{\Omega }=3.0/1.5=2`$, corresponding to a set of parameters where no hysteresis is obtained, and $`I_{ac}/\mathrm{\Omega }=4.05/1.5=2.7`$, corresponding to the hysteretic regime. In Fig. 8 we plot the $`IV`$ curves and maximum Lyapunov exponents for $`\mathrm{\Omega }=1.5`$ and $`I_{ac}=3.0`$ (a-b) and $`I_{ac}=4.05`$ (c-d). The exponents are estimated from $`\lambda \lambda (t)`$ after a finite time $`t=1024T`$, with $`T=2\pi /\mathrm{\Omega }`$. For $`I_{ac}=3.0`$ we show a range in $`I_{dc}`$ where a wide transverse phase locking step exists (Fig. 8a), and the corresponding maximum Lyapunov exponent, $`\lambda `$, is shown in Fig. 8b. We see that within the step we have $`\lambda <0`$, with the most negative value of $`\lambda `$ at the center of the step. Outside the steps, the Lyapunov exponent is nearly equal to zero, $`\lambda 0`$, corresponding to quasiperiodic behavior. A different behavior is obtained for $`I_{ac}=4.05`$, shown in Fig. 8c, where we see that the step disappears and thus there is no transverse phase locking in the same $`I_{dc}`$ range where we find a step in Fig. 8a. The maximum Lyapunov exponent, plotted in Fig. 8d, is small but positive for the $`I_{dc}`$ range around the expected location of the step, the smallness of $`\lambda `$ implies that the dynamics can be either chaotic ($`\lambda >0`$) or quasiperiodic ($`\lambda =0`$) in this case. The $`IV`$ curves near the corresponding critical currents are shown as insets. We see that the absence of hysteresis in critical current is associated with the occurrence of transverse phase locking. Inversely, hysteresis in critical current is obtained for approximately the same parameters for which transverse phase locking is absent. This is in agreement with the above analysis that indicates the critical current hysteresis is present in the vicinity of $`J_0\left(\frac{I_{ac}}{\mathrm{\Omega }}\right)=0`$, which is also the location of the dynamical instabilities of the locked phase of the $`\mathrm{\Delta }S_1`$ step.

In summary, around the transverse phase locking step region we can distinguish three different voltage responses: A, B and C, which are indicated in Fig. 8. We now calculate the voltage power spectrum and Poincaré sections to distinguish these three types of dynamical behavior. This way to characterize dynamical behaviors was used before in capacitive rf-biased JJA. kautz ; thomas  We analyze both transverse and longitudinal voltage power spectra. From the instantaneous transverse voltage we obtain the transverse voltage power spectrum:

$$S_y(\omega )=\left|\frac{1}{T_t}_0^{T_t}𝑑tV_y(t)\mathrm{exp}(i\omega t)\right|^2,$$
(42)
where $`T_t=N_t\mathrm{\Delta }t`$. From the instantaneous longitudinal voltage, $`V_x`$, we obtain the longitudinal voltage power spectrum:

$$S_x(\omega )=\left|\frac{1}{T_t}_0^{T_t}𝑑tV_x(t)\mathrm{exp}(i\omega t)\right|^2.$$
(43)
For studying the nature of the attractor in the different regimes it is useful to consider a Poincaré section of the phase-space trajectories. strog  We consider the stroboscopic Poincaré section of the trajectories in the $`d\mathrm{\Phi }_x/dt`$ vs. $`\mathrm{sin}\mathrm{\Phi }_x`$ plane, recording the values taken by these variables each period of the ac drive. In Fig. 9 we show the power spectra and Poincare sections for $`\mathrm{\Omega }=1.5`$ and for the $`I_{ac}`$ and $`I_{dc}`$ values corresponding to the A, B and C regimes. For each case we show the longitudinal, $`S_x(\omega )`$, and transverse, $`S_y(\omega )`$, voltage power spectra as a function of $`\omega /\mathrm{\Omega }`$ and their corresponding Poincaré sections. Let us first discuss the case corresponding to the regime B, in which there is transverse phase-locking. This is shown in Fig. 9b for $`I_{dc}=1.66`$ and $`I_{ac}=3.0`$, which corresponds to the step with mean voltage $`v_x=V_x=\mathrm{\Omega }`$ (see Fig. 8a, regime B). We see that the longitudinal power spectrum, $`S_x(\omega )`$, presents a delta-like peak for $`\omega =2\mathrm{\Omega }`$. Thereby, the first harmonic of longitudinal voltage fluctuations is locked to $`2\mathrm{\Omega }`$, as expected for this step, since it corresponds to $`n=2`$ in the phase-locking condition $`\omega _0=n\mathrm{\Omega }`$. The phase-locking with a double frequency corresponds to the case when the vortex lattice oscillates in full synchrony with the transverse ac-current, and the ground state repeats itself after one period of the ac drive. In the transverse voltage power spectrum, $`S_y`$, there is a sharp peak at $`\mathrm{\Omega }`$. This is characteristic of transverse phase-locking: the dynamics in the transverse direction locks at half the frequency than the dynamics of the longitudinal direction. This is so because in a single period of the ac-drive, $`T=2\pi /\mathrm{\Omega }`$, the longitudinal component moves forward $`n`$ steps in the lattice period $`a`$, while the transverse component completes only the first half of its oscillation. In Fig. 9e, we show the Poincaré section corresponding to this case in the regime B. The figure shows a very localized Poincaré section since the trajectory always comes back approximately to the same location in phase space in each ac cycle, since the trajectory is periodic (closed orbit).

Now we analyze the case corresponding to the A regime, which is for a current outside the step, $`I_{dc}=1.5`$, see Fig. 8a. In this case we see again in $`S_x(\omega )`$ a peak at $`2\mathrm{\Omega }`$ and in $`S_y(\omega )`$ a peak at $`\mathrm{\Omega }`$. However, the peaks now have a small broadening, and small amplitude satellite peaks have appeared at neighboring frequencies. This is evidence of another kind of long-term behavior, namely quasiperiodic dynamics. We can corroborate this with the corresponding Poincaré section shown in Fig. 9d. It consists now on a closed one dimensional curve, which means that trajectories wind around on a torus, never intersecting itself and yet never quite closing, typical of a quasiperiodic orbit. We have also looked at the time dependent estimates of the Lyapunov exponent, $`\lambda (t)`$. We find that $`\lambda (t)<0`$ for finite $`t`$, but its absolute value tends to zero for long times as $`1/t`$, consistent with quasiperiodic behavior.

Let us now study the last case, corresponding to regime C. This is done for $`I_{dc}=1.66`$ and $`I_{ac}=4.05`$ in Fig. 9c. We see that there are broad peaks in the spectrum in both directions, $`S_x(\omega )`$ and $`S_y(\omega )`$, and that there is a marked increase in the power spectra for low frequencies. Moreover, in Fig. 9f we show the corresponding Poincaré section which consists on successive points jumping from one region of phase space to another and forming a complex curve, which does not seem to close on itself. It is rather difficult to decide from this plot if it corresponds to a quasiperiodic orbit or to a low dimensional attractor of a weakly chaotic orbit. We have obtained also the time dependent estimate of the Lyapunov exponent also for this case. We find that $`\lambda (t)>0`$ for all $`t`$, but its magnitude is decreasing with time as $`1/t`$ as far as we have been able to observe. The fact that $`\lambda (t)`$ is possitive for finite $`t`$ means that there is a dynamical instability that causes a seemingly chaotic behavior at intermediate times and a large noise as seen in the low frequency power spectrum. However, for long times it is very likely that the system will settle in a quasiperiodic dynamics with $`\lambda (t)0`$. In any case, the regime C is very different from the regime A, as can be seen by comparing the power spectra of Fig. 9a and Fig. 9c.

Another view of the dynamics can be obtained by looking at the behavior of the Lyapunov exponent and the noise in the region, where a step is expected, as a function of $`I_{ac}/\mathrm{\Omega }`$. We proceed as follows: for a given $`I_{ac},\mathrm{\Omega }`$, we compute the set of values of Lyapunov exponents $`\lambda `$ and low frequency longitudinal noise $`S_0=lim_{\omega 0}S_x(\omega )`$ that correspond to currents $`I_{dc}`$ in the region of voltage where a step is expected. (We look at $`I_{dc}`$ values for which $`\mathrm{\Omega }ϵ<V_x<\mathrm{\Omega }+ϵ`$, we consider $`ϵ=0.005`$). We plot the resulting set of values of $`\lambda `$ as a function of $`I_{ac}/\mathrm{\Omega }`$ in Fig. 10a and the values of $`S_0`$ as a functions of $`I_{ac}/\mathrm{\Omega }`$ in Fig. 10b. The vertical lines in the plot correspond to the zeroes of $`J_0(I_{ac}/\mathrm{\Omega })`$. We see clearly that near these values there are windows of dynamical instability where $`\lambda 0`$ and where the noise $`S_0`$ is large. In the regions of phase-locking we find a couple of interesting results that are worth mentioning. (i) The most negative value of the Lyapunov exponent occurs in the middle of the phase-locked step and its magnitude is proportional to the step width $`\mathrm{\Delta }S_1`$, as given by Eq. (31). (ii) The largest value of the noise $`S_0`$ occurs at the edge of the phase-locked step; its magnitude is also proportional to the step width $`\mathrm{\Delta }S_1`$, as given by Eq. (31).

### IV.4 Results for Large JJA

We will now consider the quality of the simple $`2\times 2`$ model as representing the dynamics of large $`N\times N`$ JJAs. It is known that collective effects at high currents may come into play. At high currents, the $`Z_2`$ symmetry of the ground state can be broken because a driving current can induce domain walls. mon ; niels ; marconi2  Simulations of $`IV`$ curves with the RSJ model and free boundary conditions, for $`f=1/2`$ and $`T=0`$ have reported a chaotic regime at $`I>I_c`$ related to the motion of domain walls. falo  It has been shown that open boundary conditions nucleate domain walls leading to a critical current lower than de analytic value $`I_c=0.35<\sqrt{2}1`$ at $`T=0`$.minn  Moreover, Ciria and Giovanella ciria  have shown microscopically that different dynamical states are possible for the longitudinal Shapiro steps. Besides the checkerboard ground state configuration, other stable solutions with domain-walls are possible. Then, depending on dc current value and history, domain walls can appear, which are not permitted in the four plaquette model. Therefore, in order to evaluate to what extent the four plaquette model is valid in the transverse ac driven case, we have calculated numerically $`IV`$ curves for $`N\times N`$ arrays, for $`N=8,16,32,64`$, with the full RSJ model used before in Ref. marconi, ; marconi2, . We use periodic boundary conditions in both directions in the presence of an external dc current, $`I_{dc}`$, plus a perpendicular ac current, $`I_{ac}\mathrm{sin}(\mathrm{\Omega }t)`$. We solve the dynamical equations with time step $`\mathrm{\Delta }t=0.1\tau _J`$ ($`\tau _J=2\pi cR_NI_0/\mathrm{\Phi }_0`$) and total integration time $`t_{\mathrm{int}}=2^{15}\mathrm{\Delta }t`$ after a transient $`t_{\mathrm{int}}/2`$. We calculate $`IV`$ curves as a function of $`I_{ac}`$ and $`\mathrm{\Omega }`$, increasing dc current, $`I_{dc}^{}`$, from checkerboard ground state at $`I_{dc}=0`$ and then decreasing dc current, $`I_{dc}^{}`$, from the phase configuration obtained at high current. We use a dc current step $`\mathrm{\Delta }I_{dc}=0.01`$ to obtain $`I_c`$ and $`\mathrm{\Delta }I_{dc}=0.0001`$ to calculate the step width.

One of the relevant results with the four plaquette model is the dependence of the critical current with $`I_{ac}/\mathrm{\Omega }`$ for high frequencies, as shown in Fig. 4a and Fig. 6a. We have also calculated $`I_c`$ as a function of $`I_{ac}/\mathrm{\Omega }`$ for high $`\mathrm{\Omega }`$ in large JJA arrays. In Fig. 10a we show the case for a particular high frequency value in a $`32\times 32`$ array. We see that it has the same behavior as observed in the four plaquette model: $`I_c(I_{ac},\mathrm{\Omega })I_c(0,0)`$, ranges of $`I_{ac}/\mathrm{\Omega }`$ around the maxima of $`I_c`$ where there is hysteresis, and a quadratic increase with $`I_{ac}`$ for $`I_{ac}/\mathrm{\Omega }1`$. We compare with the analytical results expected for $`I_c^{}`$, Eq. (21) and $`I_c^{}`$, Eq. (37), which are represented by dot-dashed lines. We see that $`I_c^{}`$ obtanained numerically for a large array is in excellent agreement with the analytical result for the $`2\times 2`$ system. This is quite reasonable, since $`I_c^{}`$ corresponds to the limit of stability of the checkerboard ground state, which is well represented by the $`2\times 2`$ model. On the other hand, the $`I_c^{}`$ shows some small deviation from the $`2\times 2`$ result, $`I_c^{}(N\times N)I_c(2\times 2)`$. Also the range in $`I_{ac}/\mathrm{\Omega }`$ where there is hysteresis is bigger in a large system. The current $`I_c^{}`$ corresponds to the low current limit of stability of the moving (non-zero voltage) state. In large systems, the moving state can have domain walls, as was found in Ref. niels, ; marconi2, , and the presence of domain walls can lead to a lower $`I_c^{}`$.

In order to analyze more quantitatively in which $`I_{ac}/\mathrm{\Omega }`$ ranges the collective effects could be more relevant, we focus on two cases: case $`a`$ corresponding to values that do not show hysteresis in the critical current in a small system, but are close to the edge of the $`I_{ac}/\mathrm{\Omega }`$ range of hysteresis, and case $`b`$ corresponding to values that show hysteresis in the $`2\times 2`$ system. For each case we calculate the critical current both increasing and decreasing the dc drive, and therefore they correspond to $`a`$, $`a`$, $`b`$ and $`b`$ in Fig. 10a. We show in the inset of Fig. 10a the critical currents obtained for all these cases as a function of system size, $`N`$. In case $`a`$, corresponding to the non-hysteretic region, we see that there is no size effect in $`a`$ up to $`N=64`$. Also we see that $`a=a`$ for $`N32`$, while for $`N=64`$ we find that hysteresis has appeared and $`a<a`$. In the hysteretic region, case $`b`$, size dependent critical currents are obtained for $`b`$, while $`b`$ is size independent. Moreover, the amplitude of the hysteresis, $`bb`$ weakly increases with system size.

In Fig. 10b we show the range of phase locking $`\mathrm{\Delta }S_1`$ as a function of $`I_{ac}/\mathrm{\Omega }`$ for the $`32\times 32`$ array. We find that, when there is phase-locking, the numerically obtained $`\mathrm{\Delta }S_1`$ is very accurately described by the analytical result of Eq. (31) for the $`2\times 2`$ model. In the inset of Fig. 10b we show the size dependence of $`\mathrm{\Delta }S_1`$ for the case marked as $`c`$ in the plot (it corresponds to the same $`I_{ac}/\mathrm{\Omega }`$ of case $`a`$ of Fig. 10a). There is no appreciable size dependence. As observed in the simulations of the $`2\times 2`$ system, we also find here that the phase-locking is lost near the zeros of $`J_0(I_{ac}/\mathrm{\Omega })`$ due to the presence of dynamical instabilities. Also we observe that the presence of hysteresis in the critical current is nearly coincident with the absence of phase locking. We find that with increasing system size these regimes of dynamical instability are amplified in their extension both in their $`I_{dc}`$ dependence and in their range of $`I_{ac}/\mathrm{\Omega }`$ around the zeros of $`J_0(I_{ac}/\mathrm{\Omega })`$. This means that the dynamical instabilities detected in the four plaquette system can lead to an increased spatiotemporal chaos in larger systems where collective effects are important.

## V CONCLUSIONS

It is important to point out that there are no trivial connections between vortex dynamics in the fully frustrated JJA and that of a commensurate vortex lattice moving in a rectangular pinning potential in a bulk superconductor. reichhardt2  This is so because the London model with vortices interacting through pair potentials apply to JJAs only in the limit of very low vortex density, such that $`f=Ha^2/\mathrm{\Phi }_01`$. marconi  The fully frustrated case represents, in this last respect, an interesting limit for studying, where the complete phase field, rather than just the positions of vortices, should be taken into account to describe the dynamics.

We have found transverse phase locking steps in fully frustrated JJA. This new type of (fractional) giant phase-locking steps presents marked differences with the well known longitudinal fractional giant Shapiro steps. Particularly, the presence of the transverse ac force increases the critical depinning current with respect to the case without ac drive (or with a longitudinal ac drive). We have analyzed both analytically and numerically the behavior of the steps as a function of ac amplitude $`I_{ac}`$ and frequency $`\mathrm{\Omega }`$. For $`I_{ac}/\mathrm{\Omega }1`$, the depinning critical current and the phase-locked step width $`\mathrm{\Delta }S_1`$ for $`V=\mathrm{\Omega }\mathrm{}/2e`$, increase quadratically with $`I_{ac}`$. For $`I_{ac}/\mathrm{\Omega }>1`$ we have found windows of $`I_{ac}/\mathrm{\Omega }`$ where depinning is hysteretic and phase locking is destroyed due to dynamical instabilities. The emergence of a weakly chaotic behavior at zero temperature, in a system with non capacitive junctions, is another particular characteristic of transverse phase locking which is absent in longitudinal phase locking in overdamped JJA. Comparing with the behavior of large fully frustrated arrays we have found that transverse phase locking can be well described by an effective four plaquette model, and that collective effects become more important close to the regions of dynamical instability of the four-plaquette model. Our results could be observed experimentally in JJA. In particular, the enhancement of the critical depinning current with a transverse ac drive could be an interesting experimental consequence of the phenomena reported here.

###### Acknowledgements.

We acknowledge discussions with G. Carneiro, H. Pastoriza, C. Reichhardt and D. Shalom. Parts of this work were supported by the Director, Office of Advanced Scientific Computing Research, Division of Mathematical, Information and Computational Sciences of the U.S. Department of Energy under contract number DE-AC03-76SF00098. The work in Argentina was supported by CNEA, Conicet and ANPCyT (PICT99-03-06343).