format: text

text:
# Melting and transverse depinning of driven vortex lattices in the periodic pinning of Josephson junction arrays

## I INTRODUCTION

The behavior of superconducting vortices in the presence of periodic pinning shows very rich static and dynamic phenomena. The competition between the repulsive vortex-vortex interaction and the attractive periodic pinning potential results in novel vortex structures at low temperatures.  The equilibrium phase transitions of these vortex structures and their various dynamical regimes when driven out of equilibrium are of great interest both experimentally and theoretically. Several techniques have been developed to fabricate in superconducting samples an artificial periodic pinning structure: thickness modulated superconducting films, superconducting wire networks, Josephson junction arrays, magnetic dot arrays, submicron hole lattices, and pinning induced by Bitter decoration. The ground states of these systems can be either commensurate or incommensurate vortex structures depending on the vortex density (i.e the magnetic field). In the commensurate case, a “matching” field is defined when the number of vortices $`N_v`$ is an integer multiple of the number of pinning sites $`N_p`$: $`N_v=nN_p`$. A “submatching” or “fractional” field is defined when $`N_v`$ is a rational multiple of $`N_p`$: $`N_v=fN_p`$ with $`f=p/q`$. One of the main properties of periodic pinning is that there are enhanced critical currents and resistance minima both for fractional and for matching magnetic fields, for which the vortex lattice is strongly pinned.

In the case of Josephson junction arrays (JJA) the discrete lattice structure of Josephson junctions induces an effective periodic pinning potential (the so-called “egg-carton” potential) which at low temperatures confines the vortices at the centers of the unit cells of the network. There are strong commensurability effects for submatching fields $`f=p/q`$, for which the vortices arrange in an ordered $`q\times q`$ superlattice that is commensurate with the underlying array of junctions. The transition temperature $`T_c(f)`$ and the critical current $`I_c(f)`$ have maxima for rational $`f=p/q`$, which have been observed experimentally. Moreover, as we will discuss here, the model that describes the physics of the JJA can be thought as a discrete lattice London model for thin film superconductors with periodic arrays of holes. However, this comparison can be valid only for low submatching fields since it can not describe the effects of interstitial vortices.

The equilibrium phase transitions at finite temperatures of two dimensional systems with periodic pinning have been studied in the past.  It is possible to have a depinning phase transition of the commensurate ground states at a temperature $`T_p`$ and a melting transition of the vortex lattice at a temperature $`T_M`$.  Franz and Teitel have studied this problem for the case of submatching fields. For $`0<T<T_p`$ there is a pinned phase in which the vortex lattice (VL) is pinned commensurably to the periodic potential and has long-range order. For $`T_p<T<T_M`$ there is a floating VL which is depinned and has quasi-long-range order. For high submatching fields ($`f1/30`$ ) both transitions coincide, $`T_p=T_M`$, while for low submatching fields ($`f1/30`$) both transitions are different with $`T_p<T_M`$.

The non-equilibrium dynamics of driven vortex lattices interacting either with random or periodic pinning shows an interesting variety of behavior.  Many recent studies have concentrated in the problem of the driven VL in the presence of random pinning.  When there is a large driving current the effect of the pinning potential is reduced, and the nature of the fastly moving vortex structure has been under active discussion. The moving vortex phase has been proposed to be either a crystalline structure, a moving glass, a moving smectic or a moving transverse glass.  These moving phases have been studied both experimentally  and numerically.  Motivated by these results, the dynamical regimes of the moving VL in the presence of periodic pinning has also become a subject of interest. At zero temperature, the dynamical phases of vortices driven by an external current with a periodic array of pinning sites has been studied in very detail by F. Nori and coworkers. A complex variety of regimes has been found, particularly for $`N_v>N_p`$ where the motion of interstitial vortices leads to several interesting dynamical phases.

In this paper we study the dynamical regimes of a moving VL in the periodic pinning of a Josephson junction array (JJA) for finite temperatures in the case of a submatching field ($`f=1/25`$). We obtain a phase diagram as a function of the driving current $`I`$ and temperature $`T`$, which is shown in Fig. 1. We find that when the VL is driven by a low current, the depinning and melting transitions can become separated even for a field for which they coincide in equilibrium. Moreover, we can distinguish between the depinning of the VL in the direction of the current drive, and the transverse depinning in the direction perpendicular to the drive. This later case corresponds to the vanishing of the transverse critical current in a moving VL at a given temperature $`T_{tr}`$, or equivalently, to the vanishing of the transverse superconducting coherence. We obtain three distinct regimes at low temperatures: (i)Pinned vortex lattice: for $`0<T<T_p(I)`$ there is an ordered VL which has crystaline long-range order, superconducting coherence (i.e., a finite helicity modulus) and zero resistance both in the longitudinal and transverse directions. (ii)Transversely pinned vortex lattice: for $`T_p(I)<T<T_{tr}(I)`$ there is a moving VL which has anisotropic Bragg peaks, quasi-long range order, transverse superconducting coherence and zero transverse resistivity. There is a finite transverse critical current. This regime also has strong orientational pinning effects in the  and  lattice directions. (iii)Floating vortex lattice: for $`T_{tr}(I)<T<T_M(I)`$ there is a moving VL which is unpinned in both directions and it has quasi-long range crystalline order with a strong anisotropy. Some of the results we discuss here were reported by us previously in a shorter and less detailed version. After our work in Ref. , further studies of thermal effects in a moving vortex lattice with a periodic pinning array have been reported for $`N_v=N_p`$  and for $`N_v>N_p`$.  Some of these results are similar to ours.

The remainder of our paper is organized as follows. In Sec. II we introduce the theoretical model used for the dynamics of the JJA. We also discuss how this model can be mapped to a superconducting film with a periodic array of holes. In Sec. III we discuss the melting transition near equilibrium when the JJA is driven by a very small current. In Sec. IV we present our results for the transport properties, analyzing the temperature dependence of the current-voltage curves and the onset of resistivity at $`T_p`$ for different currents. In Sec. V we present our results for the transverse depinning transitions at different temperatures and driving currents. We study the transverse current-voltage characteristics both at $`T=0`$ and at finite $`T`$. We define the transverse depinning temperature $`T_{tr}`$ from the behavior of the transverse resistivity. In Sec. VI we present in detail our results for the different dynamical regimes. For various fixed driving currents, we analyze the behavior of structure factor, longitudinal and transverse resistance and helicity modulus through the different regimes as a function of temperature. These results constitute the central core of this paper, from which the phase diagram of Fig. 1 was obtained. In Sec. VII we discuss the orientational pinning effects, which are observed when the direction of the driving current is rotated with respect to the JJA square lattice. Finally in Sec. VIII we present our discussion and conclusions. We also add an Appendix at the end, where we provide a detailed definition of the adequate periodic boundary conditions for a JJA with an external magnetic field and an external driving current, as well as the algorithm used for the numerical simulation.

## II MODEL AND DYNAMICS

### A RSJ dynamics

We consider a square Josephson junction array with $`L\times L`$ superconducting nodes. The nodes are in the lattice sites $`𝐧=(n_x,n_y)`$ and their superconducting phases are $`\theta (𝐧)`$. We study the dynamics of JJA using the resistively shunted junction (RSJ) model for the junctions of the square network.\[32, 33, 34, 35, kim93vorstroud, 38\] In this case, the current flowing in the junction between two superconducting nodes in the JJA is modeled as the sum of the Josephson supercurrent and the normal current:

$$I_\mu (𝐧)=I_0\mathrm{sin}\theta _\mu (𝐧)+\frac{\mathrm{}}{2eR_N}\frac{d\theta _\mu (𝐧)}{dt}+\eta _\mu (𝐧,t)$$
(1)
where $`I_0`$ is the critical current of the junction between the sites $`𝐧`$ and $`𝐧+\mu `$, ($`\mu =\widehat{𝐱},\widehat{𝐲}`$), $`R_N`$ is the normal state resistance and

$$\theta _\mu (𝐧)=\theta (𝐧+\mu )\theta (𝐧)A_\mu (𝐧)=\mathrm{\Delta }_\mu \theta (𝐧)A_\mu (𝐧)$$
(2)
is the gauge invariant phase difference with

$$A_\mu (𝐧)=\frac{2\pi }{\mathrm{\Phi }_0}_{𝐧a}^{(𝐧+\mu )a}𝐀𝑑𝐥.$$
(3)
The thermal noise fluctuations $`\eta _\mu `$ have correlations

$$\eta _\mu (𝐧,t)\eta _\mu ^{}(𝐧^{},t^{})=\frac{2kT}{R_N}\delta _{\mu ,\mu ^{}}\delta _{𝐧,𝐧^{}}\delta (tt^{})$$
(4)
In the presence of an external magnetic field $`H`$ we have

$`\mathrm{\Delta }_\mu \times A_\mu (𝐧)`$ $`=`$ $`A_x(𝐧)A_x(𝐧+𝐲)+A_y(𝐧+𝐱)A_y(𝐧)`$ (5)
$`=`$ $`2\pi f,`$ (6)
$`f=Ha^2/\mathrm{\Phi }_0`$ and $`a`$ is the array lattice spacing. We take periodic boundary conditions (p.b.c) in both directions in the presence of an external current $`𝐈=I\widehat{𝐲}`$ in arrays with $`L\times L`$ junctions. (See the Appendix). The vector potential is taken as

$$A_\mu (𝐧,t)=A_\mu ^0(𝐧)\alpha _\mu (t)$$
(7)
where in the Landau gauge $`A_x^0(𝐧)=2\pi fn_y`$, $`A_y^0(𝐧)=0`$ and $`\alpha _\mu (t)`$ allows for total voltage fluctuations under periodic boundary conditions. In this gauge the p.b.c. for the phases are:

$`\theta (n_x+L,n_y)`$ $`=`$ $`\theta (n_x,n_y)`$ (8)
$`\theta (n_x,n_y+L)`$ $`=`$ $`\theta (n_x,n_y)2\pi fLn_x.`$ (9)
We also consider local conservation of current,

$$\mathrm{\Delta }_\mu I_\mu (𝐧)=\underset{\mu }{}I_\mu (𝐧)I_\mu (𝐧\mu )=0.$$
(10)
After Eqs. (1,10) we obtain the following equations for the phases,

