{'fraction_non_alphanumeric': 0.021180350046503763, 'fraction_numerical': 0.0061300414306248415, 'mean_word_length': 4.60279488394126, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
What limits supercurrents in high temperature superconductors? A microscopic model of cuprate grain boundaries

The interface properties of high-temperature cuprate superconductors have been of interest for many years, and play an essential role in Josephson junctions, superconducting cables, and microwave electronics. In particular, the maximum critical current achievable in high-Tc wires and tapes is well known to be limited by the presence of grain boundaries, regions of mismatch between crystallites with misoriented crystalline axes. In studies of single, artificially fabricated grain boundaries the striking observation has been made that the critical current Jc of a grain boundary junction depends exponentially on the misorientation angle. Until now microscopic understanding of this apparently universal behavior has been lacking. We present here the results of a microscopic evaluation based on a construction of fully 3D YBCO grain boundaries by molecular dynamics. With these structures, we calculate an effective tight-binding Hamiltonian for the d-wave superconductor with a grain boundary. The critical current is then shown to follow an exponential suppression with grain boundary angle. We identify the buildup of charge inhomogeneities as the dominant mechanism for the suppression of the supercurrent.

current achievable in high-T c wires and tapes is well known to be limited by the presence of grain boundaries 1 , regions of mismatch between crystallites with misoriented crystalline axes. In studies of single, artificially fabricated grain boundaries the striking observation has been made that the critical current J c of a grain boundary junction depends exponentially on the misorientation angle 2 . Until now microscopic understanding of this apparently universal behavior has been lacking. We present here the results of a microscopic evaluation based on a construction of fully 3D YBa 2 Cu 3 O 7−δ grain boundaries by molecular dynamics. With these structures, we calculate an effective tight-binding Hamiltonian for the d-wave superconductor with a grain boundary. The critical current is then shown to follow an exponential suppression with grain boundary angle. We identify the buildup of charge inhomogeneities as the dominant mechanism for the suppression of the supercurrent.
To explain the exponential dependence of the critical current J c on the misorientation angle α Chaudhari and collaborators 3 introduced several effects which can influence the critical current that are particular to high-temperature superconductor (HTS) grain boundaries (GB). First, a variation with angle can arise from the relative orientation of the d-wave order parameters pinned to the crystal lattices on either side of the boundary. This scenario was investigated in detail by Sigrist and Rice 4 . However, such a modelling cannot explain the exponential suppression of the critical current over the full range of misorientation angles.
Secondly, dislocation cores, whose density grows with increasing angle, can suppress the total current. A model assuming insulating dislocation cores which nucleate antiferromagnetic regions and destroy superconducting order was studied by Gurevich and Pashitskii 5 . However, for grain boundary angles beyond approximately 10 • when the cores start to overlap this model fails. Finally, variations of the stoichiometry in the grain boundary region, such as in the oxygen concentration, may affect the scattering of carriers and consequently the critical current. Stolbov and collaborators 6 as well as Pennycook and collaborators 7 have examined the bond length distribution near the grain boundary and calculated the change in the density of states at the Fermi level, or the change in the Cu valence, respectively. In the latter work the authors used the reduced valences to define an effective barrier near the boundary whose width grows linearly with misorientation angle.
A critical examination of the existing models shows that the difficulty of the longstanding HTS "grain boundary problem" arises from the multiple length scales involved: atomic scale reconstruction of the interface, the electrostatic screening length, the antiferromagnetic correlation length, and the coherence length of the superconductor. Thus it seems likely that only a multiscale approach to the problem can succeed.
Our goal in this paper is to simulate, in the most realistic way possible, the nature of the actual grain boundary in a cuprate HTS system, in order to characterize the multiple scales which cause the exponential suppression of the angle dependent critical current. To achieve this goal we proceed in a stepwise fashion, first simulating the atomic structure of realistic YBCO grain boundaries and assuring ourselves that our simulations are robust and duplicate the systematics of actual grain boundaries. Subsequently we construct an effective disordered tight-binding model, including d-wave pairing, whose parameters depend on the structures of the simulated grain boundaries in a well-defined way. Thus for any angle it will be possible to calculate the critical current; then, for a given pairing amplitude (reasonably well known from experiments on bulk systems) the form of J c (α) and its absolute magnitude is calculated.
We simulate YBCO grain boundaries by a molecular dynamics (MD) procedure which has been shown to reproduce the correct structure and lattice parameters of the bulk YBa 2 Cu 3 O 7−δ crystal 9 , and adapt techniques which were successfully applied to twist grain boundaries in monocomponent solids 10 . We only sketch the procedure here, and postpone its details to the supplementary information and a longer publication. The method uses an energy functional within the canonical ensemble where m i is the mass of the ion and U(r ij ) is the effective potential between ions, taken to be of the form U(r ij ) = Φ(r ij ) + V (r ij ). Here V (r) is the screened Coulomb interaction V (r) = ±e −κr Z 2 e 2 /(4πǫ 0 r) with screening length κ −1 , and Φ(r) is a short range Buckingham potential Φ(r) = A exp(−r/ρ) − C/r 6 . We take the parameters A, ρ and C from Ref. 9, and the initial lattice constants are a = 3.82Å, b = 3.89Å and c = 11.68Å. To construct a grain boundary with well defined misorientation angle, we must fix the ion positions on both sides far from the boundary, and ensure that we have a periodic lattice structure along the grain boundary and also along the c-axis direction. Therefore we apply periodic boundary  Table   1 of Ref. 7 (for notation see this reference).
conditions in the molecular dynamics procedure in the direction parallel to the GB and also in the c-axis direction. In the direction perpendicular to the GB only atomic positions with a distance smaller than six lattice constants from the GB are reconstructed. This method restricts our consideration to grain boundary angles that allow a commensurate structure parallel to the GB, e.g. angles α = 2 · arctan[N 1 a/(N 2 b)], N 1 , N 2 ∈ N 11 . An important step to be taken before starting the reconstruction of the GB is the initial preparation of the GB. Since we use a fixed number of ions at the GB within the molecular dynamics algorithm we have to initialize it with the correct number of ions. If we start with two perfect but rotated crystals on both sides of the interface, in which all ions in the half space behind the imaginary boundary line are cut away, we find that several ions are unfavorably close to each other. Here we have to use a set of selection rules to replace two ions by a single one at the grain boundary. These selection rules have to be carefully chosen for each type of ion since they determine how many ions of each type are present in the grain boundary region; the rules are detailed in the supplementary information. Ultimately they should be confirmed by a grand canonical MD procedure. However, such a procedure for a complex multicomponent system as the YBCO GB is technically still not feasible.
While the selection rules are ad hoc in nature, we emphasize that they are independent of misorientation angle, and reproduce very well the experimental TEM structures 7 .
With the initial conditions established, the MD equations of motion associated with Eq. 1 are iterated until all atoms are in equilibrium. As an example, we show in Fig. 2 the reconstructed positions of Y, Cu, and O ions for a (410) boundary. We emphasize that the MD simulation is performed for all the ions in the YBCO full 3D unit cells of two misoriented crystals except for a narrow "frame" consisting of ions that are fixed to preserve the crystalline order far from the grain boundary. The sequence of typical structural units identified in the experiments 7 is also indicated and we find excellent agreement.
We next proceed to construct an effective tight-binding model which is restricted to the Cu sites, the positions of which were determined through the algorithm described. We cal- probabilities along the boundary induce scattering processes of the charge carriers and consequently contribute to an "effective barrier" at the grain boundary. We note that the angular dependence of the critical current J c (α) is not directly related to the changes in averaged hopping parameters observed for different misalignment angles. In fact we found that the variation of the hopping probabilities with boundary angle cannot account by itself for the exponential dependence of J c (α) over the whole range of misalignment angles.
The structural imperfection at the grain boundary will necessarily lead to charge inhomogeneities that will contribute-in a similar way as the reduced hopping probabilities-to the effective barrier that blocks the superconducting current over the grain boundary. We include these charge inhomogeneities into the calculations by considering them in the effec-  Finally we have to translate the charge on the Cu sites (or better CuO 4 squares) into effective on-site lattice potentials. The values of screening lengths in the cuprates, and particularly near grain boundaries, are not precisely known, but near optimal doping they are of order of a lattice spacing or less. We adopt the simplest approach and assume a 3D Yukawa-type screening with phenomenological length parameter ℓ of this order, and find a potential integrated over a unit cell V 0 of whereq is the charge in units of the elementary charge, a 0 is the Bohr radius, and ℓ is taken to be 2Å while a = 4Å. Thus we find a surplus charge of a single elementary charge integrated over a unit cell to produce an effective potential of around 10 eV. In the following we will use an effective potential of either 6 or 10 eV, reflecting the uncertainty in this parameter.
The value of the effective potential will affect the scale of the final critical current.
To calculate the order parameter profile and the current across the grain boundary we solve the Bogoliubov-de Gennes mean field equations of inhomogeneous superconductivity self-consistently. The Hamiltonian iŝ where the effective hopping parameters t ij are determined for a given grain boundary by the procedure described above and the onsite energies ǫ i = u i − µ are a sum of the effective charge potentials u i and the chemical potential µ. Performing a Bogoliubov transformation, we find equations for the particle and hole amplitudes u n and v n The self-consistency equation for the d-wave order parameter is then where we use N sc supercells in the direction parallel to the GB and k y is the corresponding Bloch wave vector. We adjust the chemical potential µ to ensure a fixed carrier density in the superconducting leads corresponding to 15% hole doping. The definition of the d-wave pair potential V ij in the vicinity of the grain boundary is not straightforward since the bonds connecting a Cu to its neighbors are not exactly oriented perpendicular to each other. We With these preliminaries, the current itself can finally be calculated by imposing a phase gradient across the sample (see, e.g. Ref. 19 for details) from the eigenfunctions u n , v n and eigenvalues E n of the BdG equations for the grain boundary, The critical current J c is defined as the maximum value of the current as a function of the phase. In the tunnelling limit, when the barrier is large, this relationship is sinusoidal so the maximum current occurs at phase π/2. However, for very low angle grain boundaries we observe deviations from the tunnelling limit, i.e. higher transparency grain boundaries with non sinusoidal current-phase characteristics, although for the parameters studied here these deviations are rather small.
It is instructive to examine the spatial pattern of supercurrent flow across a grain boundary, which is far from simple, as illustrated in To account for the difference of the calculated gap magnitude from its experimental value we multiply the current by a factor of ∆ exp /∆ 0 , where ∆ exp is the experimentally measured order parameter and ∆ 0 is its self-consistently determined bulk value. We have checked that at low temperatures T ≪ T c an approximately linear scaling of the critical current as a function of the order parameter holds true for all grain boundary angles.
In Fig. 6, the critical current is plotted as a function of misorientation angle for a set of grain boundary junctions from (710)

The initial setup of the grain boundary
The simplest way of modelling a tilt grain boundary is to stitch two perfect crystals together that are rotated into one another around an axis in the plane of the GB perpendicular to a lattice plane; atoms which cross into the region of space initially occupied by the other crystal are then eliminated (see Fig. 7 a). If we examine the structures constructed in this way we find on the one hand that some ions are left unphysically close to each other, while on the other hand large "void" regions may also remain. Since a molecular dynamics scheme employing an energy functional in the canonical ensemble does not allow for the creation and annihilation of ions in the GB region these faults in the setup will persist during the reconstruction process. To improve the initial structures we develop a set of "selection rules": (i) We introduce an overlap of the two crystals, extending them into a region behind the virtual GB plane to prevent the creation of "void" regions. (ii) We replace two ions that are too close to each other by a single one at the GB. Due to the different ionic radii, the different multiplicities and the different internal positions these "selection rules" have to be adjusted for every ion in the unit cell in order to reproduce the correct stoichiometry found in TEM experiments 7 . In Table I we show the maximum overlap of ionic positions behind the GB plane as well as the minimum distance between two ions allowed before replacing them by a single ion at the GB. These values have proven to reproduce the experimentally observed GB structures for all misorientation angles examined within this work. An example of such a constructed GB is shown in Fig. 7 b).

