Algorithms and analysis for simulation of deterministic and stochastic Ginzburg-Landau models

(1)

ALGORITHMS AND ANALYSIS FOR SIMULATION OF DETERMINISTIC

AND STOCHASTIC GINZBURG-LANDAU MODELS

by Ty Thompson

(2)

A thesis submitted to the Faculty and the Board of Trustees of the Colorado School of Mines in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Mathematical and Computer Sciences).

Golden, Colorado Date Signed: Ty Thompson Signed: Dr. Mahadevan Ganesh Thesis Advisor Golden, Colorado Date Signed: Dr. Willy Hereman Professor and Department Head Applied Mathematics and Statistics

(3)

ABSTRACT

We focus on the efficient simulation of nondeterministic critical phenomena in the Ginzburg-Landau (GL) model for superconductivity. Deterministic GL is widely used to study the formation of vortex configurations in thin superconductors. When such phenomena are essentially nondeterministic, the Langevin version of the time dependent GL problem applies, and simulation presents a significant computational challenge.

To investigate nondeterministic dynamics, we work on 2-manifolds having rotational symmetry about the axis of a constant magnetic field, and consider an ideal (superconductor or normal) initial state. Using highly efficient deterministic schemes, we first motivate the stochastic extension, and demonstrate an ad hoc ap-proach for simulating the evolution of dense, locally stable vortex configurations. For stochastic simulations, we identify the opportunity to use a spectral Galerkin approach. Building upon a linearized Crank-Nicolson scheme, we use generalized spectral decompositions to reliably achieve a reduction in dimensionality. The efficiency of the method is examined through comparisons with Monte Carlo estimators, and our efforts at qualifying the perturbation approach are discussed.

(4)

LIST OF FIGURES

2.1 A comparison of FEM approximations obtained for the L = 9, j = 10 eigen-pairing on the planar disk S1with corresponding reference approximations is shown. A good agreement with

our optimal accuracy estimates is clearly apparent. . . 83 2.2 Selected error measurements taken from the results of an r = 6, λmax= 25 FEM computation

when S is a sphere of raduis 4 are shown. . . 86

3.1 Contour plots depicting GIE estimates of the modulus squared of a Galerkin solution on the planar disk S1at time t = 1 are shown. . . 130

3.2 Contour plots of the error in GIE and GCN simulations of a Galerkin solution on the Bessel surface S3 are shown. . . 136

3.3 A graphical representation of the results included in Table 3.3 is shown. A large improvement in execution times, while retaining the accuracy properties of the GCN scheme, was clearly achieved with the LGCN implementation. . . 148 3.4 Surface contour plots of the error in GIE, GCN and LGCN estimates of a Galerkin solution

on the hollow sphere S7 are shown. With execution times fixed at around 3 hours, LGCN

clearly outperformed both GIE and GCN. . . 149 3.5 Data plots of the product j_M∗ L∗_M versus M for three choices of S are shown. The validity of

the quadrature complexity approximation O(M ) is clear. . . 175 3.6 Execution times from a single step of the GIE algorithm with LAPACK on a hollow sphere are

shown. Cubic complexity in M is observed, and the approximation of the Jacobian consistently requires more time than the serial system solver. . . 178 3.7 Speedup in Jacobian construction times at several mode levels are shown. Performance is

better sustained with increasing p at the highest mode levels. . . 179 3.8 Output data from the first step of the GIE algorithm with LAPACK, taken with M = 1500

on 256 processes, is shown. Small communication times and quadratic Newton convergence are apparent. . . 179 3.9 The experimentally determined serial fraction corresponding to the data given in Figure 3.7

is shown for p ≥ 16. Speedup degradation is due to serial computations preceding the use of DGESV, rather than parallel communications overhead. . . 180 3.10 Execution time ratios between two LGCN implementations, measured during the “corrector”

step at various choices of M and p, are shown. The benefits of using Jacobian symmetry rules with DGESV are clear, as is the need for larger problem sizes when using SCALAPACK. . . 180 3.11 An execution time and speedup comparison for two implementations of the predictor-corrector

LGCN step at several mode levels is shown. Our LAPACK/DGESV approach was faster in all cases, but the optimal performance of the SCALAPACK/PDGESV technique is steadily evolving. . . 181 3.12 Speedup measurements for the LGCN/PDGESV technique at elevated choices for M are

shown. Overall optimality of the predictor-corrector step is achieved, and optimal PDGESV performance begins to emerge. . . 182

(7)

3.13 Modulus squared contour plots of the computed solutions obtained from nonadaptive LGCN simulations on S4, at several time steps and choices for nθ, are shown. The impact of

unmit-igated computational errors is clearly apparent. . . 183 3.14 Rowstacked representations of the results in Figure 3.13 are shown. The violation of E0

invariance for n ≥ 400 in each computed solution, corresponding to their observed spatial symmetry loss, is clearly apparent. . . 185 3.15 Relative error approximations, computed from the results of a collection of nonadaptive LGCN

simulations over the projected Galerkin space P2402E0, are shown. The nearly optimal

per-formance at low Galerkin dimension indicates an excellent level of spatial accuracy. . . 186 3.16 Relative error approximations, computed from adaptive and nonadaptive LGCN results for the

Meissner initial condition on the hollow sphere S4, are shown. Saturation effects appropriately

diminish with increasing λmax, and nearly optimal performance is apparent in the largest

Galerkin spaces. . . 187 3.17 A simulated vortex ring nucleation, obtained with adaptive LGCN from a critical perturbation

of the global limit approximation for the Meissner initial state on S4, is shown. The simulation

achieved an apparent equilibrium approximation, but stability computations suggest that another transition may occur. . . 189 3.18 Modulus squared contour plots generated from adaptive LGCN results for a pair of critical

perturbations of the 15 vortex ring state are shown. The innermost vortex formations closely resemble an Abrikosov lattice. . . 195 3.19 Modulus squared contour plots generated by the adaptive LGCN implementation for a pair of

critical perturbations of the 15 vortex ring state are shown. The trajectories themselves are clearly different, but both eventually achieve the same locally stable configuration. . . 197 3.20 Contour plots generated from the results of an adaptive LGCN simulation for a critical

per-turbation of the 15 vortex ring state are shown. The corresponding equilibrium approximation that results is highly symmetric, but slightly unstable. . . 198 3.21 Contour plots generated from adaptive LGCN for a critical perturbation of the 15 vortex

pentagonal state found by a previous simulation are shown. The equilibrium approximation that results eventually regains the symmetry of the initial state, while also achieving local stability. . . 198 3.22 GL energy trajectories, computed from LGCN simulation results produced by the

perturba-tion technique, are shown. The locally stable equilibrium approximaperturba-tions found have similar energy levels, but exhibit a variety of spatial symmetry. Evolutions toward an Abrikosov-type arrangement achieved the lowest energy levels. . . 199

4.1 Approximations for the first moment of R₀Tku(t, ·)k2dt, obtained from LGCN-GSD-SI sim-ulation results for a stochastic TDGL problem on a planar disk, are shown. Evidence of approximately optimal spectral Galerkin performance is absent at R = 4 and R = 8, but clearly emerges at the reduction level R = 16. . . 236 4.2 Computed ˜µ values obtained from LGCN-GSD-SI simulations over a wide temporal

discretiza-tion range are shown. At a reducdiscretiza-tion order of 16, more than a 7-fold increase in the Galerkin dimension affected only minor changes in the results; at m = 512, the observed change was less than 0.003%. . . 239

(8)

4.3 A comparison of moment approximations obtained from MC and LGCN-GSD-SI results for a coarsened Langevin GL model is shown. A good agreement between the mean and standard deviation approximations produced by both methods at m = 32, M = 421 can be clearly observed. . . 241 4.4 A comparison of µ approximations by the MC and LGCN-GSD-SI methods for a coarsened