$$\frac{\mathrm{}}{2eR_N}\mathrm{\Delta }_\mu ^2\frac{d\theta (𝐧)}{dt}=\mathrm{\Delta }_\mu [S_\mu (𝐧)+\eta _\mu (𝐧,t)]$$
(11)
where

$$S_\mu (𝐧)=I_0\mathrm{sin}[\mathrm{\Delta }_\mu \theta (𝐧)A_\mu ^0(𝐧)\alpha _\mu ],$$
(12)
and the discrete Laplacian is

$`\mathrm{\Delta }_\mu ^2f(𝐧)`$ $`=`$ $`f(𝐧+\widehat{𝐱})+f(𝐧\widehat{𝐱})+f(𝐧+\widehat{𝐲})+f(𝐧\widehat{𝐲})`$ (14)
$`4f(𝐧).`$
The Laplacian can be inverted with the square lattice Green’s function $`G_{𝐧,𝐧^{}}`$:

$$\mathrm{\Delta }_\mu ^2G_{𝐧,𝐧^{}}=\delta _{𝐧,𝐧^{}}.$$
(15)
Since we take periodic boundary conditions (see Appendix), the total current has to be fixed by:

$`I_x`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}\left[{\displaystyle \underset{𝐧}{}}I_0\mathrm{sin}\theta _x(𝐧)+\eta _x(𝐧,t)\right]+{\displaystyle \frac{\mathrm{}}{2eR_N}}{\displaystyle \frac{d\alpha _x}{dt}},`$ (16)
$`I_y`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}\left[{\displaystyle \underset{𝐧}{}}I_0\mathrm{sin}\theta _y(𝐧)+\eta _y(𝐧,t)\right]+{\displaystyle \frac{\mathrm{}}{2eR_N}}{\displaystyle \frac{d\alpha _y}{dt}},`$ (18)
These equations determine the dynamics of $`\alpha _\mu (t)`$.  For this case we take $`I_x=0`$ and $`I_y=I`$. After Eqs. (11,15,LABEL:tot) we obtain the following set of dynamical equations,

$`{\displaystyle \frac{d\theta (𝐧)}{dt}}`$ $`=`$ $`{\displaystyle \underset{𝐧^{}}{}}G_{𝐧,𝐧^{}}\mathrm{\Delta }_\mu \left[S_\mu (𝐧^{})+\eta _\mu (𝐧^{},t)\right],`$ (19)
$`{\displaystyle \frac{d\alpha _\mu }{dt}}`$ $`=`$ $`I_\mu {\displaystyle \frac{1}{L^2}}{\displaystyle \underset{𝐧}{}}S_\mu (𝐧)+\eta _\mu (𝐧,t),`$ (20)
where we have normalized currents by $`I_0`$, time by $`\tau _J=2eR_NI_0/\mathrm{}`$, and temperature by $`I_0\mathrm{\Phi }_0/2\pi k_B`$.

### B Comparison with thin film with a periodic array of holes

Let us consider a superconducting thin film with a square array of holes, which act as pinning sites for vortices. There are $`N_p=L^2`$ pinning sites separated by a distance $`a`$. The current density in the superconducting film is given by the sum of the supercurrent and the normal current:

$`𝐉`$ $`=`$ $`𝐉_S+𝐉_N`$ (21)
$`𝐉`$ $`=`$ $`{\displaystyle \frac{ie\mathrm{}}{m^{}}}\left[\mathrm{\Psi }^{}𝐃\mathrm{\Psi }(𝐃\mathrm{\Psi })^{}\mathrm{\Psi }\right]+{\displaystyle \frac{\sigma \mathrm{\Phi }_0}{2\pi c}}{\displaystyle \frac{}{t}}\left(\theta {\displaystyle \frac{2\pi }{\mathrm{\Phi }_0}}𝐀\right)`$ (22)
with $`𝐃=+i\frac{2\pi }{\mathrm{\Phi }_0}𝐀`$, $`\mathrm{\Psi }(𝐫)=|\mathrm{\Psi }(𝐫)|\mathrm{exp}[i\theta (𝐫)]`$ the superconducting order parameter and $`\sigma `$ the normal state conductivity. These equations are valid everywhere in the film except in the hole regions. If the number of vortices $`N_v=BL^2a^2/\mathrm{\Phi }_0`$ is much smaller than the number of pinning sites $`N_p`$, all vortices will be centered in the holes in equilibrium. In this case we can assume that $`|\mathrm{\Psi }(𝐧)||\mathrm{\Psi }_0|`$ is homogeneous in the superconducting film. Therefore the dynamics is given by the superconducting phase $`\theta (𝐫)`$, corresponding to a London model in a sample with holes. After considering current conservation $`𝐉=0,`$ we obtain the London dynamical equations for the phases in this multiply connected geometry. Since $`N_vN_p`$, we make the approximation of solving the equations in a discrete grid of spacing $`a`$. This means that we take as the relevant dynamical variables the phases $`\theta (𝐫_𝐧)`$ defined in the sites which are dual to the pinning sites. They represent the average superconducting phase in each superconducting square defined by four pinning sites. Therefore, we take the discretization $`𝐫_𝐧=(n_xa,n_ya)=a𝐧`$. \[Pinning sites at centered at positions $`𝐫_𝐩=(n_x+1/2,n_y+1/2)a`$\]. The derivatives in the supercurrent are discretized in a gauge-invariant way as

$$D_\mu \mathrm{\Psi }(𝐫)\frac{1}{a}\left[\mathrm{\Psi }(𝐧+\mu )e^{i2\pi A_\mu (𝐧)/\mathrm{\Phi }_0}\mathrm{\Psi }(𝐧)\right].$$
(23)
After doing this, we obtain an equation analogous to (1). Now $`I_\mu (𝐧)`$ has to be interpreted as current density normalized by $`J_0=2e\mathrm{}|\mathrm{\Psi }_0|^2/ma=\mathrm{\Phi }_0/(8\pi ^2\lambda ^2a)`$, time normalized by $`\tau =c/(4\pi \sigma \lambda ^2)`$, and the fraction of vortices is $`f=N_v/N_p=Ba^2/\mathrm{\Phi }_0`$. This leads to a set of dynamical equations of the same form as Eqs.(19,20). Therefore, we expect that for $`f1`$ the model for a JJA also gives a good representation of the physics of a superconducting film with a square array of holes (meaning that effects of interstitial vortices are neglected for $`N_vN_p`$). In other words, we expect that for a low density of vortices the specific shape of the periodic pinning potential (being either an egg-carton or an array of holes) will not be physically relevant.

### C Quantities calculated and simulation parameters

The Langevin dynamical equations (19,20) are solved with a second order Runge-Kutta-Helfand-Greenside algorithm with time step $`\mathrm{\Delta }t=0.1\tau _J`$. The discrete periodic Laplacian is inverted with a fast Fourier + tridiagonalization algorithm as in Ref.. (See also the Appendix).

We calculate the following physical quantities:

(i) Transverse superconducting coherence: We obtain the helicity modulus $`\mathrm{{\rm Y}}_x`$ in the direction transverse to the current as