The molecular dynamics procedure
The GB structures determined by applying the "selection rules" described in detail in the previous section are now used to initialize the molecular dynamics process.     we take the parameters A, ρ and C from molecular dynamics studies by Zhang and Catlow 9 leading to a stable YBCO lattice with reasonable internal coordinates of each atom (see Table II). In addition we use the lattice constants a = 3.82Å, b = 3.89Å, and c = 11.68Å when setting up the initial GB structure and we fix them in the leads far away from the GB. The CuO chains in the CuO layer are directed parallel to the b-axis direction (compare Fig. 8).
To construct a GB with well defined misalignment angle we have to fix the atomic po- As a first step we calculate the hopping between the Cu-d x 2 −y 2 orbital and the O-p x orbital that are connected by the vector r = d e x and therefore l = 1 and m = n = 0. The angular dependence is introduced as In the next step we have to calcuate the distance-dependent potential of the σ-bond Again we need in addition the distance-dependent potential of the σ-bond of two Cu-d orbitals: and the total energy associated with the direct overlap of two Cu-d orbitals can finally be calculated as Now we can compare this energy to the kinetic energy describing the superexchange between two Cu sites via the O-p orbitals  In Fig. 9 we show the averaged hopping as a function of the distance to the grain boundary plane for different GBs. If we fit the suppressed hopping values in the vicinity of the grain boundary by a Gaussian form and integrate over the "effective barrier" derived in this way, we find only a linear variation with misalignment angle (see Fig. 10).  Fig. 9.

