text:
# 1 Introduction

## 1 Introduction

Complex Ginzburg-Landau (CGL) equations is a class of universal models describing pattern formation in media of various physical nature, combining nonlinearity with dissipative and dispersive linear properties. They find important applications to traveling-wave convection (see Ref.  and references therein), fiber optics , and in other fields. Nonlinear terms in a CGL equation, as well the linear ones, may be both dissipative, accounting for loss and gain, and conservative.

An important class of patterns that appear in different versions of the CGL equations are solitary pulses (SPs). In particular, a GL model with the cubic-quintic (CQ) nonlinearity was first introduced by Petviashvili and Sergeev  in the two-dimensional (2D) case, with an objective to construct stable fully localized 2D pulses. The nonlinearity of the CQ type was necessary, as the stability of a pulse requires, first of all, the zero background to be stable, which implies that the linear part of the equation could not contain gain, hence a gain term had to be placed in the cubic nonlinearity. Finally, to provide for the overall stability of the model, the cubic gain had to be capped with a quintic lossy term.

Stable SPs (“solitons”) in 2D CGL equations of various types were found in Refs. . Stable axisymmetric solitons with an internal vorticity (“spin”) have been found too in the isotropic CQ CGL equation . On the other hand, the 2D CGL equation may find its most plausible experimental realization as a model governing evolution of optical spatiotemporal pulses in a planar nonlinear waveguide (2D optical cavity); in that case, however, the equation is strongly anisotropic due to the natural difference between the propagation direction ($`z`$) and the transverse one ($`x`$) along which the electromagnetic wave diffracts. In a rather general form, the corresponding anisotropic 2D model was put forward in Ref. :

$$iu_z+\frac{1}{2}u_{xx}(\frac{i}{2}+\beta )u_{tt}+\left[iu+\left(1i\gamma _1\right)|u|^2u+i\gamma _2|u|^4u\right]=0.$$
(1)
In this equation, $`u(z,x,t)`$ is the local amplitude of the electromagnetic field, $`tTz/V_0`$ is the so-called reduced time, where $`T`$ is time proper, and $`V_0`$ is the group velocity of the carrier wave propagating in the $`z`$–direction . Further, coefficients in front of the terms accounting for the transverse diffraction ($`\frac{1}{2}u_{xx}`$), linear loss ($`iu`$), optical filtering ($`\frac{1}{2}iu_{tt}`$, which formally looks like diffusion in the $`t`$–direction), and Kerr nonlinearity ($`|u|^2u`$) are all normalized to be $`1`$. The remaining free parameters $`\beta `$, $`\gamma _1`$, and $`\gamma _2`$ control the temporal dispersion, nonlinear cubic gain, and quintic loss ($`\gamma _1`$ and $`\gamma _2`$ must be positive, while $`\beta `$ may have either sign or be equal to zero). SP solutions to Eq. (1) in the form $`u(z,x,t)=\mathrm{exp}(ikz)`$$`v(x,t)`$, with the function $`v(x,t)`$ exponentially vanishing at $`\left|x\right|\mathrm{}`$ and $`\left|t\right|\mathrm{}`$, were found in Ref. , and their stability region was identified in the corresponding parameter space.

The model based on Eq. (1) assumes that the planar waveguide carrying the field is equipped with uniformly distributed linear bandwidth-limited gain (which induces the filtering term) combined with a saturable absorber , that together give rise to the CQ nonlinearity. For applications to nonlinear optics, it is also quite interesting to consider models where wave patterns are supported, instead of the intrinsic gain, by an external pump field, that leads to driven CGL equations, in which it is sufficient to take into regard the ordinary cubic (Kerr) nonlinearity, while the linear diffraction, dispersion, and filtering terms remain the same as in Eq. (1). In fact, two different driving terms may appear in the accordingly modified equation: a parametric drive and direct forcing, see below. It is quite interesting to find 2D spatiotemporal SPs in the driven optical models, as they have a real chance to be observed in the experiment, and are also of interest in their own right as completely localized pulses in the anisotropic 2D CGL equations of the new type.

In this work, we find these pulses and demonstrate their stability. Besides that, we also study quasi-1D solitons of two different types, which represent spatial and temporal solitons in the underlying optical models. It will be demonstrated that the spatial solitons are stable in a broad parametric region (although not everywhere); in particular, they coexist as stable attractors with the stable 2D pulses, hence both models are bistable in the corresponding regions. On the contrary to that, the temporal solitons are always modulationally unstable, and simulations show that the instability splits them into an array of 2D pulses, which subsequently coalesce into two 2D solitons. A uniform nonzero state in the parametrically driven model is also subject to modulational instability, which leads to a formation of a large number of 2D SPs; later, they merge into few final 2D solitons.

The rest of the paper is organized as follows. The two versions of the driven model are introduced in section 2, and SP solutions found in them are displayed in section 3. The results of extensive numerical simulations are summarized in the form of explicit stability diagrams in the corresponding three-dimensional parameter spaces. The paper is concluded by section 4.

## 2 The driven models

### 2.1 The parametric drive

In the 1D case, the parametrically driven CGL equation takes the form which generalizes the earlier studied parametrically driven damped nonlinear Schrödinger (NLS) equation (see Ref. , in which exact SP solutions to the latter equation were found, and their stability was investigated, and further results in Refs. ). The generalization to the 2D case of the same type as in Eq. (1), which describes a planar waveguide (optical cavity) with the Kerr nonlinearity, is straightforward:

$$iu_z+\frac{1}{2}u_{xx}(\frac{i}{2}+\beta )u_{tt}+\frac{1}{2}|u|^2u=(ki)u+\gamma u^{},$$
(2)
where the asterisk stands for the complex conjugation. The gain, which is accounted for by the last term in Eq. (2) with a real coefficient $`\gamma `$ is provided by a pump wave copropagating with the $`u`$–field at the double frequency. The real coefficient $`k`$ determines a wavenumber mismatch between the signal and pump waves.

Dropping the filtering term and assuming $`\beta <0`$ (which physically corresponds to anomalous dispersion of the electromagnetic waves in the waveguide ) transform Eq. (2) into a 2D isotropic damped parametrically driven NLS equation, which is used as a phenomenological model of the Faraday ripples in hydrodynamics  , in which case the evolutional variable $`z`$ in Eq. (2) is replaced by time $`T`$, and the former temporal variable $`t`$ is replaced by the second transverse coordinate $`y`$. A more general isotropic equation, which includes the diffusion term,

$$iu_z+(\frac{1}{2}iD)^2u+\frac{1}{2}|u|^2u=(ki)u+\gamma u^{},D>0,$$
(3)
is also used as a model to describe localized 2D states, the so-called oscillons, observed in experiments with the Faraday resonance . Getting back to the damped isotropic NLS equation with the parametric drive, we note that this equation (with the evolutional variable $`z`$) finds applications to the description of spatial patterns in 3D driven optical cavities . In a recent preprint , it has been shown that axisymmetric solitons may be stable in the latter model, despite the presence of the collapse in its undamped undriven version (which is the NLS equation proper).

In the present work, the objective is to find stable SPs in the strongly anisotropic full CGL model (2), which will have the meaning of spatiotemporal optical solitons  propagating in the planar waveguide. It seems obvious that pulses may only exist if the gain parameter exceeds the loss coefficient, i.e., $`\gamma >1`$. On the other hand, an obvious necessary condition for the stability of pulses is the stability of the zero background against perturbations of the form

$$u=a\mathrm{exp}\left(\sigma z+iqxi\omega t\right)+b\mathrm{exp}\left(\sigma ^{}ziqx+i\omega t\right)$$
(4)
with infinitesimal amplitudes $`a`$ and $`b`$, where $`q`$ and $`\omega `$ are arbitrary real wavenumber and frequency of the perturbation. A dispersion relation which determines the (complex) instability growth rate $`\sigma `$ as a function of $`q`$ and $`\omega `$ can be easily found from the linearized equation (2):

$$\left(\frac{1}{2}q^2\beta \omega ^2+k\right)^2+\left[\left(\mathrm{Re}\sigma +1+\frac{1}{2}\omega ^2\right)+i\mathrm{Im}\sigma \right]^2=\gamma ^2,$$
(5)
which gives rise to two different branches of the dispersion relation, viz., $`\mathrm{Re}\sigma =\left(1+\frac{1}{2}\omega ^2\right)`$ with $`\mathrm{Im}\sigma 0`$, and the one with $`\mathrm{Im}\sigma =0`$ and $`\sigma `$ determined by the equation

$$\left(\sigma +1+\frac{1}{2}\omega ^2\right)^2=\gamma ^2\left(\frac{1}{2}q^2\beta \omega ^2+k\right)^2.$$
(6)
Obviously, the former branch always satisfies the stability condition, $`\mathrm{Re}\sigma (q,\omega )0`$. Straightforward analysis of Eq. (6), with regard to the above assumption $`\gamma >1`$, shows that, if $`\omega =0`$, the stability condition for the latter branch amounts to the well-known inequalities  $`\gamma ^21+k^2`$, $`k>0`$.

The consideration of Eq. (6) with $`\omega 0`$ demonstrates that, in this case, the most dangerous perturbations are those with $`q=0`$. If $`\beta 0`$, they do not generate any additional instability, and if $`\beta `$ is positive, a new stability condition appears, $`\gamma ^2(k+2\beta )^2/(1+4\beta ^2)`$; it is easy to check that the above-mentioned necessary condition $`\gamma ^21+k^2`$ is a straightforward corollary of this inequality. Thus, the eventual system of the conditions providing for the stability of the zero solution to Eq. (2) (provided that the system is infinitely extended in both $`x`$– and $`t`$–directions) is $`k0`$ and

$$\begin{array}{cc}\gamma ^21+k^2,\hfill & \mathrm{if}\beta 0;\hfill \\ \gamma ^2(k+2\beta )^2/(1+4\beta ^2),\hfill & \mathrm{if}\beta >0.\hfill \end{array}$$
(7)
It is also relevant to consider nonzero spatially uniform solutions to Eq. ( 2), which have the form

