text:
# Long-Range Electron Transfer and Electronic Transport Through Macromolecules

## Abstract

A theory of electrical transport through molecular wires is used to estimate the electronic factor in the intramolecular electron transfer (ET) in porphyrin-nitrobenzene supermolecules, and to analyze its structure. The chosen molecules have complex donor and acceptor configurations, and relatively simple structure of the bridge, which enables us to concentrate our studies on donor/acceptor coupling to the bridge. We present analytical and numerical results concerning the effect of donor/acceptor coupling to the bridge on the ET process in molecules with complex donor/acceptor subsystems.

Long-distance electron transfer (ET) plays a central role in many biological processes 1 . It has been established that molecular ET is essentially a combination of nuclear environment fluctuations and electron tunneling 2 . The expression for the ET rate including both electronic and nuclear factors was first proposed by Marcus 3 ; 4 ; 5 . The existing theories of long-range ET are mainly concentrated on the analysis of the electron tunneling through the bridge (Refs 6 ; 7 ; 8 ; 9 ; 10 ; 11 ; 12 ). Hitherto, less attention has been paid to systematic studies of donor/acceptor coupling to the bridge and its contribution to the electronic transmission coefficient. Usually, theoretical models of ET present both donor and acceptor as single sites with one state per site. In realistic biological systems, however, donor and acceptor subsystems are complex and include numerous sites which provide effective coupling to the bridge.

In this letter, we study the effect of the donor/acceptor subsystems on the characteristics of the long-range ET in macromolecules. We use an apparent close similarity between this process and an electron transport through molecular wire placed between metallic contacts. This enables us to apply the methods developed in Refs. 13 ; 14 ; 15 ; 16  to calculate the electron transmission function $`T(E)`$ which we can further analyze to separate out contributions from different donor/acceptor sites participating in the ET process. We also calculate current-voltage characteristics of the intramolecular electron transport.

In further studies, we simulate our donor and acceptor subsystems as sets of semi-infinite tight-binding homogeneous chains. Each chain is attached to a site of donor/acceptor which can be effectively coupled to the bridge. To simplify the following calculations we assume that the chains do not interact with each other. This model for the donor/acceptor subsystem is a generalization of that adopted in the theory of electron transport through molecular wires 13 ; 14 ; 15 ; 16  where each metallic electrode was treated as a single chain. However, in the present analysis we need this generalization to take into account the complexity of the donor/acceptor subsystems.

To proceed we consider each element of the donor-bridge-acceptor system as an array of potential wells. Each well represents one atomic orbital of the corresponding part of the molecule (donor, acceptor or bridge). These wells are sites of the donor-bridge-acceptor system. Within this approach, any atom will be represented by a set of sites corresponding to its states, and we do not distinguish between sites corresponding to the same atom, and those corresponding to different atoms of our molecular system 8 ; 17 .

In the following calculations we start from a tight-binding Hamiltonian for the bridge:

$$H=H_0+H_1$$
(1)
whose matrix elements between states $`|i>`$ and $`|j>`$ corresponding to the $`i`$-th and $`j`$-th sites are given by:

$$\left(H_0\right)_{ij}=\alpha _i\delta _{ij};\left(H_1\right)_{ij}=𝒱_{ij},$$
(2)
where $`\delta _{ij}`$ is Kronecker symbol, $`𝒱_{ij}`$= 0 when $`i=j`$ and only states associated with valence electrons are considered. The diagonal matrix elements $`\alpha _i`$ are ionization energies of electrons at sites $`i`$, while the off-diagonal matrix elements $`𝒱_{ij}=𝒱_{ji}`$ are the hopping strengths between the $`i`$th and $`j`$th sites. Both direct and exchange energy contribution are included in $`𝒱_{ij}.`$ As well as it was carried out in 15 ; 16  we also take into account self-energy corrections arising due to the coupling of the donor $`\left(\mathrm{\Sigma }_𝒟\right)`$ and acceptor $`\left(\mathrm{\Sigma }_𝒜\right)`$ to the bridge. As a result, we arrive at the effective Hamiltonian for the bridge, i.e.,

$$H_{eff}=H+H_𝒟+H_𝒜.$$
(3)
Generalizing the results of Ref. 18  for the case when we have several donor and acceptor sites coupling to the bridge, we obtain

$$(H_𝒟)_{ij}=\left(\mathrm{\Sigma }_D\right)_i\delta _{ij}\underset{k}{}\frac{𝒟_{ki}^{\mathrm{\hspace{0.17em}2}}}{Eϵ_k\sigma _k}\delta _{ij},$$
(4)
$$(H_𝒜)_{ij}=\left(\mathrm{\Sigma }_A\right)_i\delta _{ij}\underset{l}{}\frac{𝒜_{il}^{\mathrm{\hspace{0.17em}2}}}{Eϵ_l\sigma _l}\delta _{ij}.$$
(5)
Here, $`𝒟_{ki}`$ and $`𝒜_{il}`$ are, respectively, coupling strengths between the $`k`$-th donor site or the $`l`$-th acceptor site and the $`i`$-th site of the bridge, and $`\sigma _{k,l}=\frac{1}{2}\left\{\theta _{k,l}i\sqrt{4\gamma _{k,l}^2\theta _{k,l}^2}\right\}`$ are the self-energy corrections of the semiinfinite chains attached to the corresponding sites 13 . The parameters $`\theta _{k,l}=Eϵ_{k,l}`$, $`ϵ_{k,l}`$ and $`\gamma _{k,l}`$ are the ionization energies of electrons at the corresponding donor/acceptor sites, and the nearest-neighbor hopping integrals for the chains. Summation in Eqs. (4) and (5) is carried out over all donor/acceptor sites coupled to the bridge. Due to the presence of the self-energy corrections $`\mathrm{\Sigma }_{𝒟,𝒜}`$, the eigenvalues of the effective Hamiltonian (3) include imaginary parts $`\mathrm{\Gamma }_i`$ which represent broadening of the bridge energy levels $`E_i`$ originating from the coupling of the bridge to the donor and acceptor systems. The energy levels are broadened further if include scattering processes in the bridge. Here, however, for simplicity, we separate the effect of electronic dephasing, assuming that coherent electronic tunneling through a bridge to be the predominant mechanism of ET.

Treating the long-range ET processes, we can neglect the direct coupling of donor sites to acceptor sites. Then the probability amplitude for the transition between the $`k`$-th donor site and the $`l`$-th acceptor site is given by 19

$$T_{kl}=\underset{i,j}{}𝒟_{ki}𝒢_{ij}𝒜_{jl}.$$
(6)
Here, the sum is carried out over the bridge sites, $`G_{ij}`$ is the matrix element of the Green’s function corresponding to $`H_{eff}`$, i.e.,

$$𝒢_{ij}=i\left|\left(EH_{eff}\right)^1\right|j.$$
(7)
An electric tunneling current $`I`$ flowing from donor to acceptor through the bridge in the presence of a small applied voltage $`V`$ has the form 20

$$I=\frac{2\pi e}{\mathrm{}}\underset{k,l}{}f\left(E_k\right)[1f\left(E_l+eV\right)]T_{kl}^2\delta \left(E_kE_l\right)$$
(8)
where $`f(E)`$ is the Fermi function; $`e`$ is the charge of the electron. Assuming that the applied potential varies linearly in the molecule, we can use the approximation $`V=V_0(L_b/L_{mol}),`$ where $`L_b`$ and $`L_{mol}`$ are the lengths of the bridge and the whole molecule, respectively, and $`V_0`$ is a voltage applied across the whole molecule. Starting from the expression (8) and following a usual way we arrive at the standard formula 21 :

