{'fraction_non_alphanumeric': 0.035806451612903224, 'fraction_numerical': 0.014157706093189963, 'mean_word_length': 4.308409436834094, 'pattern_counts': {'":': 0, '<': 3, '<?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': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Abundance of Weyl points in semiclassical multi-terminal superconducting nanostructures

We show that the quasi-continuous gapless spectrum of Andreev bound states in multi-terminal semi-classical superconducting nanostructures exhibits a big number of topological singularities. We concentrate on Weyl points in a 4-terminal nanostructure, compute their density and correlations in 3D parameter space for a universal RMT model as well as for the concrete nanostructures described by the quantum circuit theory. We mention the opportunities for experimental observation of the effect in a quasi-continuous spectrum.

The topological properties of quantum spectra in condensed matter systems got considerable attention in the past decade and are still under active consideration [1][2][3][4]. A large research field that has been formed thereby addresses gapped phases of insulators [5] and superconductors [6] characterized by globally defined topological numbers, and the edge modes [7] at the interfaces separating such phases. In addition to this, the spectra can exhibit topological singularities in the form of level crossings where the topological charge is defined at the singularity rather than globally. The simplest example of such singularity is a Weyl point (WP) [8] corresponding to crossing of two levels in a point in 3D space of parameters. Physical realizations of WPs include special points in the bandstructure of 3D solids [9], spectra of polyatomic molecules [10] and nanomagnets [11], quantum transport systems [12].
The occurrence of WPs have been recently predicted in the spectrum of Andreev bound states (ABS) of generic 4terminal superconducting nanostructures [13] where the 3D parameter space is formed by three independent superconducting phases of the terminals. Most important WPs are the crossings at zero energy that define the topology of the ground state. These WPs in 3D give rise to 2D global Chern numbers that are directly manifested as quantized transconductances of the nanostructure. The ideal periodicity of the space of superconducting phases allows to model higher-dimensional bandstructures with the multi-terminal superconducting nanostructures (MTSN). These ideas resulted in outburst of theoretical [14][15][16][17][18][19][20] and experimental [21][22][23][24][25] activities in the field of MTSN.
A separate recent development concerned semiclassical MTSN where a big number of ABS form a quasi-continuous spectrum. It has been predicted [26] that this spectrum can be either gapped or gapless depening on speficics of the MTSN and the point in the space of the superconducting phases. A specific topology can be introduced in semiclassical MTSN. It has been discovered and confirmed experimentally [27,28] that the gapped phases are characterised by topological numbers, and the gapless phase is explained by topological protection of these numbers [29]. The protection-unprotection transition has been discussed in this context [30].
In this Letter, we analyse the gapless spectrum at the level of discrete states and reveal the abundance of zero-energy topological singularities (Fig. 1 c). In 4-terminal structures, those are isolated WPs separated by a typical distance l c (G/G Q ) −1/2 2π. (G is a typical conductance of the nanostructure, G Q ≡ e 2 /π ). The positions of WPs are random determined by details of electron interference in the structure, while their averaged density and its correlations are determined by the structure design. We relate the density of WPs to the parameter l c governing the universal parametric correlations [31,32] in random matrix ensembles, show how to compute this density for concrete nanostructures, investigate the density correlations manifested as the transconductace of the structure, and shortly discuss the opportunities of experimental detection of the WPs in semiclassical MTSN's.
Let us start with qualitative estimations. Given a 4-terminal nanostructure of a typical conductance G one expects G/G Q , G Q ≡ e 2 /π , conduction channels, and, correspondingly, G/G Q discrete Andreev bound states affected by superconductivity. This estimation is valid both for "short" nanostructures with the typical size smaller than the superconducting coherence length, where these levels are spread in energy interval ∆, ∆ being the superconducting coherence length, and "large" nanostructure where these levels are concentrated in a much smaller energy interval E Th (G/G Q )δ S , δ S being the mean level spacing in the normal state. The energies of these levels depend on the 3 superconducting phases. Owing to periodicity in phases, the spectrum is to be considered in a Brillouin zone (BZ) of the size 2π. The RMT of parametric correlations suggests that the level energies wiggle randomly. The energies change at the scale of the level spacing at a typical distance l c in the parameter space [31,32]. This distance is determined from comparison of the mean fluctuations of the derivatives of the energies with respect to the parameters and this level spacing. For our situation, the estimation l c G/G Q in the space of phases holds for both long and short nanostructures. To understand WP's we concentrate on the level that is closest to zero energy. Upon wiggling, it will reach zero at a typical distance of the order of l c . Therefore, the total number of WP's in the Brillouin zone can be estimated as N w (l c ) −3 (G/G Q ) 3/2 .
Our detailed results (See Fig. 3) indeed give for the cross-like structures with the arm conductances G where A = 0.40 for the ballistic conductor and A = 0.16 for the diffusive one. The dimensionless coefficient A < 1, this is explained by a rather small fraction of BZ volume taken by the gapless phase (25% for ballistic and 18% for diffusive cross).
As it was shown in [13] the transconductance of the structure is defined by a Chern number of a plane traversing the BZ. The difference of two Chern numbers corresponding to two different planes is given by the total charge of the WP's enclosed between the planes. A naive estimation of the variance of this difference would be the number of WP's enclosed, d l c being the separation of the planes. This estimation would hold for randomly placed uncorrelated WP's. However, they do correlate similar to ions in an electroneutral gas: a charge of a WP is screened by other points at the distance of the order of their separation, that is, of l c . Therefore, only WP's at a distance l c contribute to the fluctuation of the Chern number and . A typical transconductance is thus G Q G. Our quantitative results are obtained in the course of three activities: A. we study numerically a generic RM model to relate the density of WPs to l −3 c and quantify the correlations of WP's. B. we develop a theory to compute l −3 c (φ) for any MTSN described by the quantum circuit theory [33] and derive concrete expressions for a single-node circuit. C. We find numerically the positions of WP's in the ballistic cross junction (Fig. 3 a) to prove the consistency of the results obtained in the activities A and B. The details of all activities are given in [34].
Activity A. The studies of statistics of spectral crossings have been pioneered by Wilkinson et al. [35][36][37]. They introduced a convenient RMT model in a 3D parameter space In this model, X i , Y i are 2N × 2N random Hermitian matrices with independent normally distributed elements of variance 1/3, N 1. Since we address the WPs in superconducting stuctures at zero energy, in distinction from [35][36][37], we choose these random matrices to obey BdG mirror symmetry of the spectrum (class C [38]). Generally, l c is defined as l −3 is the mean level spacing at the corresponding energy. For the model in use, l c = π 3/2N conveniently does not depend on φ so that the WP density is uniform. We search the positions of WPs by an iterative minimization of the energy of the closest to zero level. To make sure we find all the WPs, we repeat the iteration cycle starting it from a randomly chosen point in the parameter space. We have to do this a number of times that by a factor exceeds the expected number of points. The execution time of the algorithm thus scales as N 9/2 so we cannot access very large N and work with N = 40 − 80. For the WP concentration, we compute This is lower than the concentration of the level crossings in GUE ensemble [36] (2/3) √ πl −3 c ≈ 1.18l −3 c . We reproduce this result searching for the crossings of 10th and 11th levels.
We address the correlator of charges of the WP's, Q(0)Q(r) , r being the vector distance in units of l c . To enhance the statistics, we have evaluated an equivalent correlator of Berry curvatures of the closest to zero level. The results of 10 5 runs per point are presented in Fig. 2 and can be fitted with Since the charge density of WPs is given by the divergence of Berry curvature ( [39,40]),  D ≈ 0.5 from our calculations, dS being an area element of the surface. Activity B. While there are no perturbative methods to compute the density of WPs directly, they are available for the mesoscopic parametric correlations [41,42]. With those, one can compute l −3 c for any system characterized by electronic Green functions. We make use of the quantum circuit theory [33] that is a powerful finite-element technique for electronic Green functions. In quantum circuit theory, the structure is subdivided into reservoirs and nodes, the network is formed by connectors of various kinds, for instance, ballistic, tunnel or diffusive. The Green functions are presented by the ma-tricesĜ a ,Ĝ 2 = 1, TrĜ = 0 defined in the nodes and and the reservoirs. The semiclassical solution is obtained by minimization of an action with respect toĜ in the reservoirs at fixedĜ in the nodes.
The mesoscopic parametric correlations for a general circuit theory have been derived in [43]. For this, one substitutes to the actionĜ of double dimension, two diagonal blocks corresponding to the parameter sets 1, 2. Near the minimum, the action can be expanded up to quadratic terms with respect to non-diagonal deviations ofǦ,M being the matrix characterizing the quadratic expansion. The correlator of mesoscopic fluctuations of the action values at two parameter sets reads [43] S 1 S 2 = ln det M where 'prime' excludes the zero eigenvalues ofM from the determinant.
We implement this general technique to compute l −3 c for concrete superconducting nanostructures. It is known [44] that the energies of Andreev bound states are expressed in terms of the effective scattering matrix SS * , S being the electron scattering matrix in the space of all channels coming to the nanostructure, S * being the hole scattering matrix,the superconducting phases included. The circuit-theory action at the imaginary energy ∆ sin θ (see e.g. [45]) before the averaging over the mesoscopic fluctuations can be expressed in terms of eigenvalues SS * → −e iλ of the effective scattering matrix, these eigenvalues coming in ± pairs the last equality holding for close to zero energies/eigenvalues. The correlator of the velocities in this limit is related to the correlator of the action values, The action can be represented in a quantum circuit theory of 2 × 2 matrices, and the correlator is to be computed with the aid of Eq. 7. In Supplementary Material [34], we derive an explicit expression of l −3 c for a single-node structure. We concentrate on two simplest example MTSN (Fig. 3  a,b): a chaotic cavity connected to four leads by ballistic conductors of the same conductance G, ballistic cross, and the corresponding diffusive structure, diffusive cross. With the expressions obtained, we compute the distribution of l −3 c , and, consequently, the WP density, over the phase space, by evaluating l −3 c in random points and collecting the data into histograms: this gives a fraction of the phase space volume dV /d(l −3 c ), at a given l −3 c . The histograms for these two examples are qualitatively similar but distinct. Summing up the histograms and employing the result (1) gives the already mentioned estimations of the number of WP's, Eq. 3.
Activity C. We explicitly compute the WP positions for random chaotic cavities. For this, we pick up the electron scattering matrix S from the circular orthogonal ensemble, and find all phase settings at which SS * has an eigenvalue −1 [44]. We find 75-95 WPs for N = G/G Q = 50 conform to the results of the activities A,B and verify the scaling of the number of points with N . We also perform a more thorough check evaluating l −3 c in the random positions found and collecting the data to the histogram. The resulting estimation of dV /d(l −3 c ) that involves 2686 WP's coincides with the results of activity B within the statistical error (Fig. 3 c). In Fig. 3 e, we plot the positions of WP's found for a realization of S at G/G Q = 50. We choose the coordinate system in the space of phases to be consistent with the symmetry of the structure, χ 2,3 are defined by the above relation with cyclically permuted indexes. In these coordinates, the BZ is the truncated octahedron, as for a fcc lattice. The gapped region is in the centre of the BZ, the gapless region encloses its boundary [46]. The special points where the gapless region becomes infinitesimally thin [30] are located in the centres of the squares and hexagons, and, as seen in the Figure, the WPs are mostly concentrated in the corners of the BZ. Let us shortly discuss the methods of experimental detection of WPs in MTSN. For sufficiently large level splitting G G Q , the WPs can be found spectroscopically as the zeros of the lowest Andreev state. For G G Q where the level splitting is small not exceeding k B T , the detection is more challenging. For this case, we envisage the following detection methods: i. (Telegraph) noise measurements of the inductive or Berry curvature response of the MTSN. These responses diverge for a single discrete state at WP position. While the averaging over the states with thermal Boltzmann weights cancels the divergence, it is manifested in the noise at the time scale of the order of the time of switching between the states. ii. Transconductance (noise) measurements. We predict a transconductance GG Q . While in the presence of thermal averaging this transconductance is not quantized, its value will exhibit fluctuations as a function of the control phase [13] that can be used for scanning the WP positions. iii. Transport spectroscopy. If the MTSN is in a weak tunnel contact with a normal lead, the differential conductance of this tunnel junction exhibits low-voltage anomalies at the WP positions [47]. There is also a WP signature persisting at high voltage bias [47].
The source code and the raw data can be found at [48]. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 694272).

SUPPLEMENTAL MATERIAL
In this Supplemental Material, we present additional details about the activites A, B, and C.

ACTIVITY A: DENSITY AND CHARGE CORRELATIONS OF WP IN A UNIFORM PARAMETER SPACE
The goal of this activity is to relate the actual density of WP's with the universal correlation parameter l c . We implement a variation of Wilkinson model with random class C 2N × 2N matrices (Eq. 1 of the main text) where l c ≡ π 3/2N does not depend on the position in the 3d space of the phases. N=60 lc = 0.496 N=80 lc = 0.430 To find the WP's, we implement an iterative optimization procedure in the 3D parameter space. The optimization function is the smallest in modulus eigenvalue of the matrix. Since it involves the matrix diagonalization, the computation time scales as N 3 . The initial position is chosen randomly. The coordinates of a WP are found after several tens of iterations. To find all WP's, we repeat the procedure again and again, keeping the list of WP's found to exclude the duplicates. We learned from the experience that the procedure has to be repeated five times the expected number of WP's: further runs do not deliver new points. Since the expected number of points scales as N 3/2 , the total computation time scales as N 9/2 and really big N are not accessible for practical calculations. The following table summarizes our concrete results for the number of W P averaged for a number of runs : In the course of calculations, we have accumulated significant statistics of WP coordinates and their charges. We hoped that these statistics suffices to compute the charge-charge correlations of WP distribution. However, this did not work. The histograms approximating the charge density at a given distance from a WP exhibited significant fluctuations at relatively large distances. Our attempts to smooth these fluctuations considering the Laplace transform of the charge-charge correlator initially led us to an erroneous conclusion of a power-like tail in this correlator.
Fortunately, we checked these conclusions with an alternative method. We have computed the correlator of Berry curvatures at given distances. The most general form of the correlator of two vector quantities B α in the dimensionless coordinates r reads Since Q(r) = (4π) −1 ∂ α B α (r), the charge-charge correlator is then expressed as We compute the correlator at each r separately accumulating the statistics of Berry curvatures B α 1,2 in two random points 1, 2 separated by the distance r. The computation is relatively fast so for each point we can accumulate 10 4 samples for N = 40 and 10 5 samples for N = 20. It may seem that the two independent functions in the correlator are just given by the average products However, the evaluation is not so simple. As has been mentioned in [39], the distribution of B α has long power-law tails resulting in an infinite variance. Owing to this, the accuracy of computed averages does not increase with the number of samples in the statistics.
The universal prescription to evaluate the averages in this situation is to disregard the large values. We implement it in the following fashion: we rescale the accumulated values of B α to decrease it if large, Here, B 0 is a parameter and the modulus of the rescaledB α never exceeds B 0 . The variance ofB α is thus finite, and usual statistical considerations do apply. The value of B 0 should not be taken too small since the averages would not approximate the correlator, nor too large since the large values of B α would not be suppressed. In practice, we plot the averages versus B 0 and pick up the value of the average that persist in a large interval of B 0 . Within statistical error, B 1 (r) = 0 for all r. We can prove this analytically for r 1. With this, AS stated in the main text, the correlator of Berry curvatures can be approximated with This has no trace of long-distance power-law correlations. This expression also proves the electro-neutrality of the WP distribution, that is, the presence of a WP with the charge +1 at a point results in a depletion of average charge density around the point, the charge depleted being −1.
After completion of these calculation, we became aware of a similar calculation of the Berry curvature correlations. [40]. The authors address the correlator in a general GUE ensemble that is similar but distinct from near-zero energy correlator of interest. However, they use the same fitting function and end up with similar coefficients. In our notations, they give rB(r) ≈ 7.42e −3.56x−2.03x 2 . (S9)

ACTIVITY B: FINDING lc IN QUANTUM CIRCUIT THEORY
The goal of this activity is to find the scale governing universal parametric correlations for concrete nanostructures that can modelled with quantum circuit theory [33]. For a 3D parameter space, this scale is defined as where v i ≡ ∂E/∂φ i is the "velocity" of an energy level in the spectrum and δ S is the mean level spacing. The parameter l −3 c depends on energy as well as on a point in the parameter space.
The action, mean level spacing, and the velocity correlator At quantum level, the nanostructure is characterized by an electron scattering matrix S in the space of of the quantum channels incoming from the leads where the superconducting phase is incorporated with a factor e iφi ascribed to a channel coming from the lead i. One can derive (see e.g. [45]) an action for imaginary-time Green functions characterizing the nanostructure, being the imaginary energy, E ≡ √ 2 + ∆ 2 . It is known [44] that the Andreev bound states are related to the eigenvalues of SS * . For these eigenvalues, we will use The eigenvalues λ come in ± pairs. The energies of the bound states correspond to zeros of the expression under the log. We introduce the notation /E = sin θ, −π/2 < θ < π/2, and rewrite the action as S(θ) = − λ ln 1 − cos 2 θ cos 2 (λ/2) (S13) Since we are interested in characteristics of the spectrum near zero energy, we can expand in small θ, λ so the action becomes: ln θ 2 + λ 2 /4 (S14) In this limit, each λ gives a bound state at energy ∆λ. Let us compute the derivative of the action with respect to θ: where we made a semiclassical approximation replacing the summation over the discrete λ with integration over their continuous density ρ λ . The semiclassical action has therefore a cusp at θ = 0, the value of the cusp determines the density of the eigenvalues, that is directly related to the mean level spacing. Let us look at the random velocities of the levels, the velocity with respect to a parameter α being v α ≡ ∂ α λ The derivative of the action then reads The velocities correlate at the same level only. For the correlator of the derivatives, this gives This implies that if we know S(θ, φ)S(θ , φ) , we can evaluate the velocity correlators and l −3 c

Semiclassics: saddle point
In quantum circuit theory approach, the same action is expressed as a functional of the matricesǦ,Ǧ 2 = 1, TrǦ = 0 that is defined in the nodes and reservoirs of the nanostructure. The situation in hand we can describe with 2 × 2 matrices. These matrices are fixed in the leadsǦ i = sin θ cos θe −iφi cos θe iφi − sin θ . (S18) The matrices in the nodes are obtained by minimization of the action. In case of a short nanosctructure, the action is a sum of contributions of each connector, The function F c is proportional to the conductance of the connector and depends on the type of the connector. For instance, for a quantum point contact of conductance G F(x) = −(G/G Q ) ln((x + 1)/2). Generally, a connector is characterized by a set of transmission coefficients T p , and Long nanostructures can be described with addition of "leakage" connectors [33], yet we do not need this extension since the density of Weyl points is determined by the spectrum properties at zero energy that do not depend on the size of the nanostructure.
We will restrict ourselves to the simplest situation with a single node in the nanostructure and any number of the reservoirs. The connectors can be labelled with the lead index i, and whereǦ is the matrix in the node. We will make use of the fact thatǦ iǦ +ǦǦ i is a number rather than a matrix for any 2 × 2 matrices. It is constructive to map the 2 × 2 matrices on the corresponding 3D vectors, while the node matrixǦ → (cos Θ cos Φ, cos Θ sin Φ, sin Θ) This gives the following inner products for each connector, so the action reads To findǦ, we minimize with respect to Θ, φ. This gives two conditions: In the limit of θ → 0, this becomes We can extract level spacing from the cusp of the action at small θ. By virtue of optimization, sgnΘ = sgnθ. Therefore, at small values of θ s i → s i + |θ|| sin Θ|, and the cusp part of the action reads The density of the eigenvalues is then extracted with the aid of Eq. S15.

Semiclassics: correlations
To compute the correlations of the action at two different parameter settings, (θ, φ i ), and (θ , φ i ), we have to double the dimension of the matrices. So we consider 4 × 4 matrices. We need to do this separately for diffusion and Cooperon channels [43]. For the reservoirs, these matrices are made from two diagonal blocks, each corresponding to a setting of the corresponding reservoir. We will distinguish the settings marking or not marking then with a prime, For the diffusion channel, G i is just given by Eq. S18. For the Cooperon channel, G i is transposed.
The optimization of the action results in the block-diagonal matrix in the nodê To compute the correlations, we have to derive the quadratic expansion of the action near this optimum. With the quadratic accuracy,Ĝ is given byǦ We need to substitute this to the action and expand it to the terms quadratic inǧ. The first-order terms cancel sinceǦ 0 corresponds to the minimum of the action. This calculation is made most efficiently in the basis where G 0 is a diagonal matrix, In this basis, a matrix of a reservoir readĜ where s, u posess the index i, s 2 + |u| 2 = 1, s being inner product introduced earlier. The minimization equation in these basis reads: Let us specify those more explicitly. We choose a basis in the space of 3d vectors with the z-axis in the direction of G 0 , the angle φ is counted from Φ, x = (− sin Θ, 0, cos Θ) (S37) y = (0, 1, 0) (S38) z = (cos Θ, 0, sin(Θ)) (S39) From this, At θ → 0, Let us proceed with the expansion. We need to choose theǧ in the non-diagonal blocks and guarantee that it anti-commutes withǦ -otherwise, it will not modifyǦ and the value of the action will not change. The most general matrix of this kind can be parametrized as A straightforward but lengthy calculation results in the following form: This expression is for the diffusion channel, for the Cooperon channel we need to replace u → u * .

Close points
For our task, we need to analyse the quadratic form in close points, φ i → φ i . The form of the correlator given by Eq. S17 suggests that there is an eivenvalue of this matrix that goes to 0 at φ i → φ i and θ → 0, and the parametric dependence of this eigenvalue defines the correlations. The calculation shows no such eigenvalue in the diffusion channel, so from now on we concentrate on the Cooperon one. In this case, At θ → 0, i F (s i )s i → 0 and the determinant A 2 − |B 2 | vanishes indicating the small eigenvalue expected. The determinant can be presented as where A − ReB goes to zero in close points at θ → 0, A + ReB does not, and ImB goes to zero in close points irrespective of θ.
Let us compute A − B in conciding points. Here, the summation over i is implied. We know that i F (s i )u i = 0 at any θ. We note that s i = sin θ sin Θ + cos θ cos Θ cos(φ i − Φ) (S48) Re(u i ) = − sin Θ cos θ cos(φ i − Φ) + sin θ cos Θ (S49) To this end, we evaluate Generaling to small θ and the same points in the phase space, we obtain Let us compute A + B. We can neglect θ and the difference between the points to obtain Let us compute the terms in A − ReB that are proportional to the squares of the differences between the points. We start with Let us represent u = 1 − s 2 e iµ ; u = 1 − s 2 e iµ (S54) With this, the difference becomes There are two contributions to the difference. One comes from δµ i and reads For another one, one can set δµ = 0.
To simplify, we may add F (s)s + F (s )s that is zero under sum provided θ = 0. This gives With this, we get for the difference We also have to inspect ImB.
With all this, ln det = ln(C + D) (S63) The latter part can be presented as D = 1 2 D αβ δφ α δφ β (S67) α, β labelling the independent phases. For this, we need to express δs i , δµ i in terms of δφ i . The corresponding formulas are straightforward but rather cumbersome. In fact, we do not use those in numerical calculations, but rather compute δs i , δµ i in terms of δφ i to evaluate the quadratic form D αβ . So we do not give these formulas here.

Resulting relation
We recall that the density in the phase space can be expressed as We have derived that