Langevin GL model is shown. Improved precision relative to p in the LGCN-GSD-SI-based approximations is apparent when either R or icnt is increased. . . 242 4.5 Approximations for the error in moment computations produced by MC and LGCN-GSD-SI

at m = 32, M = 421 are shown. To a reasonable approximation, reducing the error in the MC estimator to that achieved by the LGCN-GSD-SI method at R = 16, icnt = 15 requires 12 million LGCN simulations. . . 244 4.6 Eigenspace-projected moments computed from LGCN-GSD-SI simulation results taken at

m = 512, M = 421, p = 3 and R = 16 are shown. The moments are ordered according to the size of their associated eigenvalues, and no significant growth along directions orthogonal to the forcing domain is expected. . . 246

(9)

LIST OF TABLES

2.1 The algorithm used to compute approximations for the pairings ˜λj,L, ˜χj,L, where L ∈ Z and

˜

λj,L∈ [0, λmax), is shown. . . 82

2.2 The algorithm used to obtain an orthonormal basis for the space spanned by consecutive FEM eigenfunction estimates having nearly equal eigenvalue approximations is shown. . . 87 2.3 FEM approximation data for L = 0 on the sphere of radius 4 is shown, both before and after

the GS post-processing orthogonalization procedure was used to update results having nearly equal eigenvalue approximations. . . 87 2.4 Measured relative eigenvalue differences and their corresponding approximations for r,

ob-tained from the results of FEM simulations on the conical surface S11, are tabulated. Optimal

performance is observable over a wide range of discretization choices. . . 88 2.5 Eigenfunction approximation differences, along with the experimental values for r they yield,

are shown. The data was obtained from the same FEM simulations used to construct Table 2.4, and is again in agreement with the optimal accuracy predictions. . . 89

3.1 Error and rate of convergence estimates computed using results obtained from an implemen-tation of the GIE scheme on a parallel computing system are shown. The data indicates, as expected, the validity of the GIE error estimate O(∆t) when u(t, ·) ∈ G_2M, but we observe large execution times when relative error is below 1%. . . 129 3.2 Error and rate of convergence estimates computed from the results of GIE and GCN

simula-tions are shown. The validity of the GCN error estimate O(∆t2_{) when u(t, ·) ∈ G}

2M is clearly

apparent. . . 136 3.3 Error and rate of convergence estimates computed from the results of GIE, GCN and LGCN

simulations are shown; large reductions in execution time and second order in time accuracy was achieved by the LGCN implementation. . . 148 3.4 The algorithm used for the GIE and GCN schemes with Newton’s method is shown. . . 158 3.5 The algorithm for efficiently estimating one of the terms in the Jacobian matrix required by

Newton’s method in the GIE and GCN schemes is shown. . . 160 3.6 The algorithm for efficiently estimating one of the terms in both the Newton’s method Jacobian

matrix for the GIE and GCN schemes, and matrices in the LGCN scheme, is shown. . . 161 3.7 The algorithm used for the LGCN scheme is shown. A considerable simplification over the

GIE and GCN algorithms with Newton’s method is apparent in the main loop, as no iterative technique at step n is required. . . 164 3.8 A computational implementation for the two dimensional quadrature rule called when

approx-imating inner products during the GIE, GCN and LGCN algorithms is shown. The complexity estimate O(nsnθ) can be easily seen. . . 167

3.9 The procedure for approximating a Galerkin function on the quadrature grid Sns,nθ is shown.

The algorithm allows for parallelization by Galerkin mode, resulting in speedup even when all processes await the results from a single input vector. . . 168

(10)

3.10 The parallel algorithm for approximating the Galerkin component vector for f ∈ L2_{(S; C)} is shown. The partitioning scheme used by the procedure is naturally compatible with the algorithms it calls, and for efficiency a maximum iteration count terminates the procedure. . 170 3.11 Data collected from an implementation of the Galerkin component approximation algorithm

in Table 3.10 on the torus S9is shown. The observed results reasonably support the a priori

quadrature order rule ns> jmax/2. . . 171

3.12 The adaptive algorithm for computing the three types of Galerkin component vectors required by the GIE, GCN and LGCN schemes is shown. The procedure may be called to construct rows in matrices of the type Ax _{or B}y_{. . . 173}

3.13 The algorithm used to compute a critical perturbation of an unstable GL equilibrium approx-imation is shown. A maximally dissipative direction, from those whose gradient projection onto eigenspace is parallel, is selected. . . 193 3.14 Energy and stability computations produced from adaptive LGCN simulation results are shown.194

4.1 The recursive algorithm used to specify stochastic Galerkin space for a choice of total degree p and sample space dimension N is shown. . . 210 4.2 Example output of the C routine used to execute the algorithm in Table 4.1 when N = 4 and

p = 3 is shown. . . 211 4.3 The recursive algorithm used to generate estimates for integrals of the form

R

RHp1Hp2Hp3Hp4dP (x) is shown. . . 213

4.4 The algorithm used to generate an estimate for E[ξi1ξi2Hi3Hi4] given the stochastic Galerkin

coefficient vectors for ξi1, ξi2 ∈ Pp R

N_{is shown. Both symmetry and orthogonality are used}

to reduce computational costs. . . 213 4.5 The algorithm used to generate an approximation for E[ξi1ξi2ξi3ξi4] given the stochastic

Galerkin coefficient vectors for ξi1, ξi2, ξi3, ξi4 ∈ Pp R

N_{is shown. . . 214}

4.6 The LGCN-GSD fully discrete spectral Galerkin scheme for the stochastic TDGL problem is shown. Each computation requires that a term in our approximation sequence simultaneously satisfies a pair of natural orthogonality criteria. . . 218 4.7 The LGCN-GSD-SI computational procedure for approximating the GSD sequence

corre-sponding to a choice of u0, S, T, m, M, M and reduction order R is shown. . . 226

4.8 The parallel algorithm used to construct the first system matrix A1(˜y, ˜W ), given ˜W and

approximations for its deterministic and stochastic building blocks, is shown. . . 229 4.9 The parallel algorithm used to construct the second system matrix A2(˜y, ˜Y ), given ˜Y and

approximations for its deterministic and stochastic building blocks, is shown. . . 230 4.10 The parallel algorithm used for approximating b1 y, ˜˜ W , ˜Wn−1, ˜Yn−1, dn−

1

2 is shown. The

main loop iterates only ≈ √p times, so the reduction approach used will not compromise efficiency. . . 232 4.11 The parallel procedure for approximating b2 y, ˜˜ Y , ˜Wn−1, ˜Yn−1, dn−

1

2 from appropriate input

and the preliminary SI data is shown. As was the case in the second system matrix algorithm, additional expectation matrices are required to resolve nonlinear contributions to the result. . 233 4.12 Four of the ˜µ values depicted in Figure 4.2, corresponding to the reduction order R = 16, are

shown. Neither a doubling of the temporal order, nor a 7-fold increase in Galerkin dimension, altered our approximation by more than 0.02%. . . 237

(11)

4.13 Moment approximations computed from MC and LGCN-GSD-SI simulation results for a coars-ened Langevin GL model are shown. Error approximations, and the computational costs associated with each underlying simulation, are included. . . 243

(12)

ACKNOWLEDGEMENTS

First and foremost, I would like to thank my parents and family for their unwavering love and support upon which all of my endeavours have been founded. Next, I wish to extend heartfelt appreciation to my advisor Ganesh, whose unflinching support and expertise guided me through the most challenging aspects of this work. I would also like to thank Dr. James W. Swift from Northern Arizona University, whose knowledge of symmetry phenomena in PDEs aided my approach and understanding immeasurably. This work required the use of extensive computing resources, which were provided by the Golden Energy Computing Organization (GECO). The tireless support of GECO, and the reliable availability of their computing cluster RA, is gratefully acknowledged.

(13)

CHAPTER 1 INTRODUCTION