The bond valence analysis
The structural imperfection at the grain boundary will necessarily lead to charge inhomogeneities that will contribute -in a similar way as the reduced hopping -to the effective barrier that blocks the superconducting current over the GB. We can include these charge inhomogeneities in our calculations by "translating" them into on-site potentials on the Cu sites. It is evident that we also have to include the charges on the O sites in our considerations although the O sites themselves have already been integrated out in the effective one-band tight-binding model. We will start our considerations by assigning every atom of the perfect crystal a formal integer ionic charge: Y 3+ , Ba 2+ , O 2− . The requirement of charge neutrality leaves us with 7 positive charges to be distributed on the 3 Cu atoms: Thus we will have two Cu 2+ , and one Cu 3+ ion per unit cell. For a crystal with bonds that are neither strictly covalent nor strictly ionic it is convenient to introduce a fractional valence for each ion, that is determined with respect to its ionic environment in the unit cell of the crystal.
Here we calculate the bond valence of a cation by where the sum is over all neighboring anions, in our case the neighboring negatively charged O and r ij is the cation-anion distance. Here we take B = 0.37Å following Ref. 16 On the other hand, this should be equal to the sum of a fraction of bond valences V

Solving this set of equations for ξ
(3) i allows us to determine the fraction of Cu ions with the valence 3+: In a similar way we can proceed assuming that we have a Cu ion that could be in a 1+ or a 2+ oxidation state. Here we will use r 0 (Cu 2+ ) = 1.679Å and r 0 (Cu 1+ ) = 1.6Å and we In the case that Cu 2+ and Cu 3+ are the most probable oxidation states we will find that ξ . This method ensures that we end up with a charge neutral crystal.
In the bulk system we derive fractional valences close to the formal integer oxidation states. In the vicinity of the GB we find however deviations from the bulk values due to missing or displaced neighboring ions.

