UFDC Home  myUFDC Home  Help 



Full Text  
STRUCTURAL TRANSITIONS OF THE VORTEX LATTICE IN ANISOTROPIC SUPERCONDUCTORS AND FINGERING INSTABILITY OF ELECTRON DROPLETS IN AN INHOMOGENEOUS MAGNETIC FIELD By ALEXIOS KLIRONOMOS A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2003 ACKNOWLEDGMENTS First and foremost I would like to thank Professor Alan T. Dorsey for his kind and motivating mentorship. His dedication, ability and encouragement have been inspiring in difficult times. I would also like to extend my gratitude to my immediate family and to my close and dear friends for their support and tolerance. TABLE OF CONTENTS ACKNOWLEDGMENTS .. ............. . ii ABSTRACT . . . . . . . .. v CHAPTERS 1 SUPERCONDUCTIVITY . .. 1.1 Conventional Superconductivity .......... 1.2 GinzburgLandau Theory .............. 1.3 London M odel .................... 1.4 Anisotropy and Structural Phase Transitions . 1.5 Elasticity and Melting of the Vortex Lattice . 2 GINZBURGLANDAU THEORY . ..... 2.1 Anisotropic GinzburgLandau Theory . . 2.2 Virial Theorem . .............. 3 MEAN FIELD THEORY . . 3.1 Solution of the First GinzburgLandau Equation 3.1.1 Landau Gauge . .......... 3.1.2 Symmetric Gauge . . . 3.2 Structure of the Vortex Lattice and Minimization c . . . I 2 . . . 7 .. . 13 . . 15 . . 18 . 20 .. . 20 .. . 25 .. . 28 . . . 28 .. . 29 . . . 3 1 f the Free Energy 35 4 ELASTICITY THEORY FOR THE VORTEX LATTICE .. ..... 5 STRUCTURAL PHASE TRANSITIONS . ..... 5.1 Nondispersive Elastic Moduli . ........... 5.2 Fluctuations ............... ....... 6 FINGERING INSTABILITY OF ELECTRON DROPLETS . 6.1 Manybody Wavefunction . ............. 6.2 Conformal Mapping Method ................ 6.3 Monte Carlo Simulation . .............. 7 CONCLUSION . . . . . . A EXTENSION OF THE VIRIAL THEOREM . .... A.1 A Useful Identity . ................ A.2 Generalized Abrikosov Identities .............. A.3 Free Energy, Magnetization, Gibbs Energy . ... . 47 . . 47 . . 50 . . 56 . . 58 . . 61 . . 64 . 77 . . 79 . . 79 . . 81 B CALCULATION OF /A, ( X *(l n )2'), \ 2*(II 2 I)22). . C M OM ENTS . . . . . . . C.1 Conformal Map z = f(w) = aw bwM ............... C.I.1 Exterior M om ents ........................ C.1.2 Interior M om ents ........................ C.2 Conformal Map z = f(w) .= aw M ... C.3 Conformal Map z = f(w) = rw Q/(w wo) + Q/(w wo) . C.4 Connection with Magnetic Field Inhomogeneity .. ......... C.4.1 Exterior M moments .. .................... C.4.2 Potential V (z) . .................... D CODE FOR THE NUMERICAL CALCULATIONS . D.1 driverl.m. .................. .. D.2 energy.m ... .. ................ D .3 m oduli.m . . . . . . D.4 driver2.m ....... ........... ... D.5 csquash.m ..................... E MONTE CARLO CODE .. ........ REFERENCES . . . . . . . 97 . . . 98 . . . 99 . . . .. 109 . . . 120 . . . 12 1 . 160 BIOGRAPHICAL SKETCH . . . . . 167 S 83 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy STRUCTURAL TRANSITIONS OF THE VORTEX LATTICE IN ANISOTROPIC SUPERCONDUCTORS AND FINGERING INSTABILITY OF ELECTRON DROPLETS IN AN INHOMOGENEOUS MAGNETIC FIELD By Alexios Klironomos Mi., 2003 Chairman: Alan T. Dorsey lM.. ir Department: Physics I present a derivation of the nondispersive elastic moduli for the vortex lattice within the anisotropic GinzburgLandau model. I derive an extension of the virial theorem for superconductivity for anisotropic superconductors, with the anisotropy arising from sd mixing or an anisotropic Fermi surface. The structural transition from rhombic to square vortex lattice is studied within this model along with the effects of thermal fluctuations on the structural transition. The reentrant transition from square to rhombic vortex lattice for high fields and the instability with respect to rigid rotations of the vortex lattice, predicted by calculations within the nonlocal London model, are also present in the anisotropic GinzburgLandau model. I also study the fingering of an electron droplet in a special Quantum Hall regime, where electrostatic forces are weak. Performing Monte Carlo simulations I study the growth and fingering of the electron droplet in an inhomogeneous magnetic field as the number of electrons is increased. I expand on recent theoretical results and find excellent agreement between my simulations and the theoretical predictions. CHAPTER 1 SUPERCONDUCTIVITY The fascinating phenomenon of superconductivity, although discovered in 1911 by H. Kamerlingh Onnes, has in the last two decades begun to play a very impor tant role in everyday life through technology. This technological revolution was even more apparent to scientists and laymen alike after the discovery of high temperature superconductors (HTSCs) by Bednorz and Miller in 1986. The theory of conventional superconductivity is well understood today. There exists an accepted microscopic theory, the so called BCS theory, which explains ad equately the behavior of conventional or isotropic, low temperature superconductors like Pb, Al, Nb, Nb3Sn. However, HTSCs have been proven difficult beasts to tame. No microscopic theory with the degree of success of BCS in describing the relevant phenomena exists as of today. There are many phenomenological models that describe the phenomena in both kinds of superconductors with various degrees of success. One very complicated and extremely important issue for technological applica tions is the morphology and properties of the vortex lattice in typeII superconductors. The vortex lattice is an arrangement of magnetic filaments inside the bulk of the su perconductor which enter under certain conditions when the superconductor is placed inside a magnetic field. The distinction between typeI and typeII superconductors will be clarified later in this introductory survey. In this work I will use one particular phenomenological model, the Ginzburg Landau theory, to study the elastic properties of the vortex lattice in anisotropic superconductorsin other words, the response of the vortex lattice to deformations. I also address the issue of structural phase transitions of the vortex lattice and the effects of thermal fluctuations. 2 The survey of superconductivity that follows will attempt to define and ex plain the basic concepts required and thus provide the theoretical background for the understanding of the main body of the research that is reported here. For a more detailed introduction to superconductivity, the reader can refer to various sources on the subject [13]. 1.1 Conventional Superconductivity The two distinctive phenomena associated with superconductivity are the dis appearance of resistance at some critical temperature T, (perfect conductivity) and the expulsion of magnetic fields from the bulk of a metal once it enters the supercon ducting phase (perfect diamagnetism or Meissner effect). Conversely, the application of a sufficiently strong magnetic field is found to destroy superconductivity. There exists a classification scheme for superconductors according to their be havior when an external magnetic field is applied. For typeI superconductors, the magnetic field cannot enter the bulk of the material until the socalled thermody namic critical field H,(T) is applied. The thermodynamic critical magnetic field is related to the free energy difference between the normal and the superconduct ing phases. For typeII superconductors the magnetic field penetration begins at a stronger field, called H2(T), and persists down to a smaller critical field, H, i(T), with H,2 > H, > H,1. The characteristic behavior of each type of superconductor when a magnetic field is applied is illustrated in figure (1.1). TypeII superconductors and especially the intermediate phase, also called the mixed phase or Shubnikov phase, are of interest to me. In the mixed phase the magnetic field enters the bulk of the material in the form of flux tubes, each carrying a quantum of magnetic flux he o= = 2.07 x 107 gausscm2. (1.1) 2e Type II Type I 0 /HI H, Hc2 H Figure 1.1: Flux penetration behavior for typeI and typeII superconductors with the same thermodynamic critical field H,. For most of the region between the two critical fields H~i(T) and Hc2(T), the flux tubes form a regular array that in a sense resembles the ( 'I.1l lattice of ordinary ( v I .,l The major difference is that the vortex lattice is essentially two dimensional as the flux tubes can be considered as straight lines which, however, can bend and deform. The existence of magnetic fields inside the superconductor and superconductiv ity may at first seem to be two mutually exclusive events. However, the magnetic tubes have a vortex of supercurrent swirling around them which screens the magnetic field from the rest of the bulk. Superconductivity is completely destroyed in the core of the magnetic flux tube where the magnetic field attains its largest value. The interactions between the vortex lines are repulsive. This interaction, bal anced by the external magnetic field pressure, stabilizes the vortex lattice. Alexei Abrikosov, in a seminal work [4], first studied and predicted the existence of the vortex lattice in conventional, typeII superconductors. There are two very important length scales associated with superconductivity. The penetration length A determines the electromagnetic response of the supercon ductor providing a measure of the spatial extent of the magnetic field in the bulk of the material. The coherence length C provides a measure of the size of the normal core. The coherence length, C is also the characteristic lengthscale associated with variations of the superconducting phase in more general situations. It can be shown that H,1 (o/47iA2 and Hc2 I 0/27r2. In addition, K = A/C is a dimensionless q(111.nl ilt approximately independent of temperature, called the GinzburgLandau pa rameter. For a typical elemental superconductor A 500 A and C 3000 A which gives a GinzburgLandau parameter K < 1. The value of the GinzburgLandau parameter facilitates the characterization of a material as typeI or typeII. Actually, c = 1//2 serves as the boundary value between the two regimes. For smaller values one has a typeI superconductor while for larger values a typeII. Typical HTSCs are extreme typeII superconductors with K of the order of 102. There are several methods that one can use to observe and study the static properties and the morphology of the vortex lattice. Small angle neutron scattering (SANS) [514], Bitter decoration [15, 16], scanning tunneling microscopy (STM) [17, 18], magnetooptical methods [19] and muon spin rotation methods (p/SR) [20] are the ones widely used. Of all these methods, small angle neutron scattering stands out as the most useful to make bulk measurements as the majority of the mentioned methods are limited to imaging the flux lines at the surface of the superconductor.1 In short, the SANS method takes advantage of the weak interaction of neutrons with mat 1One can probe the bulk of the superconductor using pSR. However, sophisticated modeling is necessary to interpret the data. ter through short range nuclear forces and dipole interactions with spatially varying magnetic fields. The neutrons are scattered from the bulk of the material revealing information about ( ',I.,l structure, the vortex lattice and magnetic ordering. More specifically, from the diffraction pattern one can deduce the symmetry of the vortex lattice. Measuring the scattered intensity enables one to determine the characteristic length scales A and and the field profile around the flux lines, which in turn allows the determination of the critical fields. Also, from the width of the diffraction peaks one can estimate the vortex lattice correlation lengths. The existence of this new superconducting phase and the phase transition which accompanies it leads naturally to the idea of the order parameter. Although proposed by Ginzburg and Landau [21] before the discovery of the BCS theory, the order pa rameter can be thought of as the wavefunction for the center of mass motion of the Cooper pairs. I will devote a separate section to the study of the GinzburgLandau theory. For typeII superconductors the mean field HT phase diagram is rather simple, as illustrated in figure (1.2). The characterization min .111 field" implies that fluctua tion effects are not considered, which is a drastic simplification. Although this dissertation deals with mainly static properties of the vortex lattice, I will address here some of the dynamic properties of the vortex system to unveil the complexity that emerges once one chooses to depart from the simple (but nonetheless important) approach of mean field theory. For an excellent extensive review the interested reader might refer to the classic work of Blatter and others [22]. For a more pedagogical introduction two other nice reviews are available [23, 24] When an external current density is applied to the vortex system, the flux lines will start to move under the influence of the Lorentz force F = j x B/c. In a homogeneous system the only resistance to the motion is provided by the friction force, which is proportional to the steady state velocity v of the vortex system. The H H,2(T) Normal Metal Mixed Phase H, (T) Meissner Phase 0 T T Figure 1.2: The HT phase diagram of a typeII superconductor. At high fields and temperatures the material is in the normal metallic phase. At low temperatures and fields the material is in the Meissner phase while for intermediate temperatures and fields the material is in the mixed phase. steady state velocity will then be: v = j x B/plC, where pt is a generalized friction coefficient. However, as the flux tubes move, an electric field E = B x v/c appears parallel to the external current density. With the current density j and the electric E in parallel there is power dissipated in the system which implies that the fundamental property of dissipationfree current flow in a superconductor is lost. Disorder comes to the rescue, pinning the flux lines and restoring dissipationfree current flow. Of course there is a critical depinning current density je which is always bounded above by the depairing current density jo. The implications in a technological context are obvious. Another kind of disorder is thermal fluctuations. Quenched disorder as de scribed above is static while thermal fluctuations are dynamic. Considering thermal fluctuations opens up the possibility for melting of the vortex solid into a vortex liq uid or the emergence of a glass phase. It is apparent that quenched disorder and thermal fluctuations are competing forms of disorder. Thermally activated jumps over the pinning barriers can cause the vortices to move, a phenomenon called flux creep, while thermal fluctuations of a single vortex line lead to a dynamical sampling and averaging of the confining potential arising from disorder, which in turn reduce the critical current density j,. All these effects modify the mean field phase diagram considerably and at the same time enrich and complicate the 1ir, i. involved. In high temperature superconductors the effects of the fluctuations are more pro nounced, due to weaker pinning of the vortices and enhanced thermal fluctuations. This becomes immediately apparent since the typical elemental superconductor has a critical temperature T,c 10 K compared to T,c 100 K of a HTSC. As a conse quence larger parts of the phase diagram contain melted or depinned vortex lattice, and also new exotic phases, such as entangled or disentangled vortex liquids, along with a zoo of vortex glasses, are possible. As the reader can infer from this limited introduction, the 1.li,i. associated with the study of superconductivity and the multitude of effects and phenomena it encompasses is difficult but the research is extremely fascinating and rewarding at the same time. 1.2 GinzburgLandau Theory The GinzburgLandau theory for superconductivity [21] is an extension of Lan dau's theory of second order phase transitions to a spatially varying complex order parameter 1(r). Originally proposed well before the discovery of BCS theory it has been proven an indispensable tool for the phenomenological treatment of a variety of problems in many fields. Although conceptually simple, it is an elegant theory of considerable predictive power. The GinzburgLandau theory has initially a somewhat limited region of applica bility. It is valid near the transition temperature T,. However, this restriction can be relaxed using various arguments or techniques, or even ignored for the sake of making qualitative predictions. Two fields are necessary for the formulation of the theory. One is the complex order parameter T(r) = f(r)e4(r) and the other is the vector potential A(r), related to the local magnetic field B(r) = V x A(r). The order parameter is zero at the transition point and has a nonzero value when the material is in the superconducting state. Based on the continuity of the change of state in a phase transition of the second kind, an expansion of the free energy in terms of the order parameter is feasible. The GinzburgLandau free energy functional F[I, A] is built using symmetry arguments. Odd powers of T are immediately excluded as the free energy density (an observable) should be invariant under a transformation of the form ei "(x)4. Another simplification one makes is the assumption that the crystal has cubic sym metry. The spatial variation of the order parameter then dictates the inclusion of gradient terms in the free energy functional. First order derivatives and second order derivatives of the form 02/XzizXk are excluded because these are essentially surface terms. The lowest order acceptable gradient term is of the form IV I'2, which has to be made gauge invariant by combining it with the vector potential resulting in Infl 2 = (iV e*A/hc)l2. Taking into account the above considerations, the GinzburgLandau free energy functional can be written now in the following form: F[I, A] = F, + dxL a12 1ip4 2 h+ B12 (1.2) 2 2m* 87 J Here F, denotes the free energy of the normal state and b is a positive coefficient without any temperature dependence. On the other hand a is a function of the temperature given by a = ao(T T,), (1.3) with the constant ao being positive. The constants e* and m* can be identified with the charge and the mass of a Cooper pair, namely e* = 2e and m* = 2m, where e and m are respectively the charge and the mass of the bare electron. To elucidate a few fundamental concepts I consider first the case of a homoge neous superconductor with no external field. The order parameter is then independent of the coordinates. The minimum of the free energy density occurs when l (T,~ T). (1.4) b One sees that the square of the order parameter, which is proportional to the super conducting electron density, decreases linearly to zero at the transition point. The difference in the free energies of the superconductor and normal metal is found to be 2 F F =V (T T)2. (1.5) The difference in the energies can be equated to VH,2/87, which is the energy density for the critical magnetic field. This immediately yields the temperature dependence of the critical field near the transition point H, (T T). (1.6) Differentiating the difference of the free energies with respect to the temperature one obtains the difference in the entropies and finally the jump of the specific heat at the transition point Cs c, = v (1.7) b The reader should be cautioned that such trivial manipulations are only possible in the simplest of cases. The solution of the GinzburgLandau equations, the equations of motion derived from the GinzburgLandau free energy, is a difficult task as will be made apparent shortly. An interesting observation is that only the modulus of the order parameter f(r) = I(r)l and the gauge invariant supervelocity Q = (hMo/27m*)Vq(r) e*A(r)/m*c enter the GinzburgLandau theory, but not T and A separately [25]. There is a completely equivalent formulation of the GinzburgLandau in terms of Q and f, implying that since only real and gauge invariant quantities are involved, any calculations and results (even approximations) should contain gauge invariant equations. One can derive the equations of motion for the two fields T and A minimizing the GinzburgLandau free energy functional with respect to I* and A. One obtains h i2 2 20 iV ?A + a+ b112 = 0 (1.8) 2m* ch ) ieh =h(f*Vf IV\*). (1.9) 2m* In the process of deriving the equations of motion a boundary condition has to be imposed to ensure that the surface integrals in the variation 6F are zero. With n being the normal vector to the surface of the superconductor the boundary condition has the form n iV A 0. (1.10) \hC ) An immediate result is that the normal component of the supercurrent vanishes at the boundary as expected: n j = 0, since one assumes that no current is fed into the material. It can be shown that the characteristic lengths and A have the following forms within the GinzburgLandau theory: (T) (T T), (1.11) F m*c2b 1 A(T) M b (T T) (1.12) [47c*2e* OI In addition, quite generally, one can include an electric field and time dependence through the assumption ( ie* DF with D being the scalar potential and 7 being a constant [2628]. One then obtains the so called time dependent GinzburgLandau equations. At this point I have to stress that the validity criteria for the time dependent version are different from the simple GinzburgLandau equations. The interested reader should consult the appropriate literature on the subject. The theory can also be extended to superconductors with crystal symmetry other than cubic with the introduction of an anisotropic effective mass tensor that multiplies the gradient term in the free energy density [29, 30]. This subject will be discussed extensively along with other ways to account for anisotropy later in this work. The full solution of the GinzburgLandau equations, except in trivial cases, is possible only with the use of numerical methods [3134] since it is a nonlinear differ ential equation. However, useful approximations can be made which produce valuable and accurate results. One frequently "'aplvid is the assumption that the nonlinear term T112 in the first GinzburgLandau equation is small enough to be neglected. This assumption works well in many cases but should be used with appropriate care. One case for which this term requires special treatment is in the calculation of the dispersive elastic moduli of the vortex lattice, cll and c44, where it introduces un pll, .i .,l divergences that are removed as soon as this term is included and properly treated [35,36]. A well behaved solution, in the symmetric gauge, that represents the vortex lattice near Hc2 with one flux quantum per flux line was found by Brandt [3540]: '(x, y) (constant) exp [ 2 (x2 2+ )] 2 [( x) + iy )]. (1.14) Here the product is over the (periodic) flux line positions (x,, y,). Notice the similarity with the Laughlin wavefunction for a filled (v = 1) Landau level within the context of the Quantum Hall effect (see Chapter 6 later in this thesis). This particular solution can be used to calculate the energy of structural defects in the vortex lattice ;', 11] and also the elastic moduli [25, 3538, 42] of the vortex lattice. One needs to define one very important parameter here, the Abrikosov param eter PA. It is the ratio PA = (1.15) (12) 2' where the brackets denote spatial averaging. By means of the Schwartz inequality the range of the Abrikosov parameter is [1, oo). In the mixed phase the Abrikosov param eter is a geometrical constant, attaining different values for different types of vortex lattices. It is easily shown that for an isotropic superconductor the minimization of the GinzburgLandau free energy is equivalent to the minimization of the Abrikosov parameter AA. For example, for a triangular lattice PA 1.1659, while for a square lattice PA 1.1803, thus making the triangular lattice the lowest energy configura tion of the vortices in an isotropic superconductor. One notable detail is the small difference of the two values. Initially Abrikosov obtained the incorrect result, namely the square lattice being the equilibrium configuration. This small difference makes the existence of structural transitions very plausible. For an anisotropic superconductor the situation is more complicated and more interesting because the anisotropy opens up the possibility of phase transitions from the triangular phase to even more general rhombic phases and eventually to the square phase under certain conditions. 1.3 London Model The GinzburgLandau theory, being an expansion in powers of the order param eter which is assumed to be of small magnitude, is valid near H,2. Near Hj1 on the other hand, where the vortex core overlap is not significant one can put I1(r)l = 1 outside the vortex core. This is the starting point for the formulation of the London theory. Either starting from the simplified GinzburgLandau free energy density, or from first principles one obtains the free energy functional of the London theory F[B] = I dx B2 +2(V x B)2 (1.16) 87i J Here A = [m*c2/4 nse*2] is the London penetration depth, with ns being the density of the superconducting electrons. In this section B = B(r) is the local magnetic field. Minimizing (1.16) with respect to the local magnetic field B(r) using the fun damental relation V B = 0 and also adding the contribution of the vortex core in the form of 6 function singularities one obtains the London equation for an isolated vortex A2V2B + B = o6(r). (1.17) Here 4o is a unit vector along the flux line direction. The 6 function approximates the core region which is of finite extent (of the order of ). The solution of the equation is easy to obtain for a field directed along the z axis, (Do B ( K (1.18) where Ko(r/A) is the zeroth order modified Bessel function of an imaginary argument. The ..i mptotic forms of the solution are more useful: B =2( In +(constant) ,for < r < A, I ro irA _ A SB C/A for r >>A. 2 A2 2r The same method can be used for two parallel vortices directed along the z axis. One simply adds one more singularity in the right hand side of equation (1.17) and the solution is obtained in a similar manner as before. One obtains for the magnetic field of the ith vortex Go r T Bi 272 Ko (1.19) where ri is the position of the ith vortex. Substituting back in the London free energy density (1.16) after a short calculation one finds for the interaction energy per length of the two vortices Ui, = 4oB12/4r, where B12 is given by #o ri r2 B12= 2 Ko (l (1.20) Thus the interaction between two vortices is repulsive. It decreases exponentially at large distances and diverges logarithmically at short distances. The London theory has been successfully used in the investigation of the prop erties of the vortex lattice. Actually it is relatively easier to obtain the elastic moduli in the London limit than it is using the GinzburgLandau theory. What is usually done is that some property is calculated in the high field limit and in the low field limit and finally the appropriate matching is done in the intermediate region in which both London and GinzburgLandau theory fail. Although in this work the GinzburgLandau theory is the main tool used in the investigation, I use many results obtained from the London theory. 1.4 Anisotropy and Structural Phase Transitions There are many sources of anisotropy; all are found to modify fundamental properties such as the shape of the vortex cores, the interactions between the vortices, the elastic response and also the symmetry of the vortex lattice. The mechanisms are complicated and sometimes lead to significant deviations from the behavior of the isotropic superconductor. The effect that this work is mainly focused on is the structural transitions of the vortex lattice. The familiar arrangement of the vortices in a hexagonal lattice is not unique, as experiments have shown. Distorted rhombic and square lattices are pos sible. The possibility for alternate arrangements becomes apparent once one realizes that the reason behind the formation of the hexagonal lattice 2 is the isotropic inter action between the vortices. An additional indication would be the small difference in energy between the square and hexagonal lattice. One source of anisotropy is that which arises from the existence of an under lying crystal lattice of different symmetry than the one usually assumed, namely cubic. This condition is met in the cuprate family, which are layered, extreme typeII HTSCs. They exhibit uniaxial ( ,i I.,1 symmetry. A way to account for that is to introduce a different effective mass along one direction (uniaxial symmetry). In that case the superconductor is characterized by two penetration lengths, Aab for the ab (or basal) plane and A, along the caxis, and also two coherence lengths ,ab and &,.3 The anisotropy ratio F = Ac/Xab = Cab/Cc quantifies the strength of the anisotropy. 2For a fixed area (density), vortices on a hexagonal lattice are placed the farthest apart possible from each other. 3In many materials the situations is even more complicated, with Ana A'b. Figure 1.3: Vortex profiles for applied field along the z and b axis respectively. Typical values for a HTSC are F > 1 and Aab 1000A. The GinzburgLandau pa rameter is defined as K = Xb/ab. An immediate consequence is that the profiles of the vortices are different for fields along the z axis and along the ab plane, see figure (1.3). Also, the critical fields Hc1 and Hc2 become functions of the angle between the applied field and the caxis of the material. The analysis in general is difficult but there exists a scaling scheme that reduces the anisotropic problem to an isotropic one at the initial level of GinzburgLandau theory or London model [43]. Using that method, one essentially generalizes isotropic results to the anisotropic case. The flux line lattice in anisotropic materials can be distorted and depending on the presence of pinning, defects and other factors one may see the formation of domains of different flux line configuration and also transitions from one lattice sym metry to another. Structural phase transitions can be studied with the same method, by promoting the mass m* to a mass tensor, thus incorporating the anisotropy in the model. However, a significant failure of this approach is that cubic materials (e.g. VaSi) would not exhibit structural phase transitions, a fact which contradicts experi ments. The results of the research presented in this work can easily be generalized to an anisotropic mass tensor using the scaling scheme [43]. Another source of anisotropy is an anisotropic Fermi surface. This is typical in the borocarbide family (RNi2B2C, with R Y, Lu, Tm, Er, Ho, Dy), which are most likely conventional but highly anisotropic swave superconductors [44]. The anisotropic Fermi surface affects strongly the electromagnetic interactions in the su perconductor [45]. To account for an anisotropic Fermi surface, the London model has to be amended in order to include nonlocal corrections via the relation between the cur rent density j and the vector potential A. For tetragonal superconductors like the borocarbides the corrections (in Fourier space) have the following form [46]: 47 1 3 ,2 (m % A2rijlmklkm)aj. (1.21) Here aj = Aj + (Io0O/27, with 0 the phase of the order parameter and A the vector potential. The inverse mass tensor m?1 and the fourth rank tensor nijim which cou ples supercurrents with the crystal are evaluated as particular averages of the Fermi velocity over the Fermi surface and in that way are strongly dependent on the shape of the Fermi surface. More importantly, the coupling persists even for cubic materi als [47]. A qualitative explanation of the importance of nonlocality is that the current response of the superconductor must be averaged over the finite size of the Cooper pair, which is of order in contrast to the usual local electrodynamics. The London model is both successful and convenient, due to its simplicity. Among the goals of this work is to show that one obtains identical results employ ing the GinzburgLandau theory with the inclusion of an anisotropic term. For the description of conventional superconductors with anisotropic Fermi surface the addi tional term is purely phenomenological. For the anisotropy that arises from an un conventional coupling the additional term will be derived from a more general model in subhead 2.1. The anisotropy that is of particular interest to me is the one arising from the existence of an unconventional pairing or a multicomponent order parameter. It has been established that the pairing mechanism in the cuprates is of d wave nature [48 50], which results in a fourfold anisotropy of the gap and anisotropic vortex cores. The first direct observation of a well ordered square vortex lattice in a HTSC was achieved very recently [51]. Even more exotic systems have been investigated such as Sr2RuO4 with an un conventional two component order parameter [52]. Coexistence of phases of different pairing symmetry and related phase transitions have been observed for the heavy Fermion superconductors UPt3 and U1_xThxBei3. The interested reader can consult the nice review by Sigrist and Ueda [53] and also the review by Joynt and Taillefer [54]. A very recent work by Agterberg and Dodgson [55] generalizes the London model to unconventional superconductors with multiple superconducting phase transitions. Although the main body of this dissertation deals with the investigation of the effects of d,2_y2 pairing with an small admixture of s wave pairs in the condensate, it is reasonable to expect that similar results can be obtained with the same methodology for other pairing symmetries. In each case it would be necessary to find the form of the appropriate term that couples the coexisting order parameters. 1.5 Elasticity and Melting of the Vortex Lattice The structural phase transitions of the vortex lattice are closely connected to its elastic properties. The elastic response of the vortex lattice, quantified by the appro priate elastic moduli, is of great significance for the phase diagram of the superconduc tor. The theory of elasticity for the vortex lattice in an isotropic superconductor was derived by Brandt [3538] using the GinzburgLandau theory for the highinduction regime and the London model for low inductions. The theory was later confirmed by the results of Larkin and Ovchinnikov [56], who used the microscopic theory as their starting point. The theory has been generalized to anisotropic superconductors with much effort [42, 5761]. In these works the anisotropy was incorporated in the mass tensor. The elasticity and the structural phase transitions are relevant for the extremely interesting subject of the melting of the vortex lattice. Especially in HTSCs, with the enhanced role of the thermal fluctuations and the softness of the elastic moduli, the region over which a melted vortex lattice is expected is significant. One can define melting as the loss of long range correlations between the vor tices. The melted vortex lattice consists of highly mobile and flexible vortex lines. A fundamental difference in the behavior of HTSCs compared to conventional, lowT, superconductors is that at high fields, the well defined and sharp superconducting normal transition is replaced by a smooth crossover between the vortex liquid phase and the normal metallic phase. One criterion used to define the onset of melting is the socalled Lindemann criterion, according to which melting of the vortex lat tice is expected once the mean squared amplitude of the thermal displacement of the vortex lattice becomes comparable to the intervortex spacing (the lattice constant). Although this criterion is not rigorous, it provides a simple and adequate tool for the investigation of melting and the determination of the melting curve Bm(T). The elasticity enters in the calculation of the mean squared displacement (u2). For an excellent discussion of the issues pertaining to the melting of the vortex lattice and the relevant theory, the reader can consult the review by Blatter and others [22]. CHAPTER 2 GINZBURGLANDAU THEORY 2.1 Anisotropic GinzburgLandau Theory The isotropic GinzburgLandau theory has been used extensively in the study of isotropic, or swave superconductors. However, once it was realized that alternate pairing mechanisms could exist, new extended GinzburgLandau theories were pro posed and studied in the context of highTc superconductivity [6269]. Although the latter are built using symmetry arguments, microscopic derivations do exist [70]. For the highTc cuprates it is established that the major pairing channel is the d,2_y2 pairing. The dominant pairing mechanism is of dwave nature with a subdominant swave order parameter being induced by gradients of the dwave as will be shown promptly. The generalization of the GinzburgLandau theory to anisotropic superconductors in which the anisotropy is due to s and d wave mixing contains two fields, or order parameters. The corresponding GinzburgLandau free energy can be written in the following form, which includes an anisotropic term with the appropriate symmetry [54, 66, 67]: fGL = S 2 ? +1 + S 4 + d2 4 3 S 212 4 Sd2+ S2d*2) 8 (2.1) + 7(ns)(ns)* + d(nd)(nd)* + 7,[(ns)*(nd) (n,s)*(n,d) + cc] , where the brackets denote spatial averaging. Also I have defined the parameters 7, = h2/2mn, i = s, d, v and I = V/i e*A/hc is the covariant derivative as usual. From here on I set h = c = 1. The critical parameter is ad and it is generally assumed that Ts < Td, meaning that d is the dominant order parameter. The later treatment of s as a small correction arising from gradients of d is based on this assumption. The EulerLagrange equations of motion are obtained from equation (2.1) by minimizing with respect to d* and s*. One obtains (7di2 + ad)d + 2/d2 d + 31 2d + 2 2d* + 7, (f I s 0, (2.2) (7yI2 t+ ,)s + 2 11s 2S + 3d2 21 4 + 7(II _ In)d 0. (2.3) The boundary conditions obtained during the derivation of the first and second equa tion of motion are given by the formulas n {yd(Hd) + 7,[y(Iys) x(IIs)]} = 0, (2.4) n {%((Hs) + . [v(nIId) x(IId)]} 0, (2.5) where n is the unit vector normal to the surface of the superconductor. The second GinzburgLandau equation, which is the equation for the current, is obtained by minimizing the GinzburgLandau free energy with respect to the vector potential A. After a short calculation one obtains j e [s*(ns) + s(Hs)*] + 7d[d*(Hd) + d(Hd)*] + yY[s*(fd) + d(Iys)*] 1 (2.6) xY[s*(n,d) + d(n,s)*] + cc . The solution of the GinzburgLandau equations is a very complicated problem even for the isotropic superconductor which can be described in terms of only one order parameter. The additional complications of having two order parameters can be avoided with a drastic simplification: the reduction of the twofield theory to an effective onefield theory. In order to derive the single component free energy I will follow Affleck and others [71] and integrate out the s field. This task is achieved by making the substitution s > (',/as) (UI I ) d. (2.7) This approximate equation is derived from the equation of motion for the s field upon neglecting higher order terms. It is a reasonable approximation, based on the smallness in magnitude of the s component when compared to the d component. This approach is useful, for one is interested primarily in the linearized equations which are easier to solve. It will be shown that the anisotropic term added to the first Ginzburg Landau equation by this procedure is sufficient to capture the essential 'lli, i. of the problem. After all the substitutions, the linearized first GinzburgLandau equation takes the form (I2 + d (I ly)2d = 0, (2.8) as which is the central equation and the starting point in my study. For the sake of completeness the second GinzburgLandau equation is = 7 [d*(nd) + d(nd)*] 1 y [(nd)[(n )d] + d[,(n nl)d] a, (2.9) x [(Id)[( nI)d]* + d[In,(nI n )d]* + cc. The boundary condition for the dwave order parameter (2.4) can now be written as n < [d7(Hd) l [y1 (II lf2)d xii,(I2 f12)d = 0, (2.10) .. .. na)d n (2.10) which enforces the vanishing of the normal component of the supercurrent on the surface of the superconductor. This boundary condition would be correct for a superconductorinsulator interface, but for an interface between a superconductor and a normal metal it has to be modified. The modification can be found by imposing the condition n j = 0 on the current equation (2.9). It is convenient to switch to dimensionless quantities using the unit conventions that are frequently used in the literature r Ar, (2.11) d d dod, do 2 Id /\2, (2.12) A (hc/e*)A, (2.13) H (hc/Ae*)H. (2.14) The above equations define the natural lengthscales of the various quantities involved. The GinzburgLandau parameter c is defined c = A// with the penetration depth A and the coherence length C given by A2 (m*c2*47r*2) 2/ d), (2.15) C2 h2/2m* ad. (2.16) The dimensionless Ginzburg Landau equations obtain a much simpler form using these unit conventions. The first equation is simply i2d d f2 )2d 0, (2.17) while the second GinzburgLandau equation, using the vector potential instead of the supercurrent, has the form V x V x A Ad2 +  (dVd* d*Vd) dx1 [ )(12 12)d* 2K 2 I x (2.18) +d[II,(I2 I 2)(] y (IIyd)[(II )d]* + d[IIy(II II)d]* + cc where e is the dimensionless anisotropy parameter, defined as e = ( ?l.'y/as). Setting the anisotropy parameter c equal to zero one can immediately obtain the linearized GinzburgLandau equations for an isotropic superconductor. From now on I switch to the usual symbol T for the order parameter instead of d since I eliminated s. I also need to rotate the vortex lattice by an arbitrary angle 0 about the direction of the applied field. This is required because the free energy is no longer invariant under rotations about the direction of the applied field, due to the presence of the symmetry breaking anisotropic term. The preferred orientation of the vortex lattice will be determined by the minimum of the free energyit will no longer be arbitrary. This crucial point was missed in some of the previous theoretical treatments of this problem [72]. The anisotropic term provides the effective coupling to the crystal lattice, thereby making structural transitions possible. Without this c"..ilii. the vortex repulsion is isotropic for cubic or tetragonal materials resulting always in a hexagonal vortex lattice. Although the anisotropic term was derived from the two component GinzburgLandau theory appropriate for highTc superconductors with d12_y2 pairing it can serve as a phenomenological term for the investigation of anisotropy effects in tetragonal materials like the family of the borocarbides. A useful transformation of the first GinzburgLandau equation is effected by introducing creation and annihilation operators through the following relations II = fI, F illy, with i[II,, II] = h/v. niI(n_) is the creation (annihilation) operator. The transformation of these operators under rotations is shown below. The primed quantities correspond to quantities in the rotated system. I, = I, cos + Iy sin (2.19) II' = fI sin + Ily cos 9, (2.20) I 2 n2 i + 2i (2.21) Substituting in the first GinzburgLandau equation and subsequently dropping the primes I obtain II2p p ie + 2 2i n2 ] i 0, (2.22) an equation which can be reduced to the equation for an one dimensional harmonic oscillator with an appropriate perturbation. The same transformation of coordinates for the second GinzburgLandau equa tion yields a complicated result which after some manipulation can be cast in the following form i3 1 + ( 1 V x V x A = A42 *V ) + o 2 (zx VA) 2K (2.23) +(z x r) r2 b2 where To is the zeroth order in c approximation to the solution of equation (2.22), b = B/Hc2 is the reduced magnetic induction and the function A is defined as A(x, y, q) = cos (4))(x4 + y4 62y2) + 4 sin (4 ,,(x2 y2). (2.24) 2.2 Virial Theorem In this section I will present an alternative method to derive a compact and use ful form of the free energy for the anisotropic GinzburgLandau model. The method is based on the virial theorem for superconductivity discovered recently by Doria and others [73]. In Appendix A I derive in detail the generalization of the isotropic result. The virial theorem for an anisotropic superconductor has the following simple form H B = Fkinetic t 2Ffield t 2Fanisotropic, (2.25) 4x7 or the following equivalent form, with the free energy components expanded H.B In p(x)l12, + B2(X) + 2y2 ([n (x)][II(II2 f n2)T(X)]* 4x 2m* 4x as (2.26) [n UI(x)][n,(n,2 _ )(x)]* + CC ) where V x A(x) = B(x) is the microscopic magnetic field. The magnetic induction is defined B = f dx B(x)/V and the homogeneous applied field is denoted by H. Using this extension of the virial theorem, as I show in detail in Appendix A, one can cast the free energy in the following form F 22 K(I 2 _r I b (2.27) (22 1)A + 1 (2K2 1)A + 1 where b is the reduced magnetic induction, PA is the Abrikosov parameter and F is the correction to the isotropic result, which has the following form F R4 ,(2*(H _2H_)22 R( *(2 )2t2 F Re i f + A(2 1) P2) (2.28) IWPA2 12 2 2 2 It is straightforward to check that if one sets the anisotropy parameter c equal to zero one obtains the familiar isotropic result. The above expressions seem to II.w. 1 that for hard, highK, superconductors the effects of the anisotropy c are diminished. However, it will be shown in the subsequent analysis that even small deviations from isotropy produce observable and significant effects on the structure and the properties of the vortex lattice. This extension of the virial theorem allows one to easily generalize results obtained for the isotropic case. I show below the form of the Gibbs free energy G = F HB/47, the magnetization M = (B H)/4r and the magnetic induction B for the anisotropic GinzburgLandau theory (G H) (2 2 IA (2.29) (2K2_1) A 2 i > { ( )} (2.30) B =H (2 A} (2.31) B (2K2 1)A l(2 2 1) A " The generalization of the second Abrikosov identity [4] is easily obtained in the process of deriving the above relations. It reads (Kc H) I{ 2) 2Re{(4U U (1 )2 )} (1 2T 1)) {(.32 \2K Y68 ( 2 2 (2.32) + 2 f1P 12,(fl f)2,p). All the above results are invaluable for the simplification of the calculations that are going to follow. By itself, this generalization of the virial theorem can be used in numerical investigations of anisotropic superconductors presumably with the same success in reducing the overall computational complexity of the problem as the virial theorem had for conventional superconductors. The original virial theorem has inspired a significant amount of work. See for example [32, 7377], and references therein. CHAPTER 3 MEAN FIELD THEORY In order to investigate the structure of the vortex lattice one can start at the mean field theory level. Although a simplified approach, it provides useful insight about the problem and not surprisingly, many theoretical predictions which are in turn validated by experiment. 3.1 Solution of the First GinzburgLandau Equation In this section I will present the calculation of the Abrikosov parameter ?A and the calculation of the correction F to the free energy (2.27). To that effect I will solve the first GinzburgLandau equation (2.22) perturbatively in the anisotropic term. It has the form (in dimensionless units) 1 V A e(I 2 I)2 (1 I_ 2),p. (3.1) KV ) y Rotating by an angle 9 about the z axis (the direction of the applied field) and by introducing creation and annihilation operators at = ill/ 2b, a = il_/v2, with nII = IIx iFl,, one can write the first linearized GinzburgLandau equation in the following form ata_ a 2 _ e2iOa212: 0, (3.2) without specifying a particular gauge yet. The rotation about the z axis is necessary to account for the orientation of the vortex lattice which is no longer arbitrary due to the four fold symmetric term that breaks rotational symmetry about the direction of the applied magnetic field. Notice that the term 1 was grouped with the higher order term ILT121 in the right hand side of equation (3.1). This way, the linearization of the first GinzburgLandau equation is done in a consistent manner, as pointed out by Kogan [78]. To proceed with the solution it is helpful to decide on a particular gauge. Two choices are useful in this problem, the Landau gauge A = (0, bcx, 0) and the sym metric gauge A = (b y/2, brcx/2, 0). I will show the results of the solution in two separate sections. The total wavefunction will be built by superposition, using appro priate linear combinations to achieve the desired periodicity of the observable (from a quantum mechanical point of view) I~(x) 2. The wavefunction itself is neither periodic, nor gauge invariant. 3.1.1 Landau Gauge Throughout this section I will use the Landau gauge A = (0, bKx, 0), which makes the perturbation theory calculation easy, as it involves only one dimension. In fact, all the gauge sensitive calculations that are presented in this work are done with this particular choice of gauge. I use standard perturbation theory to obtain the wavefunction up to second order in the anisotropy parameter e. The relevant matrix elements are easy to calcu late realizing that the underlying theory is that of the very familiar one dimensional harmonic oscillator. I obtain _e4 2 (cb)2 F e8 I(x) = o(x) + 2H4 bx) 2 64 H8 H(b2x) (3.3) 32 32 64 + 5e4H4(b2x) + O(e3)}, (3.4) where H, is the Hermite polynomial of the nth rank and O(X) ( b exp bK } (3.5) x 2 YCL & ^t XCL Figure 3.1: The geometry of the system: the vortex lattice (VL) is rotated by an angle 0 about the caxis with respect to the crystal lattice (CL). The apex angle is defined as the angle 0 between the two basis vectors a and b. is the unperturbed solution. Notice the explicit appearance of the orientation angle c in the expansion of the order parameter. The periodic solution is constructed following Chang and others [79]. One can choose one basis vector, a, to lie on the y axis and the second basis vector b forming an angle 0 with a. The geometry is illustrated in figure (3.1). The averages (l1I2) and (1 '4), necessary for the calculation of the Abrikosov parameter IA and the correction to the free energy F, are calculated using the methodology presented in the work of Chang and others [79]. Along with the known first order results, I present below the results to second order in the anisotropy parameter (the details of the calculations can be found in Appendix B): 3A = o + ebi + (eb)2 71 72 + 51 3o (3.6) r= 8cb2 {3A( 22 1)[1+ 3cb] +3o[1 [ + 6b] i c 2b(i + 2)}. (3.7) The various quantities that appear in the above expressions are rapidly converging functions of three variables p, a, 9, which are introduced impl,ving the convenient parameterization p + i = (b/a) exp (iO) [1]. The functions 3o, 31, i7, 2 are periodic in p with period 1, and have a rich structure for a ~ 1. Their explicit forms are Ao > Anm, (3.8) nrn 1 R]L C A^4io LYA.Tn 8 2(2 n 4 67rTn2 + 3]I, (3.9) 36n2r 2) a(n2 2) 71 C eL8io A.Tn[167r 4(4 n8 1273a3n6 + 37r 2(24  4_ 45 7an2 + 105] (3.10) 72 'A. 1674 (iT4n4 127 3(3n2M2(n2 + M2) + 3 ( T2n4 + m4 16 256 where Anm = /ae2i7P(n22) 'e27a(n2 2). The prime implies that there are two summationsthe one shown, over n, m and the other with n and m replaced by (n + 1/2) and (m + 1/2) respectively. 3.1.2 Symmetric Gauge This is considered a more "natural" gauge, with no preferential treatment of one direction in space over another. It unfortunately comes with the disadvantage that for this particular problem it is difficult to construct the periodic order parameter, in sharp contrast with the Landau gauge which readily allows the construction of the order parameter with its desired periodicity. However, this gauge choice being phys ically more transparent results in an order parameter whose form has an interesting interpretation. The solution of the unperturbed first GinzburgLandau equation in the sym metric gauge reads po(x, y) =V exp (x 2 +) (3.12) lb}CV 2 bp^(+y (3.12) The correction to first order in c is easily obtained using the method of the previous section. I obtain ) (x,y) = o(x, y) 1+ e4 4 (3.13) where z = x iy. The second order correction is easy to calculate but it does not add anything to the 1ir, i. .,1 point I am trying to make. Although the periodic order parameter is difficult to construct, one can postulate an amendment of the isotropic equation (1.14) to account for the anisotropy. The form of the perturbative result leads one to propose the following form for the solution which accounts for fourfold anisotropy I(x, y) = exp [ b (x2 2) (Z z ) {1 C (z 7 Z 1)( 7 1) J (3.14) X (z + 72)(+ v 72) , with i7 = 7R+i7I = I and 72 tI+i R = I I ' and v denoting the position of the vorticesin other words z, are the zeroes of the unperturbed order parameter. The 1.li,i .,1 interpretation becomes apparent once one is familiarized with the strange form of the solution. The c term provides the four additional zeros of the in duced swave order parameter, around the central node of the d wave order parameter. The positions of the nodes are controlled by the 7 terms and they can accommodate various symmetries. I have attempted to obtain the compression (cil) and tilt (C44) elastic moduli using the solution (3.14) according to Brandt's prescription [3538] without success. The complexity unfortunately increases rapidly with each successive step. However, there are some interesting results that were obtained that are worth mentioning here. I first need to define two functions A,(x, y) and M (x, y) through A,(x, y) = X4 + Y4 6X2Y2 14 cos(4w), (3.15) M,(x, y) = 4XY(X2_ Y2)_ 4 sin(4aw), (3.16) where X = x and Y = y y, with v denoting the vth vortex. The modulus of the order parameter can be reduced to the following form f2 ftropic( + 2g), (3.17) where g is the correction due to the anisotropy and it is defined as g = [F A,. It is easy to show that the supervelocity and the quantity u = Vf/ f have the form Q = Qisotropic + Vg x z, (3.18) U = Uisotropic + Vg. (3.19) The equation Q = u x z still holds. Other useful identities that have a straightfor ward generalization in the anisotropic case are those first discovered by Brandt [35, 37] when deriving the elastic moduli of the isotropic superconductors. I obtain for the supervelocity of the Abrikosov solution (momentarily I am following Brandt's nota tion) 1 rrr QB = A z 2x i { 2 + 4e [(X3 3XY2)x + (y 3X2Y)y)] (3.20) The following central identity remains unchanged: (Vw)2 12Bw + wV2w. (3.21) Note that ci is the linearized order parameter, obtained from the Abrikosov solution for the undisplaced vortex lattice cLo, upon substitution of the displaced coordinates r, and subsequent expansion in the displacement s, to linear order. The general vortex lattice is defined by the loci of the zeros of the order parameter given by: r,(z) R, + s,(z). (3.22) Another very int rI lin_. from a ]llr,i. .,1 point of view, result is that an ar bitrary displacement s, of a vortex line and the anisotropy appear to modify the linearized order parameter in a similar manner. I obtain 1 2 wi(x,y) o(x,y) 1 + 5 (3.23) with wjo(x, y) the solution for the undisplaced vortex lattice. The generalized displace ment term q, which in the isotropic theory accounts for the displacement only of the vortex lattice, now has the form (r) 2 2 (r r ) +VW (3.24) r r,2 W = X6 + Y6 5X2Y2(X2 + y2) i4[(X2 + y2) cos(4w)) + XY sin(4w:)]. (3.25) The complexity is significant but the calculations are tractable so far. What hinders further progress is my inability to find the appropriate generalization of the free energy which serves as the starting point of the derivation of the dispersive elastic moduli cl1 and C44. In the isotropic theory it has the form [35, 37] f = + + ( )2 Q hBh) B2 (3.26) 2 4K2L, 3.2 Structure of the Vortex Lattice and Minimization of the Free Energy In order to obtain the mean field phase diagram and study the structural phase transitions formally, one has to minimize the Gibbs free energy (2.29), which is the proper free energy to minimize under the constraint of a constant external field. I found that the minimization of the Gibbs free energy 1 G = F HB/47, and the minimization of the Helmholtz free energy, equation (2.27), give identical phase diagrams. I choose to work with the Helmholtz free energy, from which one obtains the elastic moduli, as shown in Chapter 4. The particular form of the free energy (2.27) allows one to obtain the phase diagram in a direct way instead of minimizing only the Abrikosov parameter 3A as was done by other workers [72, 7981]. However, for the sake of comparison and as a consistency check to my algorithm, I minimized the Abrikosov parameter PA in order to obtain the phase diagram as a function of the anisotropy. The results are shown in figure (3.2). To first order I find the transition point at ecl = 0.0240, in agreement with Chang et al. [79] who find qci = 0.0235. The result of the numerical solution for the phase boundary of Park and Huse [72] cannot be compared with the previous results By solving for the phase boundary numerically they find qci = 0.0367.2 However, in their work the vortex lattice was not allowed to rotate, thus omitting 9 from the calculation, an essential degree of freedom. The second order correction to the Abrikosov parameter 3A moves the critical anisotropy to c2 = 0.0284 a relatively small correction which justifies the use of perturbation theory. 1After one modifies G to be a function of the reduced field b. 2The reported result, cl = 0.0734, needs to be divided by 2 due to different conventions for the anisotropy parameter. 0 75 65 F Apex angle 0 vs. anisotropy E  first order second order Orientation angle vs. anisotropy e 45 r I  first order SIsecond order '45" r  40  1 35 IEr 30 S50 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 / e 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Figure 3.2: Minimization of the Abrikosov parameter IA as a consistency check for the computer code. The apex angle 0 of the unit cell and the orientation angle 9 of the vortex lattice as a function of the anisotropy parameter c. The critical points to first and second order approximation in the anisotropy parameter c are Ec1 = 0.0235 and c1i = 0.' 1' Notice the reorientation transition in the inset panel which appears only to second order in the anisotropy parameter at c, = 0.0492. The most interesting and unanticipated fact was the observation of another phase transition at a still higher value of the anisotropy c, = 0.0492, which turns out to be connected with the vanishing of the rotational modulus as will be shown in Chapter 4. In particular, at this second transition point, the orientation of the square lattice changes continuously from the [110] direction, with 9 = r/4 now a local maximum and two degenerate minima on higher and lower angles. The vortex lattice retains its square symmetry. The minimization of the free energy shares the complications of the minimization process for the Abrikosov parameter 3A. There exist many local minima that have to be avoided in the search of the global minimum and also there are many equivalent minima which correspond to the same lattice structure, with different values for the three parameters p, a, o which determine the structure of the vortex lattice. This problem can be solved realizing that for the two subsets p E [0, 1/2], p E [1/2, 1], one obtains equivalent lattices. I chose to work in the first subset. The minimization was carried out using the fmincon function for constrained minimization included in the optimization package of Matlab. A program was written for the evaluation of the free energy, the calculation of the gradient and the Hessian. A driver was then used to automate the minimization for different values of the parameters of the problem. The interested reader will find the code in Appendix D. The results of the minimization of the free energy for e 0.11, K = 5 are shown in figure (3.3). The difference of my results compared to previous theoretical predictions [79, 80] obtained from the minimization of the Abrikosov parameter 3A is that the structure of the vortex lattice is now obtained as a function of the reduced field b, allowing the comparison with experiment. All the main features of the second order structural phase transition are there, along with the reorientation transition in the region close to the upper critical field. In order to check the results and the minimization algo rithm I calculated the correction F to second order in the anisotropy parameter. The difference between the first and second order in the anisotropy parameter results is not significant, as it can be seen in figure (3.3). For larger values of the anisotropy parameter the transition occurs at a smaller critical field, while for smaller anisotropy the transitions occurs closer to b = 1. I chose the particular value to show the reorientation transition clearly. The predictions of the anisotropic GinzburgLandau model are at variance with the predictions of the nonlocal London model, with regard to the reorientation of the vortex lattice at low fields. The nonlocal London model predicts [46, 82] that the reorientation happens at a low field H1, with the orientation angle being 4/4 Apex angle 0  I first order Ssecond order 45  40 35 30 25 0.2 Orientation angle ) 0.4 0.6 0.8 b 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 b Figure 3.3: The apex angle 0 of the unit cell, the orientation angle 9, the Abrikosov parameter IA and the dimensionless Helmholtz free energy F/i2 as a function of the reduced field b for e 0.11 and K = 5. In the first panel on the top I include the second order result for the apex angle. The reorientation transition is evident in the second panel on the top. throughout the rhombic to square transition. This has been recently verified by experiment on annealed LuNi2B2C crystals [15]. In my model the reorientation of the vortex lattice is continuous and is completed at the transition point. This is not strange, as any expectation of a quantitative agreement between the two theories at low fields (the natural region of applicability of the London theory) would be unreasonable. ' ' CHAPTER 4 ELASTICITY THEORY FOR THE VORTEX LATTICE In this section I will derive the elastic moduli of the vortex lattice. This deriva tion is independent of the particular model (GinzburgLandau theory or London model) that will be used to evaluate the elastic moduli. As is done for the stan dard theory of elasticity for crystals, one can introduce the displacement vector u of a vortex line from its equilibrium position and upon expanding the elastic energy in terms of derivatives of the displacement ui,k akUi, one obtains the various elastic moduli. In principle the derivation is the same, however one has to consider carefully the fundamental differences between the usual ( Il.1l and the vortex lattice. Before proceeding with the derivation a few comments are in order. The full calculation of the elastic response of the vortex lattice is a complicated task. I will consider only the simplest part of that, namely straight and parallel vortices, thereby reducing the complexity of the problem from three dimensional to the simpler two dimensional case. A further simplification is the so called incompressible limit, in which one con siders deformations that conserve the density of the vortices. This puts the following constraint on the displacement field: V u ui, = 0. For the shear moduli calculated in this fashion it has been proven that they are essentially nondispersive [25]. One then has to consider only uniform deformations. The vortices are inhabiting the superconductor and are thus influenced by the symmetry of the underlying crystal lattice. Moreover, rigid rotations of the vortex lattice with respect to the (I, I .,l lattice have an energy cost associated with them, a strange property never encountered in ordinary solids which live in an isotropic space. The immediate consequence is that the usual decomposition of the elastic energy in terms of strains and rotations is no longer convenient. A more useful form for the elastic energy has been proposed by Miranovi6 and Kogan [83]. I will adopt that form and present here, in detail, the derivation of the elastic moduli for the vortex lattice. The form of the elastic energy density in terms of a general fourth rank strain tensor can be written 1 E iklmUi,kUl,m (4.1) 2 The tensor Tiklm is not symmetric under the exchange of i + k or 1 + m, which would permit a decomposition into a symmetric and an antisymmetric part. In fact the only symmetries that Tikim possesses are those of the ( 1I..1 lattice and of the equilibrium vortex lattice along with the obvious Yiklm lmik I will deal with the more symmetric case first, the square vortex lattice. The energy should be invariant under an exchange of x and y coordinates. Thus, x and y appear in even numbers in the elastic energy. Under the constraint of incompressibility the elastic energy reads 1 Es [21s,2 + t 72s(Uy + u2t ) + 2m3sUx,yUy1], (4.2) with s = 2(,, ,, ), 72s = , and , Similar considerations for the hexagonal vortex lattice lead to an elastic energy of the form 1 2 2 U2 Eh = [/lh U, + 72h ,y + 73h y, 274hUx,yUy,], (4.3) with lh = (Y + r 2.,, ), 72h = 73h = r and 74h r, In order to make the elastic moduli more transparent, it helps to consider various deformations and to study the form of the elastic energy for each choice. There are four important deformations: Squash deformation: u, = p(xx yy), * Simple shear along x: u  pyx, Simple shear along y: uy = pxy, Rigid rotation about z: u, = p(xy yx). Under each of these deformations, the elastic energy assumes a particularly simple form which allows one to chose a convenient parameterization with elastic moduli Csq, C66x, C66y and cr, that are readily connected with a simple 1ir, i, .1l situation. More specifically I have for the square lattice first Es [uq] = 752, E [u,] = "72s2, E~[uy] = 172s2 E [ur] = 2(y72s 7s)2. I can then identify Csq = 'y7, C661 = C66y 72s and c, = 2(72s 73s) as the relevant elastic moduli for the square vortex lattice. In a similar fashion I obtain for the hexagonal lattice Eh[Usq] 1= /al[2 Eh [ux] = 72h 2 Eh[Uy] = T73h 2 Eh[Ur,] 1(2h + 73h 274 )2. Thus, the relevant elastic moduli for the hexagonal lattice are Csq = l1h, C66x 72h, C66y 72h and cr = (72h + 73h 24h). In figure (4.5) I present the action of the four particular displacement fields on the basis vectors of an arbitrary lattice to clarify the 1.li, i, .l1 picture. Squash deformation Shear along x / ///l/a a  /////( *''' ///////r  // / // 0 '0 II I I I l i  \^ \\ , \1\\\\s _o \\\\ss.>. .' tttt\\\\ ... 1s \\\\\ 1111rr / Rigid rotation Shear along y Figure 4.5: The action of the four basic displacement fields, represented by the small arrows, on the basis vectors of an arbitrary lattice, represented by the large arrows. ////////I I f tPII \\\\\ //'/// / t I t t 11 \ \ I.''' I I # t f t I \ N N \ N \\ m f / I I t I I I k N N N N N S 5 ^ r fc4 O.N 1 /I I I % I I I N/ /NN 1\\\\ 1 111 % / ~~~~~~~\\\\\\\\//1 J I If 0 V le e e e/J L~~~~\WIII/A/ 105         . . .. . . . 1 85 4. I  . . .  .7   I .1 t t t o* t tt t I'n 1 t s t t t I, 1 1 t ,t tt J, 1 i t I t t 4 I U t t t t ', ', t t t , ~ t t JSia ... u*lttt 1 t t I r I t ttf 11 3tt. .t t 4, 1 1 t1 I t t tt There is another deformation, although not independent from the above. It is the shear along the [110] direction for which one can easily prove that C66() = Csq + Cr. (4.4) The ratio of c66; and C66() is a useful quantity that measures the anisotropic response of the vortex lattice to shear deformations of different polarizations and will be of use later on. For a general rhombic lattice there are four independent elastic moduli Csq, c66X, c66y and cr, first proposed by Miranovib and Kogan [83]. The deformation responsi ble for the transformation of the general rhombic lattice into a square lattice is the socalled squash. The relevant modulus vanishes at the point where the structural transition takes place, signaling the instability with respect to the particular deforma tion. It will be shown that another instability occurs at higher fields. At that point the vortex lattice becomes unstable with respect to rigid rotations about the c axis. All elastic moduli will be obtained from the second derivative of the Ginzburg Landau free energy with respect to the amplitude pt of the deformation, in the limit that ft  0. I will make use of the GinzburgLandau free energy in the following form F b2 (2b2 i { F(2 i (4.5) _1)A 2 A 1 I The elastic energy can be written as E =1 2 p2C(k ) 0 202 2aF (4.6) where c(p) is the elastic modulus for each particular deformation. The amplitude ft of the deformation will enter the free energy F through the three parameters p, a and 0, which in turn depend on the lattice structure. The connection of a particular model of superconductivity and the elastic moduli is made at this stage. I chose the GinzburgLandau theory, however the evaluation of the elastic moduli is along the same lines for the nonlocal London model. It is straightforward to show that the second derivative of the GinzburgLandau free energy with respect to the amplitude of the deformation has the form (22 1)(b2F) 2 ( b)2 r S (22 1)A 1 [(22 1)A + 112 F (4.7) (2K2 1) 1 2 P 2(2K2 1)2 F(0/ A)2 (2K2 1) 1 A + [(2K2 1)/A + 1]2 Both PA and F depend on the three parameters p, a and 0. The derivatives of /A and F are thus given by the chain rule. One can define ( = (22 1)/[(22 1) A+1] and S= (1 b)2/[(22 1)A + 1]2 and obtain the following expression for each elastic modulus (denoted by c for brevity) using equations (4.6) and (4.7): c pp + + + 2Cp + 2C+2p ,0 + 2C0,7 09 (4.8) where Cij are the terms originating from the derivatives of the GinzburgLandau free energy and are given by Ci, (2 6ij) (b2 F)/A,j + ,ij (I)(A,,J +/3A,j,i +/3A, 2(/A,lA,j) }. (4.9) Notice that the differentiations in the right hand side of equation (4.9) are with respect to the three lattice parameters p, a, 9. The calculation of pp, a/, 0, is purely geometrical in nature. One needs to express the parameters of the distorted vortex lattice as a function of the magnitude of the deformation and evaluate the derivative at the limit of vanishing deformation (t  0). The geometry of the system can be seen in figure (3.1). The vortex lattice is rotated by an angle 9 with respect to the ( vI.1l lattice. The direction of the rotation is completely arbitrary, I chose the particular one for convenience. The angle 0 is the apex angle of the (generally rhombic) unit cell and a and b are the two basis vectors. I will start with the geometrical definitions of p, a and 0 which have now become functions of the amplitude of the deformation p(1 ( 2 (4.10) 2 ,iit) (4.11) cos a(f)y }, (4.12) ( )I a^ ft) I )1 I Under the influence of each deformation the two basis vectors are going to change direction and/or amplitude resulting in a deformed unit cell. The calculation of these effects on the unit cell as coded in the parameters of the vortex lattice for each individual deformation is straightforward, however one should choose carefully the branch of the inverse trigonometric functions that are involved in order to avoid sign errors. The results are Squash: (p,, (a,, ,) (2a sin(20), 2a cos(20), sin(20)), (4.13) Shear along x: (p1, (,, ,) = (a cos(20), a sin(20), cos2()), (4.14) Shear along y: (p,, o,, 0/) = (a cos(20), a sin(20), sin2 ()), (4.15) Rigid rotation: (p,, (a, 0) = (0, 0,1). (4.16) A simple substitution of the lattice parameters in equation (4.8) yields for the nondis persive elastic moduli of the vortex lattice Csq = 4Cppa2 sin2(2) + 4C,,(2 cos2 (29) + Co sin2 (29) 4CoU2 sin (49) (4.17) 4CpO, sin2(20) + 2Ca sin(40), C66, Cpp(2 Cos2(20) + C(U2 sin2 (20) + C cos4 () + C,(T2 sin (4) (4.18) 2Cpa cos2(C) cos(29) 2CO COS2 (S) sin(29), c66y = Cp2 cos2 (2)) + Cga2 sin2 (29) + Co sin4 (0) + Ca2 sin(49) (4.19) + 2Cpa sin2(C) cos(20) + 2C0,O sin2(0) sin(2 ), C, = Coo. (4.20) Up to now the elastic moduli have been independent of the particular model that is used for the study of the vortex lattice. The coefficients C, which are model dependent will be evaluated in the next chapter. CHAPTER 5 STRUCTURAL PHASE TRANSITIONS Using the results of mean field theory from Chapter 3 and the elastic moduli calculated in Chapter 4 I am in a position to calculate the elastic properties and study the phase transitions of the vortex lattice within the anisotropic GinzburgLandau model. 5.1 Nondispersive Elastic Moduli The procedure is very simple. First one has to minimize the free energy, equa tion (2.27), with respect to the three parameters p, a and 0. This will produce the equilibrium energy and the equilibrium lattice for a given value of the reduced mag netic induction b. The details of the minimization are presented in section 3.2. The second step is the calculation of the appropriate derivatives of the free energy at the equilibrium values of p, a and ( and the calculation of the four elastic moduli, given by equation (4.9), again for a given value of the reduced magnetic induction b. One has to choose values for the two parameters of the model, the anisotropy parameter e and K. The elastic moduli were obtained as a function of the reduced field for e 0.11 and c = 5. They are evaluated at each minimum of the free energy when obtaining the phase diagram in section 3. The calculated elastic moduli for c = 0.11 are shown in figure (5.1). The vanishing of the squash modulus signifies the transition from rhombic to square vortex lattice. At the same point the two shear moduli merge into one because of the higher symmetry of the square phase. At a still higher field I observed the reorientation transition signified by the vanishing of the rotation modulus c,. This instability of the vortex lattice relative to rotations was encountered in the nonlocal London model also [83]. There are experimental indications for this instability, the investigation of which is hampered so far by the high fields required for its observation [15]. 0.01 I I I I sq 66,x 0.008 C 6 66,y S ... C 2e05  I \ 0.006 le05 0.004 \ S0.85 0.9 0.95 1 b 0.002 .. 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 b Figure 5.1: The nondispersive elastic moduli as a function of the reduced field for c = 0.11 and K = 5. The first transition point at which the "ijl.ih" modulus Csq goes to zero is clearly seen in the left panel. The reorientation transition starting at the point where the rotational modulus c, vanishes is seen on the right panel. I can compare directly with the nonlocal London model results of Miranovi6 and Kogan [83]. I can see the same qualitative behavior, the vanishing of the squash modulus Csq at the first transition point b:, the merging of the two shear moduli into one after the first structural transition and also the vanishing of the rotational modulus c, at the second transition point br, which is about two times b:. However, a quantitative difference between the two results is the change in the relative magnitude of the elastic constants after the phase transition. I observe that the elastic moduli 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 b Figure 5.2: The ratio cs6,/c,.,. versus the reduced field b. in my work constantly diminish towards zero as one approaches H, (T), as expected. This behavior is attributed to the (1 b)2 factor in the free energy (2.27). I also note that the anisotropic elasticity which manifests itself with a softer shear modulus for shearing along the sides of the square lattice than along the diago nals emerges naturally in this model. This response of the vortex lattice is measured by the ratio c66r/c,.,.:, which has a maximum at the transition as can be seen in figure (5.2). 1 It has been i.'. led that this behavior explains the anisotropic orienta tional long range order observed in decoration patterns in LuNi2B2C, which manifests itself as a significant difference in the correlation lengths along the [110] and [110] directions [15]. SNote that around b = 0.95 the theory breaks down because the F term can no longer be considered a small correction. This is the origin of the discontinuity in the second peak in figure (5.2). 5.2 Fluctuations The thermal fluctuations of the vortex lattice are incorporated in the anisotropic GinzburgLandau model through the three variables p, a, 9, characterizing the struc ture of the vortex lattice. More specifically, I will replace p, a, 9, in the expressions for the free energy and the squash modulus by their thermally averaged values S= p[l (2 p)u2], (5.1) I = a[1 a(2 p)U], (5.2) S= a cot u2. (5.3) The derivation of the above equations goes along the lines of the derivation of the p,, ac and presented in Chapter 4. An arbitrary displacement field u corresponding to the thermal fluctuations is introduced with the following constraints (u(xi)U(xi)) = 0, for i / j, (5.4) (u4(Xi)Uj(X2)) = 0, (5.5) 12 (ui(x)u(x)) (u2(x)). (5.6) 2 Starting from the definitions of p, a and 9 (equations 4.10, 4.11, 4.12), under the above constraints it is straightforward to obtain after some simple algebra the thermally averaged p, a and q. One then has to find where the zeroes of the squash modulus lie in the HT plane. This procedure is analogous to the one carried out by Gurevich and Kogan [84] within the nonlocal London model, taking the simplifying approach of not calculating the phase transition line Ho(T) selfconsistently. That means that one does not take into account the effects of the fluctuations on the elastic moduli themselves. One is allowed to make this simplification based on the fact that the mean squared displace ment of the vortex u2 is finite at Ha(T) [84]. In order to show that explicitly, I will present the derivation of the mean squared displacement. One starts from the elastic energy of the square vortex lattice: 2E = CqU + c + c + c44( 2, (5.7) written in the, convenient for this case, form in terms of strains ui and rotations about the z axis Lcwy, Vuj 2 1(i + = 2 (q2uj qeqr+i + g jiij), (5.8) 1 i d2qdkqr Lxmy = (0uy 0yuX) = 2 dk eiq.r+ikz(qxuiy qyix), (5.9) 2 2f(27)3 where the u is the Fourier transformed displacement. For a more general treatment of the problem, the tilt modulus c44, related to bending deformations of the vortices, is included in the expression for the elastic energy. The condition of incompressibility takes the form Oxuy + Oyu 0 => cos f u, = sin ity, (5.10) q = q(cos 9, sin ), (5.11) from which one can obtain 2E (2) + c44k2 1+ Cs2 j ix (q)i x,(q), (5.12) (27)3 4 sin2 sin 2 with c(o) = Csq sin2(29)+cx Cos2(2)+c,. In conclusion, the thermal average u2 will be calculated from the general expression which relates the former with the components x,, and yy, of the elastic tensor (see for example Blatter et al. [22]) 7 7 7 Tf, d2qdk ( 2 8U + u =qT (2 ) ; .D(5.13) (2i7)3 D) In the case under investigation, the diagonal components of the elastic tensor are equal to the coefficient of ui(q)ui(q) in equation (5.12). I obtain / d2 qdk sin2 (5.1 )u2 = kB T 2(5.14) 1(27?)3 q2c() After performing the integral, the mean squared displacement has the form [84] U2 K Cx C,q (5.15) 66 + Cr where K is the complete elliptic integral of the first kind. It is evident that even at the instability point where Csq = 0, the mean squared displacement u2(T, B) remains finite. There are two competing parameters, the strength of the thermal fluctuations X and the anisotropy parameter c which in this model is the equivalent of the nonlo cality parameter pn of the nonlocal London model. The competition arises because thermal fluctuations of increased strength tend to smear out the fourfold symmetric character of the intervortex interaction, making the interactions isotropic and eventu ally restoring the rhombic/hexagonal arrangement as the lowest energy configuration of the vortex lattice. Following Gurevich and Kogan [84], I will take A(T) = Ao/ v t, C(T) Co/l/ t2, where t = T/TC. I then have X = Xot/ 2.2 The dimensionless mean 2Other forms of A and E with different temperature dependencies give similar results. squared displacement reads Xorl bt uq2 = (5.16) (1 t2) [(1 b)3ln(1+ 2)] ' with r ~ 3 for LuNi2B2C [84]. The anisotropy parameter c was connected to the nonlocality parameter pn via correspondence from the nonlocal London model which is derived from the extended GinzburgLandau theory [71]. One can obtain = 21 t2). (5.17) The above expression is different from what one would obtain by simply substitut ing the temperature dependence in c = ? /. I ., 'yas. Regardless of the fact that the results would be similar even with a different choice of temperature dependence, the motivation behind this particular choice was to enable the comparison with the re sults of the nonlocal London model, facilitated by the introduction of the nonlocality parameter pl. The results of the numerical solution are shown in figure (5.3). The mean field transition line without fluctuations was included in the graph also. Instead of choosing arbitrary values for all the parameters of the theory I used the same parameters used in the recent work by Gurevich and Kogan [84] to facilitate the comparison of the two theories. The transition lines are obtained for pnl = 2.5C0 and Xo = 0.0064. One can see that the results of the two theories are in reasonable agreement thus validating the elegant 1.r, i, .,1 picture concerning the competition between the thermal fluctuations and the nonlocality. The transition line for the fluctuating vortex lattice ends abruptly at some point. At that point the correction term F, equation (2.28), ceases to be small enough to Transition lines for p=2.50 xo=0.0064 0.2 0.4 0.6 0.8 Figure 5.3: The phase boundaries in the HT plane separating the square and rhombic vortex lattice obtained from GinzburgLandau theory and the nonlocal London model. The line intercepting H, (t) is the mean field GL result. be considered a correction and at this point the extended GinzburgLandau theory becomes unstable. This shortcoming of the model could be corrected by the inclusion of higher order terms or even by the inclusion of an isotropic term of the same order in the anisotropy parameter c, but that is not going to affect the qualitative behavior of the transition lines. Another important observation is that for a GinzburgLandau parameter t with value in the range appropriate for typical members of the cuprate family ( c 100) the reorientation transition does no longer appear. This seems to ii,.'. 1 that the reorientation transition is irrelevant for the cuprates. As a closing comment I would like to address a recent objection raised by Nakai and others [85] regarding the ]llii .,1 mechanism behind the reentrant structural transition. The authors perform a calculation based on the Eilenberger theory and claim to have found the true pli, i. .il origin of the phenomenon in another kind of competition. They claim that the reentrant behavior is due to intrinsic competition between superconducting gap and Fermi surface anisotropy, which provide two differ ent possible orientations of the vortex lattice, one with nearest neighbors along the gap minimum and another with nearest neighbors along the Fermi velocity minimum. Their model, as they admit in the paper, fails to predict the correct orientation of the vortex lattice in LuNi2B2C, at least for field strengths used in experiments up to now. Moreover, my most serious objection is that their work is based on assumptions about the structure of the superconducting gap in the borocarbides that remain to be experimentally verified and which are not met in the cuprates. In my opinion, the work of Gurevich and Kogan [84] along with the present work have unambiguously shown the importance of thermal fluctuations in the determination of the phase di agram, an issue that would still have to be addressed in the alternative scenario in order to achieve a complete description. CHAPTER 6 FINGERING INSTABILITY OF ELECTRON DROPLETS From here on I devote my attention to the fingering instability of electron droplets in nonuniform magnetic fields. This work was inspired by a recent discovery by Agam, Wiegmann and others [86, 87], which connects diffusion limited growth in the limit of long diffusion lengths (also called Laplacian growth) with the growth of a twodimensional electron droplet in a high, nonuniform magnetic field. The calcu lation of the harmonic moments and the development of the Monte Carlo code were my contributions to a collaboration with Alan Dorsey and Taylor Hughes [88]. Laplacian growth obtains its name from the description of a front propagating between two incompressible liquids with high viscosity contrast (e.g. water and oil). One considers a 2D flow, where the less viscous fluid (water) is supplied by a source at z = x + iy = 0 and the more viscous fluid (oil) is extracted at the same rate at infinity, resulting in the motion of the boundary between the two phases. Assuming zero surface tension, the idealized surface moves with a normal velocity vI which is proportional to the gradient of the pressure p on the surface (Darcy's law): n (6.1) On For incompressible liquids the normal velocity becomes proportional to the gradient of a harmonic field p given by: V2p = 0. (6.2) For a nice review the reader can refer to [89]. Fingered patterns form at the interface of two incompressible fluids with large viscosity contrast only when the less viscous fluid is injected into the more viscous fluid. It is essential that the system has a 2D geometry, an experimental realization of which is the HeleShaw cell. It has been found that an arbitrary deviation from a perfectly circular shape of the initial droplet will result in a pattern of growing fingers which grow as more fluid is pumped into the droplet and these fingers eventually form cusplike singularities at their tips. The same behavior has been postulated for the Quantum Hall droplet, the source of the instability in this case being the changing gradients of the applied magnetic field, producing an inhomogeneous magnetic field in the exterior of the droplet. The influence of field inhomogeneities on the shape of the electron droplet is a manifesta tion of the AharonovBohm effect. In this case also, singularities will develop in finite time. In the fluid system the singularities will be controlled by the molecular scale of the particles and also by the surface tension which is absent in the idealized Laplacian growth. In the quantum counterpart of this phenomenon the singularities will be cutoff by the only length scale in the problem, the magnetic length 1B = hc/eBo, which is the effective size of the particles involved. The fingering instability has been studied for many years since its discovery [90 94]. One can use conformal mapping methods to map the exterior of the unit circle (w plane) onto the domain external to the droplet (zplane). One fundamental property for the later connection of the two phenomena is the conservation of the external harmonic moments defined for the fluid and electron droplet respectively as tk 1 zk d2z, (6.3) 7t k J exterior tk 1 B(z)zk d2z, (6.4) 7k "7Jeterior with z x + iy, k = 1,2,... and 6B(z) being the magnetic field inhomogeneity. The moment for k = 0 is proportional to the area of the droplet and the integrals for k = 1, 2 are assumed to be properly regularized. These quantities are of central importance for the theoretical treatment of the problem since they are used to relate the magnetic field inhomogeneity to the parameters of the conformal map. This part is divided in three sections. In the first two sections I show the details of the derivation of the m.iii. body wavefunction of the 2D electrons in a perpendicular, nonuniform magnetic field and the conformal mapping method. In the last section I discuss the Monte Carlo simulation and present the comparison between theory and experiment. 6.1 Manybody Wavefunction For a 2D system of electrons in a nonuniform magnetic field perpendicular to the system, the Nparticle Hamiltonian is N H = ( V + A)2 B(rj)Uj, (6.5) j= 1 where the Coulomb interactions of the electrons have been ignored and a( is the Pauli matrix. I have also used appropriate units, namely the magnetic length 1B = hc/eBo for lengths, the average field Bo for magnetic fields, hca/2 for energies with wcL eBo/mc. If in addition one takes g = 2, equation (6.5) is the Pauli Hamiltonian and the exact solution for the manybody groundstate is known [9599]. It is interesting to point out at this moment that the single particle wavefunctions, the well known Landau levels, are derived in an almost identical method to the one presented in section 2.1. For filling factor v = 1, a filled Landau level, the ground state wave function can be written as a Slater determinant which reduces to [86, 9799] (i..., ZN) = (z z e W(z (6.6) i In the equation above, W(z) = 1z2/4 + V(z), where the first term comes from the homogeneous component of the magnetic field and V(z) is a solution of V2V 6B [B(z) Bo]/Bo. The normalization factor 1/IN!TN will not concern us in what follows. For clarity I will show explicitly the derivation. The homogeneous part of W(z) is given by the solution to V2 = Bo which becomes (if one uses Bo as the magnetic field scale) V2 = 1, which has the solution O(x, y) = Izl2/4. To calculate the solution for the presence of a small inhomogeneity, so that B = (Bo + 6B)z, one has to solve V2 = 1 6B, from which by setting = z 24 + V(z) one obtains V2V = 6B. The formal solution is given by V(z, z = f In z z' 6B(z') d2. (6.7) 27 D It is essential to assume that the field inhomogeneity is far away from the origin where the electron droplet will be located. Then, upon expanding the logarithm for z < z'l In z' z = In z' + ln(1 z/z') In z' (6.8) k= I obtain the expansion of the potential V(z) in terms of the harmonic moments of the domain D, under the constraint that the inhomogeneity vanishes in the interior of the domain D [100]: V(z, z) 2Z tk k (6.9) k=1 tk If 6Bzk d2z. (6.10) 7k D The probability density obtained from the wavefunction (6.6) is N N P2 exp 2 In Izi zj i zi2 + IV ) (6.11) i which can be interpreted as the probability density for a classical 2D plasma in a background potential of the form N N UB = z, I2 ReV(z), (6.12) j=1 j=1 with a constant Boltzmann factor / = 2. This is where the road for Monte Carlo simulations opens. The only remaining issue is the choice of V(z). One possibility is the use of thin solenoids in the form N, 6B = 6(2) (z Zn), (6.13) n=l where Ns is the number of solenoids and DT is the number of flux quanta carried by the solenoid. The potential V(z) after an integration is found to be N, V(z) = In Iz Z. (6.14) n=l This form has both advantages and disadvantages with regard to the Monte Carlo sim ulation. A discussion of these will follow in the section that deals with the simulation details. Another, by far more elegant, way to find a suitable V(z) to carry out the sim ulation is to start from a suitable chosen conformal map that describes the evolution of the boundary. First, a discussion of the conformal mapping within the context of this problem is in order. 6.2 Conformal Mapping Method The idea is to use the powerful conformal mapping methods for 2D systems to transform the highly irregular domain of the exterior of the droplet in the complex plane z onto the exterior of the unit circle wl 1 in an auxiliary complex plane w. The map is reversible by construction. I will present various conformal maps that correspond to different fingering patterns. For a given conformal map z(w,t) the evolution of the boundary in the Laplacian growth problem is governed by the following equation [91]: Re wOz (6.15) ( Ow 0t 2 For a conformal map of the form M z(w, t) = an(t)w", (6.16) n=l one obtains the following differential equations for the (real) coefficients a,(t) [91]: d, 1, (6.17) dt da n=1 Mk / The solution of these will give the time evolution of the boundary. As an example I will consider three simple conformal maps to illustrate the idea. For the conformal map z = a + bw1 one easily finds the solution z = e0 al (0) + t + ao(0), (6.19) which corresponds to a circular boundary at all times. Thus the circular boundary is stable. Considering the next order, namely M = 1, one can show that another stable boundary has the form of an ellipse. The solution is z = [ao(0) + (B + 1)A(t)] cos + i(B 1)A(t) sin 0, (6.20) with A(t) (1 B2) (6.21) B (0) (6.22) al(0) which corresponds to an ellipse with a constant eccentricity e = (B 1)/(B + 1). The mapping of order M = 2 produces the first unstable boundary, one with three cusps. The solution for initial values a2(0) = 1/40, al(0) = ao(0) = 0 and a_ = 1 is z [a2 cos(20) + a cos 0 + i[a2 sin(20) a_ sin 0], (6.23) with 2 a2 1 (6.24) 40 a_1 400 159201 800t. (6.25) The normal velocity of the boundary is given by v, = Idw/dz, which diverges at the moment the cusp is formed. This happens for any polynomial map such as the one I consider here. It can be shown that the number of cusps is equal to M+ 1 which is the number of roots of dz/dw. The divergence is caused by precisely these roots as they approach from the interior of the unit circle and eventually fall onto the boundary w = 1 at the critical time critical. The critical time for the M = 2 map and the above choice of parameters is tcritical = 199. In figure (6.1) I show the evolution of the boundary for this particular map. For higher orders in M the solution is possible only by numerical solution of the algebraic equations derived from the initial differential equations. 20 10 0 O 10 20 10 0 10 20 30 Figure 6.1: The evolution of the boundary for the conformal map with M = 2 from time t = 0 to time t = teritical = 199. An alternative and more elegant way to proceed is with use of the harmonic moments defined in the previous section, which facilitates the connection with the potential V(z). The calculation of the harmonic moments for the various confor mal maps that I used along with the connection to the potential V(z) is presented separately in Appendix C. 6.3 Monte Carlo Simulation The Monte Carlo method is a very powerful tool to simulate and investigate a system like the one I have. See for example [101]. The objective is to produce an arrangement of particles that follows the equilibrium distribution for the classical 2D plasma, equation (6.11). The standard scheme is the Metropolis algorithm which in a nutshell consists of the following steps: Starting from an initial configuration of particles, a particle is chosen at random and moved to a new position at a random direction, for a distance equal to a fraction of the magnetic length Il. If that move results in a negative energy change AE, the move is accepted. If not, it is accepted with a probability exp (2AE), since / = 2. This scheme ensures that in the limit where the number of Monte Carlo steps tends to infinity, the particles will have an equilibrium arrangement that corresponds to the minimum of the energy. In reality, one has to devise tests to ensure that the system has reached its equilibrium configuration, for the Monte Carlo steps taken. The code that I developed for the Monte Carlo simulation is presented in Ap pendix E. Here I will focus on the details and the results of the simulation. The first issue was to check the code against known analytical results for a uniform field. It can be shown that for the incompressible liquid state that the 2D plasma is in, the radius of the droplet which is formed scales as R = 2N, where N is the number of electrons in the droplet.1 The density of the droplet is uniform and equal to p = 1/27. In addition, the ;" of the droplet 2 can also be found. It has the form E N2 In (2N 1)/4. Notice that in the above formulas dimensionless 'In reality, the particles in the simulated plasma gas are not electrons, as can be inferred by their interaction energy and also the term "charge" is more appropriate than solenoids. However, to avoid confusion I retain the two terms, since the distinction is clear. 2To be precise, it is the energy of the fictitious plasma. units have been used, for which 1B 1 and B = 1. All these theoretical predictions served as consistency criteria for the simulation of the droplet in a uniform field. There was close agreement for various droplet sizes (number of particles) and many initial configurations of the droplet which indicated that the code was indeed working properly. First, the system was allowed to equilibrate and then measurements were taken. The time unit is one Monte Carlo step, the movement of one particle in the simulation. After N moves, all the particles are moved once on average. This defines the Monte Carlo sweep. Equilibration for a typical droplet size of 500 particles took about 103 sweeps, and measurements were taken typically for 4 x 104 sweeps. This ensured good statistics. The quantities measured were the total energy and the density. In the end of the measurement cycle the average density was calculated from a previously calculated histogram of the particle positions. To that end, a inning scheme was eimpl'lid with bin sizes of 0.3 x 0.3. In figure (6.6) I show the results of a simulation for a uniform field, starting with N = 100 electrons up to N = 400 electrons. The energy measured through the equilibration phase is shown for the last simulation N = 400, to illustrate the rapidity of the equilibration process. The measured energy is also shown. The density distribution is presented using pseudocolor pictures, in which blue is low density and red is high density. I find a droplet of uniform density and circular boundary. Before presenting the simulations for the inhomogeneous magnetic field a few comments are in order, on the two possible methods to implement the magnetic field inhomogeneity. Initially, I ran simulations with solenoids as the source of the inhomogeneity. The shortcomings of that method were that it was hard to predict beforehand what is the proper value for the flux and at what distance the solenoids should be placed to have the desired effect, namely drive the cusping of the droplet. These two values are interrelated and the complication was exacerbated by the fact N 100 N = 200 Equilibration energy for N=300 40000 Measuredenergy forN=300 20000 400 o 1 1016e 05 20000 60000 0000 0 10000 20000 30000 40000 40 l+weep  ........ ........ .... 00 sweeps N 300 Energy for N =300 Figure 6.6: Growth of the electron droplet in a uniform magnetic field. that for solenoids placed very close to the droplet and/or high fluxes, particles would breakoff and attach to the solenoids before the cusping occurred, an indication of an instability. However, this method proved to be invaluable in assisting in choice of the parameters of the problem and the initial understanding of the system's complicated behavior. All the above problems were alleviated with the introduction of the inhomo geneity in the magnetic field through the potential V(z) as described in Appendix C. As will be shortly shown, the agreement between the theory and the simula tions is impressive. The growth of the droplet and the shape of the boundary were 100 I 50 50  6* 50 100 ,4 I 50 0 50 100 Figure 6.7: Illustration of the breakup of the droplet for a simulation with the number of electrons N = 600 exceeding the critical number Nritical = 500. The droplet breaks up and the detached particles were accumulated at the edge of the simulation box for counting purposes. Eighty nine detached particles were counted, equally distributed among the three smaller droplets. found to be consistent with the predictions of the theory, with uniform density in the interior and a sharply defined boundary. Although the semiclassical theory breaks down at the cusping point, the simulations were carried out past that point with the same findings in each case that I tried. The droplet would break up, with the excess number of particles detaching from the tips of the cusps and moving outwards. In all cases the droplet would stabilize with a number of particles equal to the critical number Neritical. In figure (6.7) I show this behavior. In that particular simulation, an artificial barrier was placed at the edge of the simulation box in order to count the detached particles. The critical number was Neritical = 500 and the simulation 68 was run for N = 600 particles. Eighty nine detached electrons were counted, equally distributed among the three smaller droplets at the boundary. I begin with the droplet with three cusps. The appropriate conformal map is z(w) = aw + bwM, with M = 2. The appropriate potential for that situation is V(z) = t3z3/2 and t3 is related to the critical number of electrons through Ncritical 1/144t2. For the details of the derivation the reader should consult Appendix C. N =100 N = 250 N = 350 Figure 6.12: Growth of an electron droplet N 494 with three cusps. For the simulations I chose Ncitical = 500. In figure (6.12) I present the growth of the droplet for various numbers of electrons until the cusping point, showing a few representative density plots for N = 100, 250, 350 and 494. Notice that the cusping la ~s occurs for a number of electrons N slightly less than the critical number. In this case electrons started detaching from the droplet for N = 494 instead of N = 500. In figure (6.13) I show the comparison between the semiclassical curves, colored red, with the results of the Monte Carlo simulations, colored blue. 40 20 0 20 40 20 0 20 40 60 Figure 6.13: Comparison of semiclassical (red) and MC (blue) results. Although in the above simulation the appropriate potential was used, it is also possible to obtain exactly the same results using six solenoids. The solenoids have to have alternate positive and negative flux ( and should be placed at the vertices of a regular hexagon, on a ring of of radius R,. The flux ( and the distance R, are related to the relevant harmonic moment through ts = 4(./R>. Next, I turn to a simulation with two solenoids, placed on the x axis at opposite points zo, each carrying a flux ( = q; this results in a droplet with two cusps. The appropriate conformal map is the modified Joukowsky map z(w) = rw + 2Qw/(w2 N =100 N 400 Figure 6.18: Growth N 465 of an electron droplet with two cusps. w5) and the two sets of parameters are related through r zo + wo 2q r  rQ wN 1 f N Ir2 2 2Qwo 1 wo 2Q2(1 w4) (1 wo)2 4Q(1+ ) (1 w4)2 (6.26) (6.27) (6.28) For this case also, the chosen critical number of electrons was Ncritical = 500. In figure (6.18) I present the growth of the droplet for various numbers of electrons until the cusping point, showing a few representative density plots for N = 100, 250, 400 and BO 80 N = 250 465. In figure (6.19) I show the comparison between theory and experiment using the same colorcoding scheme as before. 6':4 40 4: 11 II 1 1 Figure 6.19: Comparison of semiclassical (red) and MC (blue) results. I carried out another simulation which produced beautiful pictures for which, unfortunately, very little analysis can be carried out. An annulus of 350 repulsive solenoids was placed in the middle of the simulation box with many solenoids dis tributed randomly in it. Fifty attractive solenoids were placed in the exterior of the simulation box, equivalent to a large distance, to enforce "charge" neutrality. The simulation resulted in a very complicated density plot as can be seen in figure (6.26). Finally, a simulation for an also incompressible state of the Fractional Quantum Hall Effect, with filling factor v = 1/3 is presented. I used the same potential V(z) = t3z3/2 used in the simulation for the droplet with three cusps. The results were N =100 N =180 N 220 N 260 Growth of an N = 300 a random .i i., . of solenoids. 40 30 ~o~ro Figure 6.26: electron droplet within N = 140 73 not what I initially expected. Based on theoretical grounds, the v = 1/3 state was expected to have smaller density fluctuations in the interior of the droplet but the observed behavior in the simulations was quite the contrary. As seen in figure (6.33), there is an enhancement of the density near the boundary and also, what is even more perl : ::ii _. is that the breakup of the droplet occurs in a different manner than for the = 1 case. It seems that a smaller droplet detaches from the larger initial droplet. I do not have an explanation for this behavior, although the enhancement of the density near the boundary has been encountered in previous numerical simulations of the Fractional Quantum Hall effect. 74 W W a0 a C 20 4 60 W 40 6 40 W 1 20 l N K N 60 N 120 a W 40 =2 C 20 4 60 N 40 W 40 M f 20 4 N K N 140 N 160 60 40 2O 60 a 40 M 40 2 20 4O 6 N 260 N 167 Figure 6.33: Growth of an electron droplet for v = 1/3. CHAPTER 7 CONCLUSION I have investigated the structural phase transitions of the vortex lattice within the anisotropic GinzburgLandau model, the anisotropy arising from s and d mixing or from an anisotropic Fermi surface. I have shown that the addition of the phenomeno logical higher order derivative term, is sufficient to account for the observed behaviors of the vortex lattice. I derived an extension of the virial theorem for superconduc tivity to anisotropic superconductors, which enables one to write the free energy in a compact form from which the derivation of the nondispersive elastic moduli is very easy. The particular extension of the virial theorem is shown to be a tool which facil itates the generalization of the famous Abrikosov identities along with other results obtained from the isotropic GinzburgLandau theory to the anisotropic theory. The rhombic to square structural transition was studied and its effects on the elastic moduli were analyzed. This model exhibits the same behavior at intermediate fields as the nonlocal London model, vanishing of the squash modulus cq at the rhombic to square phase transition, vanishing of the rotational modulus c, at the point where the vortex lattice exhibits rotational instability. At high fields, near H,2 the moduli vanish. I incorporated the thermal fluctuations in this model in a simplified manner, not taking into account the effect of the thermal fluctuations on H:(T). Nonetheless, this approach proves sufficient to show that the reentrant transition from square to rhombic vortex lattice, first encountered in the nonlocal London model, is present in this model also. The mechanism that causes the reentrant behavior is the thermal smearing of the anisotropy. Although the anisotropic term is not derived rigorously from a microscopic model, it succeeds in encompassing all the interesting, and at times unexpected, phe nomena pertaining to the structural transitions of the vortex lattice. I also investigated the fingering of electron droplets in the Quantum Hall regime driven by the influence of an inhomogeneous magnetic field via the AharonovBohm effect. Carrying out detailed Monte Carlo simulations of the growth of the electron droplets for various field inhomogeneities, I was able to compare the theoretical map pings I found for each inhomogeneity to the experiments. I found that the development of the boundary of the electron droplet up to the cusping point was described with remarkable accuracy by the curves obtained from the proper conformal map for each inhomogeneity. In addition I was able to predict the critical number of electrons for each particular arrangement of solenoids outside the electron droplet. Two methods were developed for the simulation of the magnetic field inhomo geneity. One in which the inhomogeneity is provided by thin solenoids distributed in the exterior of the droplet and another in which the effective potential that affects the electrons of the droplet is derived in a systematic way. Using the latter I was able to study the growth of the electron droplet after the number of electrons exceeded the critical number at which cusping occurs. I have found that the excess electrons were detached from the droplet and moved towards the boundary of the simulation box. APPENDIX A EXTENSION OF THE VIRIAL THEOREM In this appendix I show the derivation of the Virial Theorem for the anisotropic GinzburgLandau theory. As will become soon apparent, the Virial Theorem facili tates greatly the simplification of the subsequent analysis of the problem at hand. The derivation is essentially the generalization to the anisotropic case of the derivation by Doria et al. [73] for the isotropic theory. The starting point is the familiar GinzburgLandau free energy / ll 12 l) 4 I Vf (X) 12 [V x A(X)12 _,2,\ 2 F aI (x) H(x)2 (x)[Vx+A(x) + 2 2m* 8r (A.1) InT I.(x) [,(I [ '(x)][II(I ) (x)] + cc where V x A(x) = B(x) is the microscopic magnetic field. The magnetic induction is defined as B =f d3x B(x)/V and the homogeneous applied field is denoted by H; the brackets denote spatial averaging. I can also define the quantities Fkinetic W (X) 12 (A.2) /B2(X) (A.3) 87 Fanisotropic __ a, [n ( y X. (A.4) [nI'. .(x)][I,( I H (x)]*+cc lx y )] + CC , by identifying the ]lr, i, .,1 meaning of each component in the GinzburgLandau free energy. Under a scaling transformation the various components of the GinzburgLandau free energy transform in the following manner X x' =. (A.5) Y(x') (Ax'), (A.6) A,(x') AA(Ax'), (A.7) nH (A.8) B A A2B. (A.9) The next step is to take the partial derivative with respect to A and set it equal to zero. I neglect the Adependence of the order parameter and the vector field because the GinzburgLandau free energy is stationary under variations of these fields by definition. Finally, after setting A = 1, I obtain OF 2Fketic 4Ffield 4Fanisotropic + O2B = 0, (A.10) aB which can be recast in the following form which contains the applied field, using the thermodynamic relation H = 47rF/OB: HB H_ (x) 2 B2 X) 2 H + + (27, fY(x) [II,(I,2 2) )X) S2m* 47 a [T, ,./(x)][I,(I,2 _ 2)(x)]* + cc (A.11) The above expression can be simplified even further with an integration by parts. One has to use the appropriate boundary condition (2.10) ~n {q[(Yd) [yly(II 2 ] n. ..(n C n nl)d xn,(n nl ) o. (A.12) a, s L X One is thus led to the following simplified result, which will be used in the derivation of the generalized Abrikosov identities B 47 d3 H I (x) 12 B2(X) H.B 2* 4 V 2m* 47 ) [I(xk Y y +cc1 (A.13) \ Qs JJ A.1 A Useful Identity Starting from the reasonable assumption that the free energy is stationary under y + (1 + 6)y, with 6 being a constant, variation with respect to 6 leads to [3] F 2 I 2 In1"212 \2m* ) [*(I 2 n )2+C a 2 + i4 (,,2) 0, (A.14) which must hold for every 6. In this way I obtain the generalization of a well known and very useful identity \ 2m* Knfl 2 2 rn* ( ) [*( l )2 + cc] (,,2) lz (A.15) a\ 2 1 4 . The latter permits one to write the virial theorem in the following compact form, from which the derivation of the Abrikosov identities is trivial: H B a 4 a 2 J1 4 + B ( [ 2 f(n 2 . 4 7 a , (A.16) A.2 Generalized Abrikosov Identities The generalizations of the Abrikosov identities are going to be derived as was done for the isotropic superconductor by Klein and P6ttinger [74]. I begin by assuming that an expansion of the internal field can be written as B(x) = [H x(x)]z and H = [H2 A]z with A and 0 of order 11 2 and A constant. Substituting in equation (A.16) and separating terms of different orders in 12 I obtain Zeroth Order: H2 = HH. (A.17) First Order: (2A + 0)H2 47aa 1(1 12) + 2HC2A + 2HI2 +8 Rte(Y*(I2 I)2y). (A.18) Oas Averaging both sides and switching to the usual dimensionless units (see section 2.1) I obtain () (= II2 Re(*(Un n2)2y). (A.19) Here, c is the dimensionless anisotropy parameter defined as c = 21 ? I/'yas. From the relation B(x) =[H q(x)]z it follows immediately that (x) 1= 12 e n H2)2], (A.20) 1K 2 B= H y 2)+ Rey*(H F)2y). (A.21) Second Order: A2 + A0 43p( f 4) + (A2) + (2) + 2(Aq). Averaging again and using A = K H (switching once more to dimensionless units) I obtain after a short calculation H H_ 2K2 1 12 2x*(n Y) 22 k 2 U (A.22) The familiar Abrikosov identities for the isotropic superconductor are immedi ately obtained by setting the anisotropy e equal to 0 in the generalized Abrikosov identities (A.21,A.22). A very useful relation for the average magnetic field can be obtained combining the generalized Abrikosov identities and it reads K H Re(*( n2IJ R 2)2 e I 2 2 2)2 (2 2 1)A f2) (22 1)/A(1122 2 (A.23) where 3A is the Abrikosov parameter defined as 3A = (114)/(112)2. It is also conve nient to define F which is of zeroth order in  2, as S f 4 e( n2)2Y) + ReiA(2( )) 1) (A.24) (Ifp 2)2 ( I2) which simplifies the following expressions for the average magnetic field B, the cor rection A, and the applied field which in turn will be used in the derivation of the free energy: B =H (2 (A.25) S (2K2 1) A (22 1) A A (22 B) 1(2K2 1)3A + F (A.26) S(212 1) A+ H B (2 1 (22 (A.27) (22 1)A + I T2 K_ 1)A + I A.3 Free Energy, Magnetization, Gibbs Energy The free energy can be written in a very manageable form if one uses identity (A.15). I obtain F 1 F = (14) + (B2) 14 + (( A 0)2) 2 2 S(K A)(K A 2(9)) ( 1 +  ) 2). (A.28) 2 Using Abrikosov's identities again, I obtain I + (1 2) = A(B A), (A.29) 2 and now the free energy (A.28) can be written F = (B K)(r A) + BK. (A.30) Then, from (A.26) one obtains easily the generalization of the isotropic result for the free energy F(B2 B)2 F F K B2_ (1 ( 1 } (A.31) (2 1) + 1 (2 K2_ )A + 1 The magnetization M = (B H)/47 and the Gibbs free energy G = FHB/47 are calculated easily from the previous results. I find MH F (A.32) 47(2K2 1)A (2K22 1)/A (A.3 (G 1 } (A. 33) (2K2 1)A (2 2 1)3A In short, I have managed to significantly reduce the complexity of the original GinzburgLandau free energy, using the generalized Virial Theorem and at the same time obtain simple expressions for important relevant quantities. APPENDIX B CALCULATION OF /A, ( *(n2 II2)2), ( 2ip*(n I2 )21). In this section I show the details of the calculation of the necessary spatial averages of the order parameter. The calculations are done following the conventions and methodology of Chang et al. [79]. The method is essentially the same if one wants to include the second order correction. I will show the details for completeness. In order to construct the periodic solution one can form the following linear combination which can be tuned so that Il2 acquires the desired periodicity that the simple solution of the linearized first GinzburgLandau equation lacks, as S27min 27n P(x, y) c, exp ( y)(x ), (B.1) n=oo where a and b are the two basis vectors and b is the reduced field. This function is periodic in the direction of a (the y direction). One can impose periodicity in the b direction by requiring that T(x b sin 0, y + b cos 0) = (x, y) holds. This is accomplished by setting b sin 0 = 2r/bcK2a, a condition which implies that one has one vortex per unit cell. I then choose a coordinate system (X, Y) that coincides with the vortex lattice directions so that (x, y) = (Y sin 0, X Y cos 0). The spatial averages are going to be calculated integrating over one lattice cell, namely 0 < Y < b and 0 < X < a, but the final integrals will be extended over all space taking advantage of the periodicity of I2. The integrals involving I 12 are the easiest to evaluate. The integration over Y will result in a Kronecker 6 which will help in the evaluation of the first summation. The second summation will extend the integral over all space. I have The calculation of (I1 T4) is slightly more complicated. Now a product of four wavefunctions is involved and I obtain ('Fl4) bsin 0 +i 22 ib cos0 n2 m 2 + 12 q2 2 d (B.4) where j (x jb sin 0). One can choose new variables Z n m + dxq, N 1, M m q, which facilitates the evaluation of the sum under the constraint of the particular Kronecker w. The sum over i i 1, q, breaks into two sums one over Z and one for even and odd N and M, respectively. In this case it is the integration over Z that will extend the integration over one cell to an integration over all space. Before calculating the integrals, one needs to rotate the derivative term by an angle 9 about the z axis to account for the general orientation of the vortex lattice. The calculation of the integrals is tedious but straightforward. I show the results below [the prime denotes a rotated term following the notational convention introduced in (2.1)]: (2) + (cb), (B.5) 8 ( 4) =fo + (e6b)31 + (eb)2[ 1 7+/2 + 521], (B.6) (p*(n~ n_ )'f ) = 2{1 + 3(cb)}, (B.7) (2* 2 y'2 I + o 1 + (eb) 23o + 631 + 2(71 + 72) (B.8) The functions 3o, 31, 71 and 72 have the following form A = Anm, nm S3i Re e4i Anm [8722 4 'y e~e8CPAnm 161Fc7 71 Re e 8 A.n 16n4(4 n8 nffn (B.9)  6ran2 + , 12_r33 n6 2+ 3 _2( 4 (B.10) 45 2 1051 16 256J 72 Anm 4 [16rU4 nm4 12n 1 3n2M2(n2 + )+ 72( + m4 nm 45 o 2 105] + 36n2m2) (n2 + m2) + 05 16 256 (B.11) (B.12) where Anm = e2i"p(n2' 'e27(n2+2). The prime has the meaning that there are two summationsthe one shown, over n, m and the other with n and m replaced by (n + 1/2) and (m + 1/2) respectively. APPENDIX C MOMENTS In this appendix I present the calculation of the external and internal moments of various useful conformal maps. They are related to the harmonic moments through t k (C.1) 7ik I make use of one little known but important theorem from the theory of complex variables. For a function f which is analytic inside, outside and on a simple closed loop C, with the exception of a finite number of n singular points Zk inside C, with k = 1...n, it holds dz f(z) = 2 Res f(z) = 2Res f(o (C.2) C.1 Conformal Map z = f(w) = aw + bwM The conformal map z = f(w) = aw+bwM is appropriate for the investigation of the cusping of the droplet with M+ cusps, with their tips residing on a ring at regular angles m7/(M+ 1), with m = 0, 1,..., M+1. The derivative f'(w) a Mbw(M+l) has zeroes at (Mb M+I 2'ik .,= e+ with k 0 ... M, (C.3) which are in the interior of the unit circle wl 1 provided that Mb/a < 1. When these roots fall onto the unit circle, the velocity v = dw/dz diverges, signaling the formation of a singularity. C.1.1 Exterior Moments I need to evaluate the external moments in order to find the connection between the conformal map and the harmonic moments of the solenoid distribution outside the droplet. I start with the definition Cj  exterior z3 dx dy= Idwff'f . (C.4) The calculation for j = 0 is straightforward and yields co r(a2 Mb2). For the other moments I have to evaluate the integral c f dw (a2 Mb2) Mabw(M+l) + abw(M+l) jM S 2i w V) [awM+1 + b]l' II I2 13 with f (w) = wM/[awM+1 + b]j having (M + 1) poles of order j in f M+ i,(2k+1) ,c e M+I with k = 0... M. Using the result (C.2) I easily obtain for the three parts 1i, I2 and I3: 1{ j(M+I)jM1 Il = Resw=o, 00, forj(M + ) jM 1 > 0 = j > 1 and also SR o j(M+ )jM+M 0 I2 =Res w=o [=1.I 0, (C.5) (C.6) (C.7) (C.8) (C.9) for j(M + 1) M + M > 0 = j > M, along with 1( o j(M+1)jMM2 Is = Resw=o[a + =1 0, (C.10) for j(M +) jM + M 2>0 j > M + 2 with the exception of j =M + 1 for which there is a residue at w = 0 with the value a~. This is derived from the expansion w 1 qI +(i)M+lqM+l q2(M+l) (C.11) waJ(1 + qM+1 )j ""' with q = (b/a)1(M+l)w and as required q < 1. I find thus that only two moments have values other that zero co 7(a2 Mb2), (C.12) CM+1 M. (C.13) For the calculation of the critical number of electrons Ncritical at which the cusps occur I notice that when the roots of the conformal map fall onto the unit circle I have b = a/M. Solving the two equations for the exterior moments for co = 2N I obtain M 1 7 7 1 Neritical l (C.14) 2M MC t The time t in the Laplacian growth problem is equivalent to two times the number of electrons 2N in the Quantum Hall droplet growth problem. C.1.2 Interior Moments The calculation of the interior moments follows along the same line and is shown for completeness. Regardless of the closed form obtained, the results are not useful and it actually can be shown that they formally diverge. Substituting j j into the RHS of equation (C.6) I have c I dw (a2 Mb2) 1 Mabw(M+1) + bw(M+l)[aw M+ b (C.15) 12 13 Using again the theorem (C.2) I obtain for 1i, I2 and 13 1 = Reswo j+I b M+1 ' a j J The residue is different from zero for j + 1 j = m(M + 1) with m = 1,2,... I have li r (a2 Mb2)(aMb) (MM + m RS ( M k(M+)1 s,, (.16) (C.16) k(M + 1) 1 = k = j/(M + 1). So, for (C.17) 1)), and similarly for I2, 13 I obtain I2 = Ma2(aMb)" M(M + 1) mM + 1) m) m+1 m (C.18) (C.19) Combining these three results I find the final expression for the interior moments (M) (a m(M ) a2 Mb2 (20) C(M+) = (aMb) M( + l b2 (C.20) m ) +I MM+ I1 with m 1, 2.... One easily sees that these diverge for m  oo. C.2 Conformal Map z = f(w) =aw + M 1 k This is a generalization of the previous conformal map. I easily obtain following the methodology developed in the previous section the exterior moments SM co = ( a2 kb k 1 c3 = abj_1  a e (C.21) (C.22) M Skbkbk+j 7,^ 1 ) with j = 1, 2,..., (M + 1) and the implicit assumption that bo 0 and bj>M 0. In this case only Co, Cl, ..., CM+1 are nonzero. The calculation of the interior moments is forbiddingly difficult. C.3 Conformal Map z = f(w) = rw + Q/(w wo) + Q/(w + Wo) This is a modified Joukowsky map which corresponds to a droplet with two cusps on the positive and negative x axis respectively. It has the form Q Q 2Qw z = f(w) = rw + + w + w Wo W + Wo W W (C.23) The derivative of the conformal map has zeroes at f'(w) = 0, which gives w w + + 2 4( w2, 0 ?g (C.24) under the constraints that 0 < wo < 1, Q/r < 1/2 and wo > 2Qr. The external moments can now be obtained as before L dw r 2QW cyj ( WW r2 2iJw w 1 (ww) JL 2Q(w22 (2 w 2) (w2 w)2 J[( (C.25) (~w2 W2\2 2 + 2Q Although the general evaluation of the preceding integral is out of the question in this case, I can evaluate the first few nonzero external moments. With the aid of Maple for the nontrivial calculations, I obtained the following expressions for the first three nonzero external moments Co = rrr 4Q(1 w (C.26) c2 r W4)2 2QW(1 + ) (C. 27) t z z, (C.29) 2~Qw 2(1 w )2 C4 2IQwg( 0)2 I_4 w 2Qw2(+ W4( (C.28) [r(1 w4) 2Qw1]4 00 C.4 Connection with Magnetic Field Inhomogeneity In order to compare the theoretical results with the Monte Carlo simulation I need to connect the harmonic moments calculated in the previous sections with the harmonic moments of the field inhomogeneity defined as tk f 6 k d2z, (C.29) where 6B is the magnetic field inhomogeneity. This in turn will enable us to find the expansion of the potential V(z) which has the form 1 0O V(z) = tkZk. (C.30) k=1 For an arrangement of thin solenoids like the one I use in the simulations, the harmonic moments tk can be evaluated for the inhomogeneity has the form 6B = j (2) (Z n), (C.31) where (, is the flux that the nth solenoid encloses. I find that tk k (C.32) The harmonic moments tk are connected to the exterior moments through the simple relation tk (C.33) 7rk At this stage there are two ways to proceed. The objective is to compare the evolution of the conformal map in "time" (equivalent to the number N of electrons in this case) with the actual simulation. The "time" enters through the zerothorder exterior moment co = 2N. Therefore there is the obvious option of solving the equations provided by the calculation of the exterior moments which will be proven difficult to carry out in general and the not so obvious at first glance option of calculating the potential V(z) first. Both methods will be presented next. C.4.1 Exterior Moments I will focus on the solution for six solenoids on a ring of radius Rs each of which carries a flux ( with alternating signs as one travels on the ring, distributed uniformly at angles which are multiples of 7/3. The appropriate conformal map for this problem is z = f(w) = aw + bwM, with M = 2. From equations (C.32,C.33) I have co = 27N, (C.34) 127r(I ca 2 (C.35) RS which are the first nonzero exterior moments for that particular arrangement of solenoids. The solution of equations (C.12,C.13) yields after choosing the appropriate roots 2 a 2 i V/ 8coC3 (C.36) b +1 2 / SCC3 (C.37) 47T with cusps forming for coc3 = r3/8 from which I can find the critical number of electrons Ncritica ( ) (C.38) I can obtain the same equation for the critical number Nritical by direct substitution of the parameters in equation (C.14), which serves as a consistency check. The connection with the experiment has been established in this particular case as one can plot the curves obtained from the conformal map z = f(w) = aw + bw2 for a given number of electrons and compare directly with the results of the Monte Carlo simulation. Unfortunately an analytic approach of this kind is possible for the simplest of cases and fails for higher M because of the lack of analytic solutions of high order polynomials. In general, for given M one has to solve a polynomial of order 2M. There is a way around that difficulty and the method is presented in the next section. C.4.2 Potential V(z). First I need to introduce the Schwarz function for a curve C. A closed curve C in the plane is described by an equation of the form f(x, y) = 0 which in terms of complex coordinates can be written as g(z, z) = 0. For a function g which is analytic, I can solve for z and thus obtain the Schwarz function S(z), = S(z), on C. (C.39) The Schwarz function is a very interesting object, with many properties and applica tions. The interested reader can refer to a nice book by Davis [102]. For my purposes I need the exterior expansion of the Schwarz function S() (z) which acquires meaning if I analytically continue S(z) to a striplike domain which contains the curve C. It has the form [100] S(+z) = S ck1, (C.40) k1 where the external moments can be defined in terms of the Schwarz function as follows Ck I I kS(z)dz, with k = 1,2,.... (C.41) Recalling that the potential has the expansion [86] 1 V(z) Re t, (C.42) k=1 I immediately find that the potential V(z) and the external expansion of S(z) are related by V'(z) = ReS(+)(z). (C.43) 2 The objective now is to find the Schwarz function S) (z) for the conformal map that I choose. I start with the conformal map z(w) = aw + bwM. Inverting the equation and expressing w(z) as a Laurent series I obtain Z 1,if1 w(z) + O(z(2M1)) (C.44) a zM For the Schwarz function I have a b M a2 M!b2 S(z) = a + '(z) zM + 2 + O(z(M+2)). (C.45) w(z) aM z Comparing term by term with the expansion S(z) = (M + 1)tM+ilM + t + O(Z2), (C.46) z I obtain the expansion coefficients of the potential V(z), namely tM+l = baM/(M + 1) and t = a2 Mb2. The connection with the numerical experiment has again been established, as one can use the derived potential to simulate the droplet and subsequently compare with the theoretical results. What is most important, all the difficulties related to the analytical solution for this arrangement of solenoids have vanished. With regard to the simulation, the route one follows is simple. First, by the appropriate tuning of b I choose a convenient value for the critical number of electrons Ncritical, equation (C.14). Then the simulation is carried out with the appropriate tM+1. In the end, the curves obtained from the conformal map at different lim. ," are compared to the boundary of the droplet at the corresponding number of electrons. The modified Joukowsky map is also simple to treat using this method. Inverting the conformal map I obtain z x3 z2 6rQ+ 3r22 W(z) = + + (C.47) 3r 3r 3rxz where the parameter x has the form x z3 9z9zrQ zr2w + 3r z2(3OrQw2 + 6r2w 3Q2) 3z4w 3r(rw 2Q)3. (C.48) The Schwarz function is S(z) = r/w(z) 2Qw(z)/[l (w(z)wo)2]. The first term has no poles and the poles of the second term are found by solving w2 () = w2. With the help of Maple I find the simple result o + ~ (C.49) w( 1 0) 