The Ginzburg-Landau (GL) Model of Superconductivity, first introduced in 1950, remains one of the most inspiring applications of phenomenological principles. Around the time of its inception, which was formally announced in [GL50], phase transitions in matter were being explained in many contexts by working from universally appealing concepts within the framework of statistical mechanics. Landau theory, which now is regarded as one of the earliest examples of mean field theory, lead to the GL model when small scale particle to particle interactions were represented in macroscopic terms. The result provided the scientific community with the earliest known tool for meaningfully quantifying superconductivity, and its predictions regarding second order phase transitions in superconductors were soon found to be in remarkable agreement with the earliest observations of the phenomenon. To add to this early success, within the ensuing decade the first microscopic theory of superconductivity provided by Bardeen, Cooper, and Schrieffer in [BCS57] was used to derive the GL model in [Go59]. Since about 1960, the GL model has been the subject of countless numerical, theoretical, and physical studies. Perhaps its most interesting prediction, from a mathematical point of view, is that vortices may emerge in a mixed superconducting state, and laboratory experiments confirming this have only further contributed to its widespread acceptance as a meaningful tool for predicting superconductivity phenomena.

In this work, we seek to numerically simulate the evolution of vortex formations in a variety of thin superconductors by accurately approximating solutions of the Time Dependent Ginzburg-Landau (TDGL) model. This has been done in many contexts by several authors, but we take on this task in a nondeterministic setting, by working with rotationally symmetric ideal superconductors, then carefully studying the results produced from constant initial states. On one hand, there is ample evidence in the literature suggesting the emergence of vortex configurations from both the fully superconducting (Meissner), and nonsuperconducting (normal) initial states. On the other, in our context such evolutions cannot be explained in deterministic terms, since spatial symmetry loss in time would necessarily be observed.

The form of the TDGL equations we consider has been numerically investigated on a variety of planar domains, and novel techniques were demonstrated for hollow superconducting spheres in [DJ04a] and [DJ04b]. However, to our knowledge, the domain generalizations we admit in this work have never been equivalently considered. Stochastic studies of the TDGL model have also been completed in several forms, for example the work in [DDG01] considered a variety of GL noise models on superconducting squares, but we have found no work in which a spectral Galerkin approach, particularly with an accompanying reduction technique, has been implemented for a comparable problem.

The presentation that follows is divided into five chapters. Chapter 2 is devoted to theoretical and computational aspects of the Galerkin spaces we use to make spatial approximations in our numerical methods. After establishing the existence of a particularly convenient complete orthonormal system, we rigorously justify a variable separation technique, then present a highly efficient, dimensionally reduced Finite Element Method (FEM) in arbitrary order for approximating basis functions on the spatial domain of the problem. To conclude Chapter 2, we demonstrate that our chosen basis is particularly well suited for our investigations, in that it allows for the convenient and efficient representation of spatially symmetric functions. In Chapter 3, we establish a mathematical foundation for the deterministic TDGL problem, then introduce several fully discrete schemes for approximating solutions. A complete optimal analysis is provided,

(14)

and results verifying the performance we expect are presented. We describe the details of our implementation, present the results of parallel performance experiments, then close Chapter 3 with some numerical simulations that motivate our introduction of a stochastic version of the problem. In Chapter 4, a numerical simulation method for the Langevin version of the TDGL model is introduced, detailed, verified, and evaluated. We specify a fully discrete generalized spectral Galerkin scheme built from our most efficient deterministic method that reliably achieves a reduction in dimensionality. After describing our implementation for the scheme in a parallel computing environment, we present and discuss the results of a variety of numerical experiments. In the last Chapter, we summarize the study, suggest improvements, and outline several areas of the research that are worthy of continuing investigations. In the remaining sections of this chapter, we introduce the context of the study, the phenomenological foundations of the GL theory, and model realizations for rotationally symmetric manifolds. The last section outlines some notational conventions and mathematical foundations employed throughout the work.

1.1 Goals and Scope of the Study

This work is, in essence, a mathematical and computational investigation of critical phenomena in the TDGL model. We begin with a “baseline” realization of the GL equations, formed through standard deriva-tion techniques and physical idealizaderiva-tions. We then investigate and characterize the predicderiva-tions of this model, particularly those that correspond to the most crucial and physically relevant initial states. Through these efforts, we demonstrate cases in which the baseline approach is, by reasonable mathematical and com-putational standards, incomplete in its description of the global behavior of the system. This serves to motivate the need for a model generalization, which we obtain from a particular extension born from the phenomenological concepts that originally lead to the Landau theory several decades ago. The goals of the study are primarily computational. Ultimately, we seek to demonstrate how critical transitions leading to ubiquitous superconductor vortex configurations can be accurately simulated in an appropriate stochastic context. Moreover, the stochastic results should be evaluated to identify cases wherein the essence of the predictions of the extended model can be extracted from less exhaustive, nonstochastic methods.

The scope of this research is naturally limited, in the context of GL and the simulation of superconductor phenomena, by the versions of the model we selected to investigate. As mentioned, we work from two versions of the TDGL equation. Perhaps the simplest perspective is obtained when regarding the second as constitut-ing the natural phenomenological extension of the first, used to appropriately account for fluctuation effects originally discarded by the mean field approximations that underlie Landau theory itself. These effects are associated with the temperature of the system, so in this way the stochastic extension we work from is the first “nonidealization” of the baseline model one would introduce to more completely represent the effects of thermal energy in the system. Of course, other nonidealizations of the GL model have been proposed and studied, but none of these are typically used to approximate the emergence of ordered vortex formations that break the symmetry imposed by the geometry of the spatial domain. The work in [DGP92] provides an excellent, well-known summary of the applicability of the GL theory, and when particular nonidealizations of the model may be appropriate. Therein, techniques in which the anisotropy or nonhomogeneity of the superconductor itself may be incorporated into the GL equations (through a pair of temperature depen-dent material parameters) are discussed. We chose not to investigate these alternatives for this research, primarily because it seems more likely that the ordered, symmetric configurations measured and computed by other authors on rotationally symmetric domains are caused by a macroscopic effect, and are probably

(15)

not due to local material nonuniformities, that would normally vary in unpredictable ways regardless of the geometry of the superconductor. For example, in [SPD98] and [Pe00], the authors used the GL model to compute polygonal vortex ring configurations in planar superconducting disks and rings when varying the external magnetic field. These results were obtained for isotropic, homogeneous materials while using con-stant material parameters, and in some instances showed instabilities with respect to fluctuations in the GL order parameter that lead to their classification as metastable. In [Pa98] similar results are presented, and several of the approximations reported therein were obtained, and mathematically characterized, in [Th05] by way of a completely different numerical approach. In addition to these computational results, mesoscopic superconducting disk measurements are reported, along with accompanying calculations, in [KBP04]; these closely resemble the results predicted by the numerical simulations mentioned above. In the context of such experiments and measurements, our interest for this work lies in the dynamics of the formation of the states that were found, primarily as they would evolve from the Meissner initial state, when the magnetic field is introduced as an essentially fixed, normally incident quantity. The multiplicity of states found in the above citations suggests that transitions to ordered, symmetric vortex states are indeed to be expected. Appar-ently, we need to account for the fact that several possibilities may exist, but this contingency is naturally compatible with the nondeterministic mathematical expectations we have already adopted. The timing of nondeterministic GL phase transitions has not, to our knowledge, been accurately quantified on rotationally symmetric domains, but clearly this can be an invaluable tool in characterizing the metastability found in the model by other authors, in turn allowing us to distinguish it from either local stability, or the stability associated with system minimizers. Our use of the Langevin stochastic extension can be regarded as an attempt address the need for accurate approximations of such scenarios.