$$u=A_0\mathrm{exp}\left(i\varphi \right),\varphi =\frac{1}{2}\mathrm{sin}^1(1/\gamma ),A_0=\sqrt{2\left[k+\gamma \mathrm{cos}(2\varphi )\right]}.$$
(8)
This uniform solution exists for $`\gamma >1`$. The consideration of the modulational stability of the solution against small perturbations with the wavenumber $`q`$ and instability growth rate $`\sigma `$ (cf. Eq. (4)) yields a dispersion relation

$$\left(\frac{1}{2}q^2\beta \omega ^2+kA_0^2\right)^2+\left[\left(\mathrm{Re}\sigma +1+\frac{1}{2}\omega ^2\right)+i\mathrm{Im}\sigma \right]^2=\gamma ^2\gamma A_0^2\mathrm{cos}(2\varphi )+A_0^4/4.$$

As well as Eq. (5), this equation has a solution branch with $`\mathrm{Im}\sigma =0`$, then

$$\sigma =1\omega ^2/2\pm \sqrt{\gamma ^2\gamma A_0^2\mathrm{cos}(2\varphi )+A_0^4/4(q^2/2\beta \omega ^2+kA_0^2)^2}.$$
(9)
The maximum value of this expression, which is attained at $`\omega =0`$ and $`q=\sqrt{2(kA_0^2)}`$, is always positive, $`\sigma _{\mathrm{max}}=1+\sqrt{1+k^2}`$. Therefore, the uniform state (8) is always unstable against modulations in the $`x`$–direction, the same way as in the 1D damped parametrically driven NLS equation .

### 2.2 The directly forced model

The light signal in a lossy waveguide can also be directly (rather than parametrically) supported by a pump wave launched at the same (rather than double) frequency. For spatial patterns in 3D lossy optical cavities (a typical example is a domain wall ), this gives rise to the well-known isotropic directly forced 2D CGL equation, with the propagation coordinate $`z`$ playing the role of the evolutional variable, see a review  (a similar model describes a cavity with a trapped Bose-Einstein condensate in the presence of source and sink of atoms ). The optical patterns generated this way may find an application in the design of optical memory , although fully localized pulses may be subject to collapse under the action of the self-focusing Kerr nonlinearity .

The governing equation for the spatiotemporal optical signal in the directly forced nonlinear waveguide differs from Eq. (2) by the driving term:

$$iu_z+\frac{1}{2}u_{xx}(\frac{i}{2}+\beta )u_{tt}+\frac{1}{2}|u|^2u=(ki)uf,$$
(10)
where $`f`$ is the amplitude of the forcing field. Obviously, this equation can only produce localized patterns on top of a constant background $`uu_0`$, whose amplitude is determined by a cubic equation following from Eq. (10), $`(1/2)\left|u_0\right|^2u_0+\left(ik\right)u_0+f=0`$. The background can be subject to modulational instability ; however, if the amplitude $`u_0`$ is small enough, the instability will be so weak that it may be disregarded for finite values of the propagation distance that are relevant to the experiment.

## 3 Numerical results

### 3.1 Two-dimensional pulses in the parametrically driven model

Equations (2) and (10) were solved numerically by means of a pseudospectral code with $`256\times 256`$ modes and periodic boundary conditions in $`x`$ and $`t`$, with a fixed size of the system in both directions, $`L_x=L_t=20`$. As well as in Ref. , a stable 2D SP was generated for the first time, starting with an initial Gaussian pulse localized in both $`x`$ and $`t`$, with its center placed at the center of the integration domain. In order to generate stability diagrams (see below), we then varied the parameters by small steps, the initial configuration for each simulation being the stable pulse produced by the simulation at the previous step.

In Fig. 1 we display a typical example of a stable 2D completely localized pulse, which was produced by the numerical solution of Eq. (2), following the procedure outlined above. Results of extensive simulations, which were carried out for many different values of the parameters, can be presented in a compact form as a stability diagram for the 2D pulses in the three-dimensional parametric space ($`\gamma ,k,\beta `$), which is displayed in Fig. 2. We stress that this parameter space comprises all the coefficients of Eq. (2).

The stability region proper is located under the “roof” in Fig. 2. To the left of the stability border which is shown by crosses, no localized 2D pulses can be created: in this case, an initial 2D pulse configuration quickly decays into the stable zero state. On the other side, i.e., to the right of the border denoted by rhombuses, localized 2D pulses were not found either. However, the result of the evolution is different in this case: the initial 2D pulse expands along the $`t`$–axis to generate a quasi-1D pulse, which is localized in the $`x`$–direction, and delocalized in the $`t`$ –direction. Such a quasi-1D solution to Eq. (2) can be found in an exact form,

$$u(x,t)=\sqrt{2}\kappa e^{i\varphi }\mathrm{sech}(\kappa x),\varphi =\frac{1}{2}\mathrm{sin}^1(1/\gamma ),\kappa =\sqrt{2\left[k+\gamma \mathrm{cos}(2\varphi )\right]}.$$
(11)
In terms of the underlying optical models, the solution (11) corresponds to a spatial soliton. A typical example of such a soliton is shown in Fig. 3. Straightforward comparison shows that the numerically found quasi-1D pulse completely coincides with the analytical solution (11). Simulations demonstrate that this spatial soliton is stable everywhere in the stability region of the 2D pulses shown in Fig. 2, i.e., the model is bistable in this region. Beyond the right stability border of the 2D pulses, shown by rhombuses in Fig. 2, i.e., in the region where only the spatial solitons are found, they remain stable.

For comparison with these results, we have also performed direct simulations of the isotropic equation (3). Stable 2D localized pulses exist in this model too, Fig. 4 displaying an example for $`D=0.1`$. For $`D=0.1`$ and $`k=1.3`$, the stable axially symmetric 2D pulses have been found in the interval $`1.20<\gamma <1.57`$. For $`\gamma <1.20`$, the initial 2D pulse decays to the zero state, while for $`\gamma >1.57`$ the initial 2D pulse evolves to the spatially uniform state (8). In this connection, we note that a critical branch of the modulational-instability growth rate for this uniform solution to Eq. (3) can be found, as a function of the perturbation wavenumber $`q`$, in the form $`\sigma (q)=1Dq^2+\sqrt{1+k^2(q^2/2+kA_0)^2}`$, cf. Eq. (9). It is easy to check that, for $`D=0.1`$ and $`k=1.3`$, $`\sigma (q)`$ is negative for any $`q`$, hence the uniform state is stable.

### 3.2 Modulational instability of quasi-one-dimensional pulses in the parametrically driven model

In the context of the results presented above, the spatial solitons (quasi-1D pulses) were stable. However, they become unstable against $`t`$ –dependent perturbations for sufficiently large negative values of the dispersion coefficient $`\beta `$ in Eq. (2), which corresponds to anomalous dispersion in terms of nonlinear optics (in fact, modulational instability of spatial solitons against temporal perturbations in a system with anomalous dispersion is a general effect which has been known for a long time ). In particular, for $`k=1.3`$ and $`\gamma =1.3`$, this instability occurs if $`\beta <1.08`$. Figure 5 displays a result of numerical simulations in this case. The initial condition including a small perturbation is taken as

$$u(x,t,z=0)=\sqrt{2}\kappa \left[\left(1+0.05\mathrm{cos}\left(\frac{8\pi t}{L}\right)+0.01\mathrm{cos}\left(\frac{2\pi t}{L}\right)\right)\mathrm{cos}\varphi +i\mathrm{sin}\varphi \right]\mathrm{sech}\left(\kappa \left(x\frac{L}{2}\right)\right),$$
(12)
where $`\kappa `$ and $`\varphi `$ are the same constants as in Eq. (11). Figure 5(a) displays the evolution of $`|u(x=L/2,t)|`$ at $`\beta =2.5`$. While the pattern remains localized in the $`x`$–direction during the evolution, it is obvious in Fig. 5(b) that the growth of the modulation splits the original spatial soliton into four 2D spatiotemporal pulses, but three of them merge together, thus two 2D pulses survive (splitting of a spatial solitary beam into temporally localized pulses under the action of modulational instability was earlier observed in numerical  and laboratory  experiments with second-harmonic-generating optical systems). Figure 5(b) displays an eventual established state including the two 2D pulses. This modulational instability is a generic type of transverse instability for quasi-1D pulses (another generic type is the zigzag instability found in Ref.  in the context of the 2D GL equation with the CQ nonlinearity).

Another type of quasi-1D pulse solution to Eq. (2), which is localized in the $`t`$–direction and uniform in the $`x`$–direction, is also possible. Such an $`x`$–independent solution, $`u(z,t,z)=U(z,t)`$, represents a temporal soliton in the optical cavity, and satisfies the equation

$$iU_z(\frac{i}{2}+\beta )U_{tt}+\frac{1}{2}|U|^2U=(ki)U+\gamma U^{}$$
(13)
(an exact analytical solution to this equation is not available, unless $`\beta =0`$). For the same case, $`k=1.3`$ and $`\gamma =1.3`$, as considered above, this type of the quasi-1D pulse solution exists and is stable against $`t`$–dependent perturbations if $`\beta <0.812`$. Figure 6(a) displays such a pulse at $`\beta =1,k=1.3`$ and $`\gamma =1.3`$, with its center at the point $`t=L/2`$. To study its stability against $`x`$–dependent perturbations, we took an initial condition

$$u(x,t,z=0)=\left[1+0.05\mathrm{cos}(2\pi x/L)\right]\mathrm{Re}U(t)+i\mathrm{Im}U(t),$$

where $`U(t)`$ is the unperturbed $`z`$–independent SP solution to Eq. (13 ). Modulation in the $`x`$– direction grows and splits the temporal soliton into many 2D spatiotemporal pulses, as it is evident in Fig. 6(b), which displays a 3D plot of $`|u(x,t,z)|`$ at $`z=200`$. Such a temporal soliton, extended along the $`x`$–direction, seems to be always unstable against the $`x`$–dependent perturbations, breaking up into a large number of 2D pulses.