$`\mathrm{{\rm Y}}_x`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}{\displaystyle \underset{𝐧}{}}\mathrm{cos}\theta _x(𝐧){\displaystyle \frac{1}{TL^4}}\{\left[{\displaystyle \underset{𝐧}{}}\mathrm{sin}\theta _x(𝐧)\right]^2`$ (25)
$`\left[{\displaystyle \underset{𝐧}{}}\mathrm{sin}\theta _x(𝐧)\right]^2\}`$
Whenever we calculate the helicity modulus along $`x`$, we enforce strict periodicity in $`\theta `$ by fixing $`\alpha _x(t)=0`$. (See Appendix).

(ii)Transport: We calculate the transport response of the JJA from the time average of the total voltage as

$`V_x`$ $`=`$ $`v_x(t)=d\alpha _x(t)/dt`$ (26)
$`V_y`$ $`=`$ $`v_y(t)=d\alpha _y(t)/dt`$ (27)
with voltages normalized by $`R_NI_0`$.

(iii) Vortex structure: We obtain the vorticity at the plaquette $`\stackrel{~}{𝐧}=(n_x+1/2,n_y+1/2)`$ (associated to the site $`𝐧`$) as

$$b(\stackrel{~}{𝐧})=\mathrm{\Delta }_\mu \times \mathrm{nint}[\theta _\mu (𝐧)/2\pi ]$$
(28)
with $`\mathrm{nint}[x]`$ the nearest integer of $`x`$. We calculate the average vortex structure factor as

$$S(𝐤)=\left|\frac{1}{L^2}\underset{\stackrel{~}{𝐧}}{}b(\stackrel{~}{𝐧})\mathrm{exp}(i𝐤\stackrel{~}{𝐧})\right|^2.$$
(29)
## III Transition near equilibrium

We study JJA with a low magnetic field corresponding to $`f=1/25`$, for different system sizes of $`L\times L`$ junctions, with $`L=50,100,150`$. Most of the results are for $`L=100`$, except when it is explicitly specified, and for $`N_t=10^5`$ iterations after a transient of $`N_t/2`$ iterations.

The ground state vortex configuration for $`f=1/25`$ is a tilted square-like vortex lattice (VL), see Fig.2(a). We find that this state is stable for low currents and low temperatures (in fact, the structure of Fig.2(a) corresponds to $`I=0.01`$ and $`T=0.01`$). The lattice is oriented in the $`[4a,3a]`$ direction and commensurated with the underlying periodic pinning potential of the square JJA. The structure factor $`S(𝐤)`$ has the corresponding Bragg peaks at wavevectors $`𝐆`$ in the reciprocal space, as can be seen in Fig.1(b). When the temperature is increased, the VL tends to disorder and above the melting temperature $`T_M`$ a random vortex array with a liquid-like structure factor is obtained, Figs. 2(c) and 2(d).

We find a single equilibrium phase transition ($`I=0`$) at $`T_M0.050\pm 0.003`$, which is in agreement with the melting temperature obtained by Franz and Teitel and Hattel and Wheatley for $`f1/30`$.

We now apply a very low current, $`I=0.01`$, in order to study the near-equilibrium transport response simultaneously with other quantities like structure factor and helicity modulus. We find a phase transition at a temperature $`T_M(I)0.046\pm 0.001`$, which is slightly lower than the equilibrium transition. In Fig.3(a) we see that there is a large jump in the resistance $`R=V/I`$ at $`T_c`$, in good agreement with the first-order nature of the equilibrium transition.  The onset of resistivity is a signature of a depinning transition in the direction of the drive. This occurs simultaneously with a melting of the vortex lattice, corresponding to the vanishing of Bragg peaks, as shown in Fig.3(b) for the two first reciprocal lattice vectors $`𝐆_\mathrm{𝟏}=\frac{2\pi }{a}(4/25,3/25)`$ and $`𝐆_\mathrm{𝟐}=\frac{2\pi }{a}(3/25,4/25)`$. The size dependence of $`S(𝐆,T)`$ is shown in the inset of Fig.3(b). In the direction of the current drive the helicity modulus is ill-defined since total phase fluctuations are allowed. (See the Appendix). However, in the perpendicular direction to the drive the helicity modulus $`\mathrm{{\rm Y}}_x`$ can be calculated, and it is a measure of the transverse superconducting coherence. As we can see in Fig.3(c), transverse superconductivity also vanishes at $`T_M(I)`$. Above $`T_M`$, we find that the $`\mathrm{{\rm Y}}_x(T)`$ has large fluctuations around zero.

## IV Transport properties

Let us now study the transport properties for larger currents. We calculate the current-voltage (IV) characteristics for different temperatures as well as the dc resistance $`R=V/I`$ as a function of temperature for finite currents.

The zero temperature IV curve has a critical current, $`I_c(0)=0.114\pm 0.002`$, see Fig.4, which corresponds to the single vortex depinning current in square JJA,  with the typical square root depence at the onset. Similar behavior has been reported for zero temperature IV curves for low values of $`f`$. Above $`I_c(0)`$ there is an almost linear increase of voltage, corresponding to a “flux flow” regime, where there is a fastly moving VL. The structure factor of the $`T=0`$ moving VL is the same as the corresponding one of the pinned VL \[Fig.2(a)\]. The presence of periodic boundary conditions in our case prevents the occurrence of random or chaotic vortex motion near the critical current, as reported in early simulations with free boundary conditions. In what follows we will restrict our analysis for currents $`I<0.4`$, where the collective behavior of the VL is the dominant physics. (At $`I1`$ there is a sharp increase of voltage when all the junctions become normal, and $`VR_NI`$ for $`I1`$).

The IV curves for finite temperatures are shown in Fig.5. For temperatures below $`T_M`$ there is a nonlinear sharp rise in voltage which defines the apparent critical current $`I_c(T)`$. For example, we can obtain this $`I_c(T)`$ with a voltage criterion, which we choose as $`V<\frac{1}{N_t\mathrm{\Delta }t}=10^4`$. In this case, we find that $`I_c(T)`$ decreases with $`T`$, vanishing at $`T_M`$. It is interesting to point out that all the IV curves for different temperatures have a crossing point at $`I^{}=0.165`$, see Fig.5(a). A crossing in the IVs has also been reported in experiments in amorphous thin films. For temperatures $`T>T_M`$ the IV curves tend to linear resistivity for low currents. This is shown in Fig.5(b) in a log-log plot of $`R(I)=V/I`$ vs $`I`$, where we see that $`R(I)`$ tends to a low current finite value for $`T>T_M`$, while it has a strong nonlinear decrease for $`T<T_M`$.

Let us now study the dc resistance $`R=V/I`$ as a function of temperature for a given applied current in the $`y`$-direction (Fig.6). We start with the perfectly ordered VL as an initial condition at $`T=0`$ and then we slowly increase the temperature, keeping $`I`$ constant. For currents below the $`T=0`$ critical current, $`I<I_c(0)`$, the dc resistance is negligibly small at low $`T`$, and it has a steep increase at a depinning temperature $`T_p(I)`$, corresponding to the onset of vortex motion. The depinning temperature decreases for increasing currents, and the values of $`T_p(I)`$ are coincident with the apparent critical currents $`I_c(T)`$ obtained from the IV curves. For currents higher than $`I_c(0)`$, there is always a large and finite voltage for any temperature. If $`I_c(0)<I<I^{}`$, the $`R(I)`$ increases slightly with $`T`$ tending to a constant value for large $`T`$, while for $`I>I^{}`$ the $`R(I)`$ decreases with $`T`$.

## V Transverse depinning

What is its response to a small current in the transverse direction when the driven vortex lattice is moving? Is the vortex lattice still pinned in the transverse direction? Is there a transverse critical current for a moving VL? The idea of a transverse depinning current was introduced by Giamarchi and Le Doussal in Ref. for moving vortex systems with random pinning at zero temperature. The possibility of such a critical current was later questioned by Balents, Marchetti and Radzihovsky, where it was shown that this is not true for any finite temperature in random pinning; however a strong nonlinear increase of the transverse voltage was predicted at an “effective” transverse critical current. In the case of periodic pinning it is more clear that a transverse critical current will exist at $`T=0`$ since it is a commensurability effect. This has been found in the $`T=0`$ simulation work of Reichhardt et al.. It is also possible that this transverse critical current will still be non-zero at $`T0`$ in periodic pinning. In fact, we have found in our previous work  that there is a thermal transverse depinning in a periodic system, and we will now analyze this behavior in detail.

First, a high longitudinal current $`I>I_c(0)`$ is applied at zero temperature. Then, a current $`I_{tr}`$ is applied in the transverse direction. In this way, a transverse current-voltage characteristics can be obtained for each $`I`$. This is shown in Fig.7(a) for $`I=0.16`$. We clearly see that there is a finite transverse critical current $`I_{c,tr}0.12`$, which is of the order of the single vortex pinning barrier. We also show in Fig.7(a) how the longitudinal voltage $`V`$ changes when the transverse current is varied. For $`I_{tr}<I_{c,tr}`$ the longitudinal voltage $`V`$ is almost constant. At $`I_{c,tr}`$ there is a fast decay of $`V`$. When $`I_{tr}=I`$ we also have $`V=V_{tr}`$ as expected for a drive at a degree of $`\pi /4`$. Later, for $`I_{tr}>I`$ the vortex lattice becomes pinned in the other direction, since now the directions of “longitudinal” and “transverse” current are interchanged.

Another possible measurement is to study thermal transverse depinning. In this case, we start with a longitudinal current $`I`$, then a small transverse current is applied, $`I_{tr}I_{c,tr}`$, and the temperature is slowly increased. In this way, we can measure a transverse resistance $`R_{tr}=V_{tr}/I_{tr}`$. In Fig.7(b) we plot this result for $`I=0.16`$ and $`I_{tr}=0.01`$. We find that for finite low temperatures $`R_{tr}`$ is negligibly small within our numerical accuracy and it has a clear onset at a transverse critical temperature $`T_{tr}`$.

Let us now see how these results depend on the longitudinal current $`I`$ and temperature $`T`$. We have calculated the transverse IV (Tr-IV) curves for different $`I`$ and $`T`$. In Fig.8(a) we show the Tr-IV curves for $`I=0.06`$ (low current regime) and in Fig.8(b) for $`I=0.16`$ (high current regime). In both cases, there is a clear change of behavior in the Tr-IV curves when going through a characteristic transverse critical temperature $`T_{tr}(I)`$. For low temperatures there is a transverse critical current which tends to vanish when $`T`$ aproachs $`T_{tr}(I)`$ from below. In contrast, for $`T>T_{tr}`$ there is a linear resistivity behaviour.

The transverse resistivity $`R_{tr}=V_{tr}/I_{tr}`$ as a function of temperature was calculated for different longitudinal currents (Fig.9). In all cases there is an onset of transverse response at a given temperature $`T_{tr}(I)`$. At low currents, $`I<I_c(0)`$, the transverse depinning temperatures are almost constant, $`T_{tr}0.02`$ tending to increase slowly with $`I`$, see Fig.9(a). On the other hand, for $`I>I_c(0)`$ the transverse depinning temperatures increase clearly with $`I`$, see Fig.9(b).

## VI Non-equilibrium regimes

We will now study the different non-equilibrium regimes of vortex driven lattices and characterize their possible dynamical transitions. The approach we will follow in this Section is to have a fixed current applied in the system and vary the temperature. In this way, we look for the possible transitions as a function of $`T`$ in a similar way as was done near equilibrium in Sec.III.

A few similar studies were done previously in related systems. In Ref.  the melting transition of a moving vortex lattice in the three dimensional XY model was studied in this way. In this work a first order transition was found as a function of temperature in a strongly driven vortex lattice. In Ref.  a current driven two dimensional JJA at zero field was studied. The possibility of a transition as a function of temperature for finite currents below the critical current was analyzed in this case.

In the following, we will separate our study in three ranges of current: (i) low currents, $`I<0.04`$, (ii) intermediate currents $`0.04<I<I_c(0)`$, and high currents $`I>I_c(0)`$.

### A Low currents

We show the results for a low current in Fig.10 for $`I=0.03`$. The longitudinal dc response, $`V/I`$, is negligibly small at low $`T`$ and later it has a sharp increase of two orders of magnitude: this defines the depinning temperature, $`T_p`$ \[Fig.10(a)\]. At higher temperatures, $`T>T_p`$, the resistance is weakly $`T`$-dependent.

Below $`T_p`$ the vortex lattice is pinned and its structure is similar to the ground state: a vortex lattice commensurate with the underlying square array and tilted in the $`[4a,3a]`$ direction. Above $`T_p`$ the vortex lattice is moving and it has an anisotropic structural order. If we analyze the structure factor in two different reciprocal lattice directions $`S(G_1)`$ and $`S(G_2)`$, we can see this clearly (Fig.10(b)). Below $`T_p`$ the VL structure is isotropic and $`S(G_1)=S(G_2)`$. Right above $`T_p`$ the height of the peaks decreases with temperature, and the structure of the depinned VL is clearly anisotropic, $`S(G_1)S(G_2)`$. Finally, at a melting temperature $`T_M`$ the peaks vanish, and the vortex lattice melts into a vortex liquid. In Fig.11 we compare the behavior of the Bragg peaks for two system sizes $`L=50,100`$. We see that the $`S(G_{1,2})`$ are size independent below $`T_p`$ as it should be expected for a pinned phase. For large temperatures $`T>T_M`$, in the liquid phase, the value of $`S(G_{1,2})`$ is strongly size dependent, since it should go as $`L^2`$. On the other hand, for $`T_p<T<T_M`$ we find that the intensity of the Bragg peaks is weakly dependent on system size. The clear change of behavior of the size dependence gives a good criterion to determine $`T_M`$ (Fig. 11).

The helicity modulus in the direction perpendicular to the current, $`\mathrm{{\rm Y}}_x`$, decreases very slowly for $`T<T_p`$. Above $`T_p`$, $`\mathrm{{\rm Y}}_x`$ has a faster decay with important fluctuations and tends to vanish at $`T_M`$. For $`T>T_M`$, $`\mathrm{{\rm Y}}_x`$ oscillates around zero.

Therefore, for a small finite current the depinning and melting transitions become separated with $`T_p<T_M`$.

### B Intermediate currents

At intermediate currents, $`0.04<I<I_c(0)`$, a new transition appears: the transverse depinning of the moving vortex lattice. As discussed in Sec.V, one can measure transverse depinning by applying a small transverse current while the VL is driven with a fixed longitudinal current. This is shown in Fig. 12(a) for the case of $`I=0.06`$ and a small transverse current, $`I_{tr}=0.01`$. We see that there is an onset of transverse voltage at $`T_{tr}=0.02`$. We can also see that this transverse depinning temperature is above the depinning temperature $`T_p`$ for longitudinal resistance \[Fig.12(b), $`T_p=0.013`$\], and below the melting temperature $`T_M`$ for the vanishing of the Bragg peaks \[Fig.12(c), $`T_M=0.037`$\]. Therefore, this transition occurs at an intermediate temperature between the depinning and the melting transitions, $`T_p<T_{tr}<T_M`$. We also show in Fig.12(d) that the helicity modulus begins to fall down slowly at $`T_{tr}`$, while for $`T_{tr}<T<T_M`$ it has strong fluctuations, being difficult to interpret its behavior in this case.

It is interesting to study in more detail the behavior of the structure factor through all these transitions. The temperature dependence of the intensity of the Bragg peaks is shown in Fig.12(c). In Fig.13 we show examples of the structure factor $`S(𝐤)`$ at temperatures in the different regimes. In the pinned phase, the $`S(𝐤)`$ is nearly the same as in the ground state with delta-like Bragg peaks, see Fig.13(a). For the moving VL, we can see that there is less anisotropy in the transversely pinned regime $`T_p<T<T_{tr}`$ \[Fig.13(b)\] than in the floating regime $`T_{tr}<T<T_M`$ \[Fig.13(c-d)\]. Also in the inset of Fig.12(c) one can see that the intensity of the Bragg peaks has a greater dependence with size for $`T_{tr}<T<T_M`$ when compared with the $`T_p<T<T_{tr}`$ regime. Moreover, near $`T_M`$ the VL structure becomes strongly anisotropic, with the peak at $`𝐆_1`$ much larger than the peak at $`𝐆_2`$, see Fig.13(d).

The anisotropy of the Bragg peaks of the moving VL (in the regimes at $`T_p<T<T_M`$) has two characteristics: (i) the width of the peaks increases with $`T`$ in the direction of the applied current (the direction perpendicular to the vortex motion), and (ii) the height of the peaks decreases in the direction of vortex motion. This can be observed in the sequence of structure factors shown in Fig.13(b-d). The details of the broadening of the Bragg peaks are shown in Fig.14(a) for the transversely pinned regime ($`T_p<T=0.018<T_{tr}`$) and in Fig.14(b) for the floating solid regime ($`T_{tr}<T=0.023<T_M`$).

In Fig.15 we show the temperature dependence of $`S(𝐆_\mathrm{𝟏})`$ and $`S(𝐆_\mathrm{𝟐})`$ for different currents. In all the cases the three different regimes in temperatures are observed: (i) pinned regime for $`0<T<T_p`$ with isotropic and size independent Bragg peaks, (ii) transversely pinned regime for $`T_p<T<T_{tr}`$ with weakly anisotropic Bragg peaks, and (iii) the floating solid regime for $`T_{tr}<T<T_M`$ with strong anisotropy and size dependence in the Bragg peaks.

### C High currents

In the case of currents larger than $`I_c(0)`$, the vortex lattice is already depinned at $`T=0`$. As we have seen in Sec. V, this zero-temperature moving VL has a finite transverse critical current, and therefore it is pinned in the transverse direction. When we slowly increase temperature from this state, we find that the transverse resistive response is negligible for finite low temperatures. At a temperature $`T_{tr}`$ there is a jump to a finite transverse resistance $`R_{tr}=V_{tr}/I_{tr}`$. For example, this is shown for $`I=0.16`$ with a small transverse current, $`I_x=0.01`$ in Fig.16(a). The vortex lattice has an anisotropic structural order for all temperatures, i.e, $`S(𝐆_\mathrm{𝟏})S(𝐆_\mathrm{𝟐})`$, and the height of the Bragg peaks vanishes at $`T_M`$, Fig.16(b). In the inset we show $`S(𝐆_\mathrm{𝟐})`$ for different sizes, $`L=50,100`$, and we see that $`T_M`$ is size independent. Similar behaviour is found for $`S(𝐆_\mathrm{𝟏})`$. The transverse helicity modulus $`(\mathrm{{\rm Y}}_x)`$ is almost constant for low temperatures and starts to decrease at $`T_{tr}`$ presenting strong fluctuations for $`T>T_{tr}`$ , Fig16(c).

The analysis of the heights of the peaks in $`S(𝐤)`$ as a function of system size $`L`$ is a good indicator of the translational correlations in the system. This dependence is well known for two dimensional lattices.  For a pinned solid $`S(𝐆)1`$, for a floating solid, $`S(𝐆)L^{\eta _G(I,T)}`$ with $`0<\eta _G(I,T)<2`$, being this dependence a signature of quasi-long range order, and for a normal liquid $`S(𝐆)L^2`$. To assure the existence of algebraic translational correlations, we have done this scaling study for currents $`0.02<I<0.2`$ and different temperatures for system sizes of $`L=50,100,150`$. For the cases corresponding to the pinned regime ($`T<T_p`$) we found $`\eta _G0`$, as expected. In the transversely pinned regime, we show a case in Fig.17(a), we find a power law fitting with very small values of $`\eta _G(I,T)`$. In the floating regime, we find larger values of $`\eta _G`$, we show a case in Fig. 17(b). In all the cases we have obtained that $`0<\eta _G(I,T)<2`$ for $`T_p<T<T_M`$. Therefore, this finite size analysis shows the existence of quasi-long range order in the moving VL. Also, we find that $`\eta _{G_1}>\eta _{G_2}`$ for all currents and temperatures. The power-law exponent $`\eta _G`$ can be studied at a constant current, as a function of temperature. This is shown in Fig.18, for $`I=0.16`$ and different reciprocal lattice vectors. The exponent $`\eta `$ is finite for the complete temperature range. For $`T<T_{tr}`$ it has a small value $`\eta _G0.01`$. It has a fast increase near $`T_{tr}`$, where there is also a clear difference between $`\eta _{G_1}`$ and $`\eta _{G_2}`$. Finally, it reaches a value of $`\eta =2`$ near $`T=T_M`$.

## VII Orientational pinning effects

A very interesting characterization of the different regimes can be obtained by studying the effects of varying current direction.  We apply a current $`I`$ at an angle $`\varphi `$ with respect to the $`[10]`$ lattice direction,

$`I_x`$ $`=`$ $`I\mathrm{cos}\varphi `$ (30)
$`I_y`$ $`=`$ $`I\mathrm{sin}\varphi .`$ (31)
We study the voltage response when varying the orientation $`\varphi `$ of the drive while keeping fixed the amplitude $`I`$ of the current. In the parametric curves of $`V_y(\varphi )`$ vs. $`V_x(\varphi )`$ can we analyze the breaking of rotational symmetry in the different regimes of $`I`$ and $`T`$. In the case of rotational symmetry this kind of plot should give a perfect circle. However, the square symmetry of the Josephson lattice will show up in the shape of the curves. In Ref.  the motion of single vortices was studied. In this case, the “diagonal”  direction is unstable against small changes in the angle $`\varphi `$, while the  and  directions are the preferred directions for vortex motion. This shows as “horns” in parametric $`V_y`$ vs. $`V_x`$ plots, which are finite segments of points lying in the $`x`$ or the $`y`$ axis. This implies the existance of transverse pinning in these directions (thus, it corresponds to orientational pinning). Here, we perform the same analysis for the different regimes of the moving VL.

In Fig.19 we plot the voltages $`V_y`$ and $`V_x`$ when varying the orientational angle $`\varphi `$ for different current amplitudes $`I`$ and temperatures $`T`$. In Fig.19(a) we have $`I=0.06`$ and $`T=0.015`$, corresponding to the regime of a transversely pinned lattice. In this case most of the points are lying either on the axis $`V_x=0`$ or on the axis $`V_y=0`$, indicating strong orientational pinning in the symmetry lattice directions  and . When increasing $`T`$ the orientational pinning decreases and the length of the “horns” in the $`x`$ and $`y`$ axis decreases. Fig.19(b) shows results for $`I=0.06`$ and $`T=0.03`$, which correspond to $`T>T_{tr}`$ when there is a finite transverse resistance. In this case the horns have disappeared and orientational pinning is lost. However, the breaking of rotational symmetry is still present in the star-shaped curves. Also in the high current regime, Fig.19(c), for a low temperature $`T=0.015`$ ($`T<T_{tr}`$) we see that there is orientational pinning with the presence of horns, which again dissappear for $`T>T_{tr}`$ as it is shown in Fig.19(d) at $`T=0.035`$. Finally for $`TT_M`$, deep inside in the liquid phase, the stars tend to the circular shape of rotational invariance.

## VIII DISCUSSION

With all the information of the previous sections we obtain the current-temperature phase diagram which was advanced in Fig.1. For finite currents we have been able to identify three different regimes, a pinned VL for $`T<T_p(I)`$, a transversely pinned VL for $`T_p(I)<T<T_{tr}(I)`$, and a floating VL for $`T_{tr}(I)<T<T_M(I)`$. It is, however, difficult to define if the temperatures $`T_p,T_{tr},T_M`$ correspond to either phase transitions or to dynamical crossovers. The comparisons we have made of the behavior of voltages (longitudinal and transversal), structure factor and helicity modulus show that something is happening at these temperatures. Also, the comparison for different system sizes of the behavior of $`S(𝐆)`$ suggest transitions for $`T_p,T_{tr},T_M`$. The transverse helicity modulus has strong fluctuations for $`T_{tr}<T<T_M`$. These fluctuations are not reduced when increasing the simulation time in a factor of $`10`$. This could mean that actually the range of temperatures of $`T_{tr}<T<T_M`$ corresponds to a long crossover region towards a liquid state. Also, one may question if the use of the helicity modulus for these far from equilibrium states is correct, since $`\mathrm{{\rm Y}}_x`$ has been defined from the response of the equilibrium free energy to a twist in the boundary condition. Most likely, the large fluctuations in $`T_{tr}<T<T_M`$ are due to the very unstable and history-dependent steady states we find in this regime. For example the steady states obtained when increasing temperature from $`T=0`$ at fixed current differ from the steady states obtained when increasing current at fixed temperature in this regime. They differ in the degree of anisotropy and intensity of the Bragg peaks in the structure factor. We think that this reflects the fact that the VL is unpinned in all directions and the orientation of the moving state will depend on the history for this regime.

In any case, we can clearly distinguish different dynamical states which define the regimes shown in the phase diagram of Fig.1. From the experimental point of view these transitions (or crossovers) are possible to measure. The depinning temperature $`T_p`$ can be obtained from resistance measurements at finite currents. The transverse current-voltage characteristics and the transverse resistivity can be measured for different longitudinal currents and temperatures, and therefore $`T_{tr}`$ could be obtained. The helicity modulus can be measured with the two-coil technique.  It could be very interesting to see the results of such measurements at finite currents. The structure factor and melting transition can not be measured directly. However, in the presence of an external rf current, the dissappearance of Shapiro steps can indicate the melting transition,  since they are sensitive to the traslational order of the vortex lattice. Therefore, we expect that the results obtained here could motivate new experiments for fractional or submatching fields in Josephson junction arrays as well as in other superconductors with periodic pinning.

###### Acknowledgements.

We acknowledge fruitful discussions with H. Pastoriza and J. V. José. This work was supported by CONICET, CNEA, ANPCyT (PICT Nro. 03-00121-02151) and Fundación Antorchas (Proyecto Nro. 13532/1-96).

## A Periodic boundary conditions

### 1 Phases

We want to obtain the periodic boundary condition (PBC) for superconducting phases $`\theta (n_x,n_y)`$. In general we can write the PBC as:

$`\theta (𝐧+𝐋_𝐱)`$ $`=`$ $`\theta (𝐧)+u_x(𝐧)`$ (A1)
$`\theta (𝐧+𝐋_𝐲)`$ $`=`$ $`\theta (𝐧)+u_y(𝐧),`$ (A2)
where $`𝐋_𝐱=(L_x,0)`$ and $`𝐋_𝐲=(0,L_y)`$. Taking into account that all variables should be independent of the order of subsequent global translations and that the phases are defined except for an addition of $`2\pi l`$, we are lead to the consistency condition

$$u_x(𝐧)+u_y(𝐧+𝐋_𝐱)=u_x(𝐧+𝐋_𝐲)+u_y(𝐧)+2\pi l,$$
(A3)
Therefore, in order to specify the periodic boundary conditions we have to give the functions $`u_x,u_y`$, which will depend on the gauge for the vector potential $`𝐀`$.

The periodic boundary condition for the phases can be deduced by requiring that all physical quantities are invariant after a translation in the lattice size. The supercurrents $`S_\mu (𝐧)=\mathrm{sin}[\mathrm{\Delta }_\mu \theta (𝐧)A_\mu (𝐧)]`$ should satisfy:

$`S_\mu (𝐧+𝐋_𝐱)`$ $`=`$ $`S_\mu (𝐧)`$ (A4)
$`S_\mu (𝐧+𝐋_𝐲)`$ $`=`$ $`S_\mu (𝐧).`$ (A5)
This implies that the gauge invariant phase difference $`\theta _\mu (𝐧)=\mathrm{\Delta }_\mu \theta (𝐧)A_\mu (𝐧)`$ should satisfy

$`\theta _\mu (𝐧+𝐋_𝐱)`$ $`=`$ $`\theta _\mu (𝐧)+2\pi l`$ (A6)
$`\theta _\mu (𝐧+𝐋_𝐲)`$ $`=`$ $`\theta _\mu (𝐧)+2\pi l^{},`$ (A7)
with $`l,l^{}`$ any integer. This condition leads to

$`\mathrm{\Delta }_\mu u_x(𝐧)`$ $`=`$ $`2\pi l+A_\mu (𝐧+𝐋_𝐱)A_\mu (𝐧)`$ (A8)
$`\mathrm{\Delta }_\mu u_y(𝐧)`$ $`=`$ $`2\pi l^{}+A_\mu (𝐧+𝐋_𝐲)A_\mu (𝐧).`$ (A9)
We can choose the solution with $`l=l^{}=0`$. In the Landau gauge, $`A_x`$ is a linear function of $`n_y`$ and $`A_y`$ is a linear function of $`n_x`$. Taking the origin such that $`A_\mu (𝐧=0)=0`$, we obtain for the Landau gauge:

$`u_x(𝐧)`$ $`=`$ $`u_x(0)+A_y(𝐋_𝐱)n_y`$ (A10)
$`u_y(𝐧)`$ $`=`$ $`u_y(0)+A_x(𝐋_𝐲)n_x.`$ (A11)
The consistency condition (A2) requires

$$A_y(𝐋_𝐱)L_yA_x(𝐋_𝐲)L_x=2\pi l.$$
(A12)
The term in the left side is equal to the total flux $`2\pi fL_xL_y`$, therefore (A7) is equivalent to flux quantization, giving $`fL_xL_y=N_v`$, with $`N_v=l`$ the number of vortices.

If we take $`u_x(0)=u_y(0)=0`$, we obtain for the PBC

$`\theta (n_x+L_x,n_y)`$ $`=`$ $`\theta (n_x,n_y)+A_y(𝐋_𝐱)n_y`$ (A13)
$`\theta (n_x,n_y+L_y)`$ $`=`$ $`\theta (n_x,n_y)+A_x(𝐋_𝐲)n_x.`$ (A14)
A particular choice can be the gauge with $`A_x(𝐧)=2\pi fn_y`$, $`A_y(𝐧)=0`$, which leads to Eq. (7).

### 2 External currents and electric fields

In the presence of external currents or voltages the periodic boundary conditions have to be reconsidered. In this case it is possible to have $`𝐄\mathrm{𝐝𝐥}0`$ in a path that encloses all the sample either in the $`x`$ or the $`y`$ direction. Therefore, in a closed path we have to consider the Faraday’s law

$$𝐄\mathrm{𝐝𝐥}=\frac{1}{c}\frac{d\mathrm{\Phi }}{dt}.$$

The two-dimensional sample with PBC can be thought as the surface of a torus in three-dimensions. The closed paths we are considering are the two pathes that encircle the torus. The electric field is not a gradient of a potential, it is now given by $`𝐄=V\frac{1}{c}\frac{𝐀}{t}`$. One possible solution is to consider a vector potential

$$𝐀(𝐫,t)=𝐀_\mathrm{𝟎}(𝐫)\stackrel{}{\alpha }(t),$$

for which

$`𝐇_{\mathrm{ext}}`$ $`=`$ $`\times 𝐀_\mathrm{𝟎}`$ (A15)
$`𝐄_{\mathrm{total}}`$ $`=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle \frac{d\stackrel{}{\alpha }}{dt}}`$ (A16)
In our case, we take the adimensional vector potential as:

$$A_\mu (𝐧,t)=A_\mu ^0(𝐧)\alpha _\mu (t),$$
(A17)
with $`A_\mu ^0(𝐧)`$ in the Landau gauge ($`A_x^0(𝐧)=2\pi fn_y`$, $`A_y^0(𝐧)=0`$). Therefore the gauge invariant phase is:

$$\theta _\mu (𝐧,t)=\mathrm{\Delta }_\mu \theta (𝐧,t)A_\mu ^0(𝐧)+\alpha _\mu (t).$$

Then $`\alpha _\mu `$ acts as a global time-dependent phase in the $`\mu `$ direction.

In the normalized units used in this paper, the electric field in the link defined by the junction $`𝐧,\mu `$ is

$$E_\mu (𝐧)=\mathrm{\Delta }_\mu V(𝐧)\frac{dA_\mu (𝐧)}{dt},$$

where the electrostatic potential is $`V(𝐧)=\frac{d\theta (𝐧)}{dt}`$. Therefore we have

$$E_\mu (𝐧)=\mathrm{\Delta }_\mu \frac{d\theta (𝐧)}{dt}+\frac{d\alpha _\mu }{dt}.$$

The average electric field in the $`\mu `$ direction is:

$$E_\mu ^{\mathrm{av}}=\frac{1}{L_xL_y}\underset{𝐧}{}E_\mu (𝐧)=\frac{d\alpha _\mu }{dt},$$

where we have used the fact that $`_𝐧\mathrm{\Delta }_\mu \frac{d\theta (𝐧)}{dt}=0`$ (which is the discrete equivalent of $`V\mathrm{𝐝𝐥}=0`$). The current in the link $`𝐧,\mu `$ is, in normalized units,:

$$I_\mu (𝐧)=E_\mu (𝐧)+\stackrel{~}{S}_\mu (𝐧),$$

with $`\stackrel{~}{S}_\mu (𝐧)=S_\mu (𝐧)+\eta _\mu (𝐧,t)`$. Therefore the average current in the $`\mu `$ direction is

$$I_\mu ^{\mathrm{av}}=\frac{1}{L_xL_y}\underset{𝐧}{}I_\mu (𝐧)=E_\mu ^{\mathrm{av}}+S_\mu ^{\mathrm{av}}+\eta _\mu ^{\mathrm{av}}=\frac{d\alpha _\mu }{dt}+\frac{1}{L_xL_y}\underset{𝐧}{}\stackrel{~}{S}_\mu (𝐧)$$

There are two cases to consider: (i) external current source: the external current is given and the total voltage fluctuates and (ii) external voltage source: the external voltage is given and the total current fluctuates.

(i) External current source: the average current in $`\mu `$ direction is fixed by the external current: $`I_\mu ^{\mathrm{av}}=I_\mu ^{\mathrm{ext}}`$. The average electric field is a fluctuating quantity given by:

$$E_\mu ^{\mathrm{av}}(t)=\frac{d\alpha _\mu }{dt}=I_\mu ^{\mathrm{ext}}\frac{1}{L_xL_y}\underset{𝐧}{}\stackrel{~}{S}_\mu (𝐧).$$
(A18)
In this case, the $`\alpha _\mu (t)`$ is a dynamical variable, and its time evolution is given by Eq.(A10).

(ii) External voltage source: the average electric field in the $`\mu `$ direction is fixed by the external electric field: $`E_\mu ^{\mathrm{av}}=E_\mu ^{\mathrm{ext}}`$. Therefore, now the $`\alpha _\mu `$ is given by the external source:

$$\alpha _\mu (t)=E_\mu ^{\mathrm{ext}}𝑑t,$$

and $`\alpha _\mu (t)=E_\mu ^{\mathrm{ext}}t`$, if $`E_\mu ^{\mathrm{ext}}`$ is time-independent. The average current is now a fluctuating quantity given by:

$$I_\mu ^{\mathrm{av}}=E_\mu ^{\mathrm{ext}}+\frac{1}{L_xL_y}\underset{𝐧}{}\stackrel{~}{S}_\mu (𝐧).$$

Let us see how the PBC are affected by a change of gauge. The gauge transformations are the following:

$`\theta (𝐧)\theta ^{}(𝐧)`$ $`=`$ $`\theta (𝐧)+\varphi (𝐧)`$ (A19)
$`A_\mu (𝐧)A_\mu ^{}(𝐧)`$ $`=`$ $`A_\mu (𝐧)+\mathrm{\Delta }_\mu \varphi (𝐧)`$ (A20)
$`V(𝐧)V^{}(𝐧)`$ $`=`$ $`V(𝐧){\displaystyle \frac{d\varphi (𝐧)}{dt}}`$ (A21)
An interesting choice is:

$$\varphi (𝐧)=\stackrel{}{\alpha }𝐧=\underset{\mu }{}\alpha _\mu n_\mu $$

In this gauge we have $`\theta _\mu (𝐧)=\mathrm{\Delta }_\mu \theta (𝐧)A_\mu ^0(𝐧)`$ and the PBC for the phases is:

$`\theta (𝐧+𝐋_𝐱)`$ $`=`$ $`\theta (𝐧)+A_y(𝐋_𝐱)n_y+\alpha _xL_x`$ (A22)
$`\theta (𝐧+𝐋_𝐲)`$ $`=`$ $`\theta (𝐧)+A_x(𝐋_𝐲)n_x+\alpha _yL_y,`$ (A23)
and for the voltages:

$$V(𝐧+𝐋_\mu )=V(𝐧)\frac{d\alpha _\mu }{dt}L_\mu =V(𝐧)E_\mu ^{\mathrm{av}}L_\mu $$
(A24)
The equations of motion in this gauge are:

$`{\displaystyle \frac{d\theta (𝐧)}{dt}}`$ $`=`$ $`{\displaystyle \underset{\mu }{}}{\displaystyle \frac{d\alpha _\mu }{dt}}n_\mu {\displaystyle \underset{𝐧^{}}{}}G_{𝐧,𝐧^{}}\mathrm{\Delta }_\mu \stackrel{~}{S}_\mu (𝐧^{})`$ (A25)
$`{\displaystyle \frac{d\alpha _\mu }{dt}}`$ $`=`$ $`I_\mu ^{\mathrm{ext}}{\displaystyle \frac{1}{L_xL_y}}{\displaystyle \underset{𝐧}{}}\stackrel{~}{S}_\mu (𝐧)`$ (A26)
The periodic boundary conditions with a fixed external current using Eq. (A10) was used previously in Ref. for $`f=0`$ (it was called a “fluctuating twist boundary condition”) and in Ref. for $`f0`$. Also, the periodic boundary conditions in the gauge of Eqs. (A11-A12) were used in Ref.  for a time dependent $`d`$-wave Ginzburg-Landau model.

### 3 Helicity modulus

The helicitity modulus $`\mathrm{{\rm Y}}_\mu `$ expresses the “rigidity” of the system with respect to an applied “twist” in the periodic boundary conditions. The twist $`k_\mu `$ is defined as a phase change of $`L_\mu k_\mu `$ between the two opposite boundaries which are connected through the PBC in the $`\mu `$ direction

$$\theta (𝐧+𝐋_\mu )=\theta (𝐧)+L_\mu k_\mu .$$

The helicity modulus is obtained from the free energy $`F(T,k_\mu )`$ as:

$$\mathrm{{\rm Y}}_\mu =\frac{1}{L^2}\frac{^2F(T,k_\mu )}{k_\mu ^2}|_{k_\mu =0}$$
(A27)
It is clear from Eq. (A11) that $`k_\mu =\alpha _\mu `$. Then, in order to evaluate the helicity modulus, $`\alpha _\mu `$ must be set to zero. This means that the helicity modulus can not be calculated in the direction in which there is an applied current, since it gives a fluctuating twist $`k_\mu (t)`$.

### 4 Algorithm

The core of the numerical calculation is to invert the Eq. (9). This means to solve a discrete Poisson equation of the form,

$$\mathrm{\Delta }_\mu ^2f(𝐧)=d(𝐧)$$
(A28)
with the periodic boundary conditions

$`f(𝐧+𝐋_𝐱)`$ $`=`$ $`f(𝐧)`$ (A29)
$`f(𝐧+𝐋_𝐲)`$ $`=`$ $`f(𝐧).`$ (A30)
The linear system of $`L_xL_y`$ equations of (A28) is singular. Physically, the reason is that in Eq.(11), $`f(𝐧)=d\theta (𝐧)/dt`$ corresponds to a voltage, which is defined except for a constant. We choose the voltage reference such that it has zero mean:

$$\underset{𝐧}{}f(𝐧)=0,$$
(A31)
other choices are also possible (like, for example, fixing $`f(𝐧_\mathrm{𝟎})=0`$ at a given site $`𝐧_0`$).

The method we use to invert the Eq. (A28) is based on the Fourier Accelerated and Cyclic Reduction (FACR) algorithm.  In this case, we take first the discrete Fourier transform in the $`x`$ direction:

$$\stackrel{~}{d}(k_x,n_y)=\underset{n_x}{}d(n_x,n_y)e^{i\frac{2\pi k_xn_x}{L_x}},$$
(A32)
which with a fast Fourier algorithm takes a computation time of order $`L_x\mathrm{log}L_x`$. This leads to the following equation:

$$ϵ_{k_x}\stackrel{~}{f}(k_x,n_y)\stackrel{~}{f}(k_x,n_y1)\stackrel{~}{f}(k_x,n_y+1)=\stackrel{~}{d}(k_x,n_y),$$
(A33)
with $`ϵ_{k_x}=42\mathrm{cos}\frac{2\pi k_x}{L_x}`$ and with boundary condition $`\stackrel{~}{f}(k_x,L_y+1)=\stackrel{~}{f}(k_x,1)`$. This is a cyclic tridiagonal equation which can be solved with a simple LU decomposition algorithm in a computation time of order $`L_y`$. In this way the $`\stackrel{~}{f}(k_x,n_y)`$ is obtained from (A19). Finally we take the inverse Fourier transform to obtain

$$f(n_x,n_y)=\frac{1}{L_x}\underset{k_x}{}\stackrel{~}{f}(k_x,n_y)e^{i\frac{2\pi k_xn_x}{L_x}}.$$
(A34)
This algorithm takes a computation time which is of order $`L_xL_y(\mathrm{log}L_x+A)`$ with the constant $`A1`$. This is faster than the two dimensional Fourier transform method used by Eikmans and van Himbergen, which takes a computation time of order $`L_xL_y(\mathrm{log}L_x+\mathrm{log}L_y)`$.


text:
# Renormalization group flow with unstable particles

## Abstract

The renormalization group flow of an integrable two dimensional quantum field theory which contains unstable particles is investigated. The analysis is carried out for the Virasoro central charge and the conformal dimensions as a function of the renormalization group flow parameter. This allows to identify the corresponding conformal field theories together with their operator content when the unstable particles vanish from the particle spectrum. The specific model considered is the $`SU(3)_2`$-homogeneous Sine-Gordon model.

PACS numbers: 11.10Hi, 11.10Kk, 11.30Er, 05.70.Jk

The study of two-dimensional quantum field theories (2D-QFT) has turned out to be a fruitful venture since almost three decades. In particular when exploiting integrability many non-perturbative methods have been developed over the years. Besides the challenge to understand the underlying mathematical structures and the intriguing physical applications in two dimensions itself, e.g. to describe measurable quantities of carbon nanotubes , the ultimate goal is to extrapolate ones findings to higher dimensions. In particular for the celebrated c-theorem of Zamolodchikov , which originally describes the renormalization group trajectory of a function which at the renormalization group fixed point corresponds to the Virasoro central charge, various counterparts have been developed in higher dimensions, e.g. .

Fairly recently a class of massive integrable quantum field theories, the homogeneous Sine-Gordon models (HSG) , has been proposed introducing the feature of possessing unstable particles inside its particle spectrum. Despite the fact that theories containing resonances have been treated before in the context of two-dimensional massive quantum field theories, e.g. , the HSG-models are somewhat special since they constitute the first examples of theories which admit a well-defined Lagrangian description. In general the HSG models are associated to integrable perturbations of $`G`$-parafermions of level $`k`$ , i.e. WZNW-coset theories of the form $`G_k/U(1)^{\mathrm{}}`$ with $`\mathrm{}`$ being the rank of a compact Lie group $`G`$. As free parameters the model contains $`\mathrm{}`$ different mass scales and $`\mathrm{}1`$ different scales for the resonance parameter $`\sigma `$ which enter the Breit-Wigner formula . In general an unstable particle of type $`\stackrel{~}{c}`$ is described by complexifying the physical mass of a stable particle by adding a decay width $`\mathrm{\Gamma }_{\stackrel{~}{c}}`$, such that it corresponds to a pole in the S-matrix as a function Mandelstam $`s`$ at $`s=M_R^2=(M_{\stackrel{~}{c}}i\mathrm{\Gamma }_{\stackrel{~}{c}}/2)^2`$ (for a more detailed discussion see e.g. ). As mentioned in  whenever $`M_{\stackrel{~}{c}}\mathrm{\Gamma }_{\stackrel{~}{c}}`$, the quantity $`M_{\stackrel{~}{c}}`$ admits a clear cut interpretation as the physical mass. However, since this assumption is only required for interpretational reasons we will not rely on it. Transforming as usual in this context from $`s`$ to the rapidity plane and describing the scattering of two stable particles of type $`a`$and $`b`$ with masses $`m_a`$ and $`m_b`$ by an S-matrix $`S_{ab}(\theta )`$ as function of the rapidity $`\theta `$, the resonance pole is situated at $`\theta _R=\sigma i\overline{\sigma }`$. Identifying the real and imaginary parts of the pole then yields

$`M_{\stackrel{~}{c}}^2{\displaystyle \frac{\mathrm{\Gamma }_{\stackrel{~}{c}}^2}{4}}`$ $`=`$ $`m_a^2+m_b^2+2m_am_b\mathrm{cosh}\sigma \mathrm{cos}\overline{\sigma }`$ (1)
$`M_{\stackrel{~}{c}}\mathrm{\Gamma }_{\stackrel{~}{c}}`$ $`=`$ $`2m_am_b\mathrm{sinh}|\sigma |\mathrm{sin}\overline{\sigma }.`$ (2)
Eliminating the decay width from (1) and (2), we can express the mass of the unstable particles $`M_{\stackrel{~}{c}}`$ in the model as a function of the masses of the stable particles $`m_a,m_b`$ and the resonance parameter $`\sigma `$. Assuming $`\sigma `$ to be large this gives

$$M_{\stackrel{~}{c}}^2\frac{1}{2}m_am_b(1+\mathrm{cos}\overline{\sigma })e^{|\sigma |}.$$
(3)
One recognizes the occurrence of the variable $`me^{|\sigma |/2}`$, which was introduced originally in  in order to describe massless particles, i.e. one may perform safely the limit $`m0,\sigma \mathrm{}`$ and one might therefore be tempted to describe flows related to (3) as massless flows. In  the relative mass scales between the unstable and stable particles and the stable particles themselves were investigated by computing the finite size scaling function from the thermodynamic Bethe ansatz (TBA). A consistent physical picture was obtained for the overall identification of the flow between different coset models. It remained, however, an open question how to identify the operator content. In general this question is left unanswered in the context of the TBA. For theories with certain properties, it is sometimes possible to determine at least the dimension of the perturbing operators by investigating periodicities in the so-called Y-systems . Resorting to a different method, namely by appealing to sum rules which are expressible in terms of correlation functions, the major part of the operator content was successfully identified for some of the HSG models . The purpose of this manuscript is on one hand to confirm and refine the TBA results by the latter method, i.e. by investigating the renormalization group flow described by the Zamolodchikov c-function . We will precisely study the onset of the mass scale of the unstable particles and investigate how a particular coset flows to another one. On the other hand, we study in addition the flow of the operator content of one conformal field theory to another one by exploiting the flow provided by the $`\mathrm{\Delta }`$-sum rule of Delfino, Simonetti and Cardy .

Denoting by $`r`$ the radial distance and by $`t=\mathrm{ln}r^2`$ the renormalization group parameter, the functions $`c(t)`$ and $`\mathrm{\Delta }(t)`$ were defined in  and , respectively, obeying the differential equations

$`{\displaystyle \frac{dc(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{3}{4}}e^{2t}\mathrm{\Theta }(t)\mathrm{\Theta }(0)`$ (4)
$`{\displaystyle \frac{d\mathrm{\Delta }(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{𝒪(0)}}e^t\mathrm{\Theta }(t)𝒪(0).`$ (5)
The r.h.s. of these equations involve the two-point correlation functions of the trace of the energy-momentum tensor $`\mathrm{\Theta }`$ and an operator $`𝒪`$, which is a primary field in the sense of . In general these equations are integrated from $`t=\mathrm{}`$ to $`t=\mathrm{}`$ and one consequently compares the difference between the ultraviolet and the infrared fixed points. In order to exhibit the quantitative onset of the mass scale of the unstable particles we integrate these equations instead from some finite value $`t_0`$ to infinity. Restricting our attention to purely massive theories we use the fact that for those theories the infrared central charges are zero, such that

$$c(r_0)=\frac{3}{2}\underset{r_0}{\overset{\mathrm{}}{}}𝑑rr^3\mathrm{\Theta }(r)\mathrm{\Theta }(0).$$
(6)
Instead of the integral representation (6), the $`c`$-function is equivalently expressible in terms of a sum of correlators involving also other components of the energy momentum tensor . In deriving (4) these terms have been eliminate by means of the conservation law of the energy momentum tensor. We find (6) most convenient. The flow of $`c(r_0)`$ will surpass various steps: Starting with $`r_0=0`$ the theory will leave its ultraviolet fixed point and at a certain definite value, say $`r_0=r_u`$, the unstable particle will become massive such that $`c(r_0>r_u)`$ can be associated to a different conformal field theory. It appears natural to identify the mass $`M_{\stackrel{~}{c}}`$ as the point at which $`c(r_0)`$ is half the difference between the two coset values of $`c`$. As a consequence of (3) we may relate the masses of the unstable particles at different values of the resonance parameter $`\sigma `$, $`\sigma ^{}`$ and expect $`M_{\stackrel{~}{c}}(r_u,\sigma )=M_{\stackrel{~}{c}}(r_u^{},\sigma ^{})`$. We will employ the latter equality evaluated in the form (3) not only as a consistency requirement, but also as a confirmation of the fact that the renormalization group flow is indeed achieved by $`mr_0m`$. Increasing $`r_0`$ further, the energy scale of the stable particles will eventually be reached at, say at $`r_0=r_a,r_b,\mathrm{},r_n`$. Depending on the relative mass scales between the stable particles these points may coincide. Finally the flow will reach its infrared fixed point $`c(r_0=r_{ir})=0`$.

Likewise we can integrate equation (5)

$$\mathrm{\Delta }(r_0)=\frac{1}{2𝒪(0)}\underset{r_0}{\overset{\mathrm{}}{}}𝑑rr\mathrm{\Theta }(r)𝒪(0),$$
(7)
which allows to keep track of the manner the operator contents of the various conformal field theories are mapped into each other. We used that all conformal dimensions vanish in the infrared limit. Fortunately, we have $`\mathrm{\Theta }(r)𝒪(0)`$ $`𝒪(0)`$ in many applications such that the vacuum expectation value $`𝒪(0)`$ cancels often. One should note, however, that (7) is only applicable to those operators for which its two-point correlator with the trace of the energy momentum tensor is non-vanishing, such that one may not be in the position to investigate the flow of the entire operator content by means of (7).

In order to evaluate (6) and (7) we have to compute the two-point correlation functions in some way. In 2D-QFT this is probably most efficiently achieved, by expanding them in terms of n-particle form factors, i.e. the matrix elements of some local operator $`𝒪(\stackrel{}{x})`$ located at the origin between a multiparticle in-state and the vacuum denoted by $`0|𝒪(0)|V_{\mu _1}(\theta _1)V_{\mu _2}(\theta _2)\mathrm{}V_{\mu _n}(\theta _n)_{\text{in}}=:F_n^{𝒪|\mu _1\mathrm{}\mu _n}(\theta _1,\mathrm{},\theta _n)`$. Here the $`V_\mu (\theta )`$ are some vertex operators representing a particle of species $`\mu `$. Abbreviating the sum of the on-shell energies as $`E=_{i=1}^nm_{\mu _i}\mathrm{cosh}\theta _i`$, one may write

$`𝒪(r)𝒪^{}(0)={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{\mu _1\mathrm{}\mu _n}{}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d\theta _1\mathrm{}d\theta _n}{n!(2\pi )^n}}e^{rE}`$ (8)
$`\times F_n^{𝒪|\mu _1\mathrm{}\mu _n}(\theta _1,\mathrm{},\theta _n)\left(F_n^{𝒪^{}|\mu _1\mathrm{}\mu _n}(\theta _1,\mathrm{},\theta _n)\right)^{}.`$ (9)
Using this expansion we replace the correlation functions in the expression of the c-function $`c(r_0)`$ and the scaled conformal dimension $`\mathrm{\Delta }(r_0)`$ and perform the r integrations thereafter. Thus we obtain

$`c(r_0)=3{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{\mu _1\mathrm{}\mu _n}{}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d\theta _1\mathrm{}d\theta _n}{n!(2\pi )^n}}e^{r_0E}`$ (10)
$`\times \left|F_n^{\mathrm{\Theta }|\mu _1\mathrm{}\mu _n}(\theta _1,\mathrm{},\theta _n)\right|^2{\displaystyle \frac{(6+6r_0E+3r_0^2E^2+r_0^3E^3)}{2E^4}}`$ (11)
and

$`\mathrm{\Delta }(r_0)={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{\mu _1\mathrm{}\mu _n}{}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d\theta _1\mathrm{}d\theta _n}{n!(2\pi )^n}}{\displaystyle \frac{(1+r_0E)e^{r_0E}}{2E^2}}`$ (12)
$`\times F_n^{\mathrm{\Theta }|\mu _1\mathrm{}\mu _n}(\theta _1,\mathrm{},\theta _n)\left(F_n^{𝒪|\mu _1\mathrm{}\mu _n}(\theta _1,\mathrm{},\theta _n)\right)^{}.`$ (13)
We will now analyze (10), (12) and (3) for the $`SU(3)_2`$-HSG model. This model contains only two self-conjugate solitons which we denote by “$`+`$”, “$``$” and one unstable particle, which call $`\stackrel{~}{u}`$. The corresponding scattering matrix was found  to be $`S_{\pm \pm }=1`$, $`S_\pm (\theta )=\pm \mathrm{tanh}\left(\theta \pm \sigma i\pi /2\right)/2`$, which means the resonance pole is situated at $`\theta _R=\sigma i\pi /2`$. Stable bound states may not be formed. Note that for the corresponding value of $`\overline{\sigma }=\pi /2`$ and arbitrary $`\sigma `$ the condition $`M_{\stackrel{~}{u}}\mathrm{\Gamma }_{\stackrel{~}{u}}`$ is not fulfilled. However, as indicated above this condition only helps for a clearer identification of the mass parameter. For the HSG-models this condition starts to hold when the level is large, which indicates that in these type of models this interpretation is in fact a semi-classical one.

A huge class of form factors corresponding to various operators related to this model were constructed in . Labelling an operator by four quantum numbers $`\mu ,\nu ,\tau ,\tau ^{}`$ the general n-particle solution reads

$`F_{2s+\tau ,2t+\tau ^{}}^{𝒪_{\tau ,\tau ^{}}^{\mu ,\nu }|M^+M^{}}(\theta _1,\mathrm{},\theta _n)=H_{2s+\tau ,2t+\tau ^{}}^{𝒪_{\tau ,\tau ^{}}^{\mu ,\nu }|M^+M^{}}det𝒜_{2s+\tau ,2t+\tau ^{}}^{\mu ,\nu }`$ (14)
$`\left(\sigma _{2s+\tau }^+\right)^{st+\frac{\tau 1\nu }{2}}\left(\sigma _{2t+\tau ^{}}^{}\right)^{\frac{1+\tau \tau ^{}\mu }{2}t}{\displaystyle \underset{i<j}{}}\widehat{F}^{\mu _i\mu _j}(\theta _{ij}).`$ (15)
We used here a particular ordering by starting with $`2s+\tau `$ particles of the type $`\mu =+`$ followed by $`2s+\tau ^{}`$ particles of the type $`\mu =`$, collected in the sets $`M^\pm =\{\pm ,\mathrm{},\pm \}`$. Once these expressions are known, all other form factors related to it by permutations of the particles may be constructed trivially by exploiting Watson’s equations , see  for details concerning the HSG-models. The functions $`\widehat{F}^{\mu _i\mu _j}`$ for all combinations of the $`\mu `$’s are

$`\widehat{F}^{\pm \pm }(\theta )`$ $`=`$ $`i/2\mathrm{tanh}{\displaystyle \frac{\theta }{2}}\mathrm{exp}(\theta /2)`$ (16)
$`\widehat{F}^\pm (\theta )`$ $`=`$ $`2^{\frac{1}{4}}e^{{\scriptscriptstyle \frac{i\pi (11)\pm \theta }{4}}{\scriptscriptstyle \frac{G}{\pi }}_0^{\mathrm{}}{\scriptscriptstyle \frac{dt}{t}}{\scriptscriptstyle \frac{\mathrm{sin}^2\left((i\pi \theta \sigma )\frac{t}{2\pi }\right)}{\mathrm{sinh}t\mathrm{cosh}t/2}}},`$ (17)
with $`G=0.91597\mathrm{}`$ being the Catalan constant. The ($`t+s`$)$`\times `$($`t+s`$)-matrix

$$\left(𝒜_{2s+\tau ,2t+\tau ^{}}^{\mu ,\nu }\right)_{ij}=\{\genfrac{}{}{0pt}{}{\sigma _{2(ji)+\mu }^+\text{}1it}{\widehat{\sigma }_{2(ji)+2t+\nu }^{}\text{}t<is+t}$$
(18)
has as its entries elementary symmetric polynomials (see e.g.  for properties) depending on different sets of variables. We use the notation $`\sigma ^\pm `$ when they depend on the variable $`x=\mathrm{exp}\theta `$ associated to the sets $`M^\pm `$ and $`\widehat{\sigma }`$ to indicate that all variables are multiplied by a factor $`ie^\sigma `$. The overall constant was computed to

$`H_{2s+\tau ,2t+\tau ^{}}^{𝒪_{\tau ,\tau ^{}}^{\mu ,\nu }|M^+M^{}}`$ $`=`$ $`i^{s(2\tau +\tau ^{}+\nu +2)}2^{s(2s2t\tau ^{}1+2\tau )}`$ (20)
$`\times e^{s\sigma (2t+\tau ^{})/2}H_{\tau ,2t+\tau ^{}}^{𝒪_{\tau ,\tau ^{}}^{\mu ,\nu }},`$
where the value of $`H_{\tau ,2t+\tau ^{}}^{𝒪_{\tau ,\tau ^{}}^{\mu ,\nu }}`$ is fixed by the lowest non-vanishing form factor. In particular we need

$$F_{2s,2t}^\mathrm{\Theta }=\sigma _1(x_1,\mathrm{},x_n)\sigma _1(x_1^1,\mathrm{},x_n^1)F_{2s,2t}^{𝒪_{2,2}^{1,1}}.$$
(21)
Having assembled all the ingredients we can evaluate the expressions (10) and (12). We carry out the integrals by means of a Monte Carlo computation. For $`c(r_0)`$ we take contributions up to the 4-particle form factor into account and display our results in figure 1.

Figure 1: Renormalization group flow for the Virasoro central charge $`c(r_0)`$ for various values of the resonance parameter $`\sigma `$.

Following the renormalization group flow from the ultraviolet to the infrared, figure 1 illustrates the flow from the $`SU(3)_2/U(1)^2`$\- to the $`SU(2)_2/U(1)SU(2)_2/U(1)`$-coset when the unstable particle becomes massive. This confirms qualitatively the previous observation of the TBA analysis . Here we also want to compare the value of the mass of the unstable particle at different points of the resonance parameter $`\sigma `$ and $`t_0`$. Taking now the mass scales of the stable particles to be the same, i.e. $`m_+=m_{}=m`$, we compute the mass of the unstable particle according to (3), i.e. $`M_{\stackrel{~}{u}}(t_u,\sigma )m/\sqrt{2}\mathrm{exp}((\left|\sigma \right|+t_u)/2)`$. This means for different values of the resonance parameter we may still have the same value for the mass of the unstable particle when changing $`t_u`$, indeed we find

$$M_{\stackrel{~}{u}}(30.8,30)=M_{\stackrel{~}{u}}(20.8,20)=M_{\stackrel{~}{u}}(10.8,10).$$
(22)
Since the flow between the two cosets is smooth and takes place over some range of $`t_0`$, we had to select one particular point $`t_u`$. As already indicated in general, it is convenient to identify $`M_{\stackrel{~}{u}}`$ as the point at which $`c(t_0)`$ is half the difference between the two coset values of $`c`$. It is clear from figure 1, that since the overall shape of the curves between two values of $`c`$ is identical for different values of $`\sigma `$, any other value in the interval would lead to the same results in comparative considerations. This also means that when evaluating (22) the resulting value $`0.47m`$, which apparently violates the energetically necessary condition $`M_{\stackrel{~}{u}}>m_a+m_b`$, should not be taken too literally since the point $`t_u`$ is only chosen because it is easy to fix. Equations (22) confirm our general assertions outlined above.

For the evaluation of the scaled conformal dimension (12) we proceed similarly. For the solutions corresponding to the operators $`𝒪_{0,0}^{0,0}`$, $`𝒪_{0,2}^{0,1}`$ and $`𝒪_{2,0}^{1,0}`$, whose conformal dimension in the ultraviolet limit was identified  to be $`1/10`$, we take up to the 6-particle form factors into account. For the former two operators our results are presented in figures 2 and 3.

Figure 2: Renormalization group flow for the conformal dimension $`\mathrm{\Delta }(r_0)`$ of the operator $`𝒪_{0,0}^{0,0}`$ for various values of the resonance parameter $`\sigma `$.

We observe that the conformal dimension of the operator $`𝒪_{0,0}^{0,0}`$ flows to the value $`1/8`$, which is twice the conformal dimension of the disorder operator $`\mu `$ in the Ising model. The factor $`2`$ is expected from the mentioned coset structure, i.e. we find two copies of $`SU(2)_2/U(1)`$. The nature of the operator is also anticipated, since by construction $`F_n^{𝒪_{0,0}^{0,0}|M^+M^{}}`$ of the $`SU(3)_2`$-HSG model coincides precisely with $`F_n^\mu `$ of the thermally perturbed Ising model when one of the sets $`M^\pm `$ is empty. It is also clear that we could alternatively obtain (22) from the analysis of $`\mathrm{\Delta }(r_0)`$.

Figure 3: Renormalization group flow for the conformal dimension $`\mathrm{\Delta }(r_0)`$ of the operator $`𝒪_{0,2}^{0,1}`$ for various values of the resonance parameter $`\sigma `$.

Despite the fact that the explicit expressions for the form factors of $`𝒪_{0,2}^{0,1}`$ and $`𝒪_{2,0}^{1,0}`$ differ the values of $`\mathrm{\Delta }(r_0)`$ are hardly distinguishable and we therefore omit the plots for the latter case. We also note the previously observed fact , that the higher particle contributions for the latter operators are more important than for $`𝒪_{0,0}^{0,0}`$, which explains the fact that the starting point at the ultraviolet fixed point is not quite $`0.1`$. The operators also flow to the value $`1/8`$, such that the degeneracy of the $`SU(3)_2`$-HSG model disappears surjectively when the unstable particles become massive.

In comparison with other methods it would be extremely desirable to elaborate on the precise relationship between $`c(r_0)`$ and the finite size scaling function of the thermodynamic Bethe ansatz. Also the relation to the intriguing proposal in  of a renormalization group flow between Virasoro characters remains unclarified. The analogue of $`\mathrm{\Delta }(r_0)`$ still needs to be identified in the TBA as well as in the context of . In addition one may pose the question whether there exist higher dimensional counterparts of the function $`\mathrm{\Delta }(r_0)`$ in analogy to the results obtained in  for $`c(r_0)`$. Concerning the specific status of the HSG-models it remains a challenge to extend the results to other Lie groups .

Acknowledgments:  A.F. is grateful to the Deutsche Forschungsgemeinschaft (Sfb288) for financial support. O.A.C. thanks CICYT (AEN99-0589), DGICYT (PB96-0960), and the EC Commission (TMR grant FMRX-CT96-0012) for partial financial support and is also very grateful to the Institut für theoretische Physik of the Freie Universität for hospitality and for partial financial support. We are grateful to J.L. Miramontes, G. Mussardo for useful comments and A. Schilling for discussions on .