On the other hand, our research more generally represents a contribution to the mathematical and com-putational understanding of parabolic partial differential equations (PDEs) under a Langevin stochastic extension. The study of such systems has been an established branch of the mathematical sciences for some time. However, extensive numerical studies are less prevalent, mainly due to the computational expenses associated with the simulation of stochastic equations. Recent efforts to explain symmetry breaking phe-nomena through the introduction of additive noise mostly involve ordinary, coupled systems, and rely on reduction analysis to obtain relevant models. In [KC10a] and [KC10b], the authors used a Langevin noise model, along with analytical reduction techniques, to explain symmetry breaking in molecular vibrations. These studies provide a good example of a physical system in which a reasonably small number of non-deterministic phenomena can be modeled through the introduction of random, additive input. Also, they illustrate a case in which a reduction strategy produced useful insight into the nature of a stochastic problem without the use of exhaustive simulations. Another interesting study is found in [BT10], here the authors applied a “stochastic mode reduction” technique to a coupled, nonlinear Langevin system in order to make probabilistic predictions without the help of any computer simulations whatsoever. Such studies highlight the interest in reduction approaches to manage stochastic problems. Indeed, these are of the highest priority, as they typically yield both a meaningful intuitive context, as well as an alternative to taxing computations. Our approach varies somewhat, in that we seek to use expensive, but efficient, stochastic simulations to properly qualify intuitive reduction techniques. This seems justifiable in our case, because a priori reduction techniques would likely require a firm microscopic understanding of vortex phenomena in superconductors. Since progress in this area has been relatively slim, we are essentially reversing the usual approach, and instead admitting an appropriate phenomenological extension, then using advanced computer resources to seek reduced approaches through computational means. We had difficulty finding evidence of any similar

(16)

approach in the literature, likely due to reasons already discussed. However, we again submit that the con-text of our study is well suited to such a strategy, and constitutes an excellent opportunity to demonstrate how large scale simulations can lead to mathematical and scientific insight.

Numerical studies of stochastic partial differential equations have been carried out in many contexts, using a wide variety techniques. To our knowledge, the linearized spectral Galerkin technique (with generalized reduction) we implement in this work has not been applied to the Langevin TDGL equation. Our stochastic simulations began with “non intrusive” MC techniques, then we checked recent additions to the literature to find an improved method. An inherent challenge presented by the extended TDGL equation is the forcing term itself; we know little about it, other than that it is a Gaussian random function of space and time. This essentially leads us to expect a large stochastic dimension, so we knew that the additive nature of the forcing would need to be exploited in some way. Fortunately enough, the orthogonality properties of our chosen spatial basis make spectral Galerkin methods a feasible alternative. Still, “the curse of dimensionality” constitutes a primary challenge, but methods that effectively reduce the working dimensionality below that naturally imposed by a tensorization of spatial and stochastic Galerkin spaces have been demonstrated. There are several good examples in which the widely recognized Karhunen-Lo´eve (KL) expansion has been utilized to accomplish, in essence, a reduction of dimension; in [GKW08] it was used to obtain an efficient collocation technique for a coupled PDE system describing flow in porous media. However, such approaches require additional stipulations regarding the underlying spatial correlation. Having no such information in our case, we consulted the excellent overview provided in [No09] to identify an appropriate reduction generalization. The implementation of such a method, particularly in a parallel environment, is relatively complicated. However, when constructed upon our most efficient deterministic scheme, remarkably reasonable compute times can be achieved. Also, the favorable mathematical properties of our spatial basis, when coupled with the additive nature of the stochastic forcing, yield a scheme that can not only be implemented relatively easily, but also lends itself to computational parallelization in several ways.

To summarize the scope of this work, it is convenient to itemize several primary goals. These roughly correspond with the ordering of the chapters in the document. To assess the progress of our efforts, we will revisit these in the last chapter, as we identify the most promising areas for continuing research. Our primary goals are the following:

(1) Mathematically establish the existence of a complete orthonormal system of eigenfunctions of the elliptic GL Schr¨odinger operator, in a variable separated form. Implement, justify, and demonstrate an efficient technique for the numerical approximation of the eigenfunctions and eigenvalues. Investigate the compatibility of the corresponding Galerkin approximation spaces with the goals of the study.

(2) Appropriately qualify the validity of the GL model, both in the deterministic and stochastic contexts. Lay the mathematical foundation for the deterministic TDGL problem, then implement, justify, and verify several fully discrete schemes for approximating its solutions.

(3) Motivate the use of the Langevin GL extension by numerically demonstrating the presence of unstable global equilibria in the baseline TDGL system. Provide approximations for distinct, locally stable equilibria to further support the use of the probabilistic model.

(4) Implement and detail the fully discrete spectral Galerkin GSD scheme for the stochastic TDGL equa-tion. Qualify its performance through comparisons with the results of MC simulations, and illustrate its most crucial contributions to the modeling of critical GL dynamics.

(17)

(5) Process the results of stochastic simulations to determine where nonstochastic extensions of the baseline simulation techniques can be used to quantify essential characteristics the stochastic solution. Isolate any problems or challenges resulting from the method, and identify areas of future research that are likely to improve accuracy or performance.

Item (1) is achieved primarily in Chapter 2. We place emphasis here both on computational efficiency, as well as on how appropriate the choice of spatial basis is for our investigations of symmetry phenomena in the GL model. Item (2) is managed in Sections 1.2 and 1.3 in the current chapter, along with the greater part of Chapter 3. Item (3) is addressed in Section 3.3.3, where we also include visualizations of critical symmetry loss leading to dense, locally stable vortex configurations. We expect these examples to ultimately be associated with the dominant underlying solutions of the Langevin equation, and that the techniques used to realize them may be refined to recover accurate quantifications of the same. Items (4) and (5) are the subject of Chapter 4. Here, we will emphasize the details of our implementation, and the interpretation of simulation results. In the next section, we provide an overview of aspects of the GL theory that are applicable in the context of our study.

1.2 The Phenomenological Ginzburg-Landau Model

In order to establish the mathematical foundation for the problems we consider in this work, a brief overview of GL phenomenology is necessary. Over the last few decades, many approaches to GL theory have been demonstrated and detailed, but they all yield the GL energy functional as a vehicle for describing multiple phases of matter. In the scientific community, this object is generally regarded as providing an approximation for the order in a physical system near the critical temperature Tcat which a phase transition

occurs. The techniques used in obtaining it fall roughly into two categories. In the first, through a careful observation of the features shared by all mean field theories, one begins (as Landau first did) from a list of postulates, then a sequence of appropriate mathematical steps are taken to obtain a consistent result. The second includes the so-called microscopic techniques, wherein one works directly from a prescribed Hamiltonian to derive appropriate expressions for the energy functional. Both approaches use a “coarse graining” or “coarsening” process to obtain an order parameter that is defined throughout the spatial domain while approximating the partition function of the system, and both rely upon an expansion approximation near the critical temperature. It is important to keep in mind that the functional obtained (through either approach) is known to provide an incomplete description of phase transition phenomena. In particular, it was realized early on that as the system temperature approaches the critical temperature, Landau (mean field) theories lose validity due to the emergence of a large correlation length near the critical transition. Here fluctuation effects play a dominant role in the thermodynamics of the system over large distances. Since mean field theories, as the name suggests, neglect such interactions, they cannot retain their mathematical consistency in such cases. Fluctuation effects are important in other contexts as well, and in this study we are essentially investigating the use of the model usually proposed to account for them to simulate symmetry breaking vortex phenomena.

Let us briefly discuss the first method mentioned above, and how it was used to obtain an energy functional for describing the phenomena of superconductor phase transitions. An excellent discussion of Landau theory, formed through a list of postulates compatible with mean field theory (rather than by way of a system Hamiltonian) is provided in [Go92]. In the context of superconductivity, the theory is well summarized in [HT01]. In essence, Landau noticed that mean field theories applied in different contexts

(18)

shared common features, among those the fact that equations of state modeling the systems involved could all be derived from similar expressions for the Gibbs free energy. Thus, it was postulated that the Landau free energy, given by a functional of an order parameter that may also depend on the coupling constants in the system, can be found that is minimized by the state of the system near the critical transition. It can be accessed through the Landau free energy density, which for a homogeneous, isotropic system of finite volume is approximated as f = fn+ a|u|2+ b 2|u| 4 . (1.1)