One can also look for more general quasi-1D pulses of the form $`u(x,t,z)=U(xct)`$, which may be interpreted as either moving spatial solitons, or oblique temporal ones, the function $`U`$ satisfying an equation

$$\left(12c^2\beta ic^2\right)U^{\prime \prime }+|U|^2U=2(ki)U+2\gamma U^{}.$$

A detailed study of such general quasi-1D solitons is beyond the scope of this work, but one may assume that there is a critical value of the velocity $`c`$ separating completely unstable solutions and those which may be stable.

We have also simulated the development of the modulational instability of the spatially uniform state (8) in the parametrically driven model (2) (the presence of this instability was demonstrated analytically in the previous section). Formation of many 2D pulses may be expected, at sufficiently small $`\beta `$, as a result of the instability development. As an initial condition, we took the uniform state (8) randomly perturbed by an initial small disturbance. Figure 7 displays a snapshot of the pattern obtained at $`z=600`$ for $`\beta =2.5,k=1.3`$ and $`\gamma =1.3`$. At first, many 2D pulses are created; however, they subsequently merge, so that only a few 2D pulses eventually survive, which is obvious in Fig. 7.

### 3.3 The directly forced model

Numerical simulations of the directly driven CGL equation (10) also readily produce stable 2D pulses (built on top of the small-amplitude background, see above). A typical example of such a pulse is shown in Fig. 8. As well as in the parametrically driven model, in this case extensive simulations of Eq. (10) have made it possible to identify a stability region for 2D SPs in the full three-dimensional parameter space of the model, $`(f,k,\beta )`$. This region is located under the “roof” in the parameter space, see Fig. 9. Note that the parameter space in Fig. 9, as well as that in Fig. 2, incorporates all the parameters of the model.

The two boundaries of the stability region in Fig. 9 have the meaning similar to that of the stability boundaries shown in Fig. 2 for the parametrically driven model: to the left of the boundary shown by crosses, any initial pulse decays to zero, and to the right of the second boundary, shown by rhombuses, the fully localized 2D pulse evolves to a quasi-1D pulse oriented along the $`t`$–axis (the spatial soliton, in terms of the underlying optical models). As well as in the parametrically driven model, in the directly forced one the spatial solitons coexist as stable patterns with the 2D pulses (hence, the model is bistable in that region), and they remain stable across the right border of the stability region of the 2D solitons, i.e., across the surface shown by rhombuses in Fig. 9.

## 4 Conclusion

In this work, we have introduced two different two-dimensional models of the Ginzburg-Landau type with the Kerr (cubic) nonlinearity. The corresponding equations are driven either parametrically or directly by an external field. The models describe spatio-temporal dynamics in nonlinear dispersive lossy optical cavities, where the loss also includes a filtering term. Both models are strongly anisotropic, featuring only diffraction in one (spatial) direction, and a combination of effective diffusion and dispersion in the other (temporal) direction. By means of direct simulations, we have found stable two-dimensional solitary pulses in both models (in the directly driven one, the pulses are built on top of a nonzero flat background). These two-dimensional solitary pulses correspond to spatio-temporal solitons (“light bullets”) supported by the driving field in the planar nonlinear waveguide.

The main results were presented in a compact form as stability regions in the full three-dimensional parameter space of either model. The stability region is bounded by two surfaces. Beyond one of them, any initial pulse decays to zero, while stable quasi-one-dimensional pulses, extending along the temporal direction (i.e., spatial solitons, in terms of the optical cavity), were found beyond the other boundary. The spatial solitons are also stable inside the stability regions of the two-dimensional pulses, so that both models demonstrate bistability in these regions. Unlike the spatial solitons, their quasi-one-dimensional counterparts in the form of temporal solitons are always modulationally unstable. It was demonstrated that the instability splits temporal solitons into an array of spatiotemporal two-dimensional pulses, which, in turn, coalesce into two final pulses. Similarly, the instability of the uniform nonzero state in the parametrically driven model leads to formation of many two-dimensional pulses, which subsequently merge into few two-dimensional solitons.

## Figure Captions

Fig. 1. An example of the stable stationary two-dimensional solitary pulse, found as a solution of the parametrically driven complex Ginzburg-Landau equation (2). The field $`\left|u\right|`$ is shown vs. the temporal and transverse spatial coordinates $`t`$ and $`x`$. Values of the parameters are $`k=1.3,\beta =0.3,\gamma =1.30`$.
Fig. 2. The region of stable two-dimensional solitary pulses in the full parameter space of Eq. (2). The stability region is located under the “roof”. To the left of stability region, only zero solution is possible, and to the right of it, the 2D pulse evolves to a stable quasi-one-dimensional pattern, see an example in Fig. 3.
Fig. 3. An example of a stable quasi-one-dimensional pulse found in the parametrically driven complex Ginzburg-Landau equation (2). Values of the parameters are $`k=1.3,\beta =0.3,\gamma =1.36`$ (a point corresponding to these values is located slightly to the right of the stability region shown in Fig. 2).
Fig. 4 An example of a stable axisymmetric two-dimensional solitary pulse found as a solution of the parametrically driven isotropic complex Ginzburg-Landau equation (3). The field $`\left|u\right|`$ is shown vs. the temporal and transverse spatial coordinates. Values of the parameters are $`k=1.3,D=0.1`$ and $`\gamma =1.30`$.
Fig. 5 (a) Evolution of $`|u(L/2,t,z)|`$ for the parametrically driven complex Ginzburg-Landau equation (2). The initial state is a 1D localized pattern (12) (spatial soliton) including a small modulational perturbation, values of the parameters being $`k=1.3,\beta =2.5,\gamma =1.3`$ . (b) The established stationary pattern including two 2D pulses.
Fig. 6 (a) The quasi-one-dimensional pulse, which is localized in the $`t`$ –direction around $`t=L/2`$ (temporal soliton), at $`\beta =1,k=1.3`$ and $`\gamma =1.3`$, in the parametrically driven complex Ginzburg-Landau equation ( 2). (b) The result of the instability development of this soliton, shown as a plot of $`|u(x,t)|`$ at $`z=200`$.
Fig. 7 Patterns produced by the modulational instability of the uniform state in the parametrically driven complex Ginzburg-Landau equation (2) at $`z=600`$. Values of parameter are $`\beta =2.5,k=1.3`$ and $`\gamma =1.3`$.
Fig. 8. An example of the stable stationary two-dimensional solitary pulse, found on top of the nonzero flat background, in the directly forced complex Ginzburg-Landau equation (10). The field $`\left|u\right|`$ is shown vs. the temporal and transverse spatial coordinates $`t`$ and $`x`$. Values of the parameters are $`k=3,\beta =0.3,f=3.1`$.
Fig. 9. The region of stable two-dimensional solitary pulses (built on top of the nonzero flat background) in the full parameter space of the directly forced complex Ginzburg-Landau equation (10). The stability region is located under the “roof”. To the left of stability region, the 2D pulse decays to the flat-background solution, and to the right of the stability region, the 2D pulse evolves to a stable quasi-one-dimensional pattern (spatial soliton).


format: text

text:
# The ac magnetic response of mesoscopic type II superconductors

## I INTRODUCTION

The response of superconductors to an ac magnetic field has been of interest for a long time,pippard ; clem0 ; sridhar2  and particularly in the last years. clem ; brandt ; geshk ; beek ; sonin ; sridhar1 ; maeda1 ; mallozzi ; portis ; gomory ; maeda2  The microwave surface impedance $`Z_s=R_siX_s`$ and the ac magnetic susceptibility $`\stackrel{~}{\chi }=\chi ^{}+i\chi ^{\prime \prime }`$ have been extensively studied. pippard ; clem0 ; sridhar2 ; clem ; brandt ; geshk ; beek ; portis ; sonin ; sridhar1 ; maeda1 ; mallozzi ; gomory ; maeda2  Interesting behavior has been found for different values of the frequency, the magnetic field, the pinning force and the thermal fluctuations,clem ; brandt ; geshk ; beek  whereas linear and nonlinear response appear depending on the strength of the ac signal.beek  The main interest so far has been on the electrodynamics of macroscopic samples containing a large number of vortices. In this case, it is possible to use phenomenological models to describe the macroscopic behavior of the superconducting samples.clem ; brandt ; geshk ; beek

Recently there has been an increasing interest in the study of vortex physics on a mesoscopic scale, buisson ; bolech ; bolle ; geim98 ; geim00 ; peeters97 ; peeters98a ; peeters99b ; peeters00 ; peeters99c ; peeters01 ; palacios ; bruyndoncx ; chibotaru ; HD  a regime in which a small number of vortices are confined in a small sample. The behavior of mesoscopic superconductors is different from the behavior of bulk samples. In mesoscopic superconductors surface effectsbean  are very important since the interaction of vortices with the surface currents is large.HD  The magnetic properties strongly depend on the sample sizes and geometry.peeters99b ; peeters00 ; HD  The most studied geometries are the the mesoscopic disk buisson ; geim98 ; geim00 ; peeters97 ; peeters98a ; peeters99b ; peeters00 ; peeters99c ; peeters01 ; palacios  the mesoscopic slabbolech ; bolle  and the mesoscopic square. bruyndoncx ; chibotaru ; HD  In particular the mesoscopic samples can develop Abrikosov multivortex states peeters98a  and depending on the size of the sample it is possible to observe first or second order transitions.peeters97  The electric charge of vortices, peeters01  a paramagnetic effect geim98 ; palacios  and the surface barrier HD ; bean  on a mesoscopic scale have also been studied. One interesting characteristic of the magnetic properties of mesoscopic superconductors is the behavior of the dc magnetization curves. In a mesoscopic scale, the vortices that are inside the sample induce a reinforcement of the surface barrier at fields greater than the first penetration field $`H_p`$.HD  These new barriers allow for the existence of metastable states of constant vorticity as a function of magnetic field. Each metastable state becomes unstable at the $`i`$-th penetration field $`H_{pi}`$ in which vortices enter the sample and the magnetization has a discontinuous jump. These jumps have been observed numerically bolech ; peeters97 ; HD  and experimentally in mesoscopic Al disks.peeters97