The definition of the superconducting pairing interaction
To model the known momentum space structure of the superconducting order parameter in the weak coupling description of the Bogoliubov-de Gennes theory, one usually defines an interaction V ij on the bonds connecting two nearest neighbor Cu sites. Unfortunately, there is no obvious way to define an analogous pairing interaction in a strongly disordered region of the crystal, as e.g. in the vicinity of the GB, since the exact microscopic origin of this interaction is not known. Assuming that a missing or a broken Cu-O bond destroys the underlying pairing mechanism, we develop a method of tying the superconducting pairing interaction between two Cu sites to the hopping matrix elements connecting the two sites as well as to the charge imbalance between them. Hence we define the pairing interaction on a given bond as the product of a dimensionless constant V 0 , that is adjusted to reproduce the correct modulus of the gap in the bulk, and the hopping parameter t ij . In addition we impose an exponential suppression of the pairing strength with increasing charge imbalance: To avoid long range contributions to the pairing we use a threshhold of 0.2 eV for |t ij |, thus restricting the pairing interaction to the bonds between nearest neighbor Cu sites in the bulk. Here we emphasize that although we tried to model the pairing interaction in a realistic way, the exact procedure how we define the pairing interaction in the disordered region does not qualitatively change the results. The superconducting current over a weak link, as well as over a metallic or an insulating barrier is accompanied by a drastic change of the phase of the superconducting order parameter. For a true tunnel junction, for which the two superconducting regions are completely decoupled, the current-phase relation is known to be sinusoidal, and the maximum current is found for a phase jump of π/2. For a weak link, as e.g. provided by a geometrical constriction or a disordered region like a grain boundary, the phase of the superconducting order parameter has to change continously between its two bulk values in the leads on both sides of the weak link. Here one would expect deviations from the sinusoidal current-phase relation and a smooth transition of the phase over a length scale given by the dimensions of the weak link. Besides the change in the phase of the order parameter in the presence of an applied current, a weak link is also characterized by a suppression of the modulus of the or- der parameter, either due to geometric restrictions or due to strong disorder. An additional suppression of the order parameter on the length scale of the coherence length can occur due to the formation of Andreev bound states present at specifically orientated interfaces of d-wave superconductors. Here we have calculated the order parameter self-consistently from Eq. 8 in the main text and found a suppression of its modulus in the vicinity of the grain boundary that depends in its width and depth only weakly on the misalignment angle (compare Fig. 11 a,b).
For a numerical determination of the critical current -the maximum current that can be applied to a system without destroying its superconducting properties -it is convenient to enforce a certain phase difference of the order parameter between two points of the system separated by a distance d, e.g. in our example on both sides of the grain boundary region. Here we calculate the superconducting current over the grain boundary modelled by the Hamiltonian in Eq. 5 of the main text using Eqs. 9 and 10 therein and fixing the phase of the order parameter in the leads in every iteration step. The critical current for our system can then be found as the current maximum under a variation of the phase difference. If a strong perturbation limits the superconducting current -in the so-called tunneling limit -the main drop of the superconducting phase will appear in the small region of the perturbation and the change of the phase in the superconducting leads is negligible (see Fig. 11 c). However, if the perturbation is weak, the change of the phase in the superconducting leads becomes more important (see Fig. 11 d) and the current will depend on the distance d between the two points, at which the phase of the order parameter is fixed. Since we are interested in the calculation of the critical current, we have to decrease the distance, thus increasing the phase gradient, until we reach the maximum current. To calculate the current over the low angle grain boundaries, we have determined the maximum current by extrapolating it to the value expected for d = 0 (see Fig. 12).