Here, the order parameter u is zero in the disordered state, but nonzero in the ordered state, the parameters a and b generally depend on the system temperature T and coupling (material) parameters, and fn represents

the free energy of the disordered state. Actually, (1.1) represents a truncated Taylor expansion of the free energy density in powers of |u|2 _{near T}

c; the validity of this approximation is a primary stipulation of

the theory. With this at hand, the central Landau postulate amounts to the assertion that with b > 0, minimization of the free energy functional with respect to u, found itself by perhaps integrating (1.1) over a finite volume, yields a description of the state of the system. Such formalism has long been used to describe the phenomenon of superconductivity near Tc in the absence of an externally applied magnetic field, with

u = 0 (when a > 0 and T > Tc) corresponding to the normal state, and u 6= 0 (when a < 0 and T < Tc)

corresponding to the superconducting state. In this context, u is taken to be complex-valued, with nonzero values indicating the presence of superconducting charge carriers. To allow u to be spatially varying, and include effects due to the presence of a magnetic field h, the coarse graining procedure mentioned above and energy considerations motivated the inclusion of two additional terms in (1.1) by Ginzburg and Landau in [GL50]. Therein, they specified the free energy density as

f = fn+ a|u|2+ b 2|u| 4₊ 1 2m|(−i~∇ − 2eA)u| 2 +µ 2|h| 2_, _(1.2)

where m is twice the free electron mass, ~ is Planck’s constant, e is the electron charge, µ is the magnetic permeability of the material, and A is a magnetic vector potential that gives ∇ × A = µh. If we suppose that a given superconducting sample occupies the region Ω ⊆ R3_{, then}

Z

Ω

f dx, (1.3)

where f is given by (1.2), amounts to a coarse-grained approximation for the Gibbs free energy of the system. When allowing for the presence of an applied magnetic field H, the appropriate Gibbs free energy approximation becomes

Z

Ω

(f − µh · H) dx, (1.4)

and it is the minimizers of (1.4) that are taken for a description of the superconductor state near Tc. At

this point, it is conventional to subtract the constant Z Ω fn− µ 2|H|

(19)

minimizers of E(u, A) = Z Ω f − fn− µh · H + µ 2|H| 2_dx = Z Ω a|u|2+b 2|u| 4₊ 1 2m|(−i~∇ − 2eA)u| 2 +µ 2|h − H| 2_dx. (1.5)

This last step can be regarded as a mathematical convenience, used primarily to obtain favorable natural boundary conditions when applying variational methods. By no means does (1.5) underlie every GL system; even under the stipulations we have made, more complicated models have been proposed and studied. However, it does constitute a widely accepted and extensively researched version of the three dimensional GL model, and we shall henceforward work from it to realize a superconductor model applicable to the physical scenarios considered in this work.

Naturally, a nondimensionalized version of (1.5) is more convenient to work with. To obtain this, we need to first know the units associated with u. To obtain more information regarding the physical interpretation of the order parameter, the overview in [DGP92] provides a good starting point. For our purposes, it is helpful to examine the term

2e2

m |A|

2

|u|2, (1.6)

produced when calculating the the fourth term in the integrand of the right hand side of (1.5). Since (1.6) must constitute an energy density, from basic dimensional analysis we conclude that |u|−2_{has the dimensions}

of volume. In order to keep the temperature dependence isolated to the coefficients of the |u|2_{and |u|}4_terms

in the integrand of the functional, we used a somewhat unconventional set of substitutions, and retained nondimensionalized versions of both a and b. We start with the assumption that H is constant, and oriented in the z-direction as

H = H0z,ˆ (1.7)

where H0∈ (0, ∞). We introduce the substitutions

h = −H0 2 ˜ h, and H = −H0 2 ˜ H, (1.8)

then the integrand in the last term on the right hand side of (1.5) is

µ 2|h − H| 2₌ µH 2 0 8 ˜h − ˜H 2 . (1.9) Now, we identify γ = eµH0 ~ , (1.10) and take a = ~ 2_γ 2ma, b =˜ 2~2_µe2 m2 ˜b, and u = 1 2e r mγ µ u.˜ (1.11)

Using (1.11) in the integrands of the first two terms on the right hand side of (1.5), we get

a|u|2+ b 2|u| 4₌ ~2γ 2m 1 2e r mγ µ 2 ˜ a˜u 2 + 2~ 2_µe2 m2 1 2e r mγ µ 4 _˜_b 2 ˜u 4 = µH 2 0 8 ˜a ˜u 2 + ˜_b 2 ˜u 4 ! . (1.12)

(20)

We recognize that γ−1/2, where γ is as given in (1.10), has the units of length; essentially, this is the length scale we have adopted. We use this, together with the substitution

A = −~ √ γ 2e ˜ A, (1.13)

and for the third integrand on the right hand side of (1.5), get

1 2m|(−i~∇ − 2eA)u| 2 = 1 2m|−i~∇u − 2eAu| 2 = 1 2m −i~ √ γ 1 2e r mγ µ ˜ ∇˜u − 2e −~ √ γ 2e 1 2e r mγ µ ˜ A˜u 2 =µH 2 0 8 (−i ˜∇ + ˜A)˜u 2 . (1.14)

The common factor

µH2 0

8 (1.15)

appearing in (1.9), (1.12) and (1.14) has the units of energy density, as expected. Furthermore, it is straight-forward to verify that the quantities ˜a, ˜b, ˜u, ˜A, ˜h and ˜H, identified through the substitutions (1.8), (1.11) and (1.13), are dimensionless. Moreover, from (1.10) and (1.13), we have

−µH0 2 ˜ h = µh = ∇ × A = (√γ) −~ √ γ 2e ˜ ∇ × ˜A = −µH0 2 ˜ ∇ × ˜A, (1.16)

which shows that

˜

∇ × ˜A = ˜h. (1.17)

After using our length scale to obtain nondimensionalized versions of Ω and the working volume measure, which we denote as ˜Ω and d˜x, respectively, we combine (1.9), (1.12), (1.14), and (1.17) to write (1.5) as

E(u, A) = µH 2 0 8γ3/2 Z ˜ Ω ˜ a|˜u|2+ ˜ b 2|˜u| 4₊ (−i ˜∇ + ˜A)˜u 2 + | ˜∇ × ˜A − ˜H|2d˜x. (1.18) Then, with E = µH 2 0 8γ3/2E,˜ (1.19)

we henceforward restrict our attention to the nondimensionalized functional

˜ E(˜u, Ã) = Z ˜ Ω ˜ a|˜u|2+ ˜_b 2|˜u| 4₊ (−i ˜∇ + Ã)˜u 2 + | ˜∇ × Ã − ˜H|2d˜x. (1.20)

For convenience, the “ ˜ ” notation used to indicate the nondimensionalized versions of the items introduced in this section will be suppressed, unless otherwise stated. In the next section, we introduce some additional approximations that allow us to work with (1.20) on a particular class of rotationally symmetric 2-manifolds.

(21)

1.3 Model Realization on a Surface of Revolution

Henceforward, we use (R, θ, Z) to denote the usual cylindrical polar coordinates in R3_{; that is, R measures}

distance from the z-axis, Z represents the usual Euclidean z-coordinate, and θ ∈ [0, 2π) is the polar angular coordinate such that

x = R cos θ y = R sin θ z = Z,

(1.21)

for each (x, y, z) ∈ R3. Let I ⊆ R be an interval, where

inf I = 0, and sup I = s2∈ (0, ∞). (1.22)

We use s to denote an arbitrary element of I, and let the coordinates R and Z be given through functions on I, which we denote with the corresponding coordinates themselves. That is, we take

R : I → R, s 7→ R(s), (1.23)

and

Z : I → R, s 7→ Z(s). (1.24)

Then, we define S to be the range of the mapping

