Non-London electrodynamics in a multiband London model: anisotropy-induced non-localities and multiple magnetic field penetration lengths

The London model describes strongly type-2 superconductors as massive vector field theories, where the magnetic field decays exponentially at the length scale of the London penetration length. This also holds for isotropic multi-band extensions, where the presence of multiple bands merely renormalises the London penetration length. We show that, by contrast, the magnetic properties of anisotropic multi-band London models are not this simple, and the anisotropy leads to the inter-band phase differences becoming coupled to the magnetic field. This results in the magnetic field in such systems having N+1 penetration lengths, where N is the number of field components or bands. That is, in a given direction, the magnetic field decay is described by N+1 modes with different amplitudes and different decay length scales. For certain anisotropies we obtain magnetic modes with complex masses. That means that magnetic field decay is not described by a monotonic exponential increment set by a real penetration length but instead is oscillating. Some of the penetration lengths are shown to diverge away from the superconducting phase transition when the mass of the phase-difference mode vanishes. Finally the anisotropy-driven hybridization of the London mode with the Leggett modes can provide an effectively non-local magnetic response in the nominally local London model. Focusing on the two-component model, we discuss the magnetic field inversion that results from the effective non-locality, both near the surface of the superconductor and around vortices. In the regime where the magnetic field decay becomes non-monotonic, the multiband London superconductor is shown to form weakly-bound states of vortices.


I. INTRODUCTION
The work by F. London and H. London that formulated hydro-magnetostatic theory of superconductivity is one of the most influential works in condensed matter physics. It demonstrated that the magnetic field in a superconductor is described by a massive vector field theory 1 One of the direct consequences is that an externally applied magnetic field decays in a superconductor exponentially at the length scale of the London magnetic field penetration length Λ: B ∝ B 0 e −r/Λ . The superconducting current J varies at the same length scale due to its relation to the magnetic field J = ∇ × B. The work also paved the way to Anderson's demonstration of the Higgs mechanism 2 . Multi-band superconducting materials are currently of central interest. For these systems the pairing of electrons is supposed to take place in several Fermi surfaces, formed due to the overlapping of electronic bands [3][4][5][6][7][8][9] .
The range of validity of London electrodynamics is well understood in single component or single band superconductors and should not, in general, be applicable at the length scale ξ 0 associated with the superconducting carrier, namely the estimated Cooper pair size. In a weakcoupled BCS superconductor ξ 0 can exceed the magnetic field penetration length Λ, leading to markedly different electrodynamics. Such a state cannot be described by the local London model 10,11 and is termed Pippard electrodynamics 12 . In contrast, the limiting case when ξ 0 Λ is described by local London electrodynamics, which is typically applicable to extreme type-2 superconductors.
In this paper we will only be interested in multiband anisotropic materials. Furthermore we focus on the case where ξ (α) 0 Λ, where α is a band index and therefore where one can neglect the effects that lead to Pippard electrodynamics. We will study how the crystal anisotropy of superconducting materials affects their magnetic properties even at the level of London's model. In a single component case, the anisotropy effects in London electrodynamics are well studied [13][14][15] . In the single component case the anisotropy leads, in general, to an electrodynamic kernel that is characterized by two realvalued penetration lengths, corresponding to the different polarizations of the magnetic field. We show that the situation is principally different in the multiband London model due to the presence of additional massive modes, associated with the variations of the phase differences between order parameter components, known as the Leggett modes 16 . In the isotropic case the magnetic and Leggett modes are de-coupled, however we show that, in general, this coupling appears in the London model with the introduction of crystal anisotropy. On the qualitative level, this coupling arises due to differing anisotropy in the superconducting bands which enables the gradients of the inter-band phase differences to produce non-zero transverse charge currents, which generate magnetic field. The London model then exhibits N + 1 magnetic modes, where N is the number of superconducting components. Namely, the two independent components of the magnetic field and N − 1 interband phase differences, yielding in total N + 1 massive scalar fields 1745 We analyse this behaviour in detail using the minimal model of the two-component superconductor. The results can be straightforwardly generalized to the larger number of bands, e.g. the three-component model with phase frustrations yielding the intrinsically complex s + is state. This opens interesting possibilities to study the behaviour of magnetic signatures of the broken-time reversal symmetry phase transition to the s + is characterized by the soft Leggett mode which becomes massless at the transition point 18,19 . One of the interesting properties is that the masses and correspondingly the relaxation lengths can become complex under certain conditions. The unconventional behaviour of magnetic modes will be shown to result in magnetic field reversal near the boundaries and vortices leading to the formation of vortex bound states.

