Automatic exploration techniques for the numerical bifurcation study of the Ginzburg-Landau equation

TL;DR

This paper introduces a robust numerical continuation method for exploring the solution landscape of the Ginzburg-Landau equations, using bifurcation analysis to understand superconducting states.

Contribution

It develops an automatic exploration algorithm combining Lyapunov-Schmidt reduction and equivariant branching, implemented in Python for two-dimensional superconductivity models.

Findings

Successfully constructed complete solution landscapes for multiple examples.

Demonstrated robustness and effectiveness of the algorithm.

Showed applicability to a class of grids using symmetry-based methods.

Abstract

This paper considers the extreme type-II Ginzburg-Landau equations, a nonlinear PDE model that describes the states of a wide range of superconductors. For two-dimensional grids, a robust method is developed that performs a numerical continuation of the equations, automatically exploring the whole solution landscape. The strength of the applied magnetic field is used as the bifurcation parameter. Our branch switching algorithm is based on Lyapunov-Schmidt reduction, but we will show that for an important class of grids an alternative method based on the equivariant branching lemma can be applied as well. The complete algorithm has been implemented in Python and tested for multiple examples. For each example a complete solution landscape was constructed, showing the robustness of the algorithm.

Tables

Table 1. Table 1 : Approximate values for the applied magnetic fields strength ( μ 𝜇 \mu ) for the branch points displayed in figure 9 .

index	$\mu$ value	index	$\mu$ value	index	$\mu$ value	index	$\mu$ value
1	0.6989	16	1.51736	31	1.02539	46	1.34396
2	0.7319	17	1.6149	32	1.4412	47	1.3655
3	1.13777	18	1.655	33	0.58275	48	1.425396
4	1.1423	19	0.09241	34	0.551370	49	1.9115
5	1.7176	20	0.54002	35	0.80322	50	1.425179
6	1.7437	21	0.646098	36	1.0406	51	1.7430
7	1.7959	22	0.61367	37	1.3712	52	1.7063
8	0.12905	23	1.281974	38	1.41573	53	1.75378
9	0.2591	24	1.27688	39	0.5566	54	0.7922136
10	0.34209	25	1.42875	40	0.8639	55	1.10976
11	0.7594	26	1.7453	41	0.6980	56	0.761674
12	0.8029	27	1.7055	42	0.6308	57	1.312923
13	0.8635	28	1.5824	43	1.13779	58	1.517199
14	0.875147	29	0.25285	44	1.02538	59	1.557099
15	1.5048	30	1.0240	45	1.28195	60	1.67462

Equations

$$\displaystyle\begin{aligned} G(\psi,\mathbf{A})-G_{n}=\xi\frac{|\alpha|^{2}}{\beta}\int_{\Omega}\bigg{[}&-|\psi|^{2}+\frac{1}{2}|\psi|^{4}+||-\i\nabla\psi-\mathbf{A}\psi||^{2}\\ &+\kappa^{2}(\nabla\times\mathbf{A})^{2}-2\kappa^{2}(\nabla\times\mathbf{A})\cdot\mathbf{H}_{0}\bigg{]}d\Omega.\end{aligned}$$