r : I × [0, 2π) → R3, (s, θ) 7→ r(s, θ) = (R(s) cos θ, R(s) sin θ, Z(s)) . (1.25) We can visualize S geometrically as being the set formed by rotating the range of the mapping s 7→ r(s, θ0),

for any θ0 ∈ [0, 2π), about the z-axis. Suitable conditions on (1.23- 1.24) need to be imposed, in order to

ensure that S given by (1.25) has reasonable mathematical and geometric properties. We begin with the following stipulations :

(S1) R > 0 on Io.

(S2) s 7→ (R(s), Z(s)) is injective, and an element of C I; R2.

(S3) s 7→ (R(s), Z(s)) is differentiable, except perhaps at finitely many points in Io. (S4) Whenever s 7→ (R(s), Z(s)) is differentiable on J ⊆ I, we have

inf s∈J p Rs(s)2+ Zs(s)2> 0, and sup s∈J p Rs(s)2+ Zs(s)2< ∞.

(S5) If R(s0) = 0 with s0∈ ∂I, then Rs(s0) 6= 0.

(S6) If I = [0, s2], then R|∂I = 0, and Z(0) 6= Z(s2).

(S7) If I = (0, s2), then inf

s∈IR(s) > 0, and lims→0(R(s), Z(s)) 6= lims→s2

(R(s), Z(s)). (S8) If I = [0, s2), then lim s→s2 R(s) 6= 0, and lim s→s2 (R(s), Z(s)) = (R(0), Z(0)) ⇐⇒ R(0) 6= 0. (S9) If I = (0, s2], then lim s→0R(s) 6= 0, and lims→0(R(s), Z(s)) = (R(s2), Z(s2)) ⇐⇒ R(s2) 6= 0. (1.26)

In (1.26), Iodenotes the interior of the interval I, and the derivatives of R and Z have been denoted as Rsand

(22)

topological 2-manifold with respect to the subspace topology induced by the containment S ⊆ R3. By this, we mean specifically that S is locally Euclidean, and the induced subspace topology mentioned above satisfies the Hausdorff axiom. In addition, (1.26) ensures that S has favorable mathematical properties near its manifold boundary, ∂S. As we generally expect for a surface of revolution, ∂S itself is either empty, or made up of at most two circular components parallel to the xy-plane; in either case we observe that ∂S ∩S = ∅. For clarification, let us briefly comment on each stipulation in (1.26). We impose (S1) to maintain consistency with usual cylindrical coordinate system conventions, and to retain the locally Euclidean property. In (S2), the first condition is used to retain the locally Euclidean property, while preserving the invertibility of r on Io_{× [0, 2π). The second condition in (S2) ensures that S is both connected and bounded. Conditions}

(S3), (S4) and (S5) amount to minimal smoothness requirements; they allow for manifolds with limited smoothness, but disallow the occurrence of “cusps.” In particular, (S5) disallows the occurrence of cusps at any manifold points that lie in the z-axis. Through (S6), we obtain a description of how a closed version of the parameter interval I may be used. Here, S can only be a simply connected, closed (not topologically) manifold. We avoid introducing a nontrivial intersection of ∂S and S, and retain the locally Euclidean property. When I is open, (S7) applies. Both conditions are in place to avoid exceptions involving ∂S; without both, S itself might appear on either side of some subset of ∂S. Conditions (S8) and (S9) describe the use of a parameter interval that is neither open nor closed. These could be used to parameterize closed manifolds that are not simply connected, but also allow for the representation of a wide class of simply connected manifolds that have a circular boundary. In both cases, the stipulations given are adopted to avoid boundary exceptions, while retaining the locally Euclidean property. Throughout much of the sequel, we assume that strengthened forms of (S3) and (S5) from (1.26) are valid. These, we write as follows:

(S3∗) s 7→ (R(s), Z(s)) ∈ Ck _{I; R}2_{, for all k ∈ N.}

(S5∗) If R(s0) = 0 with s0∈ ∂I, then Rs(s0) 6= 0 and Zs(s0) = 0.

(1.27)

In particular, when (S3∗) and (S5∗) hold, (S4) also holds with J = I, and similar results involving higher derivatives of R and Z can be obtained. Henceforward, we refer to any parameterization of S given by (1.25), such that all of the stipulations in both (1.26) and (1.27) are satisfied, as a standard parameterization for S. In select cases, we shall carry out work on domains S whose parameterization r meets the criteria in (1.26), but not those in (1.27).

Let us suppose now that S is a manifold that is given by some standard parameterization r through (1.25), and return our attention to the GL functional in (1.20). First, we identify the (nondimensionalized) version of the region occupied by the superconductor, Ω, as being that obtained from the parametric representation ρ : (−, ) × I × [0, 2π) → R3_{given by}

ρ(ρ, s, θ) = r(s, θ) + ρ n(s, θ). (1.28)

In (1.28), n(s, θ) represents the surface normal vector given by

n = rs× rθ |rs× rθ|

, (1.29)

and > 0 is chosen sufficiently small to guarantee that each component of ∂Ω is geometrically similar to S, and always lies a distance of from S. We generally regard as being relatively small, so that Ω will then represent a thin plate or film in R3_{, perhaps with curvature, having rotational symmetry about the}

(23)

z-axis. At this stage, we introduce two additional approximations, denoted below as (T1) and (T2), that take advantage of the “thin” structure that the superconducting sample exhibits in this context. Throughout this discussion, we will comment on the mathematical validity of our steps, and indicate where other authors used similar approximations in related work. The first approximation is with regard to the magnetic vector potential. Generally speaking, A is an unknown component of the state of the superconductor; that is, a minimizer of (1.20) comes in the form of a pairing (u, A), made up of both an order parameter coordinate, and a vector potential coordinate. However, we shall henceforward use the approximation

(T1) A ≈ A0= −R ˆθ, (1.30)

where ˆθ denotes the unit coordinate vector that corresponds to the direction of increasing θ. This approxi-mation has been used extensively in the literature to obtain a two dimensional version of the GL model for sufficiently thin superconductors, and we will comment further on its justification later. For the moment, let us simply state that (1.30) leads to a model approximation that is accurate to the leading order of the thickness parameter . Our second, but essentially related, approximation is to write (1.20) instead as

(T2) E (u, A) ≈ 2 Z S a|u|2+b 2|u| 4_{+ |(−i∇ + A)u|}2 + |∇ × A − H|2dx. (1.31)

Here, any spatial variations of (u, A) with respect to the thickness of the sample have been neglected, and the underlying GL functional is now calculated over a domain of codimension 1; that is, over S ⊆ R3_{. In a certain}

sense, the GL phenomenology itself supports the use of (T2). In particular, the coarse graining process used in Landau theory to define the order parameter on a given region assumes (as the term suggests) that spatial variations can be neglected over sufficiently small subregions. Thus, for example, the approximation (1.31) is essentially redundant whenever is near the correlation length corresponding to the current temperature T . Moreover, as we shall soon show, the PDE problem obtained from an application of the variational method to (T2) corresponds to that frequently used by other researchers when investigating the GL model in similar contexts. Let us combine (T1) with (T2) by first using (1.30) to compute

∇ × A0= ∇ × − R ˆθ = −2ˆz = H, (1.32)

which shows that

Z S |∇ × A0− H|2dx = 0. (1.33) Then, we define E (u) = E(u, A0) 2 , (1.34)

and find, when combining the approximations in (1.30) and (1.31) with (1.33) and (1.34), that

E (u) =Z S a|u|2+ b 2|u| 4₊ (i∇ + R ˆθ )u 2 dx. (1.35)

To be clear, E from (1.35) represents the nondimensionalized GL free energy taken over a sample cross section S, whose minimizers we henceforward regard as describing the superconducting state corresponding to a choice of a and b.

(24)

Next, we carry out the variational method to extract a PDE problem that any sufficiently regular min-imizer of (1.35) must satisfy. In order to do this, we need to specify several vector spaces of functions defined on S that are used throughout the remainder of this work. First, we introduce the Lebesgue space of complex-valued, modulus squared-integrable functions on S, denoted as L2