II. THE MODEL
It is illustrative to view the London model as a constant density limit of the Ginzburg-Landau model for complex fields (although indeed the model is much more general and is valid at low temperatures). The simplest Ginzburg-Landau free energy density for an Ncomponent anisotropic system is given by where D = i∇ + 2πeA/c is the covariant derivative, e and c are the electron charge and light velocity respectively (hereafter we use the units with = 1), the fields ψ α = |ψ α | e iθα represent the different superconducting components. Greek indices will always be used to denote superconducting components and Latin indices will be spatial with the summation principle applied for repeated Latin indices. The inverse mass tensors λ −1 ijα represents a 3 dimensional diagonal matrix for each component, In Eq.(2) F p collects together the potential (nongradient) terms which can be any from a large range. The simplest example is the standard single band potential terms and the Josephson inter-band coupling term where ψ 0 α , γ α and η αβ are positive real constants. The second term above is the Josephson inter-band coupling, where θ αβ = θ α − θ β is the inter-band phase difference between components α and β. The simplest multicomponent system is the isotropic case λ α = I 3 without inter-band coupling η αβ = 0, which has the maximal symmetry U (1) × U (1)... = U (1) N . However the introduction of the simplest non-frustrated Josephson terms, that in the ground state locks all phase differences to zero, breaks this symmetry to a U (1) symmetry. In the isotropic case, the phase differences are neutral massive modes. In the absence of a coupling to the gauge field the phase sum gives rise to the massless Goldstone mode.
In the presence of a gauge field coupling that mode is the massive London mode.
As we are interested in strongly type 2 systems we neglect density variations, but, following Leggett, we retain the massive degree of freedom associated with the phase difference mode 1646 We consider the limit |ψ α | 2 ≈ const, leading to the following free energy (for brevity below we refer to the approximation that neglects density variations as the London limit).
Here j α are the partial superconducting currents where Φ 0 = πc/e is the flux quantum,λ k are coefficients characterizing the contribution of each band to the Meissner screening, A is the vector potential and J αβ the Josephson coupling. Parameters of the London model (5) are related to that of the Ginzburg-Landau functional Eq.(2) by J αβ = η αβ |ψ α ||ψ β | andλ −2 α = (2e|ψ α |/c) 2 λ −2 α , where the matrix indices are suppressed forλ.
Considering the total current, which is a sum of partial current contributions from each band j = α j α we get an expression for the magnetic field, which provides an extension of the London theory in multicomponent systems 20,21 where N is the number of components,λ 2 L = ( αλ −2 α ) −1 and θ αβ = θ α − θ β are the relative inter-band phases. Expression (7) shows that the phase difference gradients in the second term can generate magnetic field in anisotropic materials. We will demonstrate that the coupling between magnetic field and interband phase differences leads to the non-local magnetic response. Importantly the non-locality here has nothing to do with the Pippards non-locality 12 associated with the Cooper pair dimension 10,11 . Rather it can be obtained already within the the standard anisotropic London model, which has only "local" terms. In result of that, the magnetic response of such a system is characterized by the multiple magnetic modes with different penetration lengths which depend on the degree of anisotropy in different bands and the strength of inter-band pairing interaction. In next sections of the paper we will consider the influence of these anisotropic effects on the magnetic response of multi-band superconductors and demonstrate that they many observable physical consequences both for the Meissner and vortex states.
To connect the general model given by the Eqs. (5,6,7) with real materials we can use the multiband weakpairing pairing models and express the Josephson couplings J αβ through the microscopic parameters, such as the pairing coefficients V ij , Fermi velocities in different bands and the lengthsλ k characterising magnetic responses of each band in Eq. (6). For the generic case of a two-band superconductor there is only one Josephson energy scale E J = J 12 , which determines the Leggett mode frequency 16 . Since we are dealing with a static but spatially inhomogeneous problem and aim to describe the coupling between the phase difference and magnetic field, it is natural to construct the inverse characteristic length scale as follows k 0 = [2π 2 E J /Φ 2 0 ] 1/4 . This parameter does not depend on the anisotropy or condensate stiffness. Therefore it can be expressed through the weakcoupling pairing coefficients where ∆ i is the superconducting gap amplitude for the ith band. The density of states in the i-th band ν i can be found using the magnetic response length (λ −2 where v F α is the α-th spatial component of the Fermi velocity, and ... F S denotes the Fermi surface average. For the parameters 22,23 of the uniaxial anisotropic two-band superconductor MgB 2 we get k 0 λ 1⊥ = 1.3 for the σ-band with larger gap and k 0 λ 2⊥ = 0.9 for the π-band with smaller gap. Although the London length anisotropy in MgB 2 is quite weak 24 , it has a rather pronounced anisotropy of the partial magnetic responses in the almost cylindrical σ-band. According to the expression above it is determined by the Fermi velocity anisotropy v F ⊥ /v F z ≈ 8.6 25 . Therefore we can estimate λ 2z ≈ λ 2⊥ and λ 1z ≈ 8.6λ 1⊥ .
The other example of multiband anisotropic superconductor is Sr 2 RuO 4 . The nature of the superconducting state for this material is still highly debated. In order to make one more estimate, we consider a weak-coupling three-band model from Ref. 26 . The suggested coupling matrixV for Sr 2 RuO 4 contains two bands with strong interband interactions and the third band which has much weaker interband pairing, such that the elements V 13 and V 23 are much smaller than the others. Therefore one can use an effective two-band model to describe the lowenergy and large-scale variations of the interband phase differences θ 13 = θ 23 and assuming θ 12 = 0. Additionally, we use the same values of gaps, and lengths λ i⊥ in all bands 27 to obtain an estimation of k 0 λ ⊥i = 2.8. The overall London length anisotropy in Sr 2 RuO 4 is about λ Lz /λ L⊥ ≈ 20, although the anisotropy of each band contribution is not known. Therefore, in realistic multiband compounds k 0 can be of the order of the penetration length and the anisotropy ofλ k in each band can change in the wide limits.

III. ELECTROMAGNETIC RESPONSE
In this section we consider the system in the absence of vortices and demonstrate that the anisotropy qualitatively changes the electromagnetic modes in multiband superconductors. To obtain the equations of motion we rewrite the condensate phases in Eq.(5) as θ α = β θ αβ /N + θ Σ where θ Σ = α θ α /N . Let us assume for simplicity the same strength of Josepson coupling between all bands J αβ = E J . Then varying the free energy (5) by θ αβ and A we obtain the system of coupled equations for the phase differences and magnetic field where we have introduced the gauge invariant term p s = ∇θ Σ − 2πA/Φ 0 . In the absence of phase singularities we can choose the gauge so that the common phase is constant α ∇θ α = 0 and therefore p s = −2πA/Φ 0 . For isotropic superconductors where the total current is j ∝ p s this choice corresponds to the London gauge. In the anisotropic case the relation between p s and the current is more complicated so that the gauge in general depends on the specific choice of the anisotropy parameters.
To study the linear electromagnetic response we linearise the above equations of motion and switch to the momentum representation. Then we get the algebraic relation between the current and vector potential where Q is known as the polarisation operator and tells us how changes to the gauge field relate to the current. It is given to be, .
(12) During the derivation we used the commutatorλ −2 In the isotropic limit the kernel becomes local due to the second term in Eq.(13) becoming zero, through the London gauge choice k · A = 0. In general one can see that the non-locality scale of Q(k) is determined by the inter-band Josephson length k −1 0 which is not related to the BCS non-locality scale determined by the Cooper pair size. In the absence of coupling between θ αβ and B we obtain the usual local response, similar to the singlecomponent superconductors determined by the constant tensorλ L .
The complicated structure of the second term in Q(k) points to some unusual magnetic properties of anisotropic multi-band superconductors. In particular, that it produces multiple magnetic modes in which the magnetic field and the inter-band phase difference are coupled. This is what has lead us to calling them magnetic phase difference modes.

IV. MAGNETIC MODES AND FIELD SCREENING IN THE TWO-BAND SUPERCONDUCTOR
Let us consider the magnetic response of a two-band anisotropic superconductor. A good example of such a kind of material is the MgB 2 compound 28 which has two superconducting Fermi surfaces with the structure qualitatively similar to the one shown in Fig.(1)a. The twoband polarization operator (12) with N = 2 becomes . case the amplitude of θ 12 is small so that the condensate phases are effectively "glued together". In the second case θ 12 can be large due to negligible inter-band Josephson coupling.
The coupling between θ 12 and magnetic field also leads to especially important consequences when k 4 0 ∼ (k ·λ −2 L k). In this case the magnetic response has multiple length scales. It means that the magnetic field penetration into such a superconductor is determined by the superposition of several fundamental modes, which can be found from Eq. 9, by linearising the Josephson term near the ground state value θ 12 = 0 and searching for solutions in the form of plain waves B, θ 12 ∼ e ik·r . Thus we obtain the linear system where we denote h = 2πB/Φ 0 . This system is of the sixth order, since magnetic field has only two independent components k · h = 0. Hence in general for each direction of k there exists three different solutions with Imk > 0. Therefore the system (14) can not be solved analytically. However as we will see below, we can ascertain properties and even analytical solutions for certain symmetries, including the most realistic model of uniaxial anisotropy. Assume that the x and y axes are equivalent λ αx = λ αy and z is the anisotropy axis. In this case the general system (14,15) splits into the second order and the fourth order equations which determine the usual and unconventional magnetic modes respectively discussed below. Let us consider first the magnetic field with polarization B = B o coplanar with the anisotropy axis and the wave vector,ẑ and k as shown in Fig.(1). In this case we get k = ±iλ −1 Lx and the magnetic field is decoupled from the phase θ 12 = 0. We call this magnetic mode the ordinary one. Alternately, let us consider the extraordinary modes B = B e , which are coupled to the inter-band phase. The magnetic field is perpendicular to both k andẑ as shown in Fig.(1). In this case Eqs. (14,15) yield a finite coupling between the magnetic field and the inter-band phase difference.
To analyse the extraordinary mode in detail we parametrize the wave vector components as k z = k sin θ, k x = k cos θ sin ϕ and k y = k cos θ sin ϕ where k ⊥ = k cos θ. The direction ϕ drops out from the equations due to the rotational symmetry in xy plane. Then we get the following relation between the interband phase and magnetic field amplitudes Note that in the isotropic case with λ 1⊥ = λ 1z and λ L⊥ = λ Lz Eq.(16) yields h = 0 so that this mode becomes the non-magnetic pure phased-difference excitation. The wavenumber of the extraordinary mode is then given by the following bi-quadratic equation It has two complex solutions with Imk > 0 yielding two extraordinary magnetic modes with the polarization shown schematically in Fig Here we assume that the first band is isotropic λ 1⊥ = λ 1z = λ. The second band has either a weak anisotropy with λ 2z = 0.8λ 2⊥ (Fig.3a) or a strong one with λ 2z = 0.1λ 2⊥ (Fig.3b) and λ 2⊥ = λ. As shown in the Fig.(3) for large and small Josephson couplings, the wavenumbers of the two modes are quite different. One of them is proportional to k 2 0 and hence either diverges or goes to zero at k 0 → ∞ and k 0 → 0 respectively. At the same time the other one tends to the constant values corresponding to the local response approximation discussed above.
The general solution of Eq.(17) reads where, There are a few interesting limiting cases for k 1,2 . First, the strong inter-band coupling k 0 λ −1 i⊥ , λ −1 iz leads to b 2 4ac and hence gives purely imaginary solutions, The weak inter-band coupling k 0 λ −1 i⊥ , λ −1 iz leads again to b 2 4ac giving slightly different purely imaginary solutions, The imaginary wavenumbers obtained in the limits considered above correspond to the real-valued masses or the inverse decay length-scales of the magnetic modes. A completely different regime is possible when the masses of magnetic modes become complex, resulting in damped oscillating behaviour of the magnetic field. Indeed in the range of parameters when b 2 < 4ac the solutions (18) have finite real parts, the examples of such solutions are shown in Fig.(3) where Rek 1,2 = 0 for an interval of k 0 which expands with increasing degree of anisotropy. Generically, this regime can be realised when a strong anisotropy is applied in each band in different directions, . Then in the wide rang of we obtain the wavenumber k 1,2 = (i ± 1)k 0 / √ 2 which has the amplitudes of real and imaginary parts.

A.
Flux expulsion in the Anisotropic two-band model We now consider the problem of magnetic field screening at the surface of an anisotropic multiband superconductor. Let us consider cylindrical geometry with magnetic field applied in the y direction B = H 0ŷ parallel to the boundary of the superconducting sample as shown in Fig.(2)a. Unlike the usual isotropic superconductors, the boundary orientation with respect to the crystal axes is important and is considered by the angle θ in polar coordinates, introduced previously in considering the normal modes. The wave vectors of excited magnetic modes are directed perpendicular to the surface k = k(cos θ, 0, sin θ). We use r to indicate the coordinate orthogonal to the boundary, such that in the presence of two magnetic field penetration lengths, the magnetic field decay and phase difference follow a double-exponential law, are two magnetic field penetration lengths. We have previously made the assumption that θ 12 is small and hence linear in nature. This does not necessarily have to be the case and indeed removing this assumption could lead to more unconventional physics e.g. oscillation of the inter-band phase difference parallel to the boundary. To discover if this has an effect one would have to consider the full equations numerically. For our 1-dimensional equation, the boundary conditions are given by the value of the magnetic field on the boundary due to the external field B 1 + B 2 = H 0 and the requirement for the normal current to vanish at the boundary n · j = 0. The contribution from each of the magnetic modes (24) to the normal current can be found as n · j α = −i∇ · j α /k α = −ieE J θ 12 /k α . The relation between phase difference and magnetic field in each of the mode is given by the Eq.(16).
Then we get the following solution for the amplitudes B 1,2 in Eq. (24): Example solutions for various parameters are plotted in Figs. 5 and 6. The two key features that differ from the isotropic case is the double-exponential decay of the magnetic field and a self-induced gradient of the phase difference between the superconducting components. The origin of both effects is the hybridization of the Leggett mode with the magnetic mode. Also note that the solution matches the four-fold symmetry of the free energy, which is to be expected due to the anisotropy. Additionally, the limiting cases of k 0 → ∞ and the isotropic case decouples the phase difference and magnetic field and we are left with a single penetration length as expected.
An interesting limit to consider is where the Leggett mode becomes massless, for example the zero Josephson coupling limit (k 0 → 0). As the Josephson coupling becomes smaller one of the magnetic field penetration lengths diverges Λ 1 → ∞. The physical consequence of a diverging magnetic field penetration length in our-twoscale system is markedly different from the divergence of magnetic field penetration length at T c in a singlecomponent system; namely in our case it doesn't imply absence of magnetic screening. This is due to the amplitude of the mode vanishing as the penetration length diverges (B 1 → 0). Simultaneously the other mode be-comes the isotropic case deformed by the anisotropy, as one would expect if the anisotropy in the different components matched and could hence be rescaled.
This implies that in such a limit, most of the magnetic field's amplitude decays at the length scale Λ 2 along with an increasingly long-range penetration of a small "tail" of magnetic field. The situation should occur for example in multi-band systems close to s + is transitions where the Leggett mode becomes massless [18][19][20]29 .
Finally, and rather interestingly, the multi-mode magnetic response implies that the magnetic field is not necessarily monotonic. This non monotonic behaviour can lead to field inversion for a range of parameters. This means interactions between vortices and boundaries and also inter-vortex interactions will be non-trivial. It is likely that should the external field be increased such that vortices enter into the sample, their position will be affected by the negative magnetic field which would attract the vortices. Again it should be emphasised that the field inversion here is not related to oscillatory behaviour of the magnetic field in the non-local Pippard's model.
This field inversion can be seen most clearly in the plots of the solutions for various parameters below. In   Fig. 28 in the appendix we have plotted anisotropy in a single component for increasing strength. It is clear that the strength of the field inversion increases as the anisotropy is increased. Additionally the generated phase difference is also more pronounced. In Fig.6 we have looked at strong anisotropy in different directions in each component and in Fig.29 in the appendix we have plotted anisotropy in opposite directions for increasing strength. Again as in the other plot the field inversion and phase difference are increased in magnitude at the anisotropy amount increases. Due to the choice of parameters however, we observe a far more symmetrical dihedral solution.
Finally we have considered the negative magnetic field at θ = π 4 for various strengths of k 0 , shown in Fig.4. We can observe here the effect of decreasing the Josephson coupling strength as discussed above. This leads to one of the modes becoming weaker (and also the strength of the negative magnetic field becoming weaker) but also long range.

B. Numerical solution for the boundary problem
We now perform a numerical simulation of the Meissner state of the two-band anisotropic full Ginzburg-Landau model (2). This was performed using the FreeFem++ numerical library 30,31 , which utilises a finite element space, over which conjugate gradient flow is performed. The simulations were performed on a disc as the problem is directionally dependent. We minimise the Gibbs free energy G = R 3 F − R 3 B · H + R 2 F surf ace , where H = H z e z in an external field, applied orthogonal to our 2D system, where we have chosen a finite domain with the boundary conditions ∇ × A = H and for the current to vanish n · j = 0. We then slowly increase the external field strength |H z | in steps of 10 −3 , through the various Meissner states. When the external field is below the 1st critical value we get the Meissner state solutions shown in Fig. 7, 8 and 9. Note all simulations were run with ψ 0 α = e = = c = 1 and γ α = 10, so in the strong type 2 regime, hence our results should be comparable with the London model above.
In Fig. 7 we see a small sample in the Meissner state, with parameters λ x1 = λ y2 = 1, λ y1 = λ x2 = 10 and η αβ = 0.5. The key effect to note is the oscillation of the phase difference away from the axis: a consequence of the anisotropy driven hybridization of the Leggett mode and magnetic mode discussed above. When the angle of the boundary is away from θ = nπ/2. When looking at the magnetic field however we don't see the expected oscillation. This is due to the effect being long range, so if we increase the radius of the disc for similar parameters, as shown in Fig. 8, we can see the inversion of magnetic field.
We also have depicted the Meissner state for anisotropy in only one band in Fig. 9 for the parameters λ x1 = λ y2 = λ x2 = 1, λ y1 = 10 and η αβ = 0.5. This leads to a similar effect on the shape of the various plots, matching the symmetry of the anisotropy. However the negative component of the magnetic field is far smaller, which is as we predicted when we considered the limiting cases in the London model. Finally switching off the Josephson term removes the magnetic field inversion as was discussed in the London approximations. So the results from the full field numerics seem to qualitatively support the London model calculations and all the predictions that came from them.

V. VORTEX STRUCTURE
Let us now consider how the multiple magnetic modes in anisotropic multi-band superconductors modify the vortex states in the London model. We begin by studying vortex solutions carrying a single flux quantum, where both components have 2π phase winding around the core. We later consider a solution for fractional flux vortex i.e. the vortex that has phase winding only in one phase and carries a fraction of flux quantum (for details of flux fractionalization and energetical preference of fractional and composite vortices in multiband systems see the discussions for isotropic systems 32,33 ).

A. Integer-flux vortex
The field distribution around a single vortex line can be found in the form of a Fourier transform, where r is a coordinate vector in the plane perpendicular to the vortex line and the 2D integration is done by the corresponding momentum space cross section. The components h(k) are now determined by the nonhomogeneous system where n v is the direction of the vortex line. The anisotropy in the plane perpendicular to the vortex line can be caused by two reasons. (i) When the magnetic field is directed along the c-axis, there can be anisotropy in the ab plane either if the crystal is biaxial or as a result of the strain-induced distortions. (ii) The effective anisotropy can be caused by a misalignment be- tween the external field and the anisotropy c axis. The distribution of the magnetic field around vortices will be different in these two cases since in (i) only two of the magnetic modes are excited and in (ii) the amplitudes of all three magnetic modes are non-zero.
In this paper we consider in detail only the case (i) when the magnetic field around the vortex is h = h z z, where Here the poles k 1,2 are given by the Eq.(18). Using expressions (31,28) one can consider the asymptotics of the field far from the vortex center. We introduce the polar coordinates k = k(cos θ, sin θ, 0) and integrate first by k taking into account the symmetry h z (k) = h z (−k) which allows the integration to be extended to the domain k < 0. Using Eq.(31) we get the magnetic field distribution in the real space polar coordinates (r, ϕ) with the origin at the vortex center where h j (ϕ) = Res(h z (k, ϕ), ik j (ϕ)) is a residue of the function h z (k, ϕ) at the pole k = ik j (ϕ). For the angle integration we used an approximation k j r cos(θ − ϕ) ≈ k j r(1 − (θ − ϕ) 2 /2) which is valid provided k j r 1. In general due to the ab-plane anisotropy Eqs. (33,34)  yield the expected fourfold magnetic field profile around a vortex. But the most important point is the nonmonotonic field behaviour with field inversion at some distance from the vortex center, as with the boundary problems. For this it is necessary and sufficient to satisfy two conditions: h 2 < 0 and k 1 > k 2 , so that the mode with negative amplitude can become dominating at some distance from the vortex.
To demonstrate the possibility of field inversion we consider the parameters λ 1x = 0.7λ, λ 1y = 0.8λ, λ 2x = λ 2y = λ and k 0 = 0.7/λ which results in the correct fourfold magnetic field profile with field inversion far from the vortex core. The negative parts of field distribution B z −|B z | are shown in Fig.(10). Field inversion leads to a non-trivial four-fold non-monotonic interaction between vortices and thus bound states. This will be considered in more detail below in Sec.VI.
Finally the above modifications to vortex solutions point towards a new mechanism for vortex stripe formation in multi-band superconductors.

B. Fractional vortex solution
For fractional vortices we have a singularity appearing in just one of the components (chosen to be component 1 without loss of generality). It follows immediately that this will lead to a singularity appearing in θ 12 itself unlike in the composite case. This means the previous approach of linearising the equations and approximating sin θ 12 ≈ θ 12 is no longer valid. This is unsurprising as it is due to the presence of Josephson strings when the coupling is switched on. Here it is illustrative to consider the fractional vortex in the U (1) × U (1) model.
We follow a similar procedure by Fourier transforming the equations of motion, however the singularity now exists in ∇θ 1 , allowing us to utilise the fact that p s −∇θ 12 /2 is singularity-free. The equations of motion can be rearranged to isolate the singularity and yield in the Fourier representation the following, If the magnetic field lies in the z-direction only h = h zẑ the solution of this system is given by, where we have used k = k(cos θ, sin θ, 0) again. Using a similar method to before we can find the magnetic field to be, This illustrates that when interband Josephson coupling is zero k 0 = 0, the Leggett mode becomes non-magnetic and does not contribute to the magnetic response despite gradients of the phase difference. 47 The solution is plotted for a number of values in Fig. 11.

C. Numerical Vortex Solutions
We now consider the numerical solutions for type 2 vortices in the full anisotropic Ginzburg-Landau equations. All numerical simulations were performed using the FreeFem++ numerical library 30,31 , which utilises a finite element space over which conjugate gradient flow is performed. We take expression (2) to be our energy functional, where all the solutions were found in a type 2 superconductor with parameters γ α = 2 and ψ 0 α = e = = c = 1 for all components. We also restrict to the x-y 2-dimensional plane, where we can now refer to different solutions in terms of the number of flux quanta, We first consider the single-quantum vortex solution with weak anisotropy in one band λ 1x = λ 1y = λ 2x , λ 2y = 2λ shown in Fig. 12, and strong anisotropy in one band λ 1x = λ 1y = λ 2x , λ 2y = 10λ shown in Fig. 13. The salient point is the confirmation of magnetic field inversion for both sets of parameters, increasing in strength as the anisotropy is increased. This along with the self induced phase difference gradients confirm that the results from the above London model calculations, do transition into the full Ginzburg-Landau model as expected. Another thing to note is that the symmetry of the solutions match the broken symmetry of the energy functional as with the London model. The isotropic component appears to retain it's axial symmetry from an unbroken model, while the anisotropic components symmetry is broken to the expected 4 fold symmetry, or squashed as one might expect of a single component system.
We now consider the effect of having anisotropy exhibited in both bands in a similar direction as shown in Fig. 14. As the anisotropies in the two band approach each other the field inversion gets weaker and once they are the same the model is in many respects analogous to a single component anisotropic system. The closer the anisotropies are to each other in each of the bands the more diminished the exotic behaviour we have observed becomes and eventually vanishes as the anisotropies coincide.
For anisotropy in different directions however, we see a more pronounced effect by the anisotropy, as can be seen in Fig. 15 for parameters λ x1 = λ y2 = λ, λ y1 = λ x2 = 10λ, η 12 = 0.5. When compared to a similar solution with anisotropy in a single band in Fig. 13 the magnetic field inversion is more notable and the deformation of the shape while retaining the D 4 symmetry of the free energy is more deformed.

D. Field inversion beyond the London limit
We now consider the effect of altering the Josephson coupling strength. If we increase the coupling to become very strong it leads to a diminished disparity of the magnetic field penetration lengths. Thus the magnetic field becomes more localized and the inversion less pronounced. If we take the Josephson coupling to be of a similar scale as the covariant derivative pre-factors, then the anisotropic effects are maximal, as predicted in the London model. Finally taking the Josephson strength to be small we observe the magnetic field becoming continually longer range and the field inversion less pronounced.
In the limit of zero Josephson coupling the London model predicts that the Leggett mode decouples from the magnetic field. That is, as η 12 → 0 one of the modes becomes zero and the magnetic field decay is described by a single exponential, hence the London model predicts absence of magnetic field inversion. However this is not what is found in the full model as can be seen in Fig. 16 for the parameters γ 1 = γ 2 = 2, η 12 = 0.0, λ x1 = λ y2 = λ and λ y1 = λ x2 = 10λ. Here we see that the magnetic field still exhibits inversion, though at very long range. This result refutes the applicability of the London model for this regime.
The origin of this behaviour must be the interplay of anisotropy with an additional effect that appears beyond the London model. It has been proposed that in the isotropic Ginzburg-Landau model, the magnetic response has the form of a massive vector field coupled to Faddeev-Skyrme terms 34,35 . The Faddeev-Skyrme terms represent magnetic field contribution that is generated by cross-gradients of relative densities and relative phases of components. This leads to magnetic field inversion in fractional vortex solutions 36 , however in the standard isotropic model, for axially symmetric integer flux vortices, solutions do not have gradients of relative phases and the effect is absent. Hence in the anisotropic case, the self-induced phase differences and relative density gradients lead to magnetic field inversion even in the case of zero Josephson coupling, beyond the London limit. For finite Josephson coupling this effect coexists with the contributions discussed above in the London model.

VI. VORTEX BOUND STATES
In this section we consider inter-vortex interactions, that are likely to lead to non-trivial multi-solitons, due to the non-monotonic nature of the magnetic field and the property of field inversion. If we return to the London-Leggett model energy formulated in Eq. (5), we can expand with respect to the key terms θ 12 , p s and B, We now want to find the interaction energy of two composite vortices with winding in both components in the London model. We assume that they are well separated, such that we can write the various terms as the sums of the tail interactions of the two solitons, B = B (1) +B (2) , s . If we then separate and integrate equation 41 by parts, we can use the equations of motion to reduce the interaction energy to the simple form, Hence, by substituting in the form of the parallel field around a single composite vortex given in Eq. (33), we can calculate the total interaction energy of a given configuration. For the system in question, the only required data to represent a given configuration is then the positions of each individual composite vortex or a collection of 2n parameters, where n is the total number of quanta or winding number of the system. We can now minimise the interaction energy for a given winding number over the 2n-dim space of positions, to find the optimal configurations for type 2 composite vortices in the London model. Eq.(42) was minimised using a simulated annealing method.
We first consider the solutions for equal and opposite anisotropies for various strengths of anisotropy, some of which are plotted in Fig.s 17, 18 and 19. If we consider the solution for weaker anisotropies, presented in Fig. 17 for λ 1x = λ 2y = 0.5λ, λ 2x = λ 1y = λ and k 0 = 0.84/λ, we see that the form of the minimal energy solutions is that of polyominoes (geometric plane  Increasing the strength of the anisotropy to λ 1x = λ 2y = 0.3λ, λ 2x = λ 1y = λ and k 0 = 0.84/λ the polyominoes pattern no longer applies and the form of the minimal energy solutions are more complex, as can be seen in Fig. 18. Up to n = 4 we see similar solutions to before, but for n = 5 we see a rotated N = 4 solutions with an additional vortex in the centre of the configuration. As the winding number increases a pattern emerges, that of chains that are just off the π 4 diagonal, interlaced with each other, such that the chains are staggered.
Finally if we increase the strength of anisotropy to extremely strong values, we get a continuation of the above to more extreme behaviours. This can be seen in Fig. 19 for λ 1x = λ 2y = 0.1λ, λ 2x = λ 1y = λ and k 0 = 0.84/λ. We now see the above shifting away from the polyominoe form for n ≥ 3.
If we now consider anisotropy in a single direction we get a very different result, as shown in Fig. 20 for λ 1x = 0.3λ, λ 1y = λ 2x = λ 2y = λ and k 0 = 0.84/λ. Here we see that the form of the solutions is now of chains. As the winding number increases the chains develop kinks, looking at the solution for n = 4 this is due to the additional vortex being far enough away from the one directly above it to interact weakly, but close enough to the one at the tip of the chain to be affected by the negative magnetic field which is longer range. The chains form on a line with an angle to the x-axis determined by the form of the anisotropy. As the winding number increases further, the interlacing effect of the chains appears as before.
If we were to consider more complicated forms of anisotropy it is likely the minimal energy solution will take some hybrid of the two presented above based upon the how warped the symmetry is away from the maximal  We can now compare this to the results of the Ginzburg-Landau field theory. We naturally start with the 2 quanta N = 2 solutions, shown in Fig. 21 for anisotropy in one band and Fig. 22 for anisotropy in both bands in opposite directions, demonstrating that the bound states do exist in the full model also. The direction of separation also matches that predicted by the above London model and is based on the parameters of the model. We can find the energies of this formation by simulating the 2 quanta, 1 quanta and vacuum solution on the same grid to compare energies. For the parameters γ 1 = γ 2 = 2, η 12 = 0.5, λ x1 = λ y1 = λ x2 = λ and λ y2 = 10λ shown in Fig. 21 we get E 1 = E 1 − E 0 = 1.6416 and E 2 = E 2 − E 0 = 3.2773, where each E i is the minimal energy solution for the i quanta system. Which means the 2 quanta normalised energy is lower than the single vortex normalised energy per vortex E 2 < 2E 1 and a bound state has been formed. The binding energy of this bound state is small which is no surprise, due to the small levels of the magnetic field inversion compared to other terms in the free energy.
We now move onto the n = 3 quanta solutions, for equal anisotropies in opposite directions we are interested in stronger and weaker anisotropies and the effect on the minimal energy solutions. For weaker anisotropy, shown in Fig. 23, we observe 3 key local solutions in the form of a line, an L shape and a T shape. The first two of these are polyominoes, the second is not. If we compare the energies we get the following, E line = 4.7652, E L = 4.7678 and E T = 4.7689, which gives the line solution as the global minima, as predicted by the London model. It is instructive however to compare the energies of the other two local minima, as the London model predicts polyominoes being more favoured over interlaced solutions, thus one would expect the L solution to have lower energy than the T solution, which is the case.
For stronger anisotropy, shown in Fig. 24, we can see the same local solutions, however the energies are now given as E line = 3.1098, E L = 3.1199 and E T = 3.1110, such that the line solution is still the global minima. Here we see that Ginzburg-Landau solutions do not agree with the London model. However the energy difference has decreased between the various solutions so the trend is the same. This is likely due to the London model not allowing the form of the individual vortices to deform. Note that the minimal value when comparing the L and T solutions has now switched to the T solution. This means that the interlaced solutions are more favoured now, as predicted by the London model above.
We finally show the local solutions for N = 4 for equal anisotropies in the bands in Figs. 27 and 26. The global minima for weak anisotropy is the line however the square solution is extremely close in energy, as the London model predicted. For stronger anisotropy we observe the square being rotated and deformed, as predicted by the London model. It is now the Z solution that is close to the line solution, which is again as the London model predicted.
Finally, for anisotropy in one direction we observe sim-ilar results to the London model, with chain solutions taking the minimal energy solutions as expected. Various local and global minima solutions are plotted for N = 3, 4 in Fig. 25. Note that the energy of the L solution approaches that of the line solution, due to the deformation from the line discussed in the London model. The other local solutions have much higher energy compared to the L and line solutions. This means that for the type 2 regime, the London model is good at predicting the qualitative form for the higher quanta solutions. This is despite the fact that when one goes beyond the London limit, there appear additional terms that contribute to field inversion 36 . Hence it would be likely to continue to take the various forms predicted for higher quanta above. A different regime, that appears in the Ginzburg-Landau theory and is not captured by a London model, is where coherence lengths exceed magnetic field penetration lengths. In that regime vortex bound states form via a different mechanism (see discussion for the isotropic case in [38][39][40][41] Addition of anisotropy to these regimes also leads to anisotropic vortex cluster solutions and vortex chains 42 .

VII. CONCLUSIONS
The magnetic response of isotropic single-component superconductors can be cast in the form of a massive vector field theory characterized by the London magnetic field penetration length. However superconducting materials are often multiband and anisotropic. We have demonstrated that this leads to deviation of the magnetic properties of the model from London's hydromagnetostatics even at the level of the anisotropic multiband London model. We showed that anisotropy leads to hybridization of the Leggett and the London modes, causing the gradients of the phase difference to create transverse charge currents, generating magnetic field. This in turn leads to the existence of several magnetic field penetration lengths (in general N + 1 magnetic field penetration lengths for an N -band London model) and also to a non-local magnetic response in the nominally local London model. For example in the case of the twoband anistropic Meissner state, the magnetic field, directed along one of the crystal axes, decays according to a double-exponential law. In the general case of arbitrary directions of such two-band case, there are three magnetic modes with different penetration lengths. In the limit of vanishingly small mass for the Leggett mode (e.g. near s + is transition, one of the magnetic field penetration lengths diverged, leading to long-range, smallamplitude penetration of magnetic field even far below the superconducting phase transition.
Under certain conditions the magnetic field is described by a massive vector field theory with complex mass, which means that magnetic field decay cannot be entirely characterized by a real length scales but has oscillating behaviour. Moreover the combination of dif- ferent magnetic modes gives the overall magnetic field profile a non-monotonic form and even magnetic field inversion. This affects the nature of vortex states: the non-monotonic behaviour and field inversion leads to the formation of vortex bound states. The minimal energy bound states were shown to depend on the symmetry of the system and have forms of polyominoe vortex clusters and chains.
A number of multiband superconductors are currently the subject of detailed experimental research. The examples studied in this paper give inverted magnetic field up to 10 −3 of applied magnetic field. Such field strength makes the effect in principle measurable either by SQUIDs or in muon spin rotation experiments. The interband coupling strength is often difficult to calculate precisely. Measuring the effect that we report for samples with different boundaries, cut relative to crystaline axises, can be used as a tool to experimentally assess interband coupling strength and relative anisotropies of bands. It can additionally be used to distinguish the vortex bound states that we report from vortex clusters and chains forming for different reasons. In anisotropic multiband systems, there are at least two other mechanisms for formation of vortex bound states. Inclusion of density variations in the theory yields "type-1.5" regimes where coherence lengths are larger than magnetic field penetration lengths and vortex bound states form due to core-core interaction 42 . Compared to the core-coreinteraction-driven vortex binding, the mechanism considered on this paper yields much weaker vortex interaction forces. Thus vortex bound states considered here should be relatively easily destroyed by thermal fluctuations. Note also that one should expect vortices sticking to sample's boundaries due to the fact that those also feature inverted field. Different situation appears for strong anisotropies where integer flux vortices split into bound states of fractional vortices 43 . Those are easily distinguishable, due to the different magnetic field profile and coreless nature.