$$\{\begin{array}{ll}\nabla \times (\nabla \times \mathbf{A})=0& \text{in }\Omega ,\\ \mathbf{n}\times (\nabla \times \mathbf{A})=\mathbf{n}\times {\mathbf{H}}_0& \text{on }\partial \Omega .\end{array}$$

$$\begin{array}{rl}& \mathcal{G}\mathcal{L}:X\times \mathbb{R}\to Y,\\ & 0=\mathcal{G}\mathcal{L}(\psi ,\mu )=\{\begin{array}{ll}(-\text{i}\nabla -\mathbf{A}(\mu ){)}^2\psi -\psi (1-|\psi {|}^2)& \text{in }\Omega ,\\ \mathbf{n}\cdot (-\text{i}\nabla -\mathbf{A}(\mu ))\psi & \text{on }\partial \Omega .\end{array}\end{array}$$

$${\tau }_{\omega }:{\mathbb{R}}^2\to {\mathbb{R}}^2:\left(\begin{array}{c}x\\ y\end{array}\right)\to \left(\begin{array}{cc}\cos (\omega )& -\sin (\omega )\\ \sin (\omega )& \cos (\omega )\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right),$$

$$\sigma :{\mathbb{R}}^2\to {\mathbb{R}}^2:\left(\begin{array}{c}x\\ y\end{array}\right)\to \left(\begin{array}{c}-x\\ y\end{array}\right),$$

$${\theta }_{\eta }:X\to X:\psi \to e^{\text{i}\eta }\psi (\forall \eta \in [0,2\pi ]),$$

$${\tau }_{\omega }:X\to X:\psi (x,y)\to \psi ({\tau }_{\omega }(x,y)),$$

$$\sigma :X\to X:\psi (x,y)\to \overline{\psi (-x,y)}.$$

$$\Gamma =S^1\times D_{\infty }.$$

$${\mathcal{J}}_{\psi }(\psi ,\mu ):X\to Y:\phi \to \left((-\text{i}\nabla -\mathbf{A}{)}^2-1+2|\psi {|}^2\right)\phi +{\psi }^2\overline{\phi }.$$

$$⟨\cdot ,\cdot {⟩}_{\mathbb{R}}=\mathfrak{Re}⟨\cdot ,\cdot {⟩}_{L_{\mathbb{C}}^2(\Omega )},$$

$$\mathcal{G}\mathcal{L}({\psi }_s,{\mu }_s)=0\Rightarrow {\mathcal{J}}_{\psi }({\psi }_s,{\mu }_s)(\text{i}{\psi }_s)=0,$$

$$\mathcal{J}(\psi ,\mu ):X\times \mathbb{R}\to Y:\left(\begin{array}{c}\phi \\ \kappa \end{array}\right)\to {\mathcal{J}}_{\psi }(\psi ,\mu )\phi +\kappa {\mathcal{J}}_{\mu }(\psi ,\mu ).$$

$${\mathcal{J}}_{\psi }(\psi +\delta \psi ,\mu )\phi -{\mathcal{J}}_{\psi }(\psi ,\mu )\phi =$$

$$-\left((-i\nabla -\mathbf{A}{)}^2-1+2|\psi {|}^2\right)\phi -{\psi }^2\overline{\phi }$$

$$=$$

$${\mathcal{H}}_{\psi \psi }(\psi ,\mu ):X\times X\to Y:({\phi }_1,{\phi }_2)\to 2\left(\overline{\psi }{\phi }_1{\phi }_2+\psi \overline{{\phi }_1}{\phi }_2+\psi {\phi }_1\overline{{\phi }_2}\right).$$

$$\mathcal{H}(\psi ,\mu ):(X\times \mathbb{R})\times (X\times \mathbb{R})\to Y:$$

$$\left(\left(\begin{array}{c}{\phi }_1\\ {\kappa }_1\end{array}\right)\left(\begin{array}{c}{\phi }_2\\ {\kappa }_2\end{array}\right)\right)\to {\mathcal{H}}_{\psi \psi }(\psi ,\mu ){\phi }_1{\phi }_2+{\kappa }_1{\mathcal{H}}_{\psi \mu }(\psi ,\mu ){\phi }_2+{\kappa }_2{\mathcal{H}}_{\psi \mu }(\psi ,\mu ){\phi }_1+{\kappa }_1{\kappa }_2{\mathcal{H}}_{\mu \mu }(\psi ,\mu )$$

$$\frac{{\partial }^3\mathcal{G}\mathcal{L}}{\partial {\psi }^3}(\psi ,\mu ):$$

$$({\phi }_1,{\phi }_2,{\phi }_3)\to 2\left(\overline{{\phi }_1}{\phi }_2{\phi }_3+{\phi }_1\overline{{\phi }_2}{\phi }_3+{\phi }_1{\phi }_2\overline{{\phi }_3}\right).$$

$$\forall k\ge 4:\frac{{\partial }^k\mathcal{G}\mathcal{L}}{\partial {\psi }^k}(\psi ,\mu )=0,$$

$$\frac{{\partial }^k\mathcal{G}\mathcal{L}}{\partial {\mu }^k}(\psi ,\mu ),$$

$${\psi }^{(h)}=\left(\begin{array}{c}{\psi }_1^{(h)}\\ ⋮\\ {\psi }_n^{(h)}\end{array}\right)=\left(\begin{array}{c}\psi (x_1,y_1)\\ ⋮\\ \psi (x_n,y_n)\end{array}\right).$$

$$\forall {\psi }^{(h)},{\varphi }^{(h)}\in {\mathbb{C}}^n:\sum _{i=1}^n|V_i|\overline{{\varphi }_i^{(h)}}{\left(K^{(h)}(\mu ){\psi }^{(h)}\right)}_i=\sum _{i=1}^m\sum _{\text{edges }{e}_{j,k}\text{ of }{T}_i}$$

$$\displaystyle+\big{(}\psi_{k}^{(h)}-\overline{U_{j,k}(\mu)}\psi_{j}^{(h)}\big{)}\overline{\phi_{k}^{(h)}}\Big{)}.$$

$$\alpha_{j,k}^{(i)}=\frac{1}{2}\frac{t_{j,k}}{\sqrt{1-t_{j,k}^{2}}}\qquad\text{with }t_{j,k}=\bigg{\langle}\frac{e_{k,l}}{||e_{k,l}||_{2}},\frac{e_{l,j}}{||e_{l,j}||_{2}}\bigg{\rangle}_{2}.$$

$$U_{j,k}(\mu )=\exp \left(-\text{i}{\int }_{(x_k,y_k)}^{(x_j,y_j)}⟨e_{j,k},\mathbf{A}(\mu ,\omega ){⟩}_{2}d\omega \right)$$

$$\begin{array}{rl}& \mathcal{F}:{\mathbb{C}}^n\times \mathbb{R}\to {\mathbb{C}}^n,\\ & 0=\mathcal{F}({\psi }^{(h)},\mu )=\{\begin{array}{l}{\left(K^{(h)}(\mu ){\psi }^{(h)}\right)}_1-{\psi }_1^{(h)}\left(1-|{\psi }_1^{(h)}{|}^2\right)=0,\\ ⋮\\ {\left(K^{(h)}(\mu ){\psi }^{(h)}\right)}_n-{\psi }_n^{(h)}\left(1-|{\psi }_n^{(h)}{|}^2\right)=0.\end{array}\end{array}$$

$$⟨\cdot ,\cdot {⟩}_{\mathbb{R}}^{(h)}=\mathfrak{Re}⟨\cdot ,\cdot {⟩}_2$$

$${\tilde{\psi }}^{(1)}={\psi }^{(0)}+\Delta s{\dot{\psi }}^{(0)},$$

$${\tilde{\mu }}^{(1)}={\mu }^{(0)}+\Delta s{\dot{\mu }}^{(0)},$$

$$\{\begin{array}{ll}& \mathcal{G}\mathcal{L}({\psi }^{(1)},{\mu }^{(1)})=0,\\ & P({\psi }^{(1)},{\mu }^{(1)})=⟨{\dot{\psi }}^{(0)},{\psi }^{(1)}-{\tilde{\psi }}^{(1)}{⟩}_{\mathbb{R}}+⟨{\dot{\mu }}^{(0)},{\mu }^{(1)}-{\tilde{\mu }}^{(1)}{⟩}_{\mathbb{R}}=0\end{array}$$

$$\left(\begin{array}{cc}{\mathcal{J}}_{\psi }({\psi }^{(1)},{\mu }^{(1)})& {\mathcal{J}}_{\mu }({\psi }^{(1)},{\mu }^{(1)})\\ ⟨{\dot{\psi }}^{(0)},\cdot {⟩}_{\mathbb{R}}& ⟨{\dot{\mu }}^{(0)},\cdot {⟩}_{\mathbb{R}}\end{array}\right)\left(\begin{array}{c}\Delta \psi \\ \Delta \mu \end{array}\right)=\left(\begin{array}{c}-\mathcal{G}\mathcal{L}({\psi }^{(1)},{\mu }^{(1)})\\ -P({\psi }^{(1)},{\mu }^{(1)})\end{array}\right).$$

$${\mathcal{J}}_{\psi }(\psi ,\mu )x=y\text{ with }x\in X,y\in Y.$$

Full text

{AMS}

37M20, 37N20

1 Introduction

Superconductors are materials that, below a certain characteristic temperature ($T_{c}$), exhibit a complete loss of electrical resistivity [12] and, as a result, can generate currents that expel applied magnetic fields.

In certain types, such as bulk type II superconductors, the material is not homogenously superconducting but there are vortices through which magnetic fields can penetrate the material. These vortices organise themselves in regular patterns known as an Abrikosov lattice [1].

For small nano devices the emerging patterns of vortices depend in an intricate way on the system parameters and the geometry of the sample. Scientists and engineers are designing devices that have an improved critical field and temperatures by exploiting geometrical properties and engineering the material parameters [5]. It is important to understand the effect of the geometry on the vortex dynamics [13]. The dynamics of the superconducting materials is modelled by the Ginzburg-Landau equation, a non-linear Schrödinger equation, and the emerging patterns are steady states of this equation. Transitions between different patterns are marked by bifurcation points.

In this paper we aim to develop a numerical tool that generates automatically all connected patterns that appear in a given superconductor geometry and indentifies all the bifurcation points that define the transitions between patterns as parameters of the system change. There is wide interest from material scientists and device engineers to understand what parameters are determining the stability of superconducting states. In particular, the aim is to understand the dynamics of the transitions between these patterns and what needs to be changed in the geometry or parameter settings to prevent the system from making a transition that destroys a given pattern.

The challenge is that sparse linear algebra is required to solve the large systems that describe the Ginzburg-Landau equation. Existing tools such as Auto [9] and Matcont [8] can automatically generate bifurcation diagrams for small systems of coupled ordinary differential equations. However, these tools are based on dense linear algebra and they cannot scale to the large sparse systems that appear in the description of superconductors. In this paper we develop the required tools that can solve the large systems of equations solely based on sparse linear algebra.

Describing the state of a superconductor

The state of superconductors is described by the Ginzburg-Landau partial differential equation. We will consider samples of superconducting material that occupy an open, bounded region $\Omega$ of the two-dimensional Euclidean space, subject to an external magnetic field $\mathbf{H}_{0}$ (see Figure 1).

There are multiple types of superconductors. They are classified as type-I and type-II. A type-I superconducting material can only be in the normal or the superconducting state. In the normal state, or homogeneously non-superconducting, the material behaves like a normal conductor, i.e. the total magnetic field $\mathbf{B}$ entirely penetrates the sample (Figure 1(a)). The electrical current has a resistance. This state occurs for high temperatures ($T>T_{c}$) or for magnetic field strengths above a certain critical value $\mu_{c}$. In contrast, in the homogeneously superconducting state, the total magnetic field $\mathbf{B}$ is entirely expelled from the interior of the sample and the material exhibits zero electrical resistance (Figure 1(b)). This state occurs for low temperatures ($T<T_{c}$) provided that the external magnetic field’s strength is sufficiently low ($\mu<\mu_{c}$).

In type-II superconductors, an additional third state occurs for temperatures and magnetic field strengths between two critical values ($T_{c_{1}}<T<T_{c_{2}}$ and $\mu_{c_{1}}<\mu<\mu_{c_{2}}$). In this case, the magnetic field $\mathbf{B}$ only locally penetrates the material. This gives rise to vortices of non-superconducting material around which superconducting currents appear in the material (Figure 1(c)). This extra state that only occurs for type-II superconductors is called the mixed state [21]. It is the vortex patterns found in these type-II materials that are the focus of the current paper.

The state of a superconductor is described by two quantities: the total magnetic field $\mathbf{B}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ and the density of Cooper pairs $\rho_{C}:\Omega\cup\partial\Omega\rightarrow\mathbb{R}$. Cooper pairs are pairs of electrons that constitute superconductivity: a high density corresponds with the superconducting state. Both quantities are determined by the Ginzburg-Landau system [12], which includes a partial differential equation for an order parameter $\psi:\Omega\cup\partial\Omega\rightarrow\mathbb{C}$ with $|\psi|^{2}=\rho_{C}$. In the current paper we will only consider extreme type-II superconductors, for which the Ginzburg-Landau system decouples in a system for $\psi$ and the magnetic field is independent of $\mathbf{A}$.

Motivation for the current paper

The understanding how vortex patterns are formed and how they lose their stability is important for multiple applications of extreme type-II superconductors. This is theoretically investigated by numerically solving the steady states of the Ginzburg-Landau system. After discretization it results in a large sparse non-linear system. An efficient numerical solver for this system is described in [33], based on an AMG-preconditioned Newton-Krylov method. Unlike other popular methods (e.g. time-stepping the Ginzburg-Landau system via Gauss-Seidel iterations [2]), this solver can find both physically stable and unstable solutions. This solver is only one piece of the machinery that is required to automatically generate bifurcation diagrams.

Though the direct use of unstable steady states is limited in real situations, they still contain important information on transitions between the stable patterns. This is further motivated by example 1.1 below.

Numerical evidence in [33] shows an independency between the number of unknowns and the number of Krylov iterations required to converge, demonstrating the optimality of the presented solver. To achieve this, the solver makes full use of the sparsity properties of the system.

Example 1.1**.**

We consider a square material of side length $3$ subject to a homogeneous magnetic field with strength $\mu=1.3$. For these values, the material is in the mixed state, but can adopt two different stable patterns, shown in figure 2 [31]. Both patterns 1 and 2 are stable, their respective energies are given by $-0.23387$ and $-0.13271$. We keep the value for the magnetic field strength fixed and want to perform a transition from pattern 1 to pattern 2 through a temporary external perturbation, by supplying additional energy.

*It is, however, not sufficient to add the difference in energy between patterns 1 and 2 to the system. Instead, the minimal amount of energy that needs to be supplied is given by the energy difference between pattern 1 and the solution presented by pattern 3 (figure 2(c)). This last pattern represents an unstable solution for the same magnetic field strength, with an energy of $-0.14007$. To perform a transition between patterns 1 and 2, we first need to supply energy to the system to reach the barrier implied by pattern 3. Afterwards energy is released to reach the steady state described by pattern 2. This is further indicated in figure 3. *

To investigate how vortex patterns change when certain parameters (e.g. the strength of the applied magnetic field) are altered, a numerical continuation of the Ginzburg-Landau system is performed. This has already been the topic of some papers as well. In [31] two different sized square samples subject to a homogeneous external magnetic field are concerned and their solution landscapes are explored. In [32], square samples in the field of a magnetic disc are studied. Such bifurcation diagrams are helpful when facing problems like the one described in example 1.1 as well. This is indicated by example 1.2.

Example 1.2**.**

*We again consider the set-up of example 1.1. Figure 4 was recreated from [31] and shows the bifurcation diagram of the problem. At parameter value $\mu=1.3$ we have multiple solutions: two stable ones that lie on curves A (pattern 1 of figure 2) and curve D (pattern 2), and two unstable ones on curves B and C (pattern 3). In order to perform a transition from the solution on A to the one on D, we need to add enough energy to bypass the barrier induced by the unstable solution on curve C - the unstable solution with the lowest energy that lies on a curve connecting A and D. *

Though the solution landscapes provided by [31] and [32] already contain a lot of information, they were not automatically explored and are possibly incomplete. To fully investigate the possible steady states and their transitions, it is essential to use automatic exploration techniques in order to receive a complete solution landscape.

Besides the tools already mentioned earlier (such as AUTO [9] and MatCont [8]) there is LOCA [30] that is developed around sparse linear algebra, but is less easy to use and requires expert knowledge in HPC hardware. Furthermore it does not include a branch switching functionality.

The goal of the current paper is to extend the analysis done in [31] and [32] by providing a robust method that does not only perform a numerical continuation, but automatically explores the solution landscape as well. We will use the Python package PyNCT [11] as a basis for the implementation, and extend it with automatic exploration techniques. The solver described in [33] is an important building block for this purpose, as it fully exploits the sparsity of the equations. To the best of our knowledge, an automatic exploration of the solution landscape of the Ginzburg-Landau problem has never been provided before.

To emphasize the details of the Ginzburg-Landau problem, in the current paper we will specifically derive the methods required for automatic exploration using this problem. The algorithms are however easily extended for general problems, and a general version was implemented in PyNCT as well.

Automatic exploration techniques

Automatic exploration of the bifurcation diagram is realized in two steps. First, the bifurcation points that are encountered during the continuation of the patterns are identified. In these points the Jacobian of the system has one or multiple zero eigenvalues. A method suggested in [27] uses this property to identifiy the bifurcations; after each numerical continuation step the lowest magnitude Ritz values of the Jacobian are monitored to indicate a nearby bifurcation point. If a nearby bifurcation point is indicated an extended system of non-linear equations is solved that identifies the solution and the parameter value of the bifurcation point, simultaniously. It is solved by a Newton-Krylov solver.

The second step of this automatic exploration constructs for each bifurcation point, identified in the first step, the tangent directions of the solution curves emerging from these points. For this purpose we propose two different methods in the current paper.

The first method is based on a Lyapunov-Schmidt reduction [27]. Here, the problem is reduced to a low dimensional system of algebraic equations, that contain all of the information about the bifurcation’s behaviour and are easier to analyse. These systems define the tangent directions to the intersecting solution branches. We propose an iterative algorithm for the determination of these directions that is inpired by the analysis done in [28]. The derived algorithm can always be applied to the bifurcations that are encountered in the Ginzburg-Landau equation, but requires multiple solves of a linear system. The amount of these solves is linked to the symmetry group of the problem, and this considerably reduces the speed of the algorithm when highly symmetrical materials (e.g. a decagon) are considered.

We also propose a second method that is based on the equivariant branching lemma. Here, the lemma links the symmetry group of the problem to predict the symmetries of the branches that emerge at bifurcation points [20, 24]. In [31] the equivariant branching lemma was used for the prediction of the symmetry groups of emerging branches, but it was not used for the calculation of the tangent directions themselves. This alternative to the first approach only requires a single linear system solve and is, hence, a good alternative. However, its use is limited to the class of samples of materials with dihedral symmetry.

These two methods allow us to determine the tangent directions of the solution curves that emerge from a given branch point. For each new tangent direction, a new numerical continuation is started, and the processes of calculating branch points and constructing tangent directions are repeated. Eventually no new (connected) branch points can be found and the algorithm stops, providing us with a complete solution diagram of multiple interconnected solution curves. The Python implementation of the full algorithm was tested for multiple examples of two-dimensional superconductors, each with different symmetries. For each example a complete solution landscape was found (in the sense of interconnecting curves, the ones not in any way connected to the starting point cannot be found), emphasizing the robustness of the algorithm.

Outline

The remainder of the article is organized as follows. Section 2 reviews the Ginzburg-Landau system for extreme-type-II superconductors and its symmetries. The properties of its Jacobian and other derivatives are also discussed. Details of numerical continuation can be found in section 3: the solver described by [33] is shortly discussed and its role in our method is indicated, as well as some details of challenges we encountered for the numerical continuation itself. The main contribution of the article is contained in sections 4 and 5: while section 4 discusses the detection and determination of bifurcation points, section 5 concerns the construction of the tangent directions to new solution branches. Numerical computations of solution landscapes for multiple two-dimensional samples are included in section 6. For these examples the strength of the applied magnetic field is used as a bifurcation parameter. Section 7 concludes with a discussion of the obtained results. Appendix A contains the proof of an algorithm provided in section 5, appendix B contains some details on how to calculate certain terms in this algorithm.

Notations

For $z\in\mathbb{C}$, we denote $\mathfrak{Re}(z)$ and $\mathfrak{Im}(z)$ as respectively its real and imaginary part. $\bar{z}$ is used for complex conjugation. Similarly, for complex valued functions $\varphi:\Omega\rightarrow\mathbb{C}$ we define $\bar{\varphi}:\Omega\rightarrow\mathbb{C}$ such that $\forall x\in\Omega:\bar{\varphi}(x)=\overline{\varphi(x)}$. Adjoints of linear operators $\mathcal{A}$ are denoted by $\mathcal{A}^{*}$. For the symmetry groups under consideration, the symbol $S^{1}$ is used to denote the circle group $\{z\in\mathbb{C}:|z|=1\}$, $D_{m}$ ($m\geq 2$) for dihedral groups and $C_{m}$ ($m\geq 2$) for cyclic groups.

2 The Ginzburg-Landau system and its properties

Description of the system

We will only consider the Ginzburg-Landau problem for extreme type-II superconductors. In [31] it is shown that the Ginzburg-Landau equations decouple for this case, simplifying the problem. A short overview of the resulting equation is given below, based on the analysis done in [31].

For an open, bounded domain $\Omega\subset\mathbb{R}^{3}$, with a piecewise smooth boundary $\partial\Omega$, the Ginzburg-Landau problem is derived by minimizing the Gibbs free energy functional [12]

[TABLE]

The state $(\psi,\mathbf{A})$ is in the natural energy space such that the integral is well-defined. $\psi$ represents a scalar-valued function and is commonly referred to as the order parameter, and the magnetic vector potential corresponding to the total magnetic field is given by $\mathbf{A}$. The physical observables associated with the state $(\psi,\mathbf{A})$ are the density $\rho_{C}=|\psi|^{2}$ of the Cooper pairs and the total magnetic field $\mathbf{B}=\nabla\times\mathbf{A}$. The constant $G_{n}$ represents the energy associated with the normal (non-superconducting) state. The energy (1) depends upon the impinging magnetic field $\mathbf{H}_{0}$ and the material parameters $\alpha,\beta,\lambda,\xi\in\mathbb{R}$. It is presented in its dimensionless form, where the domain $\Omega$ is scaled in units of the coherence length $\xi$. The ratio $\kappa=\lambda/\xi$ of the penetration depth $\lambda$ and the coherence length $\xi$ completely determines the type of the superconductor: for $\kappa<1/\sqrt{2}$ it is said to be of type I, otherwise the superconductor is said to be of type II.

Using standard calculus of variations, minimization of the Gibbs free energy functional gives rise to the Ginzburg-Landau equations [12]: a boundary-value problem in the unknowns $\psi$ and $\mathbf{A}$. For extreme type-II superconductors the limit $\kappa\rightarrow\infty$ is considered, in this case the Ginzburg-Landau problem decouples for $\psi$ and $\mathbf{A}$ [31]. The magnetic vector potential $\mathbf{A}$ is completely determined by the applied magnetic field $\mathbf{H}_{0}$ through the system

[TABLE]

With the magnetic vector potential determined, the order parameter $\psi$ is derived by solving the equation

[TABLE]

The space $X$ corresponds to the natural energy space over $\Omega$ associated with the Gibbs energy (1) and $Y$ to its dual space. The dependence on the strength of the applied magnetic field ($\mu$) is explicitly given since it will be used as a bifurcation parameter. For the remainder of the paper we will write $\mathbf{A}$ instead of $\mathbf{A}(\mu)$, dropping the explicit dependency on $\mu$.

In the present paper, we will consider two-dimensional grids $\Omega\subset\mathbb{R}^{2}$ and numerically solve (LABEL:Ginzburg-Landau_vergelijking) for the order parameter. A bifurcation analysis will be carried out with the strength of the applied magnetic field ($\mu$) as the bifurcation parameter.

Symmetries

The equivariant branching lemma is crucial for one of the automatic exploration methods. To apply this lemma, the symmetry group of (LABEL:Ginzburg-Landau_vergelijking) needs to be determined. In [31] it is shown that, for square samples subject to a perpendicular, homogeneous magnetic field, this symmetry group is given by $S^{1}\times D_{4}$, with $S^{1}$ and $D_{4}$ respectively the circle group and symmetry group of the square. A more general result is given in proposition 2.1.

Proposition 2.1**.**

Let the sample $\Omega\subset\mathbb{R}^{2}$ and the applied magnetic field $\mathbf{H}_{0}$ both be invariant under the actions of a dihedral group $D_{m}=\langle\tau_{\omega},\sigma\rangle$ (with $\omega=2\pi/m$, $m\in\mathbb{N}$) defined by

[TABLE]

then (LABEL:Ginzburg-Landau_vergelijking) is invariant under the actions of $S^{1}\times D_{m}=\langle\theta_{\eta},\tau_{\omega},\sigma\rangle$, given by

[TABLE]

*A similar result holds for invariance under the actions of the cyclic group $C_{m}=\langle\tau_{\omega}\rangle$. *

Note that proposition 2.1 is only valid for 2-dimensional materials. For 3-dimensional ones an analogue result can be derived, but this is beyond the scope of the current paper.

Without any sample or magnetic field restrictions, the symmetry group of (LABEL:Ginzburg-Landau_vergelijking) is

[TABLE]

The shape of the material (and possibly its discretization) and magnetic field give rise to restrictions on this group. As example, for a triangular sample subject to a homogeneous magnetic field, the symmetry group is given by $S^{1}\times D_{3}$. A circular sample subject to a (centered) magnetic square would restrict the group to $S^{1}\times D_{4}$.

The continuous $S^{1}$ symmetry persists for every choice of material and magnetic field: every solution $\psi$ of (LABEL:Ginzburg-Landau_vergelijking) is actually a representative of a whole family of solutions $\{\theta_{\eta}\psi|\theta_{\eta}\in S^{1}\}$. Since we are mainly interested in the density of Cooper pairs $\rho_{c}=|\psi|^{2}$ and the equation $|\psi|=|\theta_{\eta}\psi|$ holds for each $\theta_{\eta}\in S^{1}$, it is sufficient to consider a single representative for each family [31].

Properties of the Jacobian operator

The solver presented in [33] is based on a Newton-Krylov algorithm. In order to apply this algorithm to (LABEL:Ginzburg-Landau_vergelijking), the (partial) Jacobian operator (to $\psi$) associated with this equation is required. In [31] an expression for this operator is derived:

[TABLE]

This operator is only linear when defined over $X$ and $Y$ as $\mathbb{R}$-vector spaces [31]. The Jacobian operator is self-adjoint with respect to the inner product

[TABLE]

which coincides with the natural inner product in $(L^{2}_{\mathbb{R}}(\Omega))^{2}$ [31]. The spectrum of (4) has been investigated in both [31] and [33]. These papers show that this spectrum is a subset of $\mathbb{R}$, but $\mathcal{J}_{\psi}(\psi,\mu)$ is generally not definite. Definiteness is completely determined by the state ($\psi,\mathbf{A}$). Solutions $\psi$ for which the Jacobian operator does not have positive eigenvalues are said to be physically stable, and physically unstable otherwise.

Another result that’s discussed in [31] and [33] is the zero eigenvalue induced by the continuous $S^{1}$ symmetry. For every $\psi_{s}\in X,\mu_{s}\in\mathbb{R}$ we have

[TABLE]

implying that the kernel of the Jacobian operator is never empty in a solution ($\psi_{s},\mu_{s}$) since $\{\i\psi_{s}\}\subset\ker(\mathcal{J}_{\psi}(\psi_{s},\mu_{s}))$.

Aside from the partial Jacobian operator to $\psi$, the one to $\mu$ is also required for numerical continuation. We will approximate the Jacobian operator $\mathcal{J}_{\mu}(\psi,\mu)$ by a second-order finite difference scheme. Finally, the automatic exploration methods require the full Jacobian of (LABEL:Ginzburg-Landau_vergelijking), defined by

[TABLE]

Other derivatives of the Ginzburg-Landau equations

The partial Hessian operators ($\mathcal{H}_{\psi\psi}(\psi,\mu)$, $\mathcal{H}_{\psi\mu}(\psi,\mu)$ and $\mathcal{H}_{\mu\mu}(\psi,\mu)$) of (LABEL:Ginzburg-Landau_vergelijking) are crucial for the determination of bifurcation points and new tangent directions. We will again apply second-order finite difference schemes to approximate the operators $\mathcal{H}_{\psi\mu}(\psi,\mu)$ and $\mathcal{H}_{\mu\mu}(\psi,\mu)$. Instead of approximating $\mathcal{H}_{\psi\psi}(\psi,\mu)$ by finite differences as well, we derive an exact expression for this operator. This is done by applying a similar technique as in [31] for the derivation of the Jacobian operator. Let $\psi,\varphi,\delta\psi\in X$ and $\mu\in\mathbb{R}$. We have

[TABLE]

The Hessian operator is obtained from this equation by neglecting the higher-order terms in $\delta\psi$:

[TABLE]

Note that this operator is independent of the strength $\mu$ of the applied magnetic field. The full Hessian operator is defined by

[TABLE]

and appears in the Lyapunov-Schmidt reduction-based method for constructing the tangent directions. Higher-order derivatives appear in this method as well. Using the same technique as before, one can show that the third partial derivative is given by

[TABLE]

This operator is independent of both the strength $\mu$ of the applied magnetic field and the order parameter $\psi$. These independencies yield

[TABLE]

Other partial derivatives are of the form

[TABLE]

These operators will be approximated by applying a second-order finite difference scheme. See [15] for more information on these schemes for general derivatives.

The discretized problem

In practice the grid and equations need to be discretized for numerical continuation to be efficiently performed. Discretization of the problem was discussed in detail in [31] and [33]. We will only provide the main results here.

Let $(x_{1},y_{1}),\dots,(x_{n},y_{n})$ be a set of discretization points of the sample $\Omega$. If possible, this set should be chosen in such a way that any symmetries of $\Omega$ are preserved. States $\psi\in X$ will be approximated by $\psi^{(h)}\in\mathbb{C}^{n}$, with

[TABLE]

Let $\{T_{i}\}_{i=1}^{m}$ be the Delauny triangulation of $(x_{1},y_{1}),\dots,(x_{n},y_{n})$ and $\{V_{j}\}_{j=1}^{n}$ its corresponding Voronoi tessellation. An edge $(x_{j},y_{j})-(x_{k},y_{k})$ of a triangle $T_{i}$ is denoted as $e_{j,k}$. We will first discretize the operator $(-\i\nabla-\mathbf{A}(\mu))^{2}$, this discretization is given by $\mathcal{K}^{(h)}(\mu)$, defined by the property [33]

[TABLE]

For a triangle $T_{i}$ consisting of the edges $e_{j,k}$, $e_{k,l}$ and $e_{l,j}$, the coefficient $\alpha_{j,k}^{(i)}\in\mathbb{R}$ is given by the formula

[TABLE]

The coefficients $\alpha_{k,l}^{(i)}$ and $\alpha_{l,j}^{(i)}$ are defined by a similar formula. The values $U_{j,k}(\mu)\in\mathbb{C}$ are given by

[TABLE]

with $\mathbf{A}(\mu,\omega)$ the magnetic vector potential evaluated at the location $\omega\in\Omega$. Using the discretization $K^{(h)}$ of $(-\i\nabla-\mathbf{A}(\mu))^{2}$, we discretize the Ginzburg-Landau equation (LABEL:Ginzburg-Landau_vergelijking) as [33]

[TABLE]

Other derivatives are discretized in a similar way. The discrete Jacobian is self-adjoint with respect to the inner product

[TABLE]

which is the discrete version of (5). In the remainder of the article we will derive methods and algorithms for the continuous problem (LABEL:Ginzburg-Landau_vergelijking). The extension of the results to the discrete case is straightforward.

3 Numerical Continuation

Pseudo-arclength continuation

To study the influence of system parameters like the magnetic field strength on the behaviour of patterns in superconductors, numerical continuation is performed. The most natural method for this purpose is considered to be pseudo-arclength continuation [3, 25]. This algorithm is a predictor-corrector method: from a given solution ($\psi^{(0)},\mu^{(0)}$) a prediction ($\tilde{\psi}^{(1)},\tilde{\mu}^{(1)}$) for the next point of the solution branch is made, which is then corrected by a Newton-Krylov algorithm. The prediction ($\tilde{\psi}^{(1)},\tilde{\mu}^{(1)}$) is made by perturbing the given solution in the direction of the tangent ($\dot{\psi}^{(0)},\dot{\mu}^{(0)}$) of the solution curve:

[TABLE]

with $\Delta s\in\mathbb{R}$ a sufficiently small step size. In practice an approximation to the tangent ($\dot{\psi}^{(0)},\dot{\mu}^{(0)}$) is used, which speeds up the algorithm. After its construction the (approximated) tangent is orthogonalized towards $\begin{pmatrix}\i\psi^{(0)}&0\end{pmatrix}^{T}$, the null vector of the (full) Jacobian induced by the continuous $S^{1}$ symmetry. Without this orthogonalization numerical continuation might eventually fail due to solution branches falsely reverting. The correction step now consists of solving the system

[TABLE]

for the unknowns $\psi^{(1)}$ and $\mu^{(1)}$ by a Newton-Krylov algorithm. Note that perpendicularity with respect to the inner product (5) is used.

A solver for the Jacobian system

Application of a Newton-Krylov algorithm to (11) requires a solver for linear systems of the form

[TABLE]

Since we want to trace solution branches for both physically stable and unstable solutions, the solver constructed in [33] will be used as a base for solving these linear systems. A summary of the construction and properties of the solver is given below, based on the analysis done in [33].

In [33] the problem

[TABLE]

is considered: a solver is constructed for linear systems with the partial Jacobian operator to $\psi$. Though this operator is singular in solutions of (LABEL:Ginzburg-Landau_vergelijking) and ill-conditioned close to such solutions, Krylov methods can be applied given that the right-hand side $y$ lies in the range of the operator $\mathcal{J}_{\psi}(\psi,\mu)$. The sparsity structure of $\mathcal{J}_{\psi}(\psi,\mu)$ promotes the use of a Krylov method as well, since only the application of an operator to given vectors needs to be known. Since the Jacobian $\mathcal{J}_{\psi}(\psi,\mu)$ is self-adjoint with respect to the inner product (5), but is generally indefinite, the Krylov method MINRES (adapted for the inner product (5)) is chosen.

Preconditioning is also discussed in [33]: approximate inverses of the operator

[TABLE]

are used as a preconditioner. Approximate inversion is realized by an algebraic multigrid (AMG) strategy. Numerical evidence of [33] shows that the approximate inversion with a single V-cycle yields the fastest solver: in this case the number of Krylov iterations is independent of the dimension of the solution space.

An implementation of the complete solver is provided in Python (package PyNosh, see [18]), and forms an important part of our automatic exploration algorithm. Though this solver is constructed for systems of the form (13), it will be combined with block elimination techniques for application to the bordered linear system (12). For more information about these techniques, see e.g. [4, 22].

4 Detection and determination of bifurcation points

Bifurcation points of (LABEL:Ginzburg-Landau_vergelijking) are points $(\psi_{b},\mu_{b})$ for which $\mathcal{J}_{\psi}(\psi_{b},\mu_{b})$ has a zero eigenvalue [3, 25] (we ignore zero eigenvalues induced by the $S^{1}$ symmetry in this definition). Both intersections of solution branches (branch points) and changes in solutions physical stability always correspond with a bifurcation, making their determination an essential part of our algorithm. In this section a condition that indicates the proximity of bifurcations is described first, then a non-linear system of equations is constructed to determine the points.

Detection

We use a method based on the analysis done in [27] to detect bifurcation points: The lowest magnitude Ritz values of $\mathcal{J}_{\psi}(\psi,\mu)$ are calculated after each numerical continuation step. When a point $(\psi,\mu)$ is close to a bifurcation, one of these Ritz values will approximate zero. Furthermore, if the eigenvalue associated with the bifurcation changes its sign after passing this point, the same happens with the corresponding Ritz value. These properties give rise to the following definition:

Definition 4.1**.**

Near bifurcation condition*

Let $\zeta_{1}^{(0)},\zeta_{2}^{(0)},\dots,\zeta_{k}^{(0)}$ and $\zeta_{1}^{(1)},\zeta_{2}^{(1)},\dots,\zeta_{k}^{(1)}$ be the $k$ lowest magnitude Ritz values of the partial Jacobian operator in respectively $(\psi^{(0)},\mu^{(0)})$ and $(\psi^{(1)},\mu^{(1)})$, two consecutive points of the solution branch calculated by pseudo-arclength continuation. Assume the Ritz values are ordered according to the same eigenvectors. A bifurcation point is close to $(\psi^{(1)},\mu^{(1)})$ if either*

[TABLE]

*In this case we say that the point $(\psi^{(1)},\mu^{(1)})$ satisfies the near bifurcation condition. *

For the calculation of the examples in section 6 we used $k=5$ Ritz values. Note that values induced by the $S^{1}$ symmetry are ignored in definition 4.1.

Except for the detection process, Ritz values are also used when determining physical stability: if the Jacobian in a point has a (high) positive Ritz value, this point is assumed to be unstable. Another use of Ritz values and Ritz vectors is deflation: multiplying the system with a certain projection (chosen to eliminate eigenvectors corresponding to low magnitude eigenvalues) accelerates Krylov algorithms [16, 17, 19].

Determination

Though the near bifurcation condition described in the previous section gives us an indication of the proximity of bifurcations, we still need to determine these points more precisely. For this purpose we construct a non-linear system of equations that’s only satisfied for bifurcation points $(\psi_{b},\mu_{b})$ of (LABEL:Ginzburg-Landau_vergelijking), this system is based on the property

[TABLE]

The system we consider is given by

[TABLE]

and is similar to the one described in [27]. When we find a point that satisfies the near bifurcation condition, we use it as an initial guess to solve (14) for $(\psi_{b},\mu_{b},\phi)$ by a Newton-Krylov algorithm. The approximate eigenvector corresponding to the near-zero Ritz value is used as an initial guess for the null vector $\phi$. Note that the guess for $\phi$ should be normalized after each Newton step, to prevent the solver from converging to the trivial solution $\phi=0$. In each Newton step the Jacobian system of (14) has to be solved, this system is given by

[TABLE]

Together with the application of block elimination techniques [4, 22] and deflation [16, 17, 19], this equation is solved by employing the solver described in section 3. Note that both (partial) Jacobians and Hessians of (LABEL:Ginzburg-Landau_vergelijking) appear in (15). As was shown in section 2, these expressions are either known, or are approximated by a finite difference scheme.

5 Construction of new tangent directions

In order to calculate new solution curves that emerge from a branch point, it is sufficient to determine the tangent directions to these curves. Given these directions, the curves are calculated by employing the pseudo-arclength continuation algorithm (see section 3). We will derive two algorithms for the construction of the tangent directions, one based on Lyapunov-Schmidt reduction, another based on the equivariant branching lemma. In the remainder of this section ($\psi_{b},\mu_{b}$) is assumed to be a bifurcation point of (LABEL:Ginzburg-Landau_vergelijking) and we denote $\mathcal{GL}^{(b)}$, $\mathcal{J}^{(b)}$, …as the evaluations of $\mathcal{GL}$, $\mathcal{J}$, …at ($\psi_{b},\mu_{b}$).

Lyapunov-Schmidt reduction

A popular approach for analysing bifurcations is the Lyapunov-Schmidt method: in a neighbourhood of the bifurcation ($\psi_{b},\mu_{b}$), the problem is reduced to low dimensional systems of algebraic equations that contain all of the information of the bifurcations behaviour. The solutions of these systems are then used to construct the tangent directions. In [28] the method is applied to general problems, leading to three equations that form the base of our first tangent construction algorithm. A short overview of the deduction of these equations (applied to the Ginzburg-Landau problem) is given below, based on the analysis done in [28].

The following assumptions are made:

$\mathcal{J}_{\psi}^{(b)}$ is a Fredholm operator of index [math] and zero is a semi-simple eigenvalue of it.

• -

Orthonormal bases for the kernels of $\mathcal{J}_{\psi}^{(b)}$ and $\mathcal{J}_{\psi}^{(b)^{*}}$ (the adjoint Jacobian with respect to (5), see [31] for more details) have been determined:

[TABLE]

for a certain $n\in\mathbb{N}$. Note that the trivial null vector $\i\psi_{b}$ induced by the $S^{1}$ symmetry is ignored in this analysis and does not count towards the dimension of these kernels. In the case of the Ginzburg-Landau equation the Jacobian $\mathcal{J}_{\psi}^{(b)}$ is self-adjoint, so we actually have $\ker\left(\mathcal{J}_{\psi}^{(b)}\right)=\ker\left(\mathcal{J}_{\psi}^{(b)^{*}}\right)$.

• -

We have

[TABLE]

Note that if this assumption does not hold and $n$ equals $1$, the bifurcation $(\psi_{b},\mu_{b})$ is in fact a turning point and no new tangent directions need to be calculated.

In practice, the kernels of $\mathcal{J}_{\psi}^{(b)}$ and $\mathcal{J}_{\psi}^{(b)^{*}}$ are approximated using numerical methods (e.g. an inverse iteration algorithm [7]).

A solution curve ($\psi(s),\mu(s)$) with $(\psi(0),\mu(0))=(\psi_{b},\mu_{b})$ is considered, where $s$ is the parametrization parameter. The fundamental theorem of linear algebra states that the space $X$ can be decomposed into $\ker\left(\mathcal{J}_{\psi}^{(b)}\right)$ and $\text{im}\left(\mathcal{J}_{\psi}^{(b)^{*}}\right)$, this implies that we can write $\psi(s)$ and $\mu(s)$ as

[TABLE]

Applying a Taylor expansion, this yields

[TABLE]

Note that the tangent ($\dot{\psi}(0),\dot{\mu}(0)$) to the solution curve is given by

[TABLE]

Equation (LABEL:Ginzburg-Landau_vergelijking) is rewritten (using again a Taylor expansion),

[TABLE]

with $R(\psi-\psi_{b},\mu-\mu_{b})$ given by

[TABLE]

The term $r_{k}$ in (18) represents the coefficient of $s^{k}$ in $R(\psi-\psi_{b},\mu-\mu_{b})$, this is explicitly given by

[TABLE]

The following three equations are now derived in [28] (for $k=1,\dots,l$):

[TABLE]

these equations will form the basis of the upcoming algorithms. For more details on Lyapunov-Schmidt reduction, see e.g. [27, 28].

An algorithm for the construction of the tangent directions

In [28] a general algorithm is deduced to construct reduced systems of equations for $\alpha_{1},\dots,\alpha_{n},\beta_{k-1}$. Together with (17) the solutions of these systems yield the tangent directions. For this method equations (20), (21) and (22) are used to write the terms $r_{k}$ and $w_{k}$ as polynomials in the variables $\alpha_{1},\dots,\alpha_{n},\beta_{1},\dots,\beta_{k-1}$. Under the assumptions

[TABLE]

it is possible to write the term $\beta_{k}$ as a polynomial in the variables $\alpha_{1},\dots,\alpha_{n}$ as well. The vector $v^{(0)}$ is implicitly defined via $w_{1}=\beta_{1}v^{(0)}$. The method described in [28] is given in algorithm 5 (rewritten to emphasize the construction of the reduced systems).

{textalgorithm}

Creation of tangent directions based on [28]

Initial:

For $k=1$, use the equation

[TABLE]

to write $w_{1}\in\text{im}\left(\mathcal{J}_{\psi}^{(b)^{*}}\right)$ as a polynomial in $\beta_{1}$.

Iteration:

For $k=2,3,\dots$ do:

Step 1: Use the polynomial expression of $w_{k-1}$ (and $\beta_{k-2}$ if $k\geq 3$) together with (19) to write $r_{k}$ as a polynomial in $\alpha_{1},\dots,\alpha_{n},\beta_{k-1}$.

Step 2: Substitute this polynomial expression in the equation

[TABLE]

This yields a reduced system of equations for $\alpha_{1},\dots,\alpha_{n},\beta_{k-1}$.

Step 3: Check whether the system has real and isolated solutions $\alpha_{1}(\beta_{k-1}),\dots,\alpha_{n}(\beta_{k-1})$ for $\beta_{k-1}$ in an open interval. Solve the system for these solutions, each one corresponds to a different tangent direction.

Step 4a: If only isolated solutions exist in the reduced system, stop the algorithm.

Step 4b: If non-isolated solutions exist as well, use the system to write $\beta_{k-1}$ as a polynomial of $\alpha_{1},\dots,\alpha_{n}$ and use this result to write $r_{k}$ as a polynomial in these variables as well. Substitute this expression for $r_{k}$ into the system

[TABLE]

and solve it to find a polynomial expression for $w_{k}$ in $\alpha_{1},\dots,\alpha_{n},\beta_{k}$.

The complexity of the polynomial expressions in algorithm 5 strongly depends on the dimension $n$ of $\ker\left(\mathcal{J}_{\psi}^{(b)}\right)$. If this kernel is one-dimensional, the algorithm only requires a single iteration: the first reduced system that is constructed does not contain non-isolated solutions. In this case the algorithm is actually equivalent to the algebraic branching equation [3, 25], the reduced system for $\alpha_{1}$ and $\beta_{1}$ is given by

[TABLE]

where $v^{(0)}$ is solved from the equation

[TABLE]

Equation (27) is solved for $\alpha_{1}$ and $\beta_{1}$, together with $w_{1}=\beta_{1}v^{(0)}$ the tangent directions are constructed by applying (17). The directions are determined up to multiplication by a constant, and should be normalized before further use.

For problems that exhibit symmetry, bifurcations arise with Jacobian kernel dimension $n\geq 2$ and the polynomials in 5 become more complex. In [28] no exact expressions for these polynomials were derived and hence only implicit reduced systems of equations are given, which limits the use of the algorithm in its current form by numerical software. Furthermore, it is not specified how these systems can be checked for real and isolated solutions. In the remainder of this section we derive exact polynomial expressions that yield explicit systems for the case $n=2$. A simple rule is used to check these for real and isolated solutions. No other cases than $n=1$ and $n=2$ were encountered for the Ginzburg-Landau equation on two-dimensional samples. Set $n=2$, assumptions (23),(24) and (25) allow us to choose $\phi_{1}^{*}$ and $\phi_{2}^{*}$ such that:

[TABLE]

A modification of algorithm 5 for the case $n=2$ is given by algorithm 5. In this last algorithm the coefficients necessary to write the different polynomial expressions are derived in order to determine explicit reduced systems of equations for $\alpha_{1},\dots,\alpha_{n},\beta_{k-1}$. A thorough derivation of algorithm 5 is provided in appendix A. {textalgorithm} Creation of tangent directions for the $n=2$ case

Initial:

Step 1: Solve the equation

[TABLE]

by applying the solver described in section 3. Note that $w_{1}$ (defined in (16)) can be written as

[TABLE]

Step 2: Calculate the following terms in the space $Y$:

[TABLE]

Note that $r_{2}$ (defined in (19)) can be written as

[TABLE]

Step 3: Calculate the scalars

[TABLE]

These values form the coefficients of the following reduced system of equations:

[TABLE]

Step 4: Check whether $b^{(2)}a_{0}^{(2,1)}=b^{(3)}a_{0}^{(2,2)}$, $a_{2}^{(2,2)}=0$ and $\forall i=1;2:b^{(2)}a_{i}^{(2,1)}=b^{(1)}a_{i-1}^{(2,2)}+b^{(3)}a_{i}^{(2,2)}$.

Step 5a: If this is the case, stop the algorithm: system (30) only contains isolated solutions, these correspond to the tangent directions.

Step 5b: If this is not the case, the only isolated solution of (30) is $(\alpha_{1},\alpha_{2},\beta_{1})=(0,0,a)$ with $a\in\mathbb{R}$. The system contains non-isolated solutions as well. Set $\kappa^{(1)}_{-1}=\kappa^{(1)}_{2}=0$ and calculate

[TABLE]

Note that $\beta_{1}$, $r_{2}$ and $x_{1}$ (defined in (16) and (19)) can respectively be written as

[TABLE]

Iteration:

For $k=3,4,\dots$ do:

Step 1: Solve the equations

[TABLE]

by applying the solver described in section 3. Note that $w_{k-1}$ (defined in (16)) can be written as

[TABLE]

Step 2: Calculate the following terms in the space $Y$:

[TABLE]

An efficient method to calculate these terms is provided in appendix B. Note that $r_{k}$ (defined in (19)) can be written as

[TABLE]

Step 3: Calculate the scalars

[TABLE]

These values form the coefficients of the following reduced system of equations:

[TABLE]

Step 4: Check whether $b^{(2)}a_{0}^{(k,1)}=b^{(3)}a_{0}^{(k,2)}$, $a_{k}^{(k,2)}=0$ and $\forall i=1,\dots,k:b^{(2)}a_{i}^{(k,1)}=b^{(1)}a_{i-1}^{(k,2)}+b^{(3)}a_{i}^{(k,2)}$.

Step 5a: If this is the case, stop the algorithm: system (36) only contains isolated solutions, these correspond to the tangent directions.

Step 5b: If this is not the case, the only isolated solution of (36) is $(\alpha_{1},\alpha_{2},\beta_{k-1})=(0,0,a)$ with $a\in\mathbb{R}$. The system contains non-isolated solutions as well. Set $\kappa^{(k-1)}_{-1}=\kappa^{(k-1)}_{k}=0$ and calculate

[TABLE]

Note that $\beta_{k-1}$, $r_{k}$ and $x_{k-1}$ (defined in (16) and (19)) can respectively be written as

[TABLE]

The initial steps of the algorithm always yield a system with isolated solution $(\alpha_{1},\alpha_{2},\beta_{1})=(0,0,a)$ ($a\in\mathbb{R}$). Application of (17) and (28) leads to a first tangent direction. If this system contains other isolated solutions as well, the corresponding tangent directions are constructed with the same formulas. In the other case the algorithm is continued, possibly yielding further systems with a sole isolated solution $(\alpha_{1},\alpha_{2},\beta_{k-1})=(0,0,a)$ ($a\in\mathbb{R}$). Relation (31) results in $\beta_{1}=0$, consequently (17) equals zero and does not actually correspond to a real tangent direction. Eventually the algorithm will find a system that only contains isolated solutions, the tangent directions are constructed from these by application of (17), (28) and (31).

Solving the reduced system of equations

Algorithm 5 eventually gives rise to a system of equations of the form

[TABLE]

with $k\geq 2$ and known coefficients $a_{0}^{(k,1)},\dots,a_{k}^{(k,1)},a_{0}^{(k,2)},\dots,a_{k}^{(k,2)},b^{(1)},b^{(2)},b^{(3)}$ such that (40) only contains isolated solutions.

To find solutions of (40) we fix $\beta_{k-1}$. This does not reduce the amount of solutions that will be found since the tangent directions are determined up to normalization. Though in some cases it is possible to solve (40) exactly (for e.g. $k=2$ solving this equation is similar to determining the intersection of two conics, for which multiple algorithms exist), we will approximate its solutions with a Newton algorithm. This requires the Jacobian of (40), given by

[TABLE]

In practice the solutions are found by applying the Newton algorithm for a high amount of different initial guesses. Together with this Newton solver, algorithm 5 allows for the construction of the tangent directions arising from bifurcations with Jacobian kernel dimension $2$.

The case $n=2$ only appears when the underlying problem exhibits symmetry and the bifurcation point is symmetric as well [27, 28]. In [27] it is proved that for a bifurcation invariant under the actions of $C_{m}$ or $D_{m}$ ($m\geq 2$), the reduced system (40) has isolated solutions for $k=m-1$ and is invariant under the actions of $C_{m}$ respectively $D_{m}$ as well. Due to the $C_{m}$ or $D_{m}$ invariance of (40), a lot of its solutions result in tangent directions to equivalent solution branches: application of the action $\tau_{\omega}$ (see proposition 2.1) on solutions of one branch lead to solutions of another. To prevent equivalent branches from being calculated, only one solution of (40) for each class of emerging branches is selected. For $m$ even/odd this results in respectively $3$ or $2$ different solutions.

Algorithm 5 was used for the results of the triangular and rotation-symmetric materials in section 6. In both examples, it yielded the required tangent directions for each branch point with Jacobian kernel dimension $2$.

The equivariant branching lemma

Though algorithm 5 can always be applied for bifurcations with Jacobian kernel dimension $2$, it becomes less practical for strongly symmetrical problems: if the bifurcation is invariant under the actions of $C_{m}$ or $D_{m}$, $k\geq m-1$ iterations might be required before a reduced system that only contains isolated solutions is obtained [28]. Before this system is found, at least $((m-1)m-2)/2$ solves of a linear equation with the Jacobian need to be performed. This is unwanted for $m\gg 1$. For $m\geq 4$ one can still apply the algorithm to determine an equation for $\beta_{1}$ and $w_{1}$ as a polynomial in the variables $\alpha_{1}$ and $\alpha_{2}$ (see (28) and (31)), this requires only one solve of the linear problem. For $D_{m}$ symmetric bifurcations the unknowns $\alpha_{1}$ and $\alpha_{2}$ can then alternatively be determined by application of the equivariant branching lemma. We will describe this second method in the remainder of this section.

We call a point $(\psi,\mu)$ $\gamma$-invariant (for $\gamma\in D_{m}$) if its transformation under the action $\gamma$ remains in the same solution family:

[TABLE]

Consider a branch point invariant under the actions of $D_{m}$ (with $m\geq 4$). If the eigenvalue corresponding to the bifurcation crosses the origin with non-zero speed, the equivariant branching lemma guarantees the existence of symmetry-breaking branches: new solution curves with symmetries corresponding to the conjugacy classes of the isotropy subgroups of $D_{m}$ will emerge from the branch point [20, 24]. These conjugacy classes are $\langle\sigma\rangle$ and $\langle\tau_{\eta_{0}}\sigma\rangle$ (with $\eta_{0}=2\pi/m$) if $m$ is even and just $\langle\sigma\rangle$ if $m$ is odd.

Using the knowledge of the symmetry $\gamma$ of a solution curve that emerges from the branch point, it is possible to calculate the remaining unknowns $\alpha_{1}$ and $\alpha_{2}$ that construct the tangent direction ($\dot{\psi},\dot{\mu}$) to this curve. First, a representative $(\psi^{\prime}_{b},\mu_{b})$ of the bifurcation point for which $\gamma\psi^{\prime}_{b}=\psi^{\prime}_{b}$ needs to be determined. In practice this is realized by applying a least-squares method: the point

[TABLE]

satisfies $\gamma\psi^{\prime}_{b}\approx\psi^{\prime}_{b}$. Let $(\psi,\mu)$ be a point of the $\gamma$-symmetric solution curve that emerges from $(\psi^{\prime}_{b},\mu_{b})$: $\gamma\psi=\psi$. If $(\psi,\mu)$ is chosen sufficiently close to $(\psi^{\prime}_{b},\mu_{b})$, we have

[TABLE]

with the tangent direction $(\theta_{\eta_{\gamma}}\dot{\psi},\dot{\mu})$ given by (see (17),(28) and (31))

[TABLE]

with unknowns $\alpha_{1}$ and $\alpha_{2}$, and $v_{0}$, $\kappa_{0}^{(1)}$ and $\kappa_{1}^{(1)}$ defined by the initial steps of algorithm 5. Since $\gamma\psi=\psi$ and $\gamma\psi^{\prime}_{b}=\psi^{\prime}_{b}$, equation (41) implies $\gamma\theta_{\eta_{\gamma}}\dot{\psi}=\theta_{\eta_{\gamma}}\dot{\psi}$. This result is combined with a least-squares method to yield the following formulas for $\alpha_{1}$ and $\alpha_{2}$:

[TABLE]

Substitution of $\alpha_{1}$ and $\alpha_{2}$ in (42) gives the desired $\gamma$-symmetric tangent.

Though the use of the equivariant branching lemma is an efficient alternative for determining the tangent directions, it can only be applied if the bifurcation point is $D_{m}$ symmetric (for a certain $m\geq 4$). Such bifurcations cannot arise if the material or the magnetic field only exhibits rotational symmetry ($C_{m}$). It is however still possible to apply algorithm 5 in this case.

Calculation of the bifurcation diagram for the square material of section 6 made use of the method based on the equivariant branching lemma. For each branch point with Jacobian kernel dimension $2$, it was able to construct the required tangent directions.

6 Numerical results

This section contains the solution landscapes for multiple shapes of materials, calculated with an extension of the Python package PyNCT [11]. We will only consider small-scale (mesoscopic) shapes, which are for example of interest when building nanoscale fluxonics devices to use in e.g. SQUIDs [26], RSFQ processors [14] and supercomputers [10, 23, 29]. The materials used in these devices are typically shaped as a triangle, square or disc. Our methods however allow for general shapes, which will be indicated by the second example.

The materials are subject to a homogeneous magnetic field, numerical continuation is performed in this fields strength ($\mu$). For each shape, the homogeneous solution $\psi\equiv 1$ in absence of a magnetic field ($\mu=0$) is taken as the starting point. The other points and bifurcations are determined using the tools discussed in this paper, yielding an automatically explored solution landscape. Note that the bifurcation diagrams are plotted in terms of the expression

[TABLE]

which is part of the Gibbs energy (1) (see [31] for more details).

Triangular material

Our first example is a material shaped as a regular triangle, the size of the edges is $6$ (scaled in units of the coherence length $\xi$, see section 2). The corresponding bifurcation diagram, together with representative solutions for each solution curve, is given in figure 5.

A total of $16$ bifurcation points were found during the continuation: $7$ turning points and $9$ branch points. For all of these branch points, the corresponding Jacobian contained a kernel with dimension $n=2$. Application of algorithm 5 yielded the tangent directions to the emerging curves. Note that only the initial steps of this algorithm had to be executed.

The trivial solution $\psi\equiv 1$, $\mu=0$ is part of solution branch A. For values of $\mu$ between [math] and approximately $1.21$, branch A is stable. For non-zero field strength, the solutions consist of zones of low supercurrent density near the edges of the domain. At $\mu\approx 1.21$ branch A destabilizes. If $\mu$ is further increased along this branch, three vortices move in from the midpoints of the edges towards the center. At field strength $\mu\approx 1.90$, the three vortices arrive at the center and form a giant vortex with multiplicity $3$. Further increase shows how this giant vortex breaks up again and three separate vortices slowly move away from the center. At $\mu\approx 2.03$ branch A restabilizes, and at a value of $\mu\approx 2.45$ the solution curve connects to the trivial zero solution $\psi\equiv 0$. Solutions between these two values consist of three separate vortices near the centre of the material.

Except for branch A, two more main solution branches (with full $D_{3}$ symmetry) were discovered: branches B and C. Branch B consists of solutions with a single vortex in the centre of the material. It is stable for values of $\mu$ between approximately $0.68$ and $1.78$. When the fields strength is increased along solution branch B, three vortices start to appear near the edges of the domain. At $\mu\approx 2.10$, these vortices start to move towards the single, central vortex. These do however not join together, as at $\mu\approx 2.23$ solution curve B connects to the trivial zero solution. Except for the bifurcations at $\mu\approx 0.68$ and $\mu\approx 1.78$, branch B contains one at $\mu\approx 2.08$ as well.

For low values of $\mu$, solution branch C also consists of solutions with a single vortex in the centre. As the field strength is increased along this branch, three more vortices, with opposite vorticity to the central one, move from the edges of the domain towards this central vortex. At $\mu\approx 1.26$, these vortices join the central one to form one giant vortex with multiplicity $2$. This vortex remains in the centre of the material during further increase of the fields strength. Solution branch C is unstable for low values of $\mu$, but stabilizes between $\mu\approx 1.41$ and $\mu\approx 2.25$. For high values of the fields strength, the edges of the domain show zones of low supercurrent density as well. The curve connects to the trivial zero solution at $\mu\approx 2.34$. Four bifurcations were found in total on solution curve C: the ones at $\mu\approx 1.41$ and $\mu\approx 2.25$, where the stability of the solutions change, and two more at $\mu\approx 0.30$ and $\mu\approx 1.04$.

Solution branch D connects the three main solution branches. Starting from the bifurcation point of branch A at $\mu\approx 1.21$ and decreasing the parameter, a single vortex moves from the edge of the material towards the centre. It arrives there at $\mu\approx 0.68$, where the curve connects with branch B. Decreasing the fields strength further shows how the vortex moves towards the opposing corner of the triangle, until a turning point at $\mu\approx 0.65$ is encountered. When $\mu$ is increased from this value, the vortex moves further towards the corner until $\mu\approx 1.12$. Around this value, a second vortex appears from the corner the original one was moving to, and both vortices start to move towards the centre of the material. At $\mu\approx 1.33$ they are joined into a giant vortex with multiplicity $2$ (this does not happen in the centre). Further increase shows how this giant vortex is split into two vortices again, one of which moves back towards the corner, the other one towards the centre. A second turning point is encountered at $\mu\approx 1.86$. Decreasing the fields strength from this value, both vortices arrive at their destination: at $\mu\approx 1.78$ solution branch D connects with branch B, this solution consists of a single vortex in the centre and zones of low supercurrent density at the edges. Decreasing $\mu$ from this bifurcation point shows how a second vortex moves from one edge of the triangle towards the centre, where it joins the central vortex into a giant one with multiplicity $2$. This corresponds with the bifurcation point at $\mu\approx 1.41$. Finally, when the parameter is decreased even further, the giant vortex splits into two vortices, which both move towards the edges of the domain. Solution curve D reconnects with itself at $\mu\approx 1.21$.

Except for branch D, branches A and C are connected through solution branch F as well. Starting from the solution of branch A at $\mu\approx 2.03$, consisting of three separate vortices near the centre of the material, the vortices move towards one of the materials corners when the field strength is decreased. A turning point is encountered at $\mu\approx 1.95$. Increase of $\mu$ from this value shows how one of the vortices completely moves to the corner while the other two return to the centre, until a turning point at $\mu\approx 2.26$ is reached. Decreasing $\mu$, zones of low supercurrent density start to appear at the edges, until a value of $\mu\approx 2.25$, where solution branch F connects with branch C. When the fields strength is decreased further, a vortex moves from one of the edges towards the centre, eventually leading back to the bifurcation at $\mu\approx 2.03$.

Solution branch E connects branches B and C as well. We start from the bifurcation of branch B at $\mu\approx 2.08$, consisting of a solution with one central vortex. Decreasing $\mu$, two vortices - one from a triangle’s corner, the other from the opposite edge - move towards the central vortex. Before the one moving from the edge reaches the centre, the other two join together into a giant vortex with multiplicity 2. Further decrease shows how this vortex splits back into two smaller ones, after which all of the vortices start to move towards the edges. Eventually two vortices remain: one in a corner, another at the opposite edge. These two vortices move towards each other, until they are joined into a giant one with vorticity 2 at the centre of the material. This happens at $\mu\approx 1.04$, where branch E connects with branch C. Decreasing $\mu$ further, this vortex splits again, with both vortices moving towards different corners. A turning point is encountered for $\mu\approx 0.99$. Increasing the fields strength from this value, the vortices continue to move towards the corners. Before these corners are reached a new vortex appears at the intermediate edge. This new vortex starts to move towards the centre, while the other two vortices now move towards the remaining edges. After encountering another turning point at $\mu\approx 2.09$, solution branch E eventually reconnects with itself at the bifurcation for $\mu\approx 2.08$.

Finally, branch G connects branch C with its symmetric counterpart for negative values of $\mu$. The solution branch starts at $\mu\approx 0.30$, a bifurcation point of branch C that shows a pattern with one central vortex connected to the edges by zones of low supercurrent density. Decreasing the fields strength, the vortex moves to one of the corners, the zone connecting it to the opposite edge stretching towards it. Increasing the parameter $\mu$ shows how the vortex moves towards one of the edges, the two zones connecting it to the other edges again stretching towards it. Almost immediately a turning point is encountered (approximately at $\mu\approx 0.30$). Decreasing $\mu$ from this point shows how the vortex first keeps moving towards the edge, later it is joined with the zones of low supercurrent density, yielding a single zone stretching from one edge to the other, bent away from the intermediate corner.

Rotation-symmetric material

The second material is shaped as a four-pointed star. It only exhibits rotational symmetry, the corresponding symmetry group of the problem is hence given by $C_{4}$. The domain is chosen such that it fits into a square with edge $6$: the distance between the central point and outer corners is given by $\sqrt{10}$, the one between the central point and inner corners by $\sqrt{2}$. The size of the short and long edges of the star itself are respectively given by $2$ and $2\sqrt{2}$. The bifurcation diagram, and representative solutions, associated with this material is given in figure 6. Several close-ups of the diagram are given in figure 7.

Multiple bifurcation points were again discovered during the continuation, this time not all of the branch points corresponded to a Jacobian with a $2$-dimensional kernel. Several ones only contained a kernel consisting of a single null vector, for these points application of algorithm 5 resulted in the (single) new tangent direction. As has been stated in section 5, for these cases the algorithm is equivalent to applying the algebraic branching equation. Other branch points did correspond to a Jacobian with a $2$-dimensional kernel, for these points a single iteration of algorithm 5 ($k=3$) was applied to construct the tangent directions. Note that these directions could not be determined by application of the equivariant branching lemma since the problem only exhibits rotational symmetry.

Contrary to the triangular material, four main branches (with full $C_{4}$ symmetry) were discovered. These are labelled A, H, J and K in figure 6. Branch A contains the trivial solution $\psi\equiv 1$, $\mu=0$. Increasing the strength of the applied magnetic field from this solution, four zones of low supercurrent density form near the inner corners of the material. It connects to the trivial zero solution $\psi\equiv 0$ at $\mu\approx 2.23$. The solution branch is stable for values of $\mu$ between [math] and approximately $1.71$. Other bifurcations are situated at $\mu\approx 1.80$, $\mu\approx 2.05$ and $\mu\approx 2.20$.

All of the solutions of branch H contain a single vortex in the centre of the material. For high values of the field strength ($\mu\gtrsim 2.15$), zones of low supercurrent density start to form near the inner corners as well. The solution branch is stable between values $\mu\approx 1.03$ and $\mu\approx 2.48$. Two more bifurcations were discovered near $\mu\approx 2.65$, immediately before the curve connects to the trivial zero solution $\psi\equiv 0$.

For low values of $\mu$, solutions that belong to branch J show a cross pattern of low superconductivity: zones of low supercurrent density appear in the centre and stretch towards the inner corners. Increasing the field strength along this branch, a bifurcation with the same pattern is discovered near $\mu\approx 0.33$. The pattern of the solutions changes for $\mu\gtrsim 0.95$: a vortex with multiplicity 2 appears in the centre of the material. Further increase of $\mu$ gives rise to more bifurcations at $\mu\approx 1.48$, $\mu\approx 1.79$ and $\mu\approx 2.01$. At this last bifurcation, branch J stabilizes. It remains stable until it connects with the trivial $\psi\equiv 0$ solution near $\mu\approx 2.90$.

The final main solution branch discovered for this problem is branch K. Similar to branch H, its solutions contain a single vortex in the centre for low values of the field strength. Increasing this strength shows how four vortices, with opposite vorticity to the central one, form near the inner corners ($\mu\gtrsim 0.20$), afterwards these vortices move towards the centre ($\mu\gtrsim 1.12$) and join with the central one into a giant vortex with multiplicity 3 ($\mu\gtrsim 2.00$). The branch eventually connects to the trivial zero solution $\psi\equiv 0$ at $\mu\approx 2.66$. A total of 6 bifurcations were discovered for solution branch K, at values of $\mu\approx 0.04$, $\mu\approx 0.47$, $\mu\approx 0.91$, $\mu\approx 2.08$, $\mu\approx 2.62$ and $\mu\approx 2.64$. Note that no stable solutions have been found for this branch.

The four main solution branches are again mutually connected by other branches. Branch A and H are connected through branches B and C. Starting from the bifurcation of branch A at $\mu\approx 1.71$ and decreasing the field strength, branch B shows how a single vortex moves from one of the inner corners towards the centre of the material, appearing there at the bifurcation $\mu\approx 1.03$ of branch H. Branch C shows the same behaviour, with the vortex appearing from one of the outer corners.

Branches D and E connect branch A with J. Decrease of the field strength from the bifurcation at $\mu\approx 1.80$ of branch A along solution branch D shows how two vortices move from opposite inner corners towards the centre. At the bifurcation point of branch J for $\mu\approx 1.79$ the vortices join together into a giant vortex with multiplicity 2. Note that branch D contains two extra bifurcations, near the ones of branch A. These will be discussed later in this section. Branch E is again similar to branch D: it connects branch A from its bifurcation at $\mu\approx 2.05$ to J at the bifurcation for $\mu\approx 1.48$. Contrary to branch D, the vortices appear at the outer corners in the solutions of branch E, and branch E does not contain other bifurcation points itself.

Branches A and K are connected through F and G. Both branches start at the bifurcation for $\mu\approx 2.20$ of A and lead to the one for $\mu\approx 2.08$ of K. Three vortices move from the corners of the material towards the centre, joining into a giant one with vorticity 3. For branch F one of these vortices appears in an outer corner, the other two at the inner ones, for branch G the opposite happens.

Solution branch H is connected to J by both L and M, both connecting these branches between the bifurcations at $\mu\approx 2.48$ (on branch H) and $\mu\approx 2.01$ (on branch J). For branch L, a vortex appears at an inner corner of the material and moves towards the central one, eventually joining into a giant vortex of multiplicity 2. The same happens in branch M, this time for a vortex appearing at one of the outer corners.

Another connection exists between solution branches H and K. Decrease of the field strength at the first bifurcation near $\mu\approx 2.65$ along branch O shows how two vortices form near the inner corners of the material and start to move towards the central one. At the bifurcation for $\mu\approx 2.64$ the three vortices join into a giant one with vorticity 3. Solution branch N is similar to O and connects branches H and K through the second bifurcation near $\mu\approx 2.65$ (of H) and the one at $\mu\approx 2.62$ (of K). For this branch the vortices appear near the outer corners of the material.

Two final connections appear between branches J and K. Starting from the bifurcation of branch J at $\mu\approx 0.33$ with the cross pattern, increase of the field strength along branch P shows how the zones of low superconductivity start to stretch towards one of the outer corners. Eventually a new vortex moves from this corner towards the centre, appearing there at the bifurcation of branch K for $\mu\approx 0.91$. Branch Q is similar, with the zones of low supercurrent density stretching towards an inner corner.

Both solution branches R and S connect K with its symmetric counterpart for negative $\mu$ values. Starting with a single vortex solution at the bifurcation for $\mu\approx 0.04$, decrease of the field strength along branch S shows how the vortex changes into a zone of low superconductivity stretching from one inner corner towards the opposite one. A similar behaviour is observed for branch R, which starts at the bifurcation for $\mu\approx 0.47$. For this branch the vortex solution changes into a zone of low supercurrent density that stretches between two opposite outer corners.

The final solution branch I is unique in that it is not a main branch, nor does it connect two main branches. It appears between the two extra bifurcation points of branch D. Starting from a solution consisting of two vortices near opposite inner corners at $\mu\approx 1.80$ and increasing $\mu$ along branch I, the vortices start to move in the same direction towards two of the the outer corners. A turning point is encountered near $\mu\approx 1.87$. Decreasing $\mu$ from this value, the vortices first continue to move towards the outer corners, but will not reach them: after a while the vortices start to move back towards each other, eventually reaching another turning point near $\mu\approx 1.68$. Increase of the field strength shows how the vortices both move back towards the centre, eventually arriving at a configuration where the two vortices lie close to each other near the centre of the material, which coincides with the bifurcation of branch D at $\mu\approx 1.79$. Note that the turning points at $\mu\approx 1.87$ and $\mu\approx 1.68$ are the only ones that appear in the bifurcation diagram of figure 6.

Square material

Our final example consists of a material shaped as a square with edges sized 5.5. This problem was also treated in [31], where a bifurcation diagram with 13 different solution branches was constructed. Application of the methods described in the current paper yields the same bifurcation diagram, with an extra 30 branches, resulting in a total of 43 solution curves. The four main branches (with full $D_{4}$ symmetry) are shown in figure 8. Several close-ups of the bifurcation diagram containing other solution curves are given in figure 10. A schematic representation of the diagram can be found in figure 9. Finally, representative solutions for each branch are given in figure 11.

A total of $60$ branch points were discovered during the continuation, their values are given in table 1. Several of these only contained a $1$-dimensional kernel for which algorithm 5 yielded the (single) new tangent direction. Other branch points corresponded to a kernel with dimension $2$. Fore these points the tangent directions were constructed by application of the equivariant branching lemma, described in the last part of section 5.

Note that except for the four main branches (A, D, F and M), solution branch B contains a stable part as well: it is stable for magnetic field strength values between approximately $0.6308$ and $0.646098$, for solutions that consist of two vortices very close to the center, on the horizontal or vertical axis. Branch B is the only non-main branch to contain a stable part, the stability of this region has been verified with a Crank-Nicolson time step method as well.

A thorough description of branches A to M can be found in [31]. Though multiple new solution branches have been discovered in the current paper, their discussion will be omitted.

7 Conclusion

This paper presents multiple tools for the automatic exploration of solution landscapes for the extreme type-II Ginzburg-Landau equation. One of the main challenges is the detection and determination of bifurcation points. We illustrate that this can be handled by tracing Ritz values that indicate a nearby bifurcation point. This triggers a precise determination by an extended system (14) that is solved with a Newton-Krylov algorithm.

The second main challenge is the construction of the tangent directions of emerging solution branches from these bifurcation points. The complexity of this problem depends on the dimension of the Jacobian’s kernel associated with the bifurcation. For $1$-dimensional kernels, algorithm 5 is applied, which is equivalent to solving the algebraic branching equation. When this kernel is $2$-dimensional, we provided an adapted version of algorithm 5, resulting in algorithm 5, to construct the tangent directions. Though this last algorithm is robust, a faster alternative based on the equivariant branching lemma is described as well. This last method can, however, only be applied for $D_{m}$ ($m\geq 4$) invariant bifurcation points.

The methods described in the current paper are implemented in Python as part of the package PyNCT. Contrary to other tools, we make use of sparse linear algebra in the algorithms. The bifurcation diagrams provided in section 6 are generated by our methods, showing their robustness. Though we specifically used the extreme type-II Ginzburg-Landau equations for the derivation of the algorithms, the generalization to other ($D_{m}$- or $C_{m}$-symmetric) problems is straightforward. For the PyNCT extension general algorithms were implemented.

In the future, we aim to apply the methods presented to more complicated superconducting systems. For example, the Ginzburg-Landau equations with space-dependent coefficients or the one that describes the states of type 1.5 superconductors. A similar preconditioner as the one presented in section 3 will likely lead to a fast solver, which can be used in combination with the current algorithms to provide automatically explored bifurcation diagrams.

Further work is required to determine robust stopping criteria for determination of bifurcation points and tangent directions. Although, for the examples presented in section 6, the desired tolerances were always achieved, they require a new tuning when applied to other samples with different discretization levels. Robust stopping criteria, that work automatically for various discretizations, would be helpful in the future.

Finally, an extension towards $3$ dimensional materials is one of the next steps. We expect that the current methods for detection and determination of bifurcation points provided are straightforward applicable. However, the construction of the tangent directions to emerging solution branches will require further efforts. An extra case will probably exist for $3$ dimensional materials, where the Jacobian associated with the bifurcation has a $3$-dimensional kernel. We believe that a similar analysis as for the derivation of algorithm 5 is possible for this case, and the equivariant branching lemma might again be applicable in some cases as well.

Appendix A Proof of algorithm 5

A sketch of the derivation of algorithm 5 is provided in this section. A more detailed one is available at author’s request. First we proof the following lemma:

Lemma A.1**.**

Consider a system of equations with unknowns $\alpha_{1}$ and $\alpha_{2}$:

[TABLE]

with $k\geq 2$, $a_{0},\dots,a_{k}$, $c_{0},\dots,c_{k}$, $e$ $\in\mathbb{R}$ and $b,d\in\mathbb{R}_{0}$. This system contains non-isolated solutions if and only if the relations

[TABLE]

*hold. Moreover, in this case the system has a single isolated solution, given by $(\alpha_{1},\alpha_{2})=(0,0)$. *

Proof A.2**.**

---

Step 1: Suppose that the system of equations contains non-isolated solutions. Since this system can be regarded as the intersection of two algebraic curves, Cramer’s theorem [6] can be applied to show that the curves must be degenerate. The degeneracy implies that*

[TABLE]

Comparing the coefficients of $\alpha_{1}$ and $\alpha_{2}$ in (45) and (46), it follows that $f^{(0)}_{0}\neq 0$. Similar comparisons of other coefficients can be combined with an induction argument to prove the statements

[TABLE]

Using (47), we get

[TABLE]

Next, (48) implies that

[TABLE]

For ease of notation, define

[TABLE]

Since $h^{(k-l)}_{k-l}=0$, we have

[TABLE]

Together with (47),(48) and (49), this implies ($\forall i=1\dots k$)

[TABLE]

Step 2*: We have yet to show that*

[TABLE]

We first rewrite the equations for the algebraic curves. The one for the first curve becomes

[TABLE]

The equation for the second curve can be rewritten as

[TABLE]

Solutions $(\alpha_{1},\alpha_{2})$ of system (43) satisfy

[TABLE]

The equation $\sum_{i=0}^{k-1}c_{i}\alpha_{1}^{i}\alpha_{2}^{k-1-i}+d=0$ again describes an algebraic curve. The solutions that satisfy this equation hence form a continuous curve and are non-isolated. The solution for which

[TABLE]

*satisfies $(\alpha_{1},\alpha_{2})=(0,0)$. This solution is isolated since $(\alpha_{1},\alpha_{2})=(0,0)$ does not belong to the curve described by $\sum_{i=0}^{k-1}c_{i}\alpha_{1}^{i}\alpha_{2}^{k-1-i}+d=0$ (because $d\neq 0$). *

The three equations

[TABLE]

with $k=1\dots l$, are crucial for the derivation of algorithm 5. These equations were derived in section 5. The terms $w_{k}$ and $\beta_{k}$ appear in the order $l$ Taylor expansion of $\psi(s)$ and $\mu(s)$ (16). $\phi_{1}$ and $\phi_{2}$ are null vectors of $\mathcal{J}_{\psi}^{(b)}$, $\phi_{1}^{*}$ and $\phi_{2}^{*}$ of $\mathcal{J}_{\psi}^{(b)^{*}}$. The terms $r_{k}$ are given by the formula

[TABLE]

The terms $v^{(0)}$, $t_{1}$, $t_{2}$, $t_{3}$, $b^{(1)}$, $b^{(2)}$, $y_{i}^{(k)}$, $a_{i}^{(k,1)}$, $a_{i}^{(k,2)}$, $\kappa_{i}^{(k)}$, $z_{i}^{(k)}$ and $q_{i}^{(k)}$ ($k=1\dots l,i=0\dots k$) that appear in the proof are defined as in algorithm 5.

Proof A.3**.**

Initial: Together with $r_{1}=0$, equation (50) implies

[TABLE]

Define $v^{(0)}$ as

[TABLE]

this leads to the following expressions for $w_{1}$ and $x_{1}$:

[TABLE]

Note that the condition $v^{(0)}\in\text{im}\left(\mathcal{J}_{\psi}^{(b)^{*}}\right)$ is necessary to satisfy (51). Formula (53) for $k=2$ yields

[TABLE]

Substitution of $x_{1}$ and bilinearity of the Hessian lead to the expression

[TABLE]

with $y_{0}^{(2)}$, $y_{1}^{(2)}$, $y_{2}^{(2)}$, $t_{1}$, $t_{2}$ and $t_{3}$ defined in algorithm 5. Substitution of this expression in equation (52) leads to the system

[TABLE]

Together with the fundamental theorem of linear algebra, assumptions (23), (24), (25) and (26) imply that we can choose $\phi_{1}^{*}$ and $\phi_{2}^{*}$ such that

[TABLE]

With these choices, (58) becomes

[TABLE]

where $b^{(1)}$, $b^{(2)}$, $b^{(3)}$, $a_{0}^{(2,1)}$, $a_{1}^{(2,1)}$, $a_{2}^{(2,1)}$, $a_{0}^{(2,2)}$, $a_{1}^{(2,2)}$ and $a_{2}^{(2,2)}$ are defined as in algorithm 5. Lemma A.1 can now be applied to check this system for non-isolated solutions. If only isolated solutions exist, the algorithm is stopped. In the other case, we have

[TABLE]

and continue by solving (59) for $\beta_{1}$:

[TABLE]

with $\kappa_{0}^{(1)}$ and $\kappa_{1}^{(1)}$ defined in algorithm 5. Substitution of $\beta_{1}$ in (57) and (56) eventually leads to the following expressions for $r_{2}$ and $x_{1}$:

[TABLE]

with $z_{0}^{(2)}$, $z_{1}^{(2)}$, $z_{2}^{(2)}$, $q_{0}^{(1)}$, $q_{1}^{(1)}$ defined in algorithm 5.

Iteration: Let $k\geq 3$ and assume terms $z_{i}^{(k-1)}\in\mathbb{C}^{n}$ (for $i=0\dots k-1$), $q_{i}^{(j)}\in\mathbb{C}^{n}$ (for $j=1\dots k-2,i=0\dots j$) and $\kappa_{i}^{(j)}\in\mathbb{R}$ (for $j=1\dots k-2,i=0\dots j$) have been derived such that

[TABLE]

Equation (50) implies

[TABLE]

Define $v_{0}^{(k-1)}$…$v_{k-1}^{(k-1)}$ as

[TABLE]

this leads to the following expression for $w_{k-1}$ ($=x_{k-1}$):

[TABLE]

Rewriting formula (53) after substitution of (61), (62) and (64), and defining the terms $y_{0}^{(k)}$,…,$y_{k}^{(k)}$ as in algorithm 5, one can show that

[TABLE]

The coefficients of $\beta_{k-1}\alpha_{1}$ and $\beta_{k-1}\alpha_{2}$ can respectively be rewritten as $t_{1}+2\kappa_{1}^{(1)}t_{3}$ and $t_{2}+2\kappa_{0}^{(1)}t_{3}$. This leads to the following expression for $r_{k}$:

[TABLE]

Substitution of this expression in equation (52) leads to the system

[TABLE]

Remember that

[TABLE]

Equation (66) hence becomes

[TABLE]

where $a_{0}^{(k,1)}$,…, $a_{k}^{(k,1)}$ and $a_{0}^{(k,2)}$,…, $a_{k}^{(k,2)}$ are defined as in algorithm 5. Lemma A.1 can again be applied to check this system for non-isolated solutions. If only isolated solutions exist, the algorithm is stopped. In the other case, we have

[TABLE]

and continue by solving (67) for $\beta_{k-1}$:

[TABLE]

with $\kappa_{0}^{(k-1)}$,…, $\kappa_{k-1}^{(k-1)}$ defined in algorithm 5. Substitution of $\beta_{k-1}$ in (65) and (64) eventually leads to the following expressions for $r_{k}$ and $x_{k-1}$:

[TABLE]

*with $z_{0}^{(k)}$, …, $z_{k}^{(k)}$ and $q_{0}^{(k-1)}$,…, $q_{k-1}^{(k-1)}$ defined in algorithm 5. With the new expressions for $r_{k}$, $x_{k-1}$ and $\beta_{k-1}$ now available, the iteration can be continued for $k\rightarrow k+1$. *

Appendix B Calculation of the $y_{i}^{(k)}$ terms in algorithm 5

An important step in the execution of algorithm 5 is the calculation of the terms

[TABLE]

for a given $k\geq 3$, $0\leq i\leq k$. The terms $q^{(m)}_{l}$, $\kappa^{(m)}_{l}$ and $v^{(k-1)}_{l}$ ($\forall m=1,\dots,k-2$, $l=0,\dots,m$) have been calculated in previous steps of the algorithm. In this appendix we derive an efficient way to calculate these $y_{i}^{(k)}$ terms. For simplicity, we define ($\forall m=1,\dots,k-2$)

[TABLE]

Given a choice for $k$ and $j$, algorithm 1 constructs a list $\mathcal{K}$ of possible index sets $K=(k_{1},\dots,k_{j})$ that satisfy

[TABLE]

The list is constructed in such a way that no two index sets are permutations of one another.

Algorithm 2 is similar to algorithm 1. Given a choice for $j$ and $i$, and a set of indices $K=(k_{1},\dots,k_{j})$, it constructs a list $\mathcal{I}$ of index sets $I=(i_{1},\dots,i_{j})$ that satisfy

[TABLE]

This time the list is constructed such that no two sets of pairs of the form ($(k_{1},i_{1}),\dots,(k_{j},i_{j})$) are permutations.

Given the choice for $k$ and $i$, algorithm 3 describes a possible way to calculate the term $y_{i}^{(k)}$. Algorithms 1 and 2 are used to construct lists of index sets required for the summation of higher-order derivative applications. It is possible to provide a bound on the order of these derivatives to speed up the process.

Acknowledgments

The authors would like to thank Jacques Tempere, Wout Van Alphen, Nico Schlömer, Przemyslaw Klosiewicz and Delphine Draelants for their fruitful discussions.