(S; C), and defined by L2_{(S; C) =} v : S → C : Z S |v|2dx < ∞ . (1.36)

The standard inner product in L2

(S; C) is the mapping L2 (S; C) × L2 (S; C) → C given by (u, v) 7→ hu, vi₂= Z S uv dx, (1.37)

where v is the complex conjugate of v. The standard norm on L2_{(S; C), k·k}₂, is that induced by h·, ·i₂, and is written as

v 7→ kvk₂= hv, vi1/2₂ . (1.38)

Although we will not formalize a proof in this work, it follows that L2

(S; C) is a Hilbert space with respect to (1.37-1.38). In many contexts, we will prefer to take the components of functions in Hilbert space as real numbers. For this purpose, it is convenient to use the real inner product h·, ·i, which we specify as

hu, vi = hu, vi₂+ hv, ui₂= 2 < hu, vi₂ (1.39) for (u, v) ∈ L2

(S; C) × L2

(S; C), where < hu, vi2 is the real part of the complex number hu, vi2. It is a

straightforward exercise to establish, for the same such u and v, that

< hu, vi₂= 1

2hu, vi , (1.40)

and

= hu, vi₂=1

2hu, ivi , (1.41)

where = hu, vi₂ is the imaginary part of the complex number hu, vi₂. Thus, the real and imaginary parts of the component of u along v measured with respect to the standard inner product may be written in terms of h·, ·i, according to (1.40) and (1.41), respectively. Moreover, when evaluating the norm of a function in L2

(S; C), we often utilize that induced by h·, ·i, given by

v 7→ kvk = hv, vi1/2. (1.42)

From the equation

kvk = hv, vi1/2=√2 hv, vi1/2₂ =√2 kvk₂, (1.43) we observe that the norms k·k₂ and k·k are equivalent, and conclude therefore that L2_{(S; C) is a (real)} Hilbert space with respect to (1.39) and (1.42). We shall also often work in the Sobolev spaces Hk

(S; C), where k ∈ N0. To define these, we take H0(S; C) = L2(S; C), and for k ∈ N write

Hk_{(S; C) =} u ∈ L2_{(S; C) :} ∂a ∂xa ∂b ∂yb ∂c ∂zcv 2 < ∞ if a, b, c ∈ N0, and a + b + c ≤ k . (1.44)

(25)

With these foundations in place, we proceed by specifying the domain ofE from (1.35) as H1_{(S; C). It is a} straightforward exercise to show thatE is well-defined on this space, and we omit the details of an exhaustive proof here. It should be clear that the translated gradient operator in the third term of the integrand on the right hand side of (1.35) makes H1

(S; C) the natural choice for the domain, and the embedding result in Theorem A.0.22 can be used to integrate the other two terms. The variational approach essentially provides another way to characterize the minima ofE as being critical (equilibrium) points of the same. By definition, these satisfy

E0_{(u; v) = 0 for all v ∈ H}1

(S; C), (1.45)

whereE0(u; v) represents the Gˆateaux derivative ofE at u, in the direction of v. As always, it is worthwhile to keep in mind that every minimizer ofE is also a critical point, but the converse is generally not valid. In the following theorem, we establish thatE is Gˆateaux differentiable.

Theorem 1.3.1. The functionalE : H1

(S; C) → R is Gˆateaux differentiable, and we have E0_{(u; v) =}D_{(i∇ + R ˆ}_{θ )u, (i∇ + R ˆ}_{θ )v}E₊_{au + b|u|}2_{u, v}

(1.46)

for all u, v ∈ H1_{(S; C).} Proof. Let u, v ∈ H1

(S; C). By definition, the Gˆateaux derivative of E at u in the v direction is given by E0_{(u; v) = lim}

τ →0

E (u + τv) − E (u)

τ . (1.47)

First, we let τ ∈ R, and find almost everywhere (a.e.) on S that

|u + τ v|2_{− |u|}2_{= (u + τ v) · (u + τ v) − |u|}2

= τ (uv + uv) + τ2|v|2_, (1.48)

and

|u + τ v|4_{− |u|}4₌_{|u + τ v|}22

− |u|4_{= |u|}2_{+ τ}2_|v|2_{+ τ (uv + uv)}2 − |u|4

= 2τ |u|2(uv + uv) + τ22|u|2|v|2_{+ (uv + uv)}2

+ 2τ3|v|2_{(uv + uv) + τ}4_|v|4_, (1.49) and (i∇ + R ˆθ ) (u + τ v) 2 − (i∇ + R ˆθ )u 2

= τ(i∇ + R ˆθ )u · (i∇ + R ˆθ )v + (i∇ + R ˆθ )u · (i∇ + R ˆθ )v+ τ2 (i∇ + R ˆθ )v 2 . (1.50)

Theorem A.0.22 can now be used with (1.48-1.50) to calculate the contribution of each term in the functional to the derivativeE0(u; v). From (1.48), we gather that

lim τ →0 R S a|u + τ v|2_{dx −}R S a|u|2_dx τ = a hu, vi + a lim τ →0τ Z S |v|2_{dx = a hu, vi .} _(1.51)

(26)

From (1.49), we get lim τ →0 R S b 2|u + τ v| 4_{dx −}R S b 2|u| 4_dx τ = b|u|2_{u, v +} b 2 lim τ →0τ Z S

2|u|2|v|2+ (uv + uv)2dx

+ blim τ →0τ 2 Z S |v|2_{(uv + uv) dx +} b 2 lim τ →0τ 3 Z S |v|4_{dx = b}_|u|2_{u, v .} (1.52) Finally, (1.50) gives lim τ →0 R S (i∇ + R ˆθ ) (u + τ v) 2 dx −R S (i∇ + R ˆθ )u 2 dx τ

=D(i∇ + R ˆθ )u, (i∇ + R ˆθ )vE+lim

τ →0τ Z S (i∇ + R ˆθ )v 2 dx

=D(i∇ + R ˆθ )u, (i∇ + R ˆθ )vE.

(1.53)

When combining (1.51-1.53), we obtain (1.46), and the proof is complete.

If u is a critical point of E that belongs to H2_{(S; C), then it must also solve the equilibrium GL problem} given by

(i∇ + R ˆθ )2u + au + b|u|2u = 0 on S ∂u

∂s = 0 on ∂S (if ∂S 6= ∅).

(1.54)

Whenever ∂S 6= ∅, the boundary condition in the second line above is imposed with the problem, and (1.54) amounts to an elliptic Boundary Value Problem (BVP) on S. If instead ∂S = ∅, then the problem presents only the PDE that appears in the first line. For notational simplicity, we shall henceforward assume that the condition on ∂S stated in any GL BVP is omitted whenever ∂S = ∅. The above assertion, along with its converse, is formalized in the next theorem. Proving it requires several results that are formally established in Section 2.1.1 of Chapter 2. By h·, ·i_∂S, we mean the real-valued inner product in L2

(∂S; C), defined whenever ∂S 6= ∅ through the mapping

(u, v) 7→ hu, vi_∂S = Z

∂S

(uv + uv) dx. (1.55)

If ∂S = ∅, any such boundary inner product is assumed to equal zero wherever it appears. Several other notations and intermediate results are used for the first time in the next proof, the reader is referred to Section 1.4 for the necessary definitions.

Theorem 1.3.2. Let u ∈ H2

(S; C). Then, we have

E0_{(u; v) = 0 for all v ∈ H}1

(S; C) (1.56)

(27)

Proof. First, let us suppose that (1.56) holds. Let v ∈ L2_{(S; C). We know that the space C}0∞(S; C) is dense

in L2

(S; C) (see Section 1.4), so we may choose a sequence {vk} from it that is convergent to v. Using (1.46)

from Theorem 1.3.1, we have

D

(i∇ + R ˆθ )u, (i∇ + R ˆθ )vk