In the last years numerical simulations of the time dependent Ginzburg Landau (TDGL) equations have been an important tool in the study of the static and dynamic properties of superconductors. bolech ; frah ; liu ; kato ; machida ; aranson1 ; aranson2 ; aranson3 ; enomoto ; alvarez ; gropp  The TDGL model has been used to study the flux growth dynamics frah ; liu , the magnetic responsebolech ; kato  and the I-V characteristics of superconducting samples.machida  In particular, Enomoto et al.enomoto  using the TDGL equations studied the temperature dependence of the ac susceptibility at a fixed frequency in the absence of a dc magnetic field for a large sample. For increasing temperature they found that $`\chi ^{}`$ exhibits a step-like change from a negative value ($`\chi ^{}=1/(4\pi )`$ to zero, while $`\chi \mathrm{"}`$ initially rises from zero, goes through a maximum, and then returns to a small value near $`T_c`$, which is qualitatively consistent with the overall behavior observed in experiments.

In this paper we perform numerical simulations of the time dependent Ginzburg-Landau equation to study the ac magnetic response of type II superconductors on a mesoscopic scale. We study the frequency and $`H_{dc}`$ dependence of the linear ac susceptibility $`\stackrel{~}{\chi }(H_{dc},\omega )=\chi ^{}(H_{dc},\omega )+i\chi ^{\prime \prime }(H_{dc},\omega )`$ of square samples with dimensions of the order of the London penetration depth $`\lambda _L`$. Our resuls were obtained in the absence of bulk pinning and in the linear regime.

The paper is organized as follows. In Sec.II we review the known results for the ac magnetic response of macroscopic superconductors, in order to later compare them with our results in mesoscopic samples. In particular, the two fluid model for the Meissner state and the Coffey and Clem model for the mixed state are reviewed. In Sec.III we describe the TDGL equations used in the simulations. In Sec.IV we present our numerical results for the ac magnetic response of mesoscopic square samples. Sec.IV.A is devoted to the study of the frequency dependence of $`\chi ^{}`$ and $`\chi \mathrm{"}`$. We find that for $`H_{dc}=0`$ the behavior of the sample is well described by the two fluid model. At $`\omega =\nu _o=c^2/(4\pi \sigma \lambda ^2)`$ a maximum in the imaginary part of $`\stackrel{~}{\chi }(\omega ,0)`$ appears. In the Meissner state for $`H_{dc}H_p`$, another maximum appears in $`\chi \mathrm{"}(\omega )`$. These behavior, characterized by the presence of two dissipation maxima is also observed in the mixed state. At the same time, for increasing frequency $`\chi ^{}(\omega ,H_{dc})`$ changes from London screening to a perfect diamagnetic state at high frequency. This transition is due to a decrease of the ac penetration length. In Sec.IV.B we study the $`H_{dc}`$ dependence of $`\chi ^{}`$ and $`\chi \mathrm{"}`$ at a fixed frequency. Depending on the sample size, we find well defined mesoscopic and macroscopic behaviors. For mesoscopic samples we observe that, in the magnetic field ranges in which the number of vortices is constant, $`\chi \mathrm{"}(H_{dc})`$ continuously increases for increasing $`H_{dc}`$. At the $`i`$-th vortex penetration field $`H_{pi}`$, the entrance of vortices produces a considerable suppression of the ac losses and $`\chi \mathrm{"}(H_{dc})`$ decreases with a discontinuous jump. On the contrary, samples of macroscopic sizes show a continuous behavior in $`\stackrel{~}{\chi }(H_{dc})`$, with $`\chi \mathrm{"}(H_{dc})`$ monotonically increasing with $`H_{dc}`$. In Sec. V we study the time evolution of the order parameter in mesoscopic superconductors. We show that at high frequencies and in mesoscopic samples, the vortices are fixed and only play a secondary role in the ac response of the sample. We show that the main dissipation mechanism in mesoscopic superconductors is due to the effect of “nascent vortices”. They cause large variations of the amplitude of the order parameter at the boundary of the sample before the entrance of vortices. Finally in Sec.VI we give a summary of our results and conclusions.

## II THE AC RESPONSE OF MACROSCOPIC SAMPLES: REVIEW OF KNOWN RESULTS

### II.1 Two fluid model and Coffey and Clem model

When an ac magnetic field is applied $`H(t)=H_{dc}+h_{ac}\mathrm{cos}(\omega t)`$, an effective complex penetration depth $`\stackrel{~}{\lambda }(\omega )`$ can be defined assuming that the fluctuating field inside a semiinfinite sample for $`x>0`$ has the form

$$\delta H=h_{ac}e^{x/\stackrel{~}{\lambda }(\omega )}e^{i\omega t}$$
(1)
or from a generalization of London’s expression for the current,

$$\times 𝐉(\omega )=\frac{c}{4\pi \stackrel{~}{\lambda }^2(\omega )}𝐁(\omega ).$$
(2)
It can be related to the frequency dependent conductivity as $`\stackrel{~}{\lambda }(\omega )=\sqrt{\frac{ic^2}{4\pi \omega \stackrel{~}{\sigma }(\omega )}}`$. In the case of a normal metal with conductivity $`\sigma _n`$, the complex penetration depth is directly related to the skin depth $`\delta _n=\sqrt{\frac{c^2}{2\pi \omega \sigma _n}}`$ as $`\stackrel{~}{\lambda }_n=(1+i)\delta _n/2`$. The complex permeability $`\stackrel{~}{\mu }(\omega )`$, or the complex susceptibility $`\stackrel{~}{\chi }(\omega )=\chi ^{}(\omega )+i\chi ^{\prime \prime }(\omega )=[\stackrel{~}{\mu }(\omega )1]/4\pi `$, depend both on $`\stackrel{~}{\lambda }(\omega )`$ and also on the specific shape of the sample, since the ac penetration of magnetic fields is a surface phenomenon. In the case of a square sample of size $`L\times L`$ an approximate expression for the permeability in the limit $`\stackrel{~}{\lambda }L`$ is

$$\stackrel{~}{\mu }_{square}\frac{4\stackrel{~}{\lambda }}{L}\left(1\frac{\stackrel{~}{\lambda }}{L}\right).$$
(3)
Exact expressions for $`\stackrel{~}{\mu }`$ for squares, cylinders and slabs can be found for example in Ref.clem, . Sometimes the complex surface impedance $`Z_s`$ is used, it can be defined when $`\stackrel{~}{\lambda }L`$ as $`Z_s=R_siX_s=i\omega \frac{4\pi }{c^2}\stackrel{~}{\lambda }(\omega ).`$

In the Meissner state the ac response has been usually described with the two fluid model.twofluid ; gittleman  In the mixed state, the magnetic flux enters in the form of quantized vortices. Therefore the dynamics of the vortices has to be included in a description of the magnetic ac response of the mixed state. clem ; brandt ; geshk ; beek  In particular, Coffey and Clem have extended the two fluid model by including the equation of motion of vortices for small displacements. clem  Similar results were obtained by Brandtbrandt  and by van der Beek et al.beek  All these models assume that there is a large number of vortices in the sample and a continuous description of the vortex lattice is used (valid for $`BH_{c1}`$).

The two fluid model consists in writing the total current as the sum of the supercurrent and the normal current,

$$𝐉=𝐉_s+𝐉_n$$
(4)
with

$$𝐉_n=\sigma _n𝐄=\frac{\sigma _n}{c}\frac{𝐀}{t}$$
(5)
and the supercurrent given by the London model

$$\times 𝐉_s=\frac{c}{4\pi \lambda _L^2}(𝐁𝐛_v).$$
(6)
Here $`\lambda _L`$ is the static London penetration depth and $`𝐛_v`$ is the local vortex magnetic field. In the absence of vortices, $`𝐛_v=0`$, the two fluid model has the characteristic time $`t_0=4\pi \lambda _L^2\sigma _n/c^2`$ for the transformation of supercurrents into normal currents, since $`𝐉_n=t_0𝐉_s/t`$. The Coffey and Clem model includes the equation of motion of vortices for small displacements $`𝐮(𝐱,t)`$ from their equilibrium positions,

$$\eta _v\dot{𝐮}+\kappa _p𝐮=\frac{1}{c}𝐉\times \mathrm{\Phi }_0\widehat{𝜶}$$
(7)
where $`\eta _v`$ is the viscous drag coefficient, $`\kappa _p`$ is the restoring force constant (Labusch parameter) of a pinning potential well and $`\widehat{𝜶}`$ the local vortex direction. Here we have neglected thermal fluctuation effects. The local vortex magnetic field depends on the vortex displacements as $`𝐛_v=\times (𝐮\times 𝐁)`$.beek

In the case of fields $`𝐇𝐡_{ac}`$ and parallel to the surface, using Eqs.(2),(4)-(7) the effective complex penetration depthclem  can be obtained as

$$\stackrel{~}{\lambda }^2=\lambda _L^2\frac{1+\left(\frac{\lambda _L^2}{\lambda _C^2}i\omega t_{ff}\right)^1}{1i\omega t_0}=\frac{\lambda _L^2+\left(\frac{1}{\lambda _C^2}i\frac{2}{\delta _{ff}^2}\right)^1}{1i\frac{2\lambda _L^2}{\delta _n^2}}$$
(8)
with the Campbell penetrationd depth $`\lambda _C^2=B\mathrm{\Phi }_0/4\pi \kappa _p`$, the flux-flow time scale $`t_{ff}=4\pi \lambda _L^2\sigma _{ff}/c^2`$, the flux-flow conductivity $`\sigma _{ff}=c^2\eta _v/B\mathrm{\Phi }_0`$ and the flux-flow skin depth $`\delta _{ff}^2=c^2/2\pi \omega \sigma _{ff}`$. Since $`\sigma _{ff}\sigma _nH_{c2}/B`$, we have $`\sigma _{ff}>\sigma _n`$ and $`t_{ff}>t_0`$. Therefore the highest characteristic frequency is $`\nu _0=1/t_0`$, which for conventional superconductors is about $`\nu _010100`$ GHz.tscales

### II.2 Meissner state

In the absence of vortices the complex penetration depth is simply $`\stackrel{~}{\lambda }^2=\lambda _L^2/(1i\omega t_0)`$. This leads to a dissipation peak in $`\chi ^{\prime \prime }`$ when $`\omega t_01`$. For large frequencies, $`\omega t_01`$ (i.e. $`\delta _n\lambda _L`$), the system behaves as a normal metal (normal state skin depth effect) with:

$$\lambda ^{}=\lambda ^{\prime \prime }\lambda _L\sqrt{\frac{1}{2\omega t_0}}$$
(9)
Therefore dissipation goes as $`\chi ^{\prime \prime }\lambda ^{\prime \prime }/L\omega ^{1/2}`$ for large $`\omega `$. For low frequencies, $`\omega t_01`$ (i.e. $`\delta _n\lambda _L`$), the system is dominated by the Meissner effect and

$$\lambda ^{}\lambda _L,\lambda ^{\prime \prime }\lambda _L\frac{1}{2}\omega t_0$$
(10)
and therefore dissipation goes as $`\chi ^{\prime \prime }\lambda ^{\prime \prime }/L\omega `$ for low $`\omega `$.

### II.3 Vortices without pinning

In the presence of vortices, the interesting frequency range is $`\omega t_01`$ since for frequencies $`\omega t_01`$ the system always behaves as in Eq.(9), corresponding to the normal state skin depth effect.

Let us first discuss the case when there is no bulk-pinning, $`\kappa _p=0`$ ($`\lambda _C=\mathrm{}`$). For frequencies such that $`\omega t_{ff}1`$, the system behaves similarly to the Meissner state, with a dissipation peak at $`\omega t_01`$, as described previously. For frequencies such that $`\omega t_{ff}1`$ the system is dominated by the “flux flow skin depth effect” due to the flux flow conductivity $`\sigma _{ff}`$. In this limit $`\stackrel{~}{\lambda }^2i\delta _{ff}^2/2=i\lambda _L^2/\omega t_{ff}`$, and therefore the susceptibility for low frequencies should diverge as $`\stackrel{~}{\chi }\omega ^{1/2}`$. This divergence is cut-off by the finite system size. At low frequencies the real part of the effective penetration depth $`\lambda ^{}`$ saturates to the system size $`L`$. This leads to a dissipation peak at a frequencygeshk

$$\omega _L\frac{c^2}{L^2\sigma _{ff}}\frac{B}{L^2}.$$
(11)
And for very low frequencies such that $`\omega \omega _L`$, dissipation diminishes linearly as

$$\chi ^{\prime \prime }\frac{\sigma _{ff}L^2}{c^2}\omega .$$
(12)
### II.4 Vortices with pinning

In the presence of bulk pinning, $`\kappa _p0`$, the relevant time scale is $`t_C=t_{ff}\lambda _C^2/\lambda _L^2=\eta _v/\kappa _p`$. Now instead of the finite size peak at $`\omega _L`$, there is a vortex dissipation peak when $`\omega t_C1`$. It is worth mentioning that this peak frequency,

$$\omega _C\frac{\kappa _p}{\eta _v},$$
(13)
is independent of magnetic field. For low frequencies, $`\omega t_C1`$, the real part of the penetration depth tends to $`\lambda ^{}\stackrel{~}{\lambda }_C=\sqrt{\lambda _C^2+\lambda _L^2}`$, while the imaginary part $`\lambda ^{\prime \prime }\omega t_C\lambda _C^2/2\stackrel{~}{\lambda }_C`$. Thus for low frequencies the dissipation diminishes linearly with $`\omega `$ as

$$\chi ^{\prime \prime }\frac{1}{2\pi }\frac{\lambda _C^2}{\stackrel{~}{\lambda }_CL}\omega t_C.$$
(14)
Therefore, in the presence of vortices there is a new dissipation peak, in addition to the “two fluid peak”, either at $`\omega _L`$ or at $`\omega _C`$ depending on the importance of pinning. At low frequencies the relevant length scale for penetration of the ac field is either $`L`$ (in the absence of pinning) or $`\lambda _C`$ (if pinning is important) and the dissipation is linear with frequency in both cases.

### II.5 Surface barrier effects

At the surface of a type II superconductor there is a potential barrier that prevents the entrance (and the exit) of vortices. This barrier has been first calculated by Bean and Livingston.bean  The surface barrier inhibits the penetration of flux at $`H_{c1}`$, where penetration is thermodinamically favorable. Instead of $`H_{c1}`$, vortices start to enter at the first penetration field $`H_p>H_{c1}`$. Above $`H_p`$, the Bean-Livinngston surface barrier can stabilize metastable states with a non-equilibrium number of vortices.clemlt13 ; ternovskii  The interval of external magnetic fields where these metastable states exist in macroscopic superconductors has been obtained in Ref. clemlt13, ; ternovskii, ; burlachkov, . It has been found that for fields $`H>H_p`$ the vortex array is separated from the boundary surface by a vortex free region because of the effect of the surface barrier. The length of this vortex free regionclemlt13 ; ternovskii ; burlachkov  is $`d_{SB}=\lambda _L\mathrm{cosh}^1(H/B)`$. The ac magnetic response has been calculated by Sonin and Traitosonin  in this case. They assume that the number of vortices is fixed, and they allow the size of the vortex free region $`d_{SB}`$ to fluctuate with the rf field $`h_{ac}e^{i\omega t}`$. A new dissipation peak is found at a frequency

$$\omega _{SB}\frac{\mathrm{tanh}^2(d_{SB}/\lambda _L)}{t_{ff}},$$
(15)
which for $`BH_p`$ is $`\omega _{SB}1/t_{ff}`$. For low frequencies $`\omega \omega _{SB}`$, the effective penetration depth goes as

$`\lambda ^{}`$ $``$ $`\lambda _L\mathrm{coth}({\displaystyle \frac{d_{SB}}{\lambda _L}})`$ (16)
$`\lambda ^{\prime \prime }`$ $``$ $`{\displaystyle \frac{\lambda _L^2}{\mathrm{sinh}^2(\frac{d_{SB}}{\lambda _L})}}\sqrt{{\displaystyle \frac{2\pi \sigma _{ff}\omega }{c^2}}}.`$
Therefore, according to this approach, the dissipation should go as $`\chi ^{\prime \prime }\lambda ^{\prime \prime }/L\omega ^{1/2}`$ if surface barrier effects are important. Experimentally, dissipation maxima attributed to surface barriers have been measured, for example in Ref.fabrega, . In platelet samples with a perpendicular magnetic field, geometrical effects enhance the surface barrier, giving place to a “geometrical barrier”. Morozov et al.morozov  have measured the effect of geometrical barriers in the ac response of a superconductor. In this case there is a dissipation maximum as a function of magnetic field for $`HH_p`$, which is independent of frequency.

## III MODEL AND DYNAMICS

Our numerical simulations are carried out using the time-dependent Ginzburg-Landau equations complemented with the appropriate Maxwell equations. In the zero-electric potential gauge we have: kato ; gropp

$`{\displaystyle \frac{\mathrm{\Psi }}{t}}={\displaystyle \frac{1}{\eta }}[(i𝐀)^2\mathrm{\Psi }+(1T)(1|\mathrm{\Psi }|^2)\mathrm{\Psi }]`$ (17)
$`{\displaystyle \frac{𝐀}{t}}=(1T)\text{Im}[\mathrm{\Psi }^{}(i𝐀)\mathrm{\Psi }]\kappa ^2\times \times 𝐀`$ (18)
where $`\mathrm{\Psi }`$ and $`𝐀`$ are the order parameter and vector potential respectively and $`T`$ is the temperature. Equations (17) and (18) are in their dimensionless form. Lengths have been scaled in units of the coherence length $`\xi (0)`$, times in units of $`t_0=4\pi \sigma _n\lambda _L^2/c^2`$, $`𝐀`$ in units of $`H_{c2}(0)\xi (0)`$ and temperatures in units of $`T_c`$. $`\eta `$ is equal to the ratio of the characteristic time $`t_0`$ for the relaxation of $`A`$ and the time $`t_{GL}`$ for the relaxation of $`\mathrm{\Psi }`$: $`\eta =t_{GL}/t_0=c^2/(4\pi \sigma _n\kappa ^2D)`$, with $`t_{GL}=\xi ^2/D`$, where $`\sigma _n`$ is the quasiparticle conductivity and $`D`$ is the electron diffusion constant. For superconductors with magnetic impurities we have $`D_{imp}=c^2/(48\pi \kappa ^2\sigma _n)`$, and therefore $`\eta =12`$ in this case.

We solve the TDGL equations following the same procedure as in our previous work in Ref.HD, , using a standard finite difference discretization scheme.gropp  We discretize space in a rectangular mesh with spacings $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }y`$ in each direction. The order parameter and vector potentials are defined at the nodes of the rectangular mesh ($`\stackrel{}{r}=(I,J)`$), and the link variables $`U_{{}_{\mu }{}^{}I,J}=\mathrm{exp}(ı\kappa h_\mu A_{\mu I,J})(\mu =x,y)`$ are introduced in order to maintain the gauge invariance under discretization. We assume that the sample has a square shape in the $`x,y`$ direction and it is infinite in the $`z`$ direction. We apply the magnetic field parallel to the $`z`$ direction, therefore the problem is reduced to two dimensions neglecting all derivatives along $`z`$. The symmetry of the problem implies for all mesh points $`A_{I,J}=(A_{xI,J},A_{yI,J},0)`$ and $`𝐁_{I,J}=(0,0,B_{zI,J})`$, where $`B_{{}_{z}{}^{}I,J}=(\times \stackrel{}{A})_z=(_xA_{{}_{y}{}^{}I,J}_yA_{{}_{x}{}^{}I,J})`$.

The dynamical equations must be complemented with the appropriate boundary conditions for both the order parameter and the vector potential. For the order parameter we have used the boundary condition:

$$(\mathrm{\Pi }\mathrm{\Psi })^{}=(i\stackrel{}{A})^{}\mathrm{\Psi }=0$$
(19)
usually known as the superconductor-insulator boundary condition (S-I) because it implies that the perpendicular component of the superconducting current is equal to zero at the surface ($`\stackrel{}{J}_s^{}=0`$). This boundary condition also minimizes the free energy at the sample surface.

The ac magnetic field is introduced in the simulation through the boundary condition for the vector potential $`A_{\mu I,J}`$. We consider the case where the ac magnetic field $`𝐡_{ac}`$ is parallel to the dc-component $`𝐇_{dc}`$, and both fields are in the $`z`$ direction:

$$B_z|_S=(\times 𝐀)_z|_S=H_{dc}+h_{ac}\mathrm{cos}(\omega t);$$

this expression is evaluated at the sample surface. We study the response of a superconductor in the linear regime, i.e. $`h_{ac}H_{dc}`$. The time dependence of the sample magnetization can be obtained from the magnetic induction averaged over the sample $`B(t)`$ through the relation:

$$4\pi M(t)=B(t)H(t)$$

The ac magnetic susceptibilities are obtained from the Fourier transform of $`M(t)`$:

$`\chi ^{}(H_{dc},\omega )={\displaystyle \frac{1}{\pi h_{ac}}}{\displaystyle _0^{2\pi }}M(t)\mathrm{cos}(\omega t)d(\omega t)`$
$`\chi \mathrm{"}(H_{dc},\omega )={\displaystyle \frac{1}{\pi h_{ac}}}{\displaystyle _0^{2\pi }}M(t)\mathrm{sin}(\omega t)d(\omega t).`$ (20)
In what follows, we have solved the TDGL equations numerically for a type II superconductor with $`\kappa =2`$, mean field temperature $`T=0.5`$ and the parameter $`\eta =12`$. We used a spatial discretization of $`\mathrm{\Delta }x=\mathrm{\Delta }y=0.5\xi `$ and, in order to make efficient calculations, we have chosen adequately the time step with values $`\mathrm{\Delta }t0.015t_0`$.

## IV THE AC SUSCEPTIBILITY IN MESOSCOPIC SUPERCONDUCTORS

In a mesoscopic sample the behavior of the ac magnetic response in the presence of a dc magnetic field can be different than in the macroscopic case. On one hand, the models typically used for macroscopic samples (see Sec.II), assume a high density of vortices (giving an almost uniform magnetic field profile inside the samplebrandt ), applied magnetic fields such that $`H_{c1}HH_{c2}`$, and very large (semiinfinite) samples. On the other hand, in mesoscopic samples there is a small number of vortices, surface barriers can not be neglected and the finite size of the sample is important. Here we perform simulations of squares samples of sizes from $`5\lambda \times 5\lambda `$ to $`20\lambda \times 20\lambda `$ using the TDGL equations.

We start by showing in Fig. 1 a typical high frequency ac magnetic response of a small sample. In Fig. 1(a) we plot the driving field $`H(t)=H_{dc}+h_{ac}\mathrm{cos}\omega t`$ and the corresponding magnetization $`M(t)`$ of the sample. We show the response of the system after a long equilibration time (typically around 100 periods of the external ac field). It can be observed that the magnetization has a sinusoidal behavior, which is a signal that we are in the linear regime and also that the sample is following the external perturbation. It can also be observed that there is a phase shift between the magnetization and the external field due to the presence of dissipation. Fig. 1(b) shows the magnetization loop obtained from the signals of Fig. 1(a). The area inside the loop is proportional to the time average of the energy dissipated in the sample and it is proportional to the imaginary part of the ac susceptibility $`\chi ^{\prime \prime }`$. The appearance of dissipation changes the relationship between $`M_{ac}(t)`$ and $`H_{ac}(t)`$ and therefore the real part of the ac susceptibility, $`\chi ^{}(H_{dc},\omega )`$, will depend on the frequency of the ac signal as well as on the bias field $`H_{dc}`$.

### IV.1 Frequency dependence

In order to understand the different dynamical regimes of the ac response, we study first the dependence of $`\stackrel{~}{\chi }=\chi ^{}+i\chi \mathrm{"}`$ with frequency $`\omega `$ for different values of the dc magnetic field $`H_{dc}`$. In Fig. 2 and Fig. 3 we show the imaginary and the real part of the susceptibility, respectively. These curves were obtained for a sample of $`20\lambda \times 20\lambda `$.

Fig. 2(a) shows the low frequency behavior of $`\chi ^{\prime \prime }`$ for several values of the magnetic field. Fig. 2(b) shows the full $`\chi ^{\prime \prime }(\omega )`$ curves obtained for fields such that there are no vortices within the sample, i.e. in the Meissner state. For fields above the first penetration field $`H_p`$, vortices enter into the sample. The full frequency dependence of $`\chi \mathrm{"}(\omega )`$ in the presence of vortices, i.e. for fields $`H_{dc}>H_p`$, is shown in Fig. 2(c).

The simplest case to understand is the behavior when there is no external field applied, $`H_{dc}=0`$. In this case the complex susceptibility agrees very well, both qualitatively and quantitatively, with the two fluid model described in Sec. IIB. For low frequencies, $`\chi ^{\prime \prime }`$ has a linear dependence with $`\omega `$ as given by Eq. (10). There is a maximum in $`\chi ^{\prime \prime }(\omega ,0)`$ at $`\omega =\nu _o=1/t_0`$. At high frequencies $`\omega \nu _0`$ the dissipation is due to the normal electrons and follows the expected dependence $`\chi ^{\prime \prime }\omega ^{1/2}`$ as given by Eq. (9), corresponding to the normal state skin depth effect.

For small magnetic fields within the Meissner state, $`0<H_{dc}H_p`$, the behavior is qualitatively similar to the $`H_{dc}=0`$ case, with a peak at $`\omega =\nu _0`$. The main difference is that the slope of the low $`\omega `$ linear dependence increases with magnetic field, as can be seen in Fig.2(a).

According to Eq. (10), the linear slope for low frequencies is proportional to $`\lambda _L`$. This implies that the London penetration depth should be dependent on magnetic field: $`\lambda _L(H)`$. Also, Fig.2(b) shows that the magnitude of $`\chi \mathrm{"}(\omega )`$ at any given $`\omega `$ increases for increasing $`H_{dc}`$. These two effects are easily understood since the Meissner screening currents deplete the magnitude of the order parameter at the surface, which leads to a larger static penetration depth $`\lambda _L(H)`$ (see also Sec.IVB) and to an increase in the normal electron dissipation. The most interesting result is that for $`H_{dc}`$ near but still below $`H_p`$ a second dissipative maximum appears in $`\chi \mathrm{"}(\omega )`$. This second peak is at a frequency $`\omega _p`$ two orders of magnitude below $`\nu _o`$.

For $`H_{dc}>H_p`$ the vortices enter into the sample. In this case we still find a low frequency linear dependence $`\chi ^{\prime \prime }\omega `$, as shown in Fig.2(a). Only for large fields close to $`H_{c2}`$ a small departure from linear dependence is observed, $`\chi ^{\prime \prime }\omega ^\alpha `$, with $`\alpha 0.8`$. The behavior for very large frequencies $`\omega \nu _0`$ is the same for all magnetic fields, since it corresponds to the normal state skin depth effect with $`\chi ^{\prime \prime }\omega ^{1/2}`$, as mentioned before. As can be observed in Fig.2(c), the main characteristic of the mixed state is that there are two dissipation maxima as a function of frequency in $`\chi ^{\prime \prime }(\omega )`$. The high frequency maximum corresponds to the “two fluid peak” at $`\omega \nu _0`$. The low frequency maximum that appears at a frequency $`\omega _p\nu _0`$, is a continuation of the maximum observed when $`H_{dc}H_p`$. The peak frequency $`\omega _p`$ shows a weak dependence with magnetic field, $`\omega _p`$ first decreaseas and then increases with increasing $`H_{dc}`$. For fields $`H_{dc}H_{c3}`$ the sample is completely in the normal state and there is a single maximum in the $`\chi ^{\prime \prime }`$ vs $`\omega `$ curve. It corresponds to the frequency at which the normal state skin deph $`\delta _n`$ equals the system size $`L`$.

Let us now analize the frequency dependence of the real part of the susceptiibility, $`\chi ^{}(\omega )`$, which is shown in Fig. 3(a) for some values of $`H_{dc}`$. From $`\chi ^{}(\omega )`$ it is possible to extract the real part of the effective ac penetration depth $`\lambda _{ac}=\lambda ^{}(\omega )`$ using Eq. (3). The length $`\lambda _{ac}`$ represents the length scale for the penetracion of the ac component of the magnetic field and its frequency dependence is shown in Fig. 3(b). In the absence of a bias field ($`H_{dc}=0`$), $`\chi ^{}(\omega )`$ is well described by the two fluid model, as given in Eqs. (9),(10). The sample changes from a London screening for $`\omega \nu _o`$, with $`\lambda _{ac}\lambda _L`$, to a perfect diamagnetic state for $`\omega \nu _o`$, with $`\lambda _{ac}0`$. At low frequencies $`\chi ^{}`$ is larger than $`1/4\pi `$ due to the penetration of the field in a region of size $`\lambda _L`$ near the sample surface, $`\chi ^{}(1/4\pi )(14\lambda _L/L)`$. At high frequency the ac penetration depth is $`\lambda _{ac}\lambda _L`$ and $`\chi ^{}(\omega )1/4\pi `$. The opposite case is for $`H_{dc}=1.9H_{c2}(T)>H_{c3}`$ with the sample mostly in the normal state (at this field the surface superconductivity disappears and there is superconductivity only at the corner of the sampleenomoto ). As can be expected, in this case the system behaves as a normal metal and is dominated by the skin-depth effect with $`\lambda _{ac}\delta _n/2`$. For low frequencies, the skin-depth $`\delta _n`$ equals the system size and thus $`\chi ^{}0`$. For magnetic fields within the Meissner state, $`0<H_{dc}H_p`$, we find that for low frequencies, $`\omega \nu _0`$, the ac penetration depth saturates to the field dependent London penetration depth $`\lambda _L(H)`$, which increseas with $`H`$. The most interesting case is for $`H_{dc}H_p`$, which is shown in Fig. 3. In this situation, one can distinguish three frequency regimes, as it was found before in the dissipative part, $`\chi ^{\prime \prime }(\omega )`$. At the frequencies for which there are maxima in the dissipation, $`\nu _0`$ and $`\omega _p`$, there is a rapid change in the value of $`\chi ^{}(\omega )`$ (or in $`\lambda _{ac}(\omega )`$). For $`\omega \nu _0`$ the normal fluid behavior of Eq. (9) is followed, as expected for all values of $`H_{dc}`$. For the frequency range $`\omega _p<\omega <\nu _0`$, there is a shoulder in $`\chi ^{}(\omega )`$ (and $`\lambda _{ac}(\omega )`$), which basically corresponds to an ac penetration depth of the order of the expected value for $`\lambda _L(H)`$ at that field. At low frequencies, $`\omega <\omega _p`$, the ac penetration depth saturates to a length $`\mathrm{}_p`$ which is smaller than the system size and larger than $`\lambda _L(H)`$.

In summary, the most interesting feature observed in $`\stackrel{~}{\chi }(\omega )`$ is the appearance of a dissipation maximum at frequencies $`\omega _p(H)\nu _0`$ in the presence of a magnetic field. A dissipation maximum in macroscopic samples is expected when $`H0`$, as we discussed in Sec. II. It can be caused by several reasons: (i) finite size, (ii) surface barriers or (iii) bulk pinning. However, the behavior predicted in any of these cases is not followed in the mesoscopic samples studied. Let us discuss each of the possibilities. (i) Finite size: in a small sample this is the first effect to analyze. The characteristic frequency for finite size effects, $`\omega _L`$, should increase linearly with $`B`$ and decrease with system size, as given by Eq. (11). However, $`\omega _p`$ has a weak and non-monotonous dependence with magnetic field, as seen in Fig. 2(c). In Fig. 4 we show $`\chi ^{\prime \prime }(\omega )`$ in a smaller sample with $`L=5\lambda `$ (and $`\eta =12`$) for $`H_{dc}H_p`$. We see that the maximum appears at the same frequency $`\omega _p`$ as in the $`L=20\lambda `$ case for the same magnetic field \[Fig. 2(b)\].

Therefore $`\omega _p`$ has no size dependence. Furthermore, in Fig. 3 we see that when $`\omega 0`$ the sample is still diamagnetic and $`\lambda _{ac}<L`$ (a finite size effect would have $`\chi ^{}(\omega 0)0`$ and $`\lambda _{ac}L`$). (ii) Surface barrier: in mesoscopic samples, surface barriers are very important.HD  In macroscopic samples, the effect of surface barriers, assuming a large number of vortices, gives a low frequency dependence for dissipation as $`\chi ^{\prime \prime }\omega ^\alpha `$ with $`\alpha =1/2`$, see Eq. (16) and Ref. sonin, . However, we find a linear frequency dependence in most cases as shown in the inset of Fig. 2(a), and even when there is a departure from linearity, it is with an exponent of $`\alpha 0.8>0.5`$. Also the characteristic frequency $`\omega _{SB}`$ of Eq. (15) has a strong dependence on magnetic field, which is not observed here. (iii) Pinning: since there is no disorder in the model, no bulk pinning is expected. However, in a very small sample one may argue that, given the smallness of the system, the sample itself acts as a confinement potential for vortices. In this case, $`\kappa _p`$ in Eq. (7) would represent the vortex elastic response for small oscillations within this confinement potential. The characteristic frequency for pinning effects, $`\omega _C`$, is independent of magnetic field in macroscopic samples, according to Eq. (13). The weak field dependence of $`\omega _p`$ may be still consistent with this result (the effective $`\kappa _p`$ may depend on the number and distribution of vortices inside the sample in this hypothetic case of “confinement potential”). Also the fact that $`\lambda _{ac}<L`$ for $`\omega \omega _p`$ in Fig. 3(b) is consistent with a finite “Campbell’s penetration depth” $`\lambda _C`$ (and then $`\mathrm{}_p\lambda _C`$). However, we also have to discard this possibility, since we find that there is a dissipation maximum at $`\omega _p`$ even when there are no vortices in the sample. As shown in Fig. 2(b) the dissipation peak already appears for magnetic fields just below the penetration field $`H_p`$. To be more precise, we have explicitly calculated the topological number of total vorticity as $`N_v=\frac{1}{\mathrm{\Phi }_o}(A+\frac{J_s}{|\mathrm{\Psi }|^2})𝑑l`$. We obtain $`N_v=0`$ at all times for the magnetic fields shown in Fig. 2(b). Therefore, all arguments based on the oscillation of vortices, as for example the Coffey and Clem modelclem  of Eq. (8), which work very well for macroscopic samples, are not enough to explain the dissipation maxima observed in these mesoscopic samples.

In Fig. 4 we show that $`\omega _p`$ depends directly on the time scale for the relaxation of the amplitude of the order parameter, $`t_{GL}=\eta t_0`$. We consider a magnetic field $`H_{dc}H_p`$ and we change the value of $`\eta `$ in Eqs.(17) and (18). The “two fluid peak” is always at the same frequency $`\nu _0=1/t_0`$, since $`t_0`$ is fixed. We find that when decreasing $`\eta `$, the frequency $`\omega _p`$ shifts to higher values and increases monotonically with $`1/t_{GL}`$, until the rather unphysical case of $`t_{GL}<t_0`$ ($`\eta <1`$), when only the two fluid peak at $`\nu _0`$ is observed. The dependence of $`\omega _p`$ with $`t_{GL}`$ shows that variations in the amplitude of the order parameter are the main mechamism of this dissipation maximum, as we will see more clearly later in Sec. V. In the inset of Fig. 4 we verify that for $`H_{dc}=0`$ the behavior is independent of $`t_{GL}`$, since this case can be simply described with the two fluid model.

The appearance of the dissipation maximum for a magnetic field below $`H_p`$, just before the penetration of vortices into the sample, leads to hystory dependent effects. The curves in Fig. 2 were generated increasing the field $`H_{dc}`$ at low frequency and then measuring $`\chi \mathrm{"}(\omega )`$ for increasing $`\omega `$. In this way, the curves reported in Fig. 2(c) have the same number of vortices at all frequencies for a given $`H_{dc}`$. If we follow the opposite procedure, i.e., we increase $`H_{dc}`$ at high frequency and then measure $`\chi \mathrm{"}(\omega )`$ for decreasing $`\omega `$, the curves can show discontinuous jumps at certain frequencies values due to vortex penetration. This occurs because increasing the field at high frequency can generate a metaestable state that can become unstable at low frequencies. In Fig. 5 we show this case, for a field near the first penetration field $`H_{dc}H_p`$. Depending on the history, the high frequency branch can have either $`N_v=12`$ vortices or $`N_v=0`$. In the case when the frequency is decreased from a high frequency state, there is a large jump in $`\chi ^{\prime \prime }(\omega )`$ at a frequency near $`\omega _p`$, from a metastable branch with $`N_v=0`$ for $`\omega >\omega _p`$ to a branch with $`N_v=12`$ for $`\omega <\omega _p`$. In contrast, we see in the inset of Fig. 5 that for a field $`H_{dc}H_p`$ the curves of $`\chi ^{\prime \prime }(\omega )`$ are independent of history.

### IV.2 Magnetic field dependence

Experimentally, the ac frequency is fixed and the dc magnetic field or the temperature can be varied. In this section we will study the ac magnetic response of mesoscopic superconductors varying the bias dc magnetic field. We will show results at different frequencies and for different sample sizes.

In the absence of the ac field, the dc magnetic behavior of mesoscopic superconductors is different from the continuos macroscopic behavior.bolech ; HD  In mesoscopic samples each vortex entrance event produces discontinuos jumps in the magnetization curve at certain $`H_{dc}`$ values. Each discontinuous jump in $`M(H)`$ corresponds to a sudden increase in the number of vortices. These jumps occur at succesive magnetic fields $`H_{pi}=H_{p1},H_{p2},H_{p3},\mathrm{}`$.HD  In the regions of $`H_{pi}<H<H_{p(i+1)}`$ the number of vortices $`N_v`$ is constant. The only penetration events occur at $`H_{pi}`$, when vortices enters into the sample. In the region of constant vorticity, $`H_{pi}<H<H_{p(i+1)}`$, one may think that no vortices enter the sample because of the effect of the surface barrier.

Fig. 6 (a) shows the behavior of $`\chi \mathrm{"}(H_{dc})`$ obtained at a fixed frequency $`\omega =0.004\nu _o<\omega _p`$ in a mesoscopic sample of size $`10\lambda \times 10\lambda `$. At small $`H_{dc}`$, in the Meissner state, $`\chi \mathrm{"}`$ increases continuously with increasing $`H_{dc}`$. The presence of the dc magnetic field in the Meissner state induces static supercurrents which deplete the order parameter at the boundary. As explained before in Sec. IVA, the depletion of the order parameter leads to an increase of dissipation for increasing $`H_{dc}`$. When the dc magnetic field is increased above the first penetration field $`H_pH_{p1}`$, the first vortices enter into the sample. Fig. 6(a) shows that there is a discontinuous jump in $`\chi \mathrm{"}`$ with a decrease of dissipation just at $`H_p`$. Two states are possible exactly at $`H_p`$: one state without vortices and high dissipation, and one sate with vortices inside the sample and low dissipation. This is in agreement with the history dependence observed in Fig. 5. Further jumps in the dissipation curve $`\chi \mathrm{"}(H_{dc})`$ are present at the other magnetic fields for vortex penetration, $`H_{p2},H_{p3},\mathrm{}`$. The jumps are followed by a later continuos increase of $`\chi \mathrm{"}(H_{dc})`$ with increasing $`H_{dc}`$ while the number of vortices $`N_v`$ remains fixed. Fig. 6(b) was obtained in a sample of the same size of Fig. 6(a) but at a higher frequency $`\omega =0.06\nu _o>\omega _p`$.

The frequencies used in Figs. 6(a) and 6(b) are at both sides of the dissipation maximum that appears at $`\omega _p`$. Both curves are similar in their qualitative features, but show a few differences. At high magnetic fields there are less jumps in $`\chi ^{\prime \prime }(H_{dc})`$ in Fig. 6(b) than in Fig. 6(a). This is because at high frequencies it is possible to remain in a metastable state with a fixed number of vortices in a wider range of magnetic field. At the same time, the higher frequency of Fig. 6(b) produces a decrease in the amplitude of the jumps. The behavior of $`\chi ^{\prime \prime }(H_{dc})`$ of Figs. 6(a) and 6(b) is different from the behavior of bulk samples. For example, Fig. 6(c) shows $`\chi \mathrm{"}(H_{dc})`$ at $`\omega =0.06\nu _o`$ obtained in a “large” sample of $`40\lambda \times 40\lambda `$. In this case there are no discontinuous jumps and the curve is almost continuous. In macroscopic samples the entrance of vortices for $`H_{dc}>H_p`$ increases the ac losses in a continuos way. In the macroscopic case vortices play a fundamental role in the ac losses, and therefore an increase on the number of vortices increases the dissipation. In contrast, the results of Figs. 6(a) and 6(b) show the opposite behavior in mesoscopic samples: there is a sudden decrease of the ac losses at the magnetic fields that correspond to an increase in the number of vortices.

Fig. 7 shows the real part of the ac susceptibility for increasing $`H_{dc}`$. We use the same parameters as in Fig. 6. The curves of $`\chi ^{}(H_{dc})`$ also have a typical mesoscopic behavior with discontinuities at each vortex penetration field. The increase of frequency from $`\omega =0.004\nu _o`$ in Fig. 7(a) to $`\omega =0.06\nu _o`$ in Fig. 7(b) shows similar differences as observed in $`\chi \mathrm{"}(H_{dc})`$. It can also be noted in Fig. 7(a) that the differential susceptibility has a paramagnetic behavior ($`\chi ^{}>0`$) near each vortex penetration field. This can be expected from the behavior of the dc magnetization curves of mesoscopic samples. bolech  In mesoscopic samples, for $`H_{dc}`$ below each vortex penetration field a magnetization maximum appears and therefore a region near $`H_{pi}`$ where the magnetization increases at increasing field. This shows in the low frequency ac response of mesoscopic samples as a differential paramagnetic behavior. Fig. 7(c) shows a continuos macroscopic dependence of $`\chi ^{}(H_{dc})`$ in a larger sample of size $`40\lambda \times 40\lambda `$. At the frequencies used in Fig. 7(c) we can see a diamagnetic behavior even in the normal state for fields greater than $`2H_{c2}(T)`$. This is simply because at high frequencies there is magnetic screening in the normal state due to the skin depth effect.

From Fig. 7 it is possible to obtain the $`H_{dc}`$ dependence of the London penetration depth. This dependence has been a subject of interest recently because of its relation with the simmetry of the order parameter. In conventional s-wave superconductors the following dependence is expected:

$$\mathrm{\Delta }\lambda (H_{dc},T)/\lambda (0,T)=\beta (T)(H_{dc}/H_{c2}(T))^2$$
(21)
where $`\mathrm{\Delta }\lambda (H_{dc},T)=\lambda (H_{dc},T)\lambda (0,T)`$. It is assumed that the sample is in the Meissner state and in the linear regime, i.e. $`|H_{dc}||h_{ac}|`$. Eq. (21) was originally obtained by Ginzburg and Landau landau  and also by Bardeen bardeen  an it has been observed experimentally in type I pippard ; sridhar2  and conventional type II superconductors maeda1 . The behavior of unconventional superconductors is expected to be different maeda2 . In the case of a d-wave superconductor Yip and Sauls yip1 ; yip2  proposed a linear field dependence $`\mathrm{\Delta }\lambda (H_{dc})H_{dc}`$. Here we use the phenomenological TDGL equation in its s-wave form. Therefore, in the Meissner state we expect a result similar to Eq.(21) for macroscopic samples. The results are shown in Fig. 8. Both the macroscopic and the mesoscopic sample show a quadratic dependence of $`\mathrm{\Delta }\lambda (H_{dc})`$ as a function of $`H_{dc}`$ in the Meissner state. The vortex penetration produces a sudden decrease of $`\mathrm{\Delta }\lambda `$ in mesoscopic samples, whereas at a fixed number of vortices $`\mathrm{\Delta }\lambda `$ increases continuously for increasing $`H_{dc}`$. For macroscopic samples the entrance of vortices at the first penetration field $`H_p`$ does not produce discontinuity in $`\mathrm{\Delta }\lambda (H_{dc})`$ but a change of slope as can be observed in the plot. For $`H_{dc}>H_p`$ the dependence changes from quadratic to approximately linear in $`H_{dc}`$.

## V NASCENT VORTICES EFFECT

From the results of the previous section, it is clear now that the ac magnetic response of a mesoscopic superconductor is different from the macroscopic behavior and it can not be explained with the known models. An interesting result is the significant decrease of dissipation each time the number of vortices increases. To understand this effect we have analyzed in detail the time variation of the order parameter and the magnetization within the sample.

Fig. 9 shows the time variation of the order parameter at different applied magnetic fields in a small sample. Each curve represents the amplitude of the order parameter along the $`x`$ direction, perpendicular to a face of the square sample, obtained at a given time. To simplify the analysis the cut was done just at the center of the face at $`y=L/2`$ (which is the preferred point for vortex penetration in a square sample). The order parameter is obtained at different times within a single period of the ac magnetic field. In particular, we show the behavior of the order parameter just before and after the two first values of the penetration fields, $`H_{p1}`$ and $`H_{p2}`$. At $`H_{dc}=0`$ the order parameter $`|\mathrm{\Psi }|`$ is equal to unity in the entire sample and the dissipation is small. In the Meissner state but for $`H_{dc}H_{p1}`$, Fig. 9 (a) shows that the magnitude of the order parameter near the boundary, and at the center of the face, has large oscillations around a small value. Fig. 9 (b) shows the time evolution of order parameter just after the first vortex entrance event, $`H_{dc}H_{p1}`$. We see that the entrance of vortices increases the order parameter at the boundary, where now $`|\mathrm{\Psi }|`$ oscillates around a value much higher than before. One can also see that the positions of the vortices are fixed inside the sample and do not oscilllate in the presence of the ac signal. Recently, in an extension of the idea of the Bean Livingston barrier to mesoscopic superconductos, we have shown that the vortex entrance produces a reinforcement of the surface barrier, due to the interaction between the vortices that are trying to enter the sample and the vortices already inside the sample.HD  The presence of this interaction decreases the surface currents. Thus, the main effect after the entrance of vortices inside a mesoscopic superconductor is that the Meissner currents decrease and the order parameter increases at the boundary. Therefore the normal state electrons decrease and then the dissipation decreases.

Figs. 9(c) and 9(d) show similar behavior for $`H_{dc}H_{p2}`$ and $`H_{dc}H_{p2}`$. After increasing $`H_{dc}`$ at a constant number of vortices, the order parameter at the boundary has decreased, as seen in Fig. 9(c). The decrease of the order parameter increases the dissipation as was observed before. If new vortices enter again into the sample, the order parameter increases at the boundary, see Fig. 9(d) and therefore the dissipation decreases again. This process is repeated in the other penetration fields $`H_{pi}`$. On the contrary, the entrance of vortices in a large sample do not produce an appreciable change in the amplitude of the order parameter at the boundary before and after the penetration event, and therefore the dissipation shows a continuous behavior \[as seen in Fig. 6(c)\].

A careful analysis of the entire profile of the order parameter for magnetic fields below but near to $`H_p`$ suggest that the depletion of the order parameter near the boundary may be attributed to “vortices” that are about to enter into the sample. These “nascent vortices” are located outside of the sample but near the surface. Even more, for higher values of $`h_{ac}`$ (beyond linear response) we find that these “nascent vortices” oscillate near the surface (i.e. that the profile of the depletion of $`|\mathrm{\Psi }|`$ near the surface oscillates with $`H(t)`$). According to this scenario, the length scale $`\mathrm{}_p`$ observed in Fig. 3(b) for low frequencies may correspond to a characteristic length scale for the nucleation of vortices from the boundary. (For large fields $`BH_p`$, the length $`\mathrm{}_p`$ could be related to the size of the vortex free region $`d_{SB}`$ of macroscopic samples). Years ago, Walton and Rosemblumwalton  suggested the idea of nascent vortices as a possible source of the high-frequency losses in superconductors. This idea had further support from the analytical work of Kramer kramer  who found static solutions of the Ginzburg-Landau equations which corresponded to vortices nucleated near the sample. The nascent vortices can oscillate under an ac field. This source of dissipation is not generally taken into account in macroscopic samples because in this case is small compared with the dissipation due to the motion of vortices inside the sample. However in a mesoscopic superconductor, with a few vortices confined inside the sample, this source of dissipation can be a significative part of the dissipation observed, as we have found here.

## VI SUMMARY AND CONCLUSION

We have studied the response of mesoscopic superconductors to an ac magnetic field. In mesoscopic superconductors vortices are confined inside the sample by the surface barrier and their dynamics play a secondary role in the dissipation. When increasing $`H_{dc}`$, discontinuities in $`\chi ^{}`$ and $`\chi \mathrm{"}`$ appear at each vortex penetration event. The dissipation is maximum before the penetration of vortices due to the effect of “nascent vortices”, which lead to large oscillations of the amplitude of the order parameter at the boundary. After the penetration of vortices, the normal currents decreases and therefore the dissipation decreases, leading to a sudden suppression of the ac losses induced by the vortex entrance.

###### Acknowledgements.

A.D.H. acknowledges Oscar Arés and C. Hart for useful comments and help and the Centro Latino-Americano de Física (CLAF) for financial support. D.D. acknowledges support from Conicet and CNEA. We acknowledge the financial support from the Argentina-Cuba cooperation project SETCIP-CITMA: CU/A00-EIII/005 and from ANPCyT (Argentina) project number PICT99-03-06343.