$$I=\frac{e}{\pi \mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑑ET(E)[f\left(E\mu _1\right)f\left(E\mu _2\right)],$$
(9)
where the chemical potentials $`\mu _1`$ and $`\mu _2`$ are determined by the equilibrium Fermi energy of the bridge $`E_F`$ and the effective voltage $`V`$ across the bridge 16 , i.e.,

$$\mu _1=E_F+\left(1\eta \right)eV;\mu _2=E_F\eta eV.$$

Here the parameter $`\eta `$ characterizes how the voltage $`V`$ is divided between the two ends of the bridge. The equilibrium Fermi energy for the bridge can be determined from the equation (See Ref. 16 ):

$$N=\frac{2}{\pi }\underset{i}{}\left(\frac{\pi }{2}+\mathrm{arctan}\frac{E_fE_i}{\mathrm{\Gamma }_i}\right),$$
(10)
where $`N`$ is the number of electrons on the bridge part of the molecule.

The electron transmission function included in the expression (9) is given by the formula

$$T(E)=2\underset{i,j}{}\mathrm{\Delta }_i|G_{ij}|^2\mathrm{\Delta }_j,$$
(11)
where the summation is carried out over the states of the bridge, and the quantities $`\mathrm{\Delta }_{i,j}`$ are defined as

$$\mathrm{\Delta }_i=\text{Im}\left(\mathrm{\Sigma }_𝒟\right)_{ii},\mathrm{\Delta }_j=\text{Im}\left(\mathrm{\Sigma }_𝒜\right)_{ii}.$$

This formula agrees with the result obtained earlier in the study of electron transport through molecular wires 15 .

We see that contributions to the transmission function $`T(E)`$ from different donor and acceptor sites are determined by the factors $`\mathrm{\Delta }_{i,j}`$. It follows from Eqs. (4) and (5) that these factors strongly depend on energy-dependent denominators. When the tunneling energy $`E`$ takes on a value close to $`ϵ_{k,l}`$, the corresponding term in Eqs. (4) and (5) can surpass all remaining terms. Thus, for different values of $`E`$, different donor and acceptor sites can predominate in the ET process. Expressions (4) and (5) enable us to determine these dominating sites for any interval of tunneling energy $`E`$. In Fig.1 we present results for the calculated I–V characteristics for the porphyrine-nitrobenzene macromolecules for different temperatures (T = 1K, 50K, 100K, 300K). In calculations we estimate $`E_F`$ using the expression (10) and we choose $`\eta =0.3.`$ Matrix elements of the bridge Hamiltonian (1) as well as the ionization energies for the donor/acceptor sites are obtained from the extended Huckel method. The parameters $`\gamma _{k,l}`$ representing hopping integrals for donor/acceptor chains are estimated as averages of matrix elements which correspond to the coupling of these sites to their nearest neighbors within the donor/acceptor part of the molecule.

Comparison of our results with those obtained for I–V curves characterizing the electric transport through molecular wires 14  shows their resemblance at room temperatures. At low temperatures we notice that our I–V curves exhibit a step-like behavior and the plateaus get narrower as the temperature raises. Associating the low-temperature I–V curves with the result for the electron transmission function $`T(E)`$ shown in part in Fig. 2 we see that the steps reflect the structure of $`T(E)`$. The latter is a set of peaks separated by the intervals where the electron transmission function takes on values close to zero. These peaks arise due to the contributions from different donor-acceptor sites as well as from the Greens function corresponding to the effective Hamiltonian for the bridge. At a given voltage $`V`$ the difference of the Fermi functions in expression (9) differs from zero only in the interior of a certain energy range including $`E_F`$. Therefore, the tunneling current $`I`$ is fully determined by the contributions of peaks of the transmission function located there. When we increase the voltage applied across the bridge the width of the relevant energy range enhances, and extra peaks of the function $`T(E)`$ contribute to the current. When a new peak appears this gives rise to a sudden enhancement of $`I`$, so that the step exhibited in the $`IV`$ curves reflect the structure of $`T(E).`$ At higher temperatures Fermi distribution functions in the integrand of (9) lose their step-like character and the plateaus are washed out. Another reason for peaks in the electron transmission function (and the steps in the I–V curves, respectively) to be eroded is the electronic phase-breaking effect which arises in complex molecules due to stochastic fluctuations of ion potential. We conjecture, however, that at low temperatures the electronic dephasing reduces to same extent, so that the structure of $`T(E)`$ can be revealed.

In summary, we have presented a theoretical model for the intramolecular ET which resembles those used to study the electric transport through molecular wires. Our model enables us to analyze in detail the donor/acceptor coupling to the molecular bridge which is important for further investigations of long-range ET through macromolecules.

Acknowledgments. We thank D. Akins for providing us the data concerning the structure of porphyrin-nitrobenzene molecules, and J. Malinsky for helpful discussions. GG acknowledges the support in part from grant # 4137308-02 from the NIH. NAZ also thanks G.M. Zimbovsky for the help with the manuscript.


format: text

text:
# Vortices in layered superconductors with columnar pins: A density functional study

## I Introduction

In the mixed phase of type-II superconductors, magnetic flux penetrates the sample as quantized vortex lines which form a special physical system known as “vortex matter”. The fascinating equilibrium and dynamical properties of vortex matter in the mixed phase of high-temperature superconductors (HTSCs) have prompted considerable experimental and theoretical attention  for more than a decade. Because of enhanced thermal fluctuations, the Abrikosov lattice in very pure samples of these highly anisotropic, layered materials is observed to undergo a first order melting transition  into a resistive vortex liquid as the temperature is increased.

The properties of the mixed phase of HTSCs are also known  to be strongly affected by the presence of pinning centers, either intrinsic to the sample or artificially generated. Understanding the effects of pinning in these systems is very important for practical applications because the presence of pinning strongly influences the value of the critical current in the mixed phase. Columnar pinning arising from damage tracks produced by heavy-ion bombardment has received much attention in this context because such extended defects parallel to the direction of the average magnetic flux are highly effective  in increasing the critical current by localizing vortex lines along their length. Columnar defects produce “strong pinning” in the sense that the pinning potential of a defect is sufficiently strong to pin a vortex line at low temperatures. Heavy-ion irradiation generally produces a random array of parallel columnar defects. The effects of such an array of extended defects on the properties of the mixed phase of HTSCs have been extensively studied experimentally , theoretically  and numerically . The thermodynamics of a collection of vortex lines in the presence of a parallel array of random columnar pins has been analyzed  by mapping the problem to the quantum mechanics of a system of two-dimensional interacting bosons in an external random potential. The main prediction of such theories is the existence of a low-temperature “Bose glass” phase , separated by a continuous phase transition from a high-temperature, entangled liquid of vortex lines. The theoretically predicted scaling behavior near the Bose glass transition has been observed . A random array of columnar pins also affects the equilibrium properties of the high-temperature vortex liquid, leading to the occurrence of anomalies  in the reversible magnetization curve near $`B=B_\varphi `$, where $`B_\varphi =\rho _p\mathrm{\Phi }_0`$ ($`\rho _p`$ is the areal density of columnar pins and $`\mathrm{\Phi }_0=hc/2e`$ is the flux quantum) is the so-called “matching field”, and $`B`$ is the magnetic induction that determines the areal density $`\rho _0`$ of vortex lines ($`\rho _0=B/\mathrm{\Phi }_0`$).

It is also possible, through the use of a variety of nanofabrication techniques , to create periodic arrays of strong (in the above mentioned sense) pinning centers in thin-film superconductors. The interplay between the lattice constant of the pin array (determined by $`B_\varphi `$) and the intervortex separation (determined by $`B`$) can lead to a variety of interesting commensurability effects in such systems. Some of these effects have been observed in recent experiments. Imaging experiments  have shown the formation of ordered structures of the vortex system at low temperatures for commensurate values of $`B/B_\varphi `$. Magnetization measurements  in the irreversible (vortex solid) regime have demonstrated the occurrence of anomalies at harmonics of $`B_\varphi `$. The effectiveness of pinning at integral values of $`B/B_\varphi `$ has been found  to produce regularly spaced sharp minima in the resistivity versus field curve. A pinning-induced reconfiguration of the vortex lattice has been observed in an experiment  on a thin-film superconductor with a rectangular array of pinning centers. Some of these effects have been studied theoretically, using analytic  and numerical  methods. Bulk HTSC samples with periodic arrays of columnar pins have not been fabricated yet, but the technology for doing this appears to be within reach .

A periodic array of strong pinning centers should have significant effects on the melting transition of the vortex lattice. We consider here the situation where $`B>B_\varphi `$, that is, when the pin array is relatively dilute. If, in addition, the value of $`B`$ is such that the melting temperature of the vortex lattice in the pure system is substantially lower than the superconducting transition temperature in zero field, then each pinning center would pin a vortex at temperatures comparable to the melting temperature of the pure vortex lattice. However, the interstitial vortices, which would be present whenever the number of pinning centers is smaller than the number of vortices (assuming that each defect can pin at most one vortex), may undergo a sharp melting transition. This would certainly be the case in the limit where the spacing of the pin array is sufficiently large. Since the vortices pinned at the defects produce a periodic potential for the interstitial ones, the melting transition of the latter is an example of a solid to liquid transition in the presence of an external periodic potential. Evidence for such melting of interstitial vortices has been found in experiments  on thin-film superconductors with periodic pinning. However, the thermodynamic behavior at the melting transition has not been characterized in these experiments. The effects of a periodic potential on the melting of two-dimensional solids have been investigated earlier using analytic  and numerical  methods. We are not aware of any theoretical study of the effects of a periodic array of columnar pins on the vortex lattice melting transition in three dimensions.

In this paper, we report the results of a study of the equilibrium properties of the mixed phase of highly anisotropic, layered superconductors in the presence of columnar pins. We consider a geometry in which both the magnetic field and the direction of the columnar defects are perpendicular to the superconducting layers. Our study is based on a model free energy functional  for a system of “pancake” vortices lying on the superconducting layers. We consider the limiting case of a vanishingly small Josephson coupling between the layers, so that the pancake vortices on different layers interact via only their electromagnetic coupling. As shown in earlier studies , this limit is appropriate for describing extremely anisotropic Bi- and Tl-based HTSCs. The Ramakrishnan-Yussouff (RY) free-energy functional  used in the present work is the same as that used in earlier studies  of vortex lattice melting in pure systems. The same free energy functional was also used, in combination with the replica method for treating quenched disorder, in a study  of the effects of random point pinning on the melting line in the $`BT`$ plane. In these earlier studies, the density distribution in the crystalline state was expressed in terms of a few “order parameters” and the free energy was minimized with respect to these parameters. In the present work, we use a different method which is more powerful and more appropriate for describing in detail pinning-induced inhomogeneities of the local density. This method, developed in our studies  of the hard-sphere system, involves direct numerical minimization of a discretized version of the free energy functional. Since both the magnetic field and the direction of the columnar pins are assumed to be perpendicular to the layers, the time-averaged local density of pancake vortices must be the same on all the layers. This simplification makes the problem effectively two-dimensional and allows a high-precision numerical investigation of the effects of columnar pins on the structure and thermodynamics of the vortex system. Furthermore, our results should apply, with little change, to thin film superconductors with strong pinning.

The model considered in our work is defined in Section II, where the method of calculation is also described. We then consider (subsection III A) the effects of an isolated columnar pin on the structure of the vortex liquid in the vicinity of the pin. This is done mainly for testing the systematics of our numerical method and also for determining appropriate values of the pinning potential to be used in subsequent calculations. We choose throughout the computations numerical parameter values appropriate for BSCCO. We determine the suitable choice of the discretization scale in order that our numerical method provides an accurate account of the density inhomogeneities produced by the trapping of a vortex at a pinning center. We also determine the range of values of the pinning potential strength for which nearly one vortex is trapped at a pinning center in the temperature range of interest. The strength of the pinning potential is kept fixed in this range in our subsequent work: pinning of multiple vortices at a pinning center is not considered because this is rarely observed in experiments. Next, in subsection III B, we consider the effects of two neighboring pinning centers on the liquid-state properties. An “effective interaction” between vortices trapped at the two pinning centers is obtained by calculating the free energy as a function of the separation between the pinning centers. This effective interaction is found to oscillate with distance. This study and the one-pin calculation mentioned above complement, in a sense, the analytic work of Ref. where the RY free-energy functional was used to analyze the structure and magnetization of a two-dimensional vortex liquid in the presence of strong pinning. However, we consider here a three-dimensional system with columnar pins, instead of a two-dimensional system with point pinning as in Ref.. Also, the numerical direct minimization method used in the present work is more accurate than the analytic variational method in the earlier study.

We next study (subsection III C) the freezing of the vortex liquid in the pure system. This is done primarily for checking the method against the results of earlier calculations . We find results in excellent agreement with those of earlier studies. Our new method also provides a very detailed and accurate account of the distribution of the density near a lattice point. We then proceed, in subsection III D, to consider the melting transition of interstitial vortices in a commensurate, triangular array of columnar pins. As discussed above, this transition provides a physical realization of three-dimensional melting in the presence of a commensurate periodic potential. Defining the concentration $`c`$ of pinning centers as $`cB_\varphi /B`$, we consider values of $`c`$ given by $`1/l^2`$ where $`l`$ is an integer. For small concentrations of pinning centers ($`l6`$), we find a first-order melting transition from a crystalline solid to an inhomogeneous liquid. As the pin concentration is increased, the transition temperature increases and the latent heat and the jump in the crystalline order parameter at the transition decrease. We find that this line of first-order transitions terminates at a critical point beyond which the thermodynamic transition is replaced by a sharp crossover. This critical point is a rare, experimentally realizable example of continuous melting in three dimensions. We show that a simple Landau theory provides a semi-quantitative understanding of most of our results. Some of our most salient results on the melting transition in the presence of a periodic pin array were summarized in a recent Letter. Here, we present many details which could not be included in that short paper. Section IV contains a summary of the main results and some concluding remarks.

## II Model and Methods

As explained in the introduction, we perform in this paper a numerical study using density functional theory, which involves, as its foundation , a model free energy functional appropriate to a system of pancake vortices in a layered superconductor. Density functional methods have long been used  with great success in the study of solidification phenomena in ordinary fluid systems. Although the theory is basically mean-field based, it works very accurately in the description of first order melting. We have ourselves performed extensive numerical studies of a density vs. disorder strength phase diagram of a hard sphere system in the presence of quenched disorder using a methodology quite similar to that employed in the present work. Density functional methods have also been successfully used  to study the melting of the vortex lattice in layered superconductors without pinning.

The starting point of our calculation is an expression for the free energy of the system, written as a functional of the time averaged local density. In our case the relevant density is $`\rho (i,𝐫)`$, the time averaged areal density of pancake vortices at point $`𝐫`$ on the $`i`$th layer. In the homogeneous vortex liquid state in the absence of pinning, this density is uniform and it is given in terms of the magnetic induction $`B`$ by $`\rho _0=B/\mathrm{\Phi }_0`$ where $`\mathrm{\Phi }_0`$ is the superconducting flux quantum. It is customary and convenient to introduce a length $`a_0`$ through the relation $`\pi a_0^2\rho _0=1`$. We will use $`a_0`$ as our standard unit of length in terms of which other lengths will be given and we will usually also normalize densities in terms of $`\rho _0`$.

We write the free energy functional in the form:

$$F[\rho ]=F_{RY}[\rho ]+F_p[\rho ].$$
(1)
The first term in the right-hand side of Eq.(1) is the free energy of the vortex system in the absence of pinning, while the second includes the pinning effects. Since the potential produced by a collection of straight columnar pins perpendicular to the layers is the same on every layer, the time-averaged density of vortices at any point $`𝐫`$ must be the same on all layers, i.e. $`\rho (i,𝐫)=\rho (𝐫)`$ for all $`i`$. Then, the free energy per layer corresponding to the first term in the right side of Eq.(1) may be written in an effectively two-dimensional form:

$`\beta F_{RY}[\rho ]={\displaystyle 𝑑𝐫\{\rho (𝐫)\mathrm{ln}(\rho (𝐫)/\rho _0)\delta \rho (𝐫)\}}`$ (2)
$`(1/2){\displaystyle 𝑑𝐫𝑑𝐫^{}\stackrel{~}{C}(|𝐫𝐫^{}|)\delta \rho (𝐫)\delta \rho (𝐫^{})},`$ (3)
where $`\beta `$ is the inverse temperature. We have defined $`\delta \rho (𝐫)\rho (𝐫)\rho _0`$ as the deviation of $`\rho (𝐫)`$ from $`\rho _0`$, the density of the uniform liquid and taken our zero of the free energy at its uniform liquid value. The function $`\stackrel{~}{C}(r)`$ is a static correlation function that contains all the required information about the interactions in the system. It is given by $`_nC(n,r)`$ where $`n`$ is the label denoting layer separation, $`r`$ is the the in-plane distance and $`C(n,r)`$ is the direct pair correlation function  of a layered liquid of pancake vortices with areal density $`\rho _0`$.

Strongly anisotropic layered superconductors can be described in terms of the Lawrence-Doniach Hamiltonian. Using that Hamiltonian as the starting point, it is relatively straightforward to calculate $`C(n,r)`$ and $`\stackrel{~}{C}(r)`$ using the hypernetted chain (HNC) approximation . Since the interaction term in Eq.(2) is of a convolution form, it is numerically most efficient to deal with it in wavevector space. This and the use of fast Fourier transform (FFT) methods reduces the computation of the interaction term in the free energy to a single sum. One begins with the expression for the Fourier transform of $`v(n,r)`$, the two-body vortex-vortex interaction, which, in the limit of vanishingly small Josephson coupling between the layers, is given by:

$$\beta v(𝐤)=\frac{2\pi \mathrm{\Gamma }\lambda ^2[k_{}^2+(4/d^2)\mathrm{sin}^2(k_zd/2)]}{k_{}^2[1+\lambda ^2k_{}^2+4(\lambda ^2/d^2)\mathrm{sin}^2(k_zd/2)]},$$
(4)
where $`k_z`$ and $`k_{}`$ are respectively the components of $`𝐤`$ perpendicular and parallel to the layer plane, $`d`$ is the layer spacing, and $`\lambda (T)`$ the penetration depth in the layer plane. The dimensionless quantity $`\mathrm{\Gamma }`$ which determines the strength of the interactions is given by:

$$\mathrm{\Gamma }=\beta d\mathrm{\Phi }_0^2/8\pi ^2\lambda ^2.$$
(5)
In coordinate space, $`v(0,r)`$ is repulsive and logarithmic in $`r`$ while $`v(n0,r)`$ is also logarithmic, but attractive and weaker than the intralayer potential by a factor of roughly $`(d/\lambda )e^{nd/\lambda }`$. Beginning with an interaction of this form, the HNC procedure of Ref. can be used to numerically compute $`C(𝐤)`$ for the appropriate values of the relevant parameters. The quantity $`\stackrel{~}{C}(k_{})`$, the two-dimensional Fourier transform of $`\stackrel{~}{C}(r)`$, is then obtained by setting $`k_z=0`$ in $`C(𝐤)`$.

The second term in Eq.(1) represents the contribution of pinning to the free energy per layer. It is given by:

$$\beta F_p[\rho ]=𝑑𝐫V_p(𝐫)\delta \rho (𝐫)$$
(6)
where $`V_p(𝐫)`$ is the dimensionless (normalized by $`k_BT`$) pinning potential at point $`𝐫`$. This quantity can be written as $`V_p(𝐫)=_jv_p(|𝐫𝐑_j|)`$, where the sum is over all pinning centers located at the points $`\{𝐑_j\}`$ on a plane, and $`v_p(r)`$ is the dimensionless form of the potential at $`𝐫`$ due to a pinning center at the origin. We take this potential to be of the truncated parabolic form:

$$v_p(r)=V_0[1(r/r_0)^2]\theta (r_0r),$$
(7)
where $`r_0`$ is the range of the pinning potential. We will write the dimensionless strength parameter $`V_0`$ as $`V_0=\alpha \mathrm{\Gamma }`$ and the quantity $`\alpha `$ will be chosen, as explained below, so that the pinning is strong enough to localize one vortex at a pinning center at the temperatures of interest, but not so strong that more than one vortex is bound to a pinning center.

In order to carry out numerical work, we have to discretize our system. We introduce for this purpose a computational triangular lattice of size $`L`$. On the sites of this lattice we define density variables $`\rho _i\rho (𝐫_i)v`$, where $`\rho (𝐫_i)`$ is the density at mesh point $`i`$ and $`v`$ the area of the unit cell in the computational lattice, proportional to the square of its lattice constant $`h`$. We have $`LNh`$, so that the computational lattice has $`N^2`$ sites. Periodic boundary conditions are used in all our calculations.

Our basic procedure is to minimize the free energy of the system given certain values of the relevant parameters and the appropriate initial conditions, that is, some initial set of values for the computational variable $`\rho _i`$. Finding the minima of the free energy is not at all trivial, since one is minimizing a function of a very large number of variables (we have used values of $`N`$ up to 2048 as discussed below) and these variables $`\rho _i`$ can take values differing by many orders of magnitude (particularly in the solid phase) and must also satisfy the nonnegativity constraint. This precludes the use of many standard minimization methods. The procedure we use here is the same as that originally employed in the hard sphere problem with the important difference that the calculations involving the interacting term are performed in wavevector space. This is for two reasons: first, it turns out to be much more efficient, since the time used by Fourier transforming back and forth using efficient FFT routines turns out to be negligible compared with the time saved by not having to perform a double sum over the large computational lattice. Second, the direct correlation function is more conveniently computed in any case in terms of the wavevector variable. The procedure we use incorporates the nonnegativity constraints and is insensitive to the large range of the variables, but requires a large number of iterations for convergence. The efficiency of the FFT method, however, still allows us to use the sizes required for the problem.

The minimization procedure finds a local minimum of the free energy. The minimum located depends on the initial values chosen for the set of variables $`\rho _i`$. The appropriate choices in each case are discussed in the next section. Generally speaking, nearly uniform initial conditions lead to liquid minima while ordered initial conditions with the proper symmetry lead to crystalline states.

## III Results

### A General and one pin

In this Section we present and discuss our numerical results. In principle, these could be given in terms of a minimal set of dimensionless parameters. However, it is more appropriate, in view of this paper’s objectives as discussed in the Introduction, to present the results in terms of physical parameters with dimensions. This is the course we have taken. The values of the material parameters that we use here have been therefore chosen as appropriate to BSCCO. These parameters are the penetration depth $`\lambda `$ and the interplane distance $`d`$, which together with the temperature and fundamental constants determine $`\mathrm{\Gamma }`$. We set $`d=15.26\AA `$, $`\lambda (T=0)=1500\AA `$, and assume a two-fluid temperature dependence of $`\lambda (T)`$ with $`T_c(H=0)=85K`$. For these parameter values, the dimensionless quantity $`\mathrm{\Gamma }`$ 2650/($`T`$ in Kelvins) at low temperatures where the $`T`$-dependence of $`\lambda `$ is negligible. We study temperatures and fields in the region where the melting transition of the vortex lattice is expected to occur. The strength of the pinning potential is described by the parameter $`\alpha V_0/\mathrm{\Gamma }`$, as introduced above. We fix the range $`r_0`$ to $`r_0=0.1a_0`$, which corresponds to about $`55\AA `$ for $`B`$ = 2kG.

To test our procedures and to find out more about the parameter range to be studied and the system sizes required, we begin by considering the simple case of an isolated pinning center in a vortex system in the liquid state. We place this pinning center in the center of the computational lattice. Since periodic boundary conditions are used, this amounts to considering a periodic array of pinning centers with spacing equal to the size $`L`$ of the computational cell. As discussed below, the values of $`L`$ used in these studies are sufficiently large, so that the behavior near a pinning center is not affected by the presence of its periodic images. We then choose the initial configuration of the variables $`\rho _i`$ as one vortex located at the pinning center and uniform density everywhere else, with the average density consistent with $`\rho _0`$. We then perform the minimization of the free energy as described above. The main issues here are the determination of the appropriate values of the pinning strength, and, from a technical standpoint, finding the system sizes, and the values of $`h/a_0`$ required. A smaller value of $`h/a_0`$ implies higher spatial resolution in describing the variation of the local density, but at constant $`N`$ this amounts to a reduction in the size of the system being studied. One must therefore strike a balance.

We have performed this procedure for fields $`B`$ = 2kG and 3kG and at several temperatures in the range of interest ($`1522`$K) at those fields. We have considered values of $`N`$ ranging from 128 to 2048 with $`a_0/h`$ from a maximum of 80 down to values of order unity. Representative results of this rather extensive study are shown in Fig. 1. In this Figure we plot the density variable $`\rho `$ (normalized by $`\rho _0`$) as a function of the dimensionless distance $`r`$ from the pin, measured in units of $`a_0`$. The data in the Figure are all taken at a temperature $`T`$ = 20.0K and a field $`B`$ = 2kG, with the parameter $`\alpha `$ set at 0.06. Results are shown for two cases: $`N`$ = 512 with $`a_0/h=20`$ and $`N`$ = 1024, $`a_0/h`$ = 40. The scaling of $`a_0/h`$ with $`N`$ ensures that we are comparing systems containing the same number of vortices. The two results are very close to each other. We have also obtained results for $`N`$ = 2048, $`a_0/h`$ = 80, which are completely indistinguishable from those at $`N`$ = 1024, so that if we had plotted them they would not be visible: the two plots at $`N`$ = 1024 and $`N`$ = 2048 would be on top of each other. This and similar data obtained at other temperatures and fields tell us the range of values of $`N`$ and $`a_0/h`$ needed to obtain high quality data. The results subsequently presented in this and the next subsection are obtained at $`N`$ = 1024 and $`a_0/h=40`$

The high peak at small $`r`$ in Fig. 1 represents the large vortex density at the pinned site. This density then decays away in an oscillatory manner, as shown in the Figure, towards its long range limit, which is unity for our normalization. As expected , the behavior of $`\rho (r)/\rho _0`$ outside the range of the pinning potential is very similar to that of the radial distribution function  of the unpinned vortex liquid at the temperature and magnetic induction being considered. Thus, the medium and long range behavior of the density is reasonable. The behavior of $`\rho (r)`$ at very short distances, inside the pinning range, is shown in Fig.2. One can see in this Figure how the results are well fit, as expected, by an exponential form $`e^{v_p(r)}`$, where the pinning potential in units of $`k_BT`$ is given by Eqn. (7). This confirms that our computational mesh is sufficiently fine even at these very small ranges.

Next, we must find the appropriate values of the pinning strength parameter $`\alpha `$, as introduced above. We wish to consider here, of course, the case of strong pinning, but nevertheless $`\alpha `$ should remain small enough so that the total amount of flux pinned at each site remains on the average below one superconducting flux quantum. To choose the appropriate value, we studied the average number of vortices pinned as a function of $`\alpha `$. Sample results are shown in Fig. 3, where the number of vortices pinned at the site (obtained by integrating the density over the pinning range) is displayed as a function of $`\alpha `$. The data shown in this Figure were taken at $`B`$ = 2kG and $`T`$ = 18.0K; data at other relevant fields and temperatures are very similar. We see that the number of pinned vortices rises very rapidly with $`\alpha `$ until a very marked kink occurs, at about the value when one vortex is pinned on the average. It is clear that one should use a value of $`\alpha `$ just below the kink in the curve, and in this study we have used the values $`\alpha `$ = 0.05, and $`\alpha `$ = 0.06 (the smaller value was used in calculations at lower temperatures). For $`T`$ in the range of interest here ($`1522`$K), these values of $`\alpha `$ correspond to $`V_079`$. This is consistent with the results of the two-dimensional study of Ref. where it was found that the average number of vortices trapped at a pinning center decreases sharply below one as the dimensionless pinning strength $`V_0`$ becomes lower than about 8.

### B Two pinning centers

Having in the previous subsection determined the properties of the density profile when the pins are, in effect, very far apart, we consider now the case where there are two pinning centers separated by a smaller distance $`d`$. With $`\alpha `$ set so that each center pins nearly one vortex, we expect to have, when the two centers are not too far apart, interactive effects as a function of $`d`$, as the density oscillations emanating from each pinning center (see Fig. 1) must be distorted to match each other. We perform this calculation by placing the two pinning centers symmetrically around the center of the computational triangular lattice, on the longer diagonal. For the initial conditions, we place one vortex on each of the two pinning sites, and a uniform density on the remaining computational sites, consistent with the average density being $`\rho _0`$. We have performed this study at fields of 2kG and 3kG and at several temperatures. Results at $`B`$ = 2kG and $`T`$ = 20K, with $`\alpha `$ = 0.06 are shown in Fig. 4. These results were again obtained at $`N`$ = 1024, $`a_0/h`$ = 40. In the main plot we give two examples of the density profile as a function of distance. Only one half of the density distribution is plotted, with the origin corresponding to the center of the lattice and the horizontal axis representing distances along the diagonal. The distribution is then symmetric about the origin, and the distance from the center of the first peak to the plot’s origin is half the interpin distance. The solid curve corresponds to the case where $`d`$ (in units of $`a_0`$) is 1.95 and the dashed curve to the case where it equals 2.60. One can see at larger values of $`r`$, away from both pins, an oscillatory decay similar to that in Fig. 1. The behavior in the interpin region near the plot origin is more complicated. At shorter interpin distances (as in the solid curve), the density is very small between the two centers, but when that distance is increased to over twice $`a_0`$ (see the dashed curve), it becomes possible to have a peak between the two pinning centers, and as $`d`$ is further increased, additional intersite peaks appear as well.

When the two centers are close enough to interact, the free energy will obviously depend on whether the oscillations in the density profiles corresponding to the two centers “lock” or not. This implies  that the free energy of the system should be an oscillatory function of interpin distance. To verify this, we have evaluated the free energy as a function of interpin distance $`d`$. Results are shown in the corner inset of Fig. 4. The dashed horizontal line is the result for the case where $`d`$ is very large, that is, twice the value for a single pin (recall that the zero of free energy is taken to be the at the uniform liquid state, see Eq.(2)). This value is $`\beta F=4.680`$ for the case plotted. The solid curve shows the behavior of $`\beta F`$ as a function of interpin distance and it clearly displays the oscillatory behavior of this quantity. As found in Ref., the free energy has minima at interpin distances approximately corresponding to the position of the maxima of the single pin density profile shown in Fig. 1. This reflects that it is easier, for those distances, to lock the oscillations corresponding to the two centers. The two pin distances corresponding to the two curves shown in the main plot of Fig. 4 were chosen so that one corresponds to a free energy minimum (solid line) while the other (dashed curve) corresponds to a nearby maximum. The higher value of the free energy in the latter case is due to greater difficulty in matching the two profiles in this case. This difficulty is reflected in the smaller height of the first peak of $`\rho (r)`$ outside the range of the pinning potential and the appearance of a small peak of $`\rho (r)`$ near $`r/a_00.6`$. The oscillatory behavior of the free energy as a function of $`d`$ implies an oscillatory dependence of the magnetization of the vortex liquid on the applied magnetic field when a periodic array of pinning centers is present. In particular, the reversible magnetization in the liquid state is expected to show minima near certain integral values of $`B/B_\varphi `$.

The integrated density inside the range of a pinning center remains close to unity when the pin separation $`d`$ is greater than $`2a_0`$. As the value of $`d`$ is decreased below $`2a_0`$, the integrated density begins to decrease and becomes substantially lower than unity for $`d<1.5a_0`$. Thus, the simultaneous occupation of two pining centers by two vortices is likely if the centers are far apart (in units of $`a_0`$), but unlikely only when the two pinning centers are separated by distances less than about $`1.5a_0`$. This is consistent with decoration experiments  which show that nearly all pinning sites are occupied by vortices when the number of pinning sites is smaller than the number of vortices ($`B>B_\varphi `$).

### C Melting in the pure system

As a preliminary step in our study of the effect of an array of pinning centers on vortex lattice melting, we have carried out calculations of the melting transition of the pure system for $`B`$ = 2kG and 3kG. This is in order to determine the “clean” limit of our subsequent results. In addition, our numerical solutions in the pure limit can be compared with those obtained from variational treatments  of the same RY free-energy functional in which the density distribution in the crystalline phase was expressed in terms of the Fourier components of the density at a few small reciprocal lattice vectors.

The computational cell used in our pure limit calculations is one triangular-lattice unit cell with lattice constant $`L`$. The spacing $`h`$ of the computational grid is chosen to have the values $`L/N`$ with $`N`$ = 16, 32, 64 and 128. Crystalline minima of the free energy are obtained by running the minimization routine with initial states that have a sharp peak of the density at the center of the computational cell. At sufficiently low temperatures such minima are found for a range of values of $`L`$. Typical results obtained for $`B`$ = 2kG, $`L=1.9884a_0`$, and two values of $`N`$ ($`N=16`$ and $`N=128`$) are shown in Fig. 5 where the dimensionless free energy $`\beta F`$ of one unit cell of the vortex crystal is plotted as a function of the temperature $`T`$. We find that the free energies of the crystal obtained for all the values of $`N`$ listed above are essentially the same, as exemplified by the data shown for two values of $`N`$, indicating that the effects of discretization are minimal provided that $`hL/160.125a_0`$.

The equilibrium value $`L_0`$ of the lattice parameter is determined by finding the value of $`L`$ that minimizes the free energy at a given $`B`$ and $`T`$. The dependence of $`\beta F`$ on the value of $`L`$ for $`B`$ = 2kG, $`T`$ = 18.5K is shown in Fig. 6. The value of $`L_0`$ is found to be close to 1.988$`a_0`$, which is slightly higher than the spacing $`\sqrt{2\pi /\sqrt{3}}a_0`$ of a triangular lattice of density $`\rho _0`$. This reflects the well-known result  that the density of a vortex lattice increases slightly at melting. The transition temperature is determined from the zero-crossing of the free energy of the crystalline state, calculated for $`L=L_0`$, as a function of $`T`$, as illustrated in Fig. 5, which shows the results of computations performed at $`L=L_0`$. For $`B`$ = 2kG, the melting temperature is then $`T_c`$ = 18.45K. This value of $`T_c`$ is slightly higher than that obtained variationally . This is expected: the free energy of the crystal obtained in our unconstrained minimization should be lower than that obtained in calculations where the free energy is minimized with respect to a few parameters only.

The discontinuity in the entropy at the crystallization transition is obtained from the numerically calculated slope of the $`\beta F`$ versus $`T`$ curve at the transition temperature. The Fourier transform of the density distribution at the crystalline minimum obtained at the transition temperature yields the value of the jump in the crystalline order parameter $`m`$, defined as the magnitude of the Fourier component of the density at the shortest reciprocal lattice vector of the triangular lattice. At $`B`$ = 2kG the entropy change $`\mathrm{\Delta }s`$ per vortex is $`0.29k_B`$, and the jump in the order parameter $`m`$ is $`\mathrm{\Delta }m=0.52`$. Very similar results are obtained for $`B`$ = 3kG: $`T_c`$ = 15.10K, $`\mathrm{\Delta }s`$ = 0.28$`k_B`$, $`\mathrm{\Delta }m`$ = 0.52, $`L_0`$ = 1.985$`a_0`$. These are in close agreement with the results of earlier studies.

Our numerical method provides detailed and accurate information about the spatial distribution of the time averaged density in a unit cell of the vortex lattice. We find that the density function near the center of the unit cell is, to a good approximation, invariant under a rotation about an axis perpendicular to the layers and passing through the center of the unit cell. Fig. 7 shows a plot of the local density $`\rho (r)`$ (angularly averaged and normalized by the average density $`\rho _0`$) at the crystalline minimum obtained for $`B`$ = 2kG, $`T`$ = 18.5K, $`L`$ = 1.988$`a_0`$, as a function of the distance $`r`$ from the center of the unit cell. As shown in the figure, a gaussian fit to the data for small $`r`$ provides a good account of the dependence of the density on $`r`$ except at larger distances. The value of the Lindemann parameter $``$ at melting may be obtained approximately from the width of the gaussian fit, or more accurately from a numerical evaluation of the root-mean-square displacement

$$r^2=\frac{r^2\rho (𝐫)𝑑𝐫}{\rho (𝐫)𝑑𝐫},$$
(8)
where the integrals are over a lattice unit cell and $`𝐫`$ is the radius vector measured from the center of the unit cell. We find that the value of $``$ at melting is 0.26 for $`B`$ = 2kG and 0.25 for $`B`$ = 3 kG. These values, similar to those found in earlier work , are substantially larger than the typical values of $``$ in simple three-dimensional solids near melting. The large value of $``$ implies that the peak of the density at the center of the unit cell is not very sharp (see Fig. 7). This helps explain why relatively coarse values of the mesh size $`h`$ (e.g. $`h=L_0/16`$) are adequate for obtaining an accurate description of the density distribution in the crystalline state.

### D Periodic array of columnar pins

Having obtained these clean limit results in the previous subsection, we proceed now with our study of the effects of a commensurate, periodic array of pins on the vortex lattice melting transition. We consider a triangular lattice of pins with spacing equal to $`lL_0`$ where $`l`$ is an integer, and $`L_0`$, as defined above, is the equilibrium value of the spacing of the pure vortex lattice at its melting point for the value of $`B`$ being considered. Thus the pin concentration is $`c1/l^2`$. The computational cell used is one unit cell of the pin lattice (which contains $`l^2`$ unit cells of the vortex lattice) with periodic boundary conditions and one pin located at the center of one of the vortex lattice unit cells. The value of $`h`$ was fixed at $`L_0/64`$ in the calculations for $`B`$ = 2kG. We also carried out a few calculations for $`B`$ = 2kG using $`h=L_0/16`$. The results obtained for this larger value of $`h`$ were found to be essentially the same as those obtained for $`h=L_0/64`$. We therefore used $`h=L_0/16`$ in our calculations for $`B`$ = 3kG. We set the pin strength parameter, as mentioned above, to $`\alpha =0.06`$ in the calculations for $`B`$ = 2kG. A slightly smaller value, $`\alpha =0.05`$, was used in the calculations for $`B`$ = 3kG because the temperature range of interest is lower in this case and the strength required to pin slightly less than one vortex is somewhat smaller.

The crystalline and liquid minima of the free energy were located from “heating” and “cooling” runs, respectively. In a heating run, the crystalline minimum was first located by minimizing the free energy starting with an initial state in which the density distribution in each of the $`l^2`$ vortex lattice unit cells contained in the computational cell was that in one unit cell of the vortex lattice of the pure system at its melting point for the same value of $`B`$. The crystalline minimum so obtained was “followed” to higher temperatures by increasing the temperature in small steps and running the minimization program at each temperature with the minimum obtained at the previous step as the initial state. In a cooling run, a liquid minimum was first obtained at a relatively high temperature by minimizing the free energy with an initial state consisting of one vortex located at the pinning center and uniform density everywhere else, so that the average density was $`\rho _0`$. This minimum was then “followed” to lower temperatures as in the heating runs, but decreasing the temperature in small steps instead of increasing it. The liquid state is not homogeneous in the presence of pinning and its free energy is nonzero.

For both values of $`B`$ studied and relatively small values of $`c`$ ($`l`$ = 10, 8, 7 and 6), we found a range of temperatures over which both crystalline and liquid minima are locally stable. The crossing of the free energies of these two minima defines a first-order transition between crystalline and liquid states. Results for the transition temperature $`T_c`$ as a function of $`c`$ for $`B`$ = 2kG and 3kG are shown in Fig. 8. The presence of columnar pins is found to increase $`T_c`$. This should be expected: columnar pins suppress the disordering effects of the lateral wandering of vortex lines, and a commensurate periodic array of such pins clearly promotes crystallization. In other words, an external potential having the same symmetry as the crystal favors the crystalline state. The results for $`B`$ = 2kG and 3kG are quite similar, with $`T_c`$ for $`B`$ = 3kG reduced by approximately 3.4K for all these values of $`c`$. The discontinuities in the entropy $`s`$ and the order parameter $`m`$ decrease as $`c`$ increases  because pinning-induced order in the liquid increases with $`c`$.

Our results yield not only bulk quantities but also very detailed information on the density distribution of the vortices. This quantity is experimentally accessible in scanning tunneling microscopy (STM) and scanning Hall probe measurements. In Fig. 9, we show the variation of the local density $`\rho `$ along a line joining two neighboring pinning centers for the crystalline and liquid minima near the transition temperature for $`c=1/49`$ and $`B`$ = 2kG. The density profile in the liquid minimum can be viewed as a superposition of liquid-like profiles near individual pins (compare with Fig. 1, noting the vertical axis scale), with small-amplitude, damped oscillations about $`\rho =\rho _0`$. In contrast, the density in the crystalline state exhibits higher, sharper and asymmetric peaks at the lattice points, with the density rising more sharply on the side closer to the pin, particularly at smaller distances from the pin site. Similar plots for $`c=1/64`$ may be found in Ref.. These plots clearly bring out the obvious differences between the structures of the coexisting crystal and inhomogeneous liquid phases.

The density distributions at the liquid minima exhibit as expected  a six-fold angular modulation. This is illustrated in Fig. 10 where we have shown the average density at the first peak of the normalized local density near a pinning center (i.e. the density averaged over the region $`1.75a_0r1.9a_0`$ where $`r`$ is the distance from the pinning center) as a function of the angle measured from the line joining the pinning center to one of its nearest neighbors. The data shown are for the liquid minimum obtained for $`B`$ = 2kG, $`T`$ = 19.6K and $`c=1/49`$. A six-fold angular modulation of the density is clearly seen in the figure. The minima of the density occur on the lines that join the pinning center to its nearest neighbors. This is different from the behavior found in the crystalline minima where density maxima occur on the lines joining neighboring pinning sites (this can be seen at a different value of $`c`$ from inspection of Fig. 12 below).

The full two-dimensional density distributions at the liquid and crystalline minima obtained near the transition temperature for $`B`$ = 2kG and $`c=1/36`$ are shown as gray scale plots in Figs. 11 and 12, respectively. The plot for the liquid minimum exhibits the usual correlation “hole” around the pinning site at the center, and concentric “rings” of alternating high and low densities with six-fold angular modulation. The angular modulation at the first ring is less pronounced (and less obvious in a gray scale plot) than that depicted in Fig. 10 for $`c=1/49`$. The plot for the crystalline minimum illustrates how the detailed structure of the periodically arranged crystalline density peaks changes with the distance from the pinning site at the center. The ability of our numerical method to provide detailed information about the density distribution in highly inhomogeneous states is clearly illustrated in these figures as well as in Figs. 9 and 10.

The degree of order in the liquid state increases with the pin concentration $`c`$, thereby decreasing the difference between the crystalline and liquid minima. This has drastic consequences for the crystallization transition. The behavior we find for $`c>1/36`$ ($`l<6`$) for both values of $`B`$ is significantly different from that described above. For $`c=1/25`$ and $`c=1/16`$, the apparent minima obtained in heating and cooling runs have almost the same free energy, but somewhat different values of the order parameter $`m`$. We have shown in Fig. 13 plots of $`m`$ versus $`T`$ obtained from heating and cooling runs for $`B`$ = 2kG and $`c=1/25`$ (circles and squares). The small difference in the heating and cooling values of $`m`$ peaks at a temperature $`TT_x`$ 21.2K. In Fig. 14 we show the two corresponding density profiles obtained at $`T`$ = 21.2K. These plots are analogous to those in Fig. 9 for $`c=1/49`$. In sharp contrast to that case, (and also to the $`c`$ = 1/36 case shown in Figs. 11 and 12) the two profiles are now very similar, with the one obtained in the heating run exhibiting only a slightly higher degree of order, consistent with the higher value of $`m`$. This leads to the suspicion that the free energy at $`c`$ = 1/25 may have only one very “flat” minimum near $`T=T_x`$ under the conditions studied. When attempting to find a minimum, our numerical routine stops when the free energy gradient becomes smaller than a certain small convergence parameter. When a minimum is very flat, it may stop at slightly different configurations when approaching it from different directions in configuration space.

If this situation occurs, then the density configuration at the true minimum of the free energy should be better approximated by a linear combination of the density configurations found in heating and cooling runs. We therefore evaluated the free energy for a set of configurations $`\{\rho _i(x)\}`$ defined by

$$\rho _i(x)=x\rho _i^{(1)}+(1x)\rho _i^{(2)},$$
(9)
where $`\{\rho _i^{(1)}\}`$ and $`\{\rho _i^{(2)}\}`$ are the density configurations, at the same temperature, at the apparent minima obtained in heating and cooling runs, respectively. The mixing parameter $`x`$ is in the range $`0x1`$. If one then plots the free energy thus obtained either as a function of $`x`$ or, equivalently, as a function of $`m(x)xm^{(1)}+(1x)m^{(2)}`$, where $`m^{(1)}`$ and $`m^{(2)}`$ are the order parameters in the two configurations, one finds that it indeed has a minimum at $`x=x_00.5`$ at temperatures near the expected transition, for all higher concentrations, $`c1/25`$. An illustration, for the case $`B`$ = 2kG, $`T`$ = 21.2K, and $`c=1/25`$ discussed above, is provided by Fig. 15. If one attempts such a procedure at lower concentrations, on the other hand, the resulting plot turns out to have a maximum, rather than a minimum, near $`x=0.5`$.

This analysis shows that the suspicions mentioned above were correct and that for $`c1/25`$, only one minimum of the free energy exists at each temperature. The value of the free energy at this minimum is lower than those found in the heating and cooling runs. Thus, there is no first-order transition at $`c=1/25`$ or higher. The line of first-order transitions found for smaller values of $`c`$ ends at a critical point which lies between $`c=1/36`$ and $`c=1/25`$ at both values of the field considered.

At $`c>1/36`$, above the critical point, a crossover rather than a sharp transition characterizes the change from liquid-like to solid-like behavior. The crossover temperature can be conveniently defined from the numerically calculated temperature derivative of the “equilibrium” value, $`m(x_0)`$, of the order parameter. Plots of both $`m(x_0)`$ and its temperature derivative are shown in Fig. 13. The temperature at which the derivative of the order parameter peaks is obviously very close to the temperature $`T_x`$ defined earlier as that at which the difference between the order parameters obtained in heating and cooling runs peaks. The crossover temperature can therefore be identified with $`T_x`$. The sharpness of the crossover suggests that $`c=1/25`$, $`T=T_x`$ 21.2K is close to the critical point for $`B`$ = 2kG, as indicated by the arrow in Fig. 8. For $`c=1/16`$ the crossover is smoother. Our results for $`B`$ = 3kG are very similar to those at the lower field, with a similar value of the critical $`c`$ but lower crossover temperatures (this is obvious from Fig. 8), with $`T_x`$ 17.6K for $`c=1/25`$.

The crossover to the crystal state at $`c>1/36`$ may be correlated with the onset of percolation of vortices which are “localized” according to a density-based criterion. Localization of vortices is, strictly speaking, a dynamical phenomenon, but some information about localization may be obtained from the distribution of the time-averaged local density. The local density near a point where a vortex is localized should be significantly higher than the average density $`\rho _0`$. Therefore, a density-based criterion for localization may be obtained by demanding that the density near a point where a vortex is localized exceed a suitably chosen cutoff value $`\rho _c`$. Our results for the density modulation around an isolated pinning center suggest an appropriate choice for this cutoff. One can see, for example, in Fig. 1, that in the temperature range of interest the local density in the neighborhood of an isolated pinning center (but outside the range of its pinning potential) does not exceed $`3\rho _0`$ if the system is in the liquid state. This suggests that values of the local density $`\rho `$ less than $`3\rho _0`$ correspond to mobile vortices. We, therefore, take $`\rho _c=3\rho _0`$. We divide the computational cell into $`l^2`$ vortex-lattice unit cells and associate a localized vortex with a unit cell if the local density exceeds $`\rho _c`$ at some point inside that cell. We then examine whether the unit cells that contain localized vortices according to this criterion percolate across the sample.

All the vortex lattice unit cells in a crystalline minimum contain localized vortices, since the maximum value of $`\rho _i`$ at the lattice sites of a crystal always exceeds $`\rho _c`$. In contrast, only the vortex lattice unit cells in which pinning centers are located and, in some cases, the nearest neighbors of such unit cells, contain localized vortices in the coexisting liquid minimum at the crystallization transition for $`c1/36`$. Thus, for $`c1/36`$, the crystallization transition trivially coincides with a percolation of unit cells containing localized vortices. For $`c1/25`$, in the crossover region, we have found that the unit cells containing localized vortices do not percolate if the temperature is higher than the crossover temperature $`T_x`$ defined above, but percolation occurs below $`T_x`$. Typical results are shown in Figs. 16 and 17. In Fig. 16, we have shown the locations of the units cells containing localized vortices in the minimum obtained in the heating run for $`B`$ = 2kG, $`T`$ = 21.2K and $`c=1/25`$. These unit cells do not percolate across the $`5\times 5`$ computational cell, while the cells containing mobile vortices do. Since the degree of localization in the minimum obtained in the heating run is higher than that in the minimum obtained in the cooling run, no percolation would be obtained at this temperature if the cooling-run minimum or the equilibrium configuration $`\rho _i(x_0)`$ were used for finding the unit cells containing localized vortices. On the other hand, as shown in Fig. 17, the unit cells containing localized vortices do percolate across the sample in the “equilibrium” configuration obtained at the slightly lower $`T`$ = 21.1K. At this temperature, the heating run shows percolation, but the cooling run does not. Thus, percolation occurs at a temperature very close to the crossover temperature $`T_x`$ 21.2K. Very similar results were obtained for $`B`$ = 2kG, $`c=1/16`$, and $`B`$ = 3kG, $`c=1/25`$, indicating that this is a general condition. This result is physically reasonable: a system in which localized vortices percolate (and consequently, the mobile ones do not percolate) should behave as a “solid” at long length scales.

It is easy to see that the occurrence of a critical point in the phase diagrams of Fig. 8 does not contradict any fundamental principles. In the presence of commensurate periodic pinning, the liquid and the crystal have the same symmetry. Since the degree of order in the liquid increases with $`c`$, it is possible for the liquid and the crystal to become indistinguishable beyond a critical value of $`c`$. One can then go from one phase to the other without crossing a sharp phase boundary.

The basic features of the phase diagrams may be understood from a simple Landau theory. From well-known symmetry arguments one can write down a Landau expansion for $`F`$:

$$\beta F=\frac{1}{2}a_2m^2\frac{1}{3}a_3m^3+\frac{1}{4}a_4m^4gm,$$
(10)
where $`m`$ is our order parameter, the constants $`a_3`$ and $`a_4`$ are positive, and $`a_2`$ decreases with $`T`$. The “ordering field” $`g`$ is proportional to the pin concentration $`c`$. A simple analysis shows that this free energy leads to a first-order transition for $`g<g_c=a_3^3/27a_4^2`$ and a critical point at $`g=g_c,a_2=a_{2c}=a_3^2/3a_4`$. The transition temperature increases with the ordering field $`g`$ in agreement with the arguments previously discussed. The latent heat and the order parameter discontinuity $`\mathrm{\Delta }m`$ vanish as $`(g_cg)^{1/2}`$ as $`g`$ approaches $`g_c`$ from below. It was shown in Ref. that our data for $`\mathrm{\Delta }s`$ and $`\mathrm{\Delta }m`$ are indeed well-described by the form $`(c_cc)^{1/2}`$ with $`c_c`$ close to 1/25.

More quantitatively, it is possible to fit our $`\beta F`$ vs. $`m`$ data to the form Eqn. (10). The best fit for $`B`$ = 2kG, $`T`$ = 21.2K, $`c=1/25`$, is shown in Fig. 15. The fitting parameter values are $`a_2=70.71`$, $`a_3=234.87`$, $`a_4=261.71`$, and $`g=7.12`$. Using the values of $`a_3`$ and $`a_4`$ obtained from the fit, we get $`a_{2c}=70.26`$, and $`g_c=7.01`$. These values are very close to, but slightly lower than the best-fit values of $`a_2`$ and $`g`$, indicating that the critical point for $`B`$ = 2kG is, as we had already stated, very close to $`c=1/25`$, $`T`$ = 21.2K. This explains the sharpness of the crossover at $`c=1/25`$. The numerical results for $`B`$ = 3kG can be analyzed in the same way. A fit for $`T`$ = 17.6K, $`c=1/25`$ yields then values of $`a_{2c}`$ and $`g_c`$ which are less than 1% lower than the best-fit values of $`a_2`$ and $`g`$, respectively. Therefore the Landau free energy gives a good semi-quantitative account of the critical behavior of our density functional computations for both values of $`B`$. This strengthens our conclusions about the existence and location of the critical point.

## IV Summary and Discussion

We have used in this paper numerical minimization of a discretized free energy functional to study the effects of columnar pinning on the structure and thermodynamics of a system of pancake vortices in the mixed phase of highly anisotropic layered superconductors. The most salient result of our study is the existence of a critical point in the vortex-lattice melting phase diagram when a commensurate, periodic array of pinning centers is present. Our results show that the line of the melting transition in the $`Tc`$ plane, which is of course first-order for small values of the concentration of pinning centers, terminates at a critical point as the pin concentration is increased. Beyond this critical point, the transition is replaced by a crossover, with smooth behavior of the order parameter and other thermodynamic quantities. To our knowledge, this is the first quantitative theoretical prediction of a continuous melting transition in a three-dimensional system. This critical point should be experimentally accessible: the pin lattice spacing for $`B`$ = 2kG, $`c`$ = 1/25 should be $``$ 0.55 $`\mu `$m, close to the spacing of the radiation-induced pin array of Ref. . The same group recently showed that columnar pins can be created in BSCCO in a highly controlled manner. We, therefore, expect that the fabrication of bulk HTSC samples with a periodic array of columnar pins with appropriate spacing is technically feasible. Fabrication of such samples and experiments to verify our theoretical predictions would be most welcome.

We have shown that most of the features of our phase diagram can be understood from a simple Landau theory. The critical point found in our study is analogous to the liquid-gas critical point in mean-field theory. Fluctuations are expected to change this correspondence because the symmetry of our order parameter is different from that for the liquid-gas transition. Theoretical studies of the universality class of this critical point would be interesting. However, the location of the first order melting line and the existence and experimental accessibility of the critical point should be quantitatively described by our work. We have also shown that the smooth crossover from liquid-like to solid-like behavior beyond the critical point might be interpreted as a percolation threshold for localized vortices.

The one-pin results reported here provide useful information about the dependence of the average number of vortices trapped at a pinning center on the temperature and the strength of the pinning potential, while our results for the interaction between two neighboring pins illustrate the occurrence of interesting effects in the liquid state arising from the commensurability of the separation of the pins with the average inter-vortex separation. Finally, our method yields very detailed results for the density distribution in the system, which is accessible through STM and scanning Hall probe measurements.

As noted in Section II, the symmetry of the three-dimensional system considered here makes the calculations effectively two-dimensional, with the function $`\stackrel{~}{C}`$ playing the role of the two-dimensional direct pair correlation function. The direct pair correlation function of a two-dimensional vortex liquid  is quite similar to the function $`\stackrel{~}{C}`$ used in the present study. We, therefore, expect that most of the results obtained here should apply, at least qualitatively, to thin-film superconductors in the presence of strong pinning centers. As noted in Section III, some of our one- and two-pin results are indeed in good quantitative agreement with those obtained in Ref. for a two-dimensional system of vortices with strong pinning. This leads us to expect that the phase diagram obtained here for the vortex lattice melting transition in the presence of a periodic array of pins should apply, with no more than fairly minor quantitative changes, to the melting transition of a two-dimensional vortex lattice in thin-film superconductors with commensurate periodic pinning. Since periodic arrays of strong pinning centers have already been fabricated  in thin-film superconductors, our predictions can be readily tested in experiments. One should, however, keep in mind that the predictions of our mean-field-like density functional calculation are less reliable in two dimensions where the effects of fluctuations are stronger. The melting transition in a pure two-dimensional system without pinning can be continuous , whereas the mean-field prediction of first-order melting is always realized in pure three-dimensional systems. Our main result about the existence of a critical point in the phase diagram should apply to thin-film superconductors with commensurate periodic pinning if the system exhibits a first-order melting transition in the absence of pinning.

The melting transition of the lattice of interstitial vortices in the presence of a commensurate, periodic array of pinning centers provides a physical example of melting in the presence of an external periodic potential. Similar melting transitions are of interest in other systems such as atoms adsorbed on crystalline substrates, and colloidal particles in interfering laser fields and arrays of optical traps. Our method and results would be of relevance to these systems also.


text:
# Faraday patterns in Bose–Einstein condensates. Amplitude equation for rolls in the parametrically driven, damped Gross–Pitaevskii equation

## I Model

We consider spontaneous pattern formation in dilute Bose–Einstein condensates (BEC) whose interatomic $`s`$wave scattering length is periodically varied in time . This modulation can be achieved by different means like, e.g., the use of magnetic , electric , or light  fields. When damping is considered  such BEC can be described by the following driven, damped Gross–Pitaevskii (GP) equation,

$$_t\psi (𝐫,t)=\left(1i\gamma \right)\left[^2\psi +V\left(𝐫\right)\psi +\left|\psi \right|^2\psi \right]+\frac{a\left(t\right)a_0}{a_0}\left|\psi \right|^2\psi .$$
(1)
which has been written in appropriate normalized variables. Parameter $`\gamma `$ accounts for damping (values on the order $`0.010.1`$ seem to be appropriate for actual BECs ), $`a\left(t\right)`$ represents the instantaneous value of the interatomic $`s`$wave scattering length and $`a_0`$ represents its mean value. We shall assume the simplest case $`a\left(t\right)=a_0\left[1+2\alpha \mathrm{cos}\left(2\omega t\right)\right]`$.

In the following we shall consider ”pancake” (2D) or ”cigar” (1D) shaped BECs in which the trapping potential $`V\left(𝐫\right)`$ strongly confines the condensate in one direction or two directions, respectively, whilst along the other direction(s) it extends sufficiently as compared with the typical wavelength of the emerging pattern. If a parabolic trapping potential $`V\left(𝐫\right)=1+\frac{1}{2}\left(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2\right)`$ is used (the arbitrary offset is set to $`1`$ for mathematical convenience) the above situation can be fulfilled whenever $`\omega _z\omega \omega _x,\omega _y`$ for the 2D case and $`\omega _y,\omega _z\omega \omega _x`$ for the 1D case. In a first approximation these inequalities allow: (i) to neglect the confined direction(s) in the descripition of the BEC dynamics, and (ii) to approximate the potential by a constant. Both approximations come from the fact that the characteristic wavelength of the selected pattern turns out to be, at the same time, much smaller than the size of the weakly confined direction and much larger than the strongly confined direction. Hence we consider along the rest of this work the following GP equation :

$$_t\psi (x,t)=\left(1i\gamma \right)\left[_x^2\psi \psi +\left|\psi \right|^2\psi \right]+2\alpha \mathrm{cos}\left(2\omega t\right)\left|\psi \right|^2\psi .$$
(2)
Eq. (2) admits the following homogeneous state

$$\psi =\mathrm{exp}\left[i\left(\alpha /\omega \right)\mathrm{sin}\left(2\omega t\right)\right],$$
(3)
which, in the absence of modulation ($`\alpha =0`$) reduces to the BEC ground state $`\psi =1`$ (the chemical potential is null in this representation because of the choice of the offset in the trapping potential).

## II Linear stability analysis: the parametric resonance

We wish to know whether the spatially homogeneous external driving is able to induce a spontaneous spatial–symmetry breaking of (3). Fot that we perform next a linear stability analysis of (3) by adding a small perturbation to that solution in the form $`\psi =\mathrm{exp}\left[i\left(\alpha /\omega \right)\mathrm{sin}\left(2\omega t\right)\right]\left[1+w\left(t\right)\mathrm{cos}\left(kx\right)\right]`$. Substitution of this expression into Eq. (2) and linearization with respect to $`w`$ leads to the following coupled equations for $`u=Rew`$, and $`v=Imw`$:

$`u^{}\left(t\right)`$ $`=`$ $`\gamma \left(2+k^2\right)u+k^2v,`$ (4)
$`v^{}\left(t\right)`$ $`=`$ $`\gamma k^2v\left[2+k^2+4\alpha \mathrm{cos}\left(2\omega t\right)\right]u,`$ (5)
which can be combined to yield

$$u^{\prime \prime }\left(t\right)+2\gamma \left(1+k^2\right)u^{}\left(t\right)+\left[\left(1+\gamma ^2\right)\mathrm{\Omega }^2\left(k\right)+4k^2\alpha \mathrm{cos}\left(2\omega t\right)\right]u\left(t\right)=0,$$
(6)
where

$$\mathrm{\Omega }\left(k\right)=k\sqrt{2+k^2},$$
(7)
is the nonlinear dispersion relation for the perturbations.

Eq. (6) is a Mathieu equation with damping, analogous to that describing parametrically driven, damped pendula, and ubiquitous in the description of parametric forcing . Its solutions are well known which, according to Floquet’s Theorem, can be written as

$$u\left(t\right)=Ref\left(t\right)e^{\mu t},$$
(8)
where $`f`$ is a periodic complex function of period $`\pi /\omega `$ and $`\mu `$ is the so called Floquet exponent (in the case of the Mathieu equation $`\mu /i`$ is known as Mathieu characteristic exponent). The BEC ground state (3) will be unstable against perturbations with wavenumber $`k`$ whenever $`Re\mu >0`$. A general property of Eq. (6) is that $`Re\mu >0`$ within a series of resonance ”tongues” very much like the parametric resonances observed in liquids vibrated vertically . For small $`\gamma `$ and $`\alpha `$ (the cases we consider here) these tongues are located, as a function of $`\mathrm{\Omega }`$, around $`\mathrm{\Omega }\left(k_n\right)=n\omega `$, $`n=1,2,3,\mathrm{}`$ Within these tongues the BEC ground state undergoes a spontaneous spatial-symmetry breaking and perturbations with wavenumber $`k_n=\sqrt{\sqrt{1+n^2\omega ^2}1}`$ amplify. The Floquet exponents can be numerically determined by using standard mathematical methods . Useful analytical information can be obtained in the limit of weak damping $`\gamma 1`$. Taking into account that for $`\gamma =\alpha =0`$ the solution to (6) is of the form $`u\left(t\right)=ReUe^{i\mathrm{\Omega }t}`$, if we consider $`\alpha \left|\omega \mathrm{\Omega }\right|\gamma `$, a perturbative expression for $`u`$ can be obtained by allowing $`U`$ to be a slowly varying function of time with the result

$`u\left(t\right)`$ $`=`$ $`Re\left[U_+\mathrm{exp}\left(\lambda _+t\right)+U_{}\mathrm{exp}\left(\lambda _+t\right)\right]e^{i\mathrm{\Omega }t},`$ (9)
$`\lambda _\pm `$ $`=`$ $`i\left(\omega \mathrm{\Omega }\right)\gamma \left(1+k^2\right)\pm \sqrt{\left({\displaystyle \frac{\alpha k^2}{\mathrm{\Omega }}}\right)^2\left(\omega \mathrm{\Omega }\right)^2}.`$ (10)
Note that (9) can be written as (8) with $`f\left(t\right)=\mathrm{exp}\left(i2\omega t\right)`$, and

$$\mu =i\omega \gamma \left(1+k^2\right)+\sqrt{\left(\frac{\alpha k^2}{\mathrm{\Omega }}\right)^2\left(\omega \mathrm{\Omega }\right)^2}.$$
(11)
Finally the condition $`Re\mu 0`$ reads

$$\alpha \frac{\mathrm{\Omega }}{k^2}\sqrt{\gamma ^2\left(1+k^2\right)^2+\left(\omega \mathrm{\Omega }\right)^2},$$
(12)
where the equality defines the boundary (neutral stability line) of the first resonance tongue. Note that (12) indicates that, for fixed $`\omega `$, the threshold for pattern formation is minimum at $`\mathrm{\Omega }=\omega `$, i.e., for wavenumbers $`k=k_1=\sqrt{\sqrt{1+\omega ^2}1}`$. The minimum parametric driving is hence predicted to be $`\alpha _{\mathrm{min}}^{n=1}\left(\omega \right)=\gamma \omega \sqrt{1+\omega ^2}/\left(\sqrt{1+\omega ^2}1\right)`$. Following a similar analysis for the second resonance tongue ($`\mathrm{\Omega }2\omega `$) in the limit $`\alpha \left|\omega \mathrm{\Omega }\right|\gamma ^{1/2}`$, it follows that $`u\left(t\right)`$ can be written as (9) with

$$\lambda _\pm =i\left(2\omega \mathrm{\Omega }\right)\gamma \left(1+k^2\right)\pm \frac{1}{3}\sqrt{20\left(\frac{\alpha ^2k^4}{\mathrm{\Omega }^3}\right)^224\frac{\alpha ^2k^4}{\mathrm{\Omega }^3}\left(2\omega \mathrm{\Omega }\right)9\left(2\omega \mathrm{\Omega }\right)^2},$$
(13)
from which the second tongue ($`Re\lambda _+0`$) runs:

$$\alpha \frac{\mathrm{\Omega }^{3/2}}{k^2}\sqrt{\frac{3}{10}}\sqrt{\sqrt{5\gamma ^2\left(1+k^2\right)^2+9\left(2\omega \mathrm{\Omega }\right)^2}+2\left(2\omega \mathrm{\Omega }\right)}.$$
(14)
The minimum driving amplitude needed in this case, $`\alpha _{\mathrm{min}}^{n=2}`$, occurs at $`\mathrm{\Omega }=2\omega `$. In this case $`\alpha _{\mathrm{min}}^{n=2}\sqrt{\gamma }`$ hence the second resonance tongue is excited at larger driving amplitudes than the first one, for which $`\alpha _{\mathrm{min}}^{n=1}\gamma `$. Note in both cases that for $`\gamma =0`$ (no damping) the threshold for both (and in fact any) resonance tongues is $`\alpha =0`$ hence all are simultaneously excited at vanishingly small values of driving amplitude. Damping hence breaks this degeneracy and selects the first resonance tongue at low drivings.

## III The roll pattern: amplitude equation

In the following we study the simplest pattern supported by the GP equation (2) under the previously described parametric instability. We consider a roll pattern in the form

$$\psi (x,t)=e^{i\left(\alpha /\omega \right)\mathrm{sin}\left(2\omega t\right)}\left[1+\epsilon w_1\left(t\right)\mathrm{cos}\left(kx\right)+\epsilon ^2w_2(x,t)+\epsilon ^3w_3(x,t)+O\left(\epsilon ^4\right)\right],$$
(15)
where $`0<\epsilon 1`$ is an auxiliary small parameter. In order to deal with a small roll component (of order $`\epsilon `$) we assume that the amplitude of the parametric driving is also small, say $`\alpha =O\left(\epsilon ^2\right)`$ (that this is the proper scaling for $`\alpha `$ is justified a posteriori by the consistency of the final result). On the other hand for the parametric excitation to be effective Eq. (12) must be fulfilled. If we want to take into account both the resonance condition and the effect of damping we must impose $`\gamma =O\left(\epsilon ^2\right)`$ and $`\left(\omega \mathrm{\Omega }\right)=O\left(\epsilon ^2\right)`$. Summarizing we consider in the following the scalings

$$\omega =\mathrm{\Omega }+\epsilon ^2\omega _2,\alpha =\epsilon ^2\alpha _2,\gamma =\epsilon ^2\gamma _2.$$
(16)
Our goal is to find an equation for the roll complex amplitude $`w_1`$. This will be done using a standard multiple timescale technique . For this purpose we introduce a slow time

$$\tau =\epsilon ^2t,$$
(17)
and allow all coefficients of the expansion to depend formally both on $`t`$ and $`\tau `$:

$`w_1\left(t\right)`$ $`=`$ $`u_{11}(t,\tau )+iv_{11}(t,\tau ),`$ (18)
$`w_j(x,t)`$ $`=`$ $`u_j(x,t,\tau )+iv_j(x,t,\tau ).`$ (19)
Finally Eqs. (15), (16), (18) and (19) are introduced into the GP equation (2). After using the chain rule for differentiation $`_t_t+\epsilon ^2_\tau `$ and equating equal powers in $`\epsilon `$ an infinite hierarchy of differential equations is obtained.

### A Order $`\epsilon ^1`$

This is the first nontrivial order and reads

$`_tu_{11}(t,\tau )k^2v_{11}(t,\tau )`$ $`=`$ $`0,`$ (20)
$`_tv_{11}(t,\tau )+\left(2+k^2\right)u_{11}(t,\tau )`$ $`=`$ $`0,`$ (21)
whose solution can be written as

$`u_{11}(t,\tau )`$ $`=`$ $`\left[r\left(\tau \right)e^{i\mathrm{\Omega }t}+\overline{r}\left(\tau \right)e^{i\mathrm{\Omega }t}\right],`$ (22)
$`v_{11}(t,\tau )`$ $`=`$ $`i{\displaystyle \frac{\mathrm{\Omega }}{k^2}}\left[r\left(\tau \right)e^{i\mathrm{\Omega }t}\overline{r}\left(\tau \right)e^{i\mathrm{\Omega }t}\right],`$ (23)
where $`r\left(\tau \right)`$ stands for a yet arbitrary complex function of the slow time, and the overbar denotes complex conjugation.

### B Order $`\epsilon ^2`$

At this order we find

$`_tu_2(x,t,\tau )+_x^2v_2(x,t,\tau )`$ $`=`$ $`\left[1+\mathrm{cos}\left(2kx\right)\right]u_{11}(t,\tau )v_{11}(t,\tau ),`$ (24)
$`_tv_2(x,t,\tau )+\left(2_x^2\right)u_2(x,t,\tau )`$ $`=`$ $`\frac{1}{2}\left[1+\mathrm{cos}\left(2kx\right)\right]\left[3u_{11}^2(t,\tau )+v_{11}^2(t,\tau )\right],`$ (25)
whose solution can be written as

$`u_2(x,t,\tau )`$ $`=`$ $`u_{20}(t,\tau )+u_{22}(t,\tau )\mathrm{cos}\left(2kx\right),`$ (26)
$`v_2(x,t,\tau )`$ $`=`$ $`v_{20}(t,\tau )+v_{22}(t,\tau )\mathrm{cos}\left(2kx\right),`$ (27)
where each of the coefficients verifies

$`_tu_{20}(t,\tau )`$ $`=`$ $`i{\displaystyle \frac{\mathrm{\Omega }}{k^2}}\left[r^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\overline{r}^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\right],`$ (28)
$`_tv_{20}(t,\tau )`$ $`=`$ $`2u_{20}(t,\tau )2\left(2+k^2\right)\left|r\left(\tau \right)\right|^2+\left(k^21\right)\left[r^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\overline{r}^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\right],`$ (29)
$`_tu_{22}(t,\tau )`$ $`=`$ $`4k^2v_{22}(t,\tau )+i{\displaystyle \frac{\mathrm{\Omega }}{k^2}}\left[r^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\overline{r}^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\right],`$ (30)
$`_tv_{22}(t,\tau )`$ $`=`$ $`\left(2+4k^2\right)u_{22}(t,\tau )2\left(2+k^2\right)\left|r\left(\tau \right)\right|^2+\left(k^21\right)\left[r^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\overline{r}^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\right],`$ (31)
which solved yield

$`u_{20}(t,\tau )`$ $`=`$ $`2\left(1+k^2\right)\left|r\left(\tau \right)\right|^2+{\displaystyle \frac{1}{2k^2}}\left[r\left(\tau \right)e^{i\mathrm{\Omega }t}+\overline{r}\left(\tau \right)e^{i\mathrm{\Omega }t}\right]^2,`$ (32)
$`u_{22}(t,\tau )`$ $`=`$ $`{\displaystyle \frac{1}{2k^2}}\left[r\left(\tau \right)e^{i\mathrm{\Omega }t}+\overline{r}\left(\tau \right)e^{i\mathrm{\Omega }t}\right]^2,`$ (33)
$`v_{20}(t,\tau )`$ $`=`$ $`v_{200}\left(\tau \right)+{\displaystyle \frac{i}{2\mathrm{\Omega }}}\left[r^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\overline{r}^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\right],`$ (34)
$`v_{22}(t,\tau )`$ $`=`$ $`{\displaystyle \frac{i\mathrm{\Omega }}{2k^4}}\left[r^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\overline{r}^2\left(\tau \right)e^{i2\mathrm{\Omega }t}\right],`$ (35)
where $`v_{200}\left(\tau \right)`$ is a yet undetermined function of the slow time, which is not fixed by the present analysis. This information should be obtainable by extending the calculation up to higher orders of the expansion. Anyway, as will be seen below, the knowledge of $`v_{200}\left(\tau \right)`$ is not relevant for our purposes.

### C Order $`\epsilon ^3`$. The amplitude equation

This is the last order we will consider. It reads

$`_tu_3(x,t,\tau )+_x^2v_3(x,t,\tau )`$ $`=`$ $`f_u(t,\tau )\mathrm{cos}\left(kx\right)+g_u(t,\tau )\mathrm{cos}\left(3kx\right),`$ (36)
$`_tv_3(x,t,\tau )+\left(2_x^2\right)u_3(x,t,\tau )`$ $`=`$ $`f_v(t,\tau )\mathrm{cos}\left(kx\right)+g_v(t,\tau )\mathrm{cos}\left(3kx\right),`$ (37)
where

$`f_u`$ $`=`$ $`_\tau u_{11}\gamma _2\left(2+k^2\right)u_{11}+\left(2v_{20}+v_{22}\right)u_{11}+\left(2u_{20}+u_{22}\right)v_{11}+\frac{3}{4}\left(u_{11}^2+v_{11}^2\right)v_{11},`$ (38)
$`f_v`$ $`=`$ $`_\tau v_{11}\gamma _2k^2v_{11}4\alpha _2\mathrm{cos}\left(2\mathrm{\Omega }t+2\omega _2\tau \right)u_{11}3\left(2u_{20}+u_{22}\right)u_{11}\left(2v_{20}+v_{22}\right)v_{11}\frac{3}{4}\left(u_{11}^2+v_{11}^2\right)u_{11}`$ (39)
$`g_u`$ $`=`$ $`u_{11}v_{22}+u_{22}v_{11}+\frac{1}{4}\left(u_{11}^2+v_{11}^2\right)v_{11},`$ (40)
$`g_v`$ $`=`$ $`3u_{11}u_{22}v_{11}v_{22}\frac{1}{4}\left(u_{11}^2+v_{11}^2\right)u_{11}.`$ (41)
The solution to Eqs. (36) and (37) can be written as

$`u_3(x,t,\tau )`$ $`=`$ $`u_{31}(t,\tau )\mathrm{cos}\left(kx\right)+u_{33}(t,\tau )\mathrm{cos}\left(3kx\right),`$ (42)
$`v_3(x,t,\tau )`$ $`=`$ $`v_{31}(t,\tau )\mathrm{cos}\left(kx\right)+v_{33}(t,\tau )\mathrm{cos}\left(3kx\right),`$ (43)
where each of the coefficients verifies

$`_tu_{31}(t,\tau )k^2v_{31}(t,\tau )`$ $`=`$ $`f_u(t,\tau ),`$ (44)
$`_tv_{31}(t,\tau )+\left(2+k^2\right)u_{31}(t,\tau )`$ $`=`$ $`f_v(t,\tau ),`$ (45)
$`_tu_{33}(t,\tau )9k^2v_{33}(t,\tau )`$ $`=`$ $`g_u(t,\tau ),`$ (46)
$`_tv_{33}(t,\tau )+\left(2+9k^2\right)u_{33}(t,\tau )`$ $`=`$ $`g_v(t,\tau ).`$ (47)
The last two equations do not give us relevant information for our purposes. On the contrary Eqs. (44) and (45) determine the evolution equation for the complex amplitude of the roll $`r`$. This comes from the fact that these equations contain resonant terms which yield divergent solutions unless a solvability condition is imposed. This is clearly seen by writing Eqs. (44) and (45) in vector form:

$$_t\left[\begin{array}{c}u_{31}(t,\tau )\\ v_{31}(t,\tau )\end{array}\right]=\left[\begin{array}{cc}0& k^2\\ \left(2+k^2\right)& 0\end{array}\right]\left[\begin{array}{c}u_{31}(t,\tau )\\ v_{31}(t,\tau )\end{array}\right]+\left[\begin{array}{c}f_u(t,\tau )\\ f_v(t,\tau )\end{array}\right].$$
(48)
Upon diagonalizing this equation we obtain

$$_t\chi (t,\tau )=i\mathrm{\Omega }\chi (t,\tau )+f(t,\tau ),$$
(49)
(and its complex conjugate) where

$`\chi (t,\tau )`$ $`=`$ $`\left(2+k^2\right)u_{31}(t,\tau )+i\mathrm{\Omega }v_{31}(t,\tau ),`$ (50)
$`f(t,\tau )`$ $`=`$ $`\left(2+k^2\right)f_u(t,\tau )+i\mathrm{\Omega }f_v(t,\tau ).`$ (51)
Eq. (49) can only be solved if the driving term $`f`$ does not contain elements oscillating as $`\mathrm{exp}\left(i\mathrm{\Omega }t\right)`$. Upon substituting Eqs. (38) and (39) into Eq. (51), and making use of Eqs. (22), (23) and (32)–(35), the solvability condition is found to be

$$\frac{dr}{d\tau }=\gamma _2\left(1+k^2\right)r+i\frac{\alpha _2k^2}{\mathrm{\Omega }}e^{2i\omega _2\tau }\overline{r}i\frac{3+5k^2}{\mathrm{\Omega }}\left|r\right|^2r.$$
(52)
Finally we define a new roll complex amplitude

$$R\left(t\right)=\epsilon r\left(\tau \right)e^{i\omega _2\tau },$$
(53)
and turn back to the original parameters by undoing the scalings (16) with the result

$$\frac{dR}{dt}=\left[\gamma \left(1+k^2\right)+i\left(\omega \mathrm{\Omega }\right)\right]R+i\frac{\alpha k^2}{\mathrm{\Omega }}\overline{R}i\frac{3+5k^2}{\mathrm{\Omega }}\left|R\right|^2R.$$
(54)
Eq. (54) is the searched roll amplitude equation. It is a Landau equation with broken phase symmetry (note the presence of the linear term proportional to $`\overline{R}`$). The roll solution (15) can be written in terms of $`R`$ making use of Eqs. (18), (22), (23) and (53). To the leading order the roll reads

$`\psi (x,t)`$ $`=`$ $`e^{i\left(\alpha /\omega \right)\mathrm{sin}\left(2\omega t\right)}\left[1+w\left(t\right)\mathrm{cos}\left(kx\right)\right],`$ (55)
$`w\left(t\right)`$ $`=`$ $`\left(1\mathrm{\Omega }/k^2\right)R\left(t\right)e^{i\omega t}+\left(1+\mathrm{\Omega }/k^2\right)\overline{R}\left(t\right)e^{i\omega t}.`$ (56)
Let us finally note that a straightforward linear stability analysis of the trivial solution $`R=0`$ (hence $`w=0`$) of Eq. (54) yields the same neutral stability curve given in (12).

Acknowledgment

I thank very useful discussions with Kestutis Staliunas (PTB, Braunschweig, Germany) and Stefano Longhi (Politécnico di Milano, Milan, Italy). This work has been supported by the Spanish DGES through project PB98-0935-C03-02.