E

+au + b|u|2u, vk = 0 (1.57)

for all k ∈ N. To (1.57), we apply the Green’s result (2.55), and find that D (i∇ + R ˆθ )2u + au + b|u|2u, vk E +D(i∇ + R ˆθ )u · ˆs, ivk E ∂S = 0 (1.58)

for all such k. If ∂S 6= ∅, we have vk|∂S = 0 for all k ∈ N, and the second term on the right hand side of

(1.58) vanishes. In any case, (1.58) therefore gives

D

(i∇ + R ˆθ )2u + au + b|u|2u, vk

E

= 0 (1.59)

for all k ∈ N, and when we let k → ∞, it follows from (1.59) that D

(i∇ + R ˆθ )2u + au + b|u|2u, vE= 0. (1.60) Since v was chosen arbitrarily in L2_{(S; C), (1.60) implies that}

(i∇ + R ˆθ )2u + au + b|u|2u = 0 on S. (1.61) Next, let us suppose that ∂S 6= ∅, and let v ∈ L2_{(∂S; C). Here, we make use of the fact that the space} C2

(∂S; C) is dense in L2

(∂S; C). Let us pick a sequence {vk} from this space convergent to v, then apply

the Trace theorem, which is included as Theorem A.0.7, to extend vk to the space H1(S; C), for each k ∈ N.

Then, just as in the above argument, we can obtain (1.58) and (1.59), and combine them to find that

D

(i∇ + R ˆθ )u · ˆs, ivk

E

∂S

= 0. (1.62)

Letting k → ∞ in (1.62) leads us to conclude that

(i∇ + R ˆθ )u · ˆs = 0 on ∂S, (1.63)

then the formulation (2.38), combined with the standard parameterization property (S4) from (1.26), allows that

∂u

∂s = 0 on ∂S. (1.64)

For the converse, let us suppose that u satisfies (1.54), and let v ∈ H1

(S; C). Here, from similar considera-tions we find that

D

(i∇ + R ˆθ )2u + au + b|u|2u, vE+D(i∇ + R ˆθ )u · ˆs, ivE

∂S

= 0, (1.65)

then (2.55) can be applied to obtainE0_{(u; v) = 0, which shows (1.56). This completes the proof.}

With our variational foundations now in place, we next introduce time dependence into our model by applying the phenomenology of linear response theory. As the term suggests, this amounts to the assumption

(28)

that the partial derivative with respect to time of some unknown is proportional to the energy derivative at that unknown, over a given time interval of interest. In the context of our study, we regard the order parameter u now as a function of time, whose values lie in the domain ofE . More specifically, we fix T > 0, then define the order parameter u on the domain [0, T ] through the mapping

u : [0, T ] → H1_{(S; C) , t 7→ u(t, ·) ∈ H}1_{(S; C) .} (1.66) Naturally, we need to specify an initial state for the order parameter from the domain of our GL functional; this we denote as u0. To preserve generality in this approach, we take the aforementioned derivative to be

the Gˆateaux derivative in (1.46). Then, the linear response theory leads us to consider the problem of finding u ∈ L2 _{0, T ; H}1

(S; C) such that D∂u

∂t, v E

+E0(u; v) = 0 a.e. on (0, T ], for all v ∈ H1_{(S; C), and} u(0, ·) = u0, (1.67) where L2 0, T ; H1_{(S; C)} = ( v : [0, T ] → H1_{(S; C) :} Z T 0 kv(t, ·)k2+ k∇v(t, ·)k2 dt < ∞ ) . (1.68)

Roughly speaking, (1.67) says that over the time interval of interest, any component of the time rate of change of the order parameter should be equal and opposite to the current energy change in the direction of that component. Now, if we instead seek a solution in the space L2 0, T ; H_∂2_{(S; C)}, we obtain from (1.67) the linear response model that results when instead using the Fr´echet derivative DE of the restriction of E to the same space as the underlying energy derivative. That is, when instead taking ∂u/∂t as being (negatively) proportional to the mapping

DE : H_∂2_{(S; C) → L}2_{(S; C) ,} (1.69)

that is given by

DE (u) = (i∇ + R ˆθ)2u + au + b|u|2u. (1.70) What results corresponds to the problem of finding u ∈ L2 _{0, T ; H}2

(S; C) such that ∂u ∂t + (i∇ + R ˆθ ) 2_{u + au + b|u|}2_{u = 0} _{on (0, T ] × S} u = u0 on {0} × S ∂u ∂s = 0 on [0, T ] × ∂S, (1.71)

where the last condition is omitted if ∂S = ∅. In Chapter 3, we will establish that (1.71) has exactly one solution, which must also satisfy (1.67), as another application of the Green’s formula in (2.55) from Theorem 2.1.1 shows. Since we have introduced another independent variable into our model, we need to explain how it is scaled to obtain the nondimensionalized version shown in (1.71). When recalling, for a moment, the scaling notations and results in Section 1.2, we recognize that the dimensionalized version of

(29)

the PDE in (1.71) must actually read as 1 2e r mγ µ ∂ ˜u ∂t + Γ~ 2_γ 4me r mγ µ i ˜∇ + ˜A0 2 ˜ u + ˜a˜u + ˜b|˜u|2u˜ = 0, (1.72)

where t has the units of time, and Γ > 0 is the linear response constant, which is measured per unit energy, per unit time. We then see that the proper time scaling is realized through the substitution

t = _2m Γ~2_γ ˜ t. (1.73)

Indeed, it is straightforward to verify that ˜t, thus specified, is a dimensionless quantity. Also, we need to know how the newly introduced initial state of the system, u0, can be appropriately chosen in the context

of our model. We will assume that the system begins in either the perfectly nonsuperconducting state, represented as u0= 0 on S, or in the perfectly superconducting state. The former corresponds to the pairing

u = 0, h = H that minimizes (1.31) whenever a > 0 (and T > Tc), and describes the normal phase of the

material under full penetration of the external magnetic field. The latter state, also called the Meissner state, is given for a < 0 (and T < Tc) by

u0=p−a/b eiφ on S, (1.74)

where φ ∈ [0, 2π) is arbitrary. Provided a < 0, the pairing of u = u0 and h = 0, where u0 is given by

(1.74), minimizes (1.31), and thus describes the ideal superconducting state characterized by the complete expulsion of the external magnetic field from the sample.

At this stage, we are motivated to assess the GL problems we have obtained so far by comparing them with those used by other researchers in similar contexts. Perhaps most recently, in [DJ04a] and [DJ04b], the problems (1.54) and (1.71) were used to simulate superconductivity in sufficiently thin superconducting hol-low spheres immersed in a constant magnetic field. Therein, the authors took S =x ∈ R3: |x| = constant , and enforced no boundary conditions on the order parameter. The equilibrium model used in [DJ04b] ex-actly corresponds to (1.54) for such S, up to the nondimensionalization scheme. The time dependent model used in [DJ04a] is similarly equivalent to (1.71), when one assumes that no external electrical voltages or currents are being applied. Essentially, the GL models we work from represent generalizations of those used in these investigations, where the superconductor sample now may be taken as occupying any sufficiently thin, rotationally symmetric region. The use of two dimensional GL models to simulate superconductivity has, of course, a more extensive history, and a variety of approaches have been applied. The often cited work in [CDG96] demonstrates the validity of leading order approximations involving externally applied fields in thin superconducting films. In [BP02], the authors used a two dimensional approximation without fixing the vector potential A to study equilibrium vortex states on a variety of thin, symmetric geometries. In [HT01], we find comprehensive information on current perspectives regarding the physical and mathematical validity of two dimensional GL models. In any case, we are confident that (1.71) represents a good approximation for the dynamics of superconductivity whenever the physical and mathematical stipulations we have described are reasonably valid.

Next, we need to make some mathematical observations that will later be required to support certain aspects of our computational approach. First, we need to properly formulate gauge invariance, as it applies