{'fraction_non_alphanumeric': 0.03675724433419777, 'fraction_numerical': 0.01080166873379186, 'mean_word_length': 4.646294881588999, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cooper pairs localization in tree-like networks of superconducting islands

We study inhomogeneous Cooper pairs distribution and localization effects in tree-like networks of superconducting islands coupled via Josephson weak links. Using a generalized Feynman's approach, reminiscent of the Bose-Hubbard model, we demonstrate that the Cooper pairs fraction which localizes on a specific network's island is limited by the network topology and, if present, by the repulsive interaction. These findings contribute to clarify the interplay between confinement effects induced by the network's topology and interaction and shed some light on recent experiments dealing with networks of Josephson junctions.

We study inhomogeneous Cooper pairs distribution and localization effects in tree-like networks of superconducting islands coupled via Josephson weak links. Using a generalized Feynman's approach, reminiscent of the Bose-Hubbard model, we demonstrate that the Cooper pairs fraction which localizes on a specific network's island is limited by the network topology and, if present, by the repulsive interaction. These findings contribute to clarify the interplay between confinement effects induced by the network's topology and interaction and shed some light on recent experiments dealing with networks of Josephson junctions.

I. INTRODUCTION
The discovery of superconductivity 1 and its theoretical explanation 2,3 , i.e. the Bardeen-Cooper-Schrieffer (BCS) theory, are among the most important achievements in physics. The BCS paradigm successfully explains a large variety of phenomena ranging from the organization of nuclear matter 4,5 to the superconductivity in metals 6 and conjectures exist about the possibility that the final quantum state of a black hole would be similar to the ground state of a BCS superconductor 7,8 . Apart from its cultural centrality, superconductivity is nowadays an essential ingredient to implement topological states of matter 9 (e.g., the Majorana's modes) and quantum computers 10 . For this reason, there is a renewed interest towards the possibility of conditioning the superconductivity by using different means. Traditionally, new superconducting materials have been identified by doping a parent compound in a controlled fashion. Even though the outcome of this procedure is hardly predictable, in fortunate cases it provides superconducting alloys with enhanced critical temperature. Extensive application of this procedure led to the discovery of high critical temperature (HTC) superconductors 11 and the impact of this success has fueled the dream of finding a room-temperature superconductor for a long time.
Recently, it has been proposed that dielectric metamaterials can be used to create effective media in which the electron-electron effective interaction (depending on the inverse of the dielectric function) can be enhanced 12,13 . The enhancement of the interaction between electrons requires the ability of engineering the dielectric function of an artificially designed system. In this way a near-zero dielectric function (ENZ) metamaterial with enhanced superconducting critical temperature can be obtained. The aforementioned theoretical framework shows a good agreement with the experimental realizations 14,15 (e.g. Al-Al 2 O 3 core-shell metamaterials) and predicts the possibility to achieve near room-temperature superconducting systems based on conventional superconductivity 16 . Interestingly, the enhancement of the superconducting critical temperature of a metamaterial based on a conventional superconductor with bulk gap ∆ (and coherence length ξ BCS ∝ ∆ −1 ) implies that the resulting superconducting gap ∆ ef f is also enhanced compared to ∆. Provided that the BCS relation between superconducting gap and coherence length is also valid for the effective medium, the above observation suggests that the coherence length of the effective medium ξ ef f is shortened compared to ξ BCS . This conclusion is not surprising since the typical distance ∼ ξ ef f ∝ ∆ −1 ef f between the paired electrons in a Cooper pair depends on the effective interaction strength g ef f ∝ ∆ ef f . A less obvious implication of this analysis is that superconducting gap and critical temperature of a synthetic superconductor can be enhanced by reducing the effective coherence length of the system. A reduced coherence length is observed in superconductors altered by the introduction of an appropriate amount of non-magnetic disorder (impurities). These systems are described under appropriate circumstances 17 by the dirty limit coherence length ξ ef f ≈ √ ℓ ξ BCS [18], which on its turn depends both on the coherence length ξ BCS of the clean system and the disorder-dependent mean free path ℓ. Assuming that ∆ ef f ∝ ξ −1 ef f , one easily get ∆ ef f ∝ ∆ ξ BCS /ℓ. Thus, a gap enhancement (i.e., ∆ ef f > ∆) would be expected when the effective coherence length ξ ef f is limited by the mean free path (i.e. when ℓ < ξ BCS ). The condition ℓ < ξ BCS is hardly reachable by using HTC superconductors because they are characterized by very short coherence length (few nanometers in cuprates 19 ). Thus, it is expected that the disorder-induced gap enhancement works well for low critical temperature superconductors. These arguments are not in contradiction with the Anderson's theorem [20] according to which nonmagnetic disorder does not affect the superconducting transition temperature in conventional superconductors. The Andreson's statement, indeed, is formulated under the restrictive condition of space-independent order parameter, which typically is not realized in disordered systems. A recent generalization of the Anderson's argument [21] confirms the possibility of enhancing the superconductivity by using disorder. Thus the above arguments appear to be plausible provided that ∆ ef f is interpreted, e.g., as space-averaged order parameter. Despite the aforementioned conclusions are based on rather naive arguments, their validity can be proven to some extent by using a microscopic theory based on the self-consistent solution of the Bogoliubov-de Gennes equations 22 . In particular, a critical temperature enhancement has been reported for disordered superconducting systems in close vicinity of the Anderson localization condition [23][24][25] . Thus, exponential localization seems to be a favorable condition to obtain a critical temperature enhancement. Exponential localization of the superconducting order parameter can be also obtained in tree-like networks of superconducting islands coupled via Josephson weak links. These systems, originally proposed as analogue model to study Bose-Einstein condensation on inhomogeneous graphs [26][27][28] , have been the object of intense and systematic experimental investigations [29][30][31][32] . Experiments confirm that network nodes with higher connectivity act as localization centres for the Cooper pairs. Moreover, indirect signatures of a gap enhancement are also reported 33 . The above evidences suggest that these artificial systems mimic the interstitial structure of the superconducting order parameter established at the phase boundary between the superconducting state and an insulating state (see Fig. 4(n) of Ref. [22]). This intuition is corroborated by a de Gennes-Alexander micronetworks formulation 34 of the problem showing a close connection between order parameter localization and critical temperature enhancement 35 . In these systems, exponential localization of the order parameter is controlled by the network's topology (instead of the impurity potentials) and for this reason the graph connectivity directly affects the superconducting phase transition 36 . In view of the relevance of the Cooper pairs localization phenomenon in graph-like structures, it would be important to verify the existence of mechanisms limiting the particles accumulation on a given system's island. In this work we demonstrate that the Cooper pairs fraction which localizes on a specific network's island is limited by the network topology and, if present, by the repulsive interaction. Treating Cooper pairs like ordinary bosons, we derive a generalized Feynman's model through which inhomogeneous Cooper pairs distribution and localization effects are systematically studied.
We will focus on tree-like structures which appear to be an ideal platform to study localization phenomena induced by the network topology. Differently from arbitrary graphs, two nodes (vertices) of a tree-like network are connected by a unique path, so that interference effects have a limited influence on the particles wavefunction. Moreover, the absence of loops limits unintentional inductive couplings and thus it is a desirable condition for the experimental implementation using Josephson junctions' networks. The work is organized as follows. In Sec. II we review the Feynman's model of a Josephson junction. Its generalization obtained by the Bose-Hubbard model is presented in Sec. III. The resulting model is a non-linear Schrödinger equation whose ground state solution allows to study the interplay between repulsive interaction and network connectivity. Selected applications of the theory are discussed both in the presence (Sec. V) or in the absence (Sec. IV) of repulsive interaction. Conclusions are given in Sec. VI, while technical details are presented in Appendix A, B and C.

II. FEYNMAN'S MODEL OF A JOSEPHSON JUNCTION
In his famous lecture 37 , Feynman proposed a simple model of Josephson junction [38][39][40] . According to the model, the dynamics of the macroscopic wavefunctions of the superconducting electrodes, i.e. ψ 1 and ψ 2 , is described by the equations: where the dot notation indicates the time derivative. The constant K describes the weak coupling of the two sides of the junction and depends on the characteristics of the tunnel barrier. The presence of an external circuit is implicitly assumed. It induces a difference between the ground state energies of the two condensates, so that E 1 − E 2 = qV represents such difference written in terms of the absolute vale of the Cooper pair charge q = 2e and the voltage drop across the junction V . The macroscopic wavefunctions of the Cooper pairs condensates can be written in the form where ρ α and θ α represent the Cooper pairs density and the phase factor pertaining to the electrode labelled by α, respectively. Using the exponential representation of the macroscopic wavefunctions in Eq.(1) and introducing the phase difference ϕ = θ 2 − θ 1 , we get: The dynamics of the phase difference ϕ is easily deduced by using Equations (3)-(4) and can be written as follows: Due to the action of the external circuit, tunneling phenomena across the junction do not change the Cooper pairs densities ρ 1 and ρ 2 , so that one can assume ρ 1 = ρ 2 = ρ 0 , being ρ 0 the time-independent particles density of the isolated electrodes. The implementation of the latter conditions is problematic within the Feynman's model and, to be consistent, one has to assume ρ 1 = ρ 0 −δρ and ρ 2 = ρ 0 + δρ with δρ/ρ 0 ≪ 1. Under these assumptions and after a first-order expansion in δρ, Eq. (2) is written as: which represents the correct current-phase relation of the Josephson junction. Under the same assumptions, Eq.(5) takes the form: the latter showing a deviation from the strict voltagefrequency relation:φ = qV .
The deviation, consisting in a cosine term in the righthand side of Eq. (7), originates from the simplifying approximations of the model. Indeed, in real systems the external circuit, neglected within the Feynman's approach, maintains δρ = 0. The latter constraint can be implemented within the semiclassical treatment proposed by Ohta 41 , where the conditionρ 1 =ρ 2 = 0 is exactly respected. Despite these limitations, the Feynman's model provides the correct constitutive equations of the Josephson effect (Equations (6) and (8)) when a zero-order expansion in δρ is considered. From the physics viewpoint, the latter expansion is justified, e.g., when the coupling constant K is so small that the Cooper pairs densities remain practically unchanged (ρ 1,2 ≈ ρ 0 ) as long as the time τ ≪ /(2K) is reached. We will show that the cosine term in Eq. (7), which is not expected in the voltage-frequency relation of a biased Josephson junction, is a necessary ingredient to get an appropriate description of isolated systems (i.e., not connected to external reservoirs). Thus, in the following, we present a generalization of the Feynman's model allowing the description of isolated networks of superconducting islands coupled via the Josephson effect.

III. FROM BOSE-HUBBARD TO THE GENERALIZED FEYNMAN'S MODEL
We are interested in studying the Cooper pairs distribution in tree-like networks of superconducting islands. Cooper pairs are composite particles obeying hard-core bosonic algebra, being the deformed bosonic algebra reminiscent of the fermionic nature of the paired electrons. Treating the Cooper pairs as genuine bosons, which is certainly an approximate picture, we can study the treelike systems mentioned above by means of the Bose-Hubbard model. This approximate picture is not too crude. For instance, quantum critical phenomena in granular superconductors have been investigated by modeling the Cooper pairs as a quantum fluid of charged bosons [42]. Similarly, quantum phase transitions of Josephson junctions arrays are often described by adopting a short-range Bose-Hubbard model [26,[43][44][45], being the latter equivalent, under appropriate conditions, to the second quantization Hamiltonian of a Josephson junctions array [43]. Second quantization Hamiltonian of the Bose-Hubbard model describing a tree-like system is written as where ǫ i is a site-dependent energy, K > 0 is a coupling constant, U > 0 is the strength of the repulsive interaction and n i = b † i b i is the number operator of the lattice site i written in terms of bosonic creation and annihilation operators. The islands connectivity is defined by the adjacency matrix of elements A ij . Adjacency matrix is real and symmetric with vanishing diagonal elements (A ii = 0). Islands labelled by the site index i and j are linked if A ij = 1, while they are disconnected if A ij = 0. Using Heisenberg equation of motion and bosonic com- Under macroscopic occupation of the islands (i.e. b † k b k = N k ≫ 1), the annihilation operator b k in Eq. (10) can be substituted with the condensate wavefunc- [46]. Moreover, due to the latter assumption, fluctuations are suppressed and the condensate wavefunction depends on the well-defined semiclassical variables N k and θ k [47]. In this way, the discrete non-linear Schrödinger equation 48 is obtained. Equation (11) preserves the total number of particles and thus the macroscopic wavefunction of the condensate can be normalized to the total number of Cooper pairs N T according to the relation k |ψ k | 2 = N T . Interestingly, Feynman's model given in Eq. (1) can be obtained from Eq. (11) by considering a two-sites system with vanishing repulsive interaction (U = 0). For these reasons, Eq. (11) can be considered a generalization of the Feynman's model. In the following analysis, we will use Eq. (11) to characterize the inhomogeneous distribution of Cooper pairs in tree-like networks of identical superconducting islands. The nonlinear term in Eq. (11) mimics, at mean-field level, the interaction between Cooper pairs 49 . The latter contribution is expected to be negligible in bulk superconductors, while it may be relevant when superconducting islands of reduced volume are considered. In view of our interest in isolated systems with identical superconducting islands, we can also set ǫ k = 0 in Eq. (11). Disorder effects (not treated in this work) can be studied by considering random site's potentials and randomly distributed hopping terms.
A. Searching for the ground state of the system: general procedure The ground state of the system can be obtained by considering the ansatz ψ k = C k exp(iθ), with C k = √ N k , which ensures vanishing Josephson currents between the superconducting islands (see Appendix A). When the trial wavefunction is substituted in Eq. (11), the following relations are obtained: A close inspection to Eqs. (12)-(13) shows thatθ cannot depend on time and thus we setθ = ω, with ω a time independent constant. The identification of the ground state requires the solution of a non-linear eigenvalues problem (see Eq. (13)) with eigenvalues of the form µ = − ω.
Before treating the general case, it is quite instructive considering the U = 0 case for which we have to solve an ordinary eigenvalues problem. In this case, the problem is described by a pure hopping Hamiltonian H = −KA which, in real space, is a real and symmetric matrix proportional to the adjacency matrix A via the constant −K. Thus, the hopping Hamiltonian and the adjacency matrix are simultaneously diagonalizable and share a common set of eigenvectors. The eigenproblem of the adjacency matrix is written as A C = λ C with λ = ω/K = −µ/K. The eigenvector C M of the adjacency matrix with greatest eigenvalue λ M also corresponds to the unique ground state of the system in view of the relation The spectral graph theory (Perron-Frobenius theorem) prescribes that the components of C M are all positive, while λ M is bounded from above by the maximum degree of the graph [50]. Moreover, there are no other positive eigenvectors of the adjacency matrix (except positive multiples of C M ) so that the theorem also implies the unicity of the bosonic ground state. Once C M is known, it provides information about the particles distribution. These arguments show the close connection between graph theory and the ground state properties of the bosonic system. The inclusion of a finite repulsive interaction modifies the aforementioned picture and the nonlinear problem has to be solved numerically. Numerical procedure is conveniently implemented by introducing eigenvectors normalized to 1 instead of the total number of Cooper pairs . Using the above notation, the nonlinear problem takes the following form: with ξ = N T U/K a small, but finite, perturbation of the linear problem. As before, the greatest eigenvalue of the nonlinear problem in Eq. (14) is associated to the ground state eigenvector. The solution of the nonlinear problem in Eq. (14) can be obtained by using an iterative procedure which is similar to that customarily used when a many-body system is studied within the Hartree-Fock approximation. The initial nonlinear problem is mapped into the linear eigenproblem 51 complemented by the self-consistency condition V kj = −ξδ kj B 2 j . The solution is iteratively determined according to the following steps. The matrix V kj is computed by using the eigenstate with the greatest eigenvalue (EGE) of the noninteracting problem (ξ = 0). Thus, eigenvectors and eigenvalues of A + V are determined. The EGE is used to recompute A + V . After diagonalization of the latter operator, a new estimate of the EGE is obtained and the procedure can be iterated. The algorithm ends when the self-consistency condition |V kj + ξδ kj B 2 j | < ε is reached with accuracy ε. The procedure typically presents fast convergence with accuracy very close to the machine precision (∼ 10 −12 in dimensionless units). Once the convergence has been reached, the solution of the nonlinear problem in Eq. (14) is obtained and the spatial distribution of the particles can be determined. In this section, we present selected applications of the theory for which analytical solutions are available. In this way, the particles distribution can be studied avoiding numerical diagonalization of the adjacency matrix. In view of their pedagogical value, fully worked examples are reported.

A. Linear chain
In the following we will explain in detail the problem of the Cooper pairs distribution on a linear chain ( Fig.   1). In particular, we consider N T non-interacting bosonic particles (i.e. the Cooper pairs) constrained to move on a linear connection of N superconducting islands. Within this assumptions, the time evolution of the bosonic wavefunctions is described by a pure hopping Hamiltonian leading to the following equation: Using ψ n = √ N n exp(iθ n ) in Eq. (16), we get the coupled equations: showing the presence of cosine terms which, as anticipated in Sec. II, are crucial for a correct description of the isolated system. When a stationary state of the closed system is considered, the number of particles residing on a single island does not change in time (Ṅ n = 0). This requirement can be obtained by setting θ n = θ(t), withθ(t) = ω a time-independent constant. Accordingly, using the notation λ = ω/K and C n = √ N n , we obtain the difference equation: complemented by the boundary conditions C 0 = C N +1 = 0. The difference equation in Eq. (19) can be presented in the form AC = λC, with A the adjacency matrix of the linear chain and C = (C 1 , ..., C N ) t . The diagonalization of the adjacency matrix provides the EGE, which is related to the particles distribution. Alternatively, we can search for solutions of Eq. (19) compatible with the boundary conditions. Such solutions present the form C (m) n ∝ sin(k m n) with k m = mπ/(N + 1). However, sign-changing solutions are not admissible and thus, after normalization, we get: with associated eigenvalue λ = 2 cos(k 1 ). C (1) also corresponds to the ground state of the system with energy eigenvalue µ 0 = −2K cos(k 1 ). The fraction of particles belonging to the n-th superconducting island is given by the following relation: showing the tendency to populate preferentially the central region of the system. This tendency originates from the fact that highly connected islands are preferentially occupied with respect to network's sites with connectivity deficiency (e.g., the system's boundaries). The aforementioned statement can be verified in simple cases since Eq. the following equations: where n ≥ 1 is a labeling of the generic network site belonging to the branch α ∈ {a, b}, while the additional relation ψ (α) 0 = ψ 0 plays the role of a boundary condition. Using the exponential notation ψ −1 = N −1 exp(iθ), where we have introduced the notation A −1 = N −1 , and λ = ω/K. Numerical diagonalization of the problem shows that the EGE presents exponential localization. The latter statement can be easily proven by using the following ansatz: being χ, A 0 and ρ < 1 parameters to be determined. Using the ansatz in Eqs. (25)- (27), we obtain the following relations: Equations (31) and (32) imply the equation λ 4 − 4λ 2 − 1 = 0, whose physically acceptable solution is given by λ = 2 + √ 5. Once λ is known, direct calculation shows that: while A 0 is fixed by imposing the wavefunction normalization condition: Using the aforementioned results and assuming N ≫ 1, the particles fraction allocated in the network sites labelled by n = 0 and n = −1 is obtained: Moreover, the particles occupation along the α branch is easily obtained by using the relation N (α) n = N 0 ρ 2n . One of the main conclusions of this section is that particles manifest the tendency to localize exponentially around the system's island presenting the maximum connectivity. Exponential localization can be further characterized by introducing the localization length being the latter solution of the equation ρ 2n = exp(−n/L). Using Eq. (33) in Eq. (38) one obtains L ≈ 2.08. An effective measure of the localization strength induced by the network connectivity can be obtained by comparing L with the localization length L v induced by a local attractive potential, namely −V , acting on a specific island of a linear chain (see Appendix B). The comparison is performed by solving the equation L = L v with respect to the variable v, being v = V /K. When this general procedure is applied to the present case, one obtains v ≈ 0.486, implying a rather strong attractive potential −0.486 K. In this way, the localizing action of the side-coupled island has been measured in terms of a sort of effective potential.

C. Star-like network
Star-like networks belong to a special class of tree-like structures (loop-free graphs) in which a central island is connected to several branches. In these systems, the degree of the central site is typically much higher than the degree of any other network's site and acts as localization center for particles. Apart from intrinsic motivations, star-like topology may describe parts of a more complex network and this observation further motivates the interest in studying these systems. The wavefunction of N T Cooper pairs confined to move on a star-like network with p branches (Fig. 3) obeys the following equations: with n ≥ 1 and j ∈ {1, ..., p}. Using exponential notation, i.e. ψ 0 = A 0 exp(iθ) and ψ n , we get the following relations: with the boundary condition A (j) 0 = A 0 . Numerical diagonalization of the problem shows that the EGE presents exponential localization. The latter statement can be easily proven by using the ansatz A (j) n = A 0 ρ n with ρ < 1 and n ≥ 1. After the substitution in Equations (41) and (42), we obtain the relations which admit the following solutions: A 0 is fixed by imposing the normalization condition: Using the aforementioned results and assuming N ≫ 1, the particles fraction allocated in the network site labelled by n = 0 is obtained: Interestingly, taking the limit p → ∞ one observes N 0 → N T /2. Moreover, the particles occupation along the j-th branch is easily obtained by using the relation N (j) n = N 0 ρ 2n . The exponential localization of the particles distribution allows to use Eq. (38) to extract the localization length in the form: while the parameter v introduced at the end of Sec. IV B takes the value: When Equations (49) and (50) are specialized to the case p = 3, one obtains L ≈ 1.44 and v ≈ 0.707. The latter estimates clearly show the relevance of connectivityinduced localization phenomena. Cooper pairs wavefunction is expected to present the same profile exhibited by the superconducting order parameter. This expectation is indeed confirmed and we have found that the wavefunction of particles confined on a three-legs star-like network closely follows the superconducting order parameter profile derived in Ref. [35] in close vicinity of the critical temperature. Following Ref. [35], the order parameter profile deviates from a pure exponential when the temperature is lowered further. This deviation has no direct correspondence with the bosonic wavefunction since temperature is not included in the present treatment. Despite the latter observation, within the Ginzburg-Landau approach proposed in Ref. [35], a temperature lowering promotes a relevant influence of the nonlinear terms affecting the order parameter equations. As will be clearly explained in Sec. V C, a similar role is here played by the repulsive interaction, whose influence increases with the dimensionless parameter ξ, which is proportional to the bosons number N T . Thus a full agreement between the distinct approaches implies that the parameter N T has to be interpreted as the average Cooper pairs number, being the latter a temperaturedependent quantity; it is vanishing above the normalsuperconductor transition temperature, while it increases when the system is progressively cooled below the critical temperature. According to Ref. [35], the order parameter localization characterizes the nucleation of the superconductivity in close vicinity of the critical temperature, while a spreading of the superconductivity to the whole system is observed when the temperature is lowered further. The mentioned behavior suggests that exponential localization of Cooper pairs characterizes the first insurgence of the superconductivity in star-like networks, while the initial localization is not able to drive the system towards a global insulating phase.

D. Flower-like network
A flower-like network 52 (Fig. 4) is obtained by linking a central island with M lateral sites, so that N = M + 1 represents the total number of islands. Despite flower-like networks can be seen as a subcase of star-like structures, they deserve a separate treatment. Indeed, results presented in Sec. IV C are here unapplicable because they are based on the assumption that the number of islands forming each branch of the star-like structure is much greater than one. The wavefunction of particles constrained to a flower-like network is obtained by solving the following equations: where ψ c = √ N c exp(iθ c ) is related to the central island, while ψ j = N j exp(iθ j ) (j ∈ {1, ..., M }) are related to the M lateral sites. Using exponential representation, Equations (51)-(52) take the following equivalent form: By imposing the stationary condition θ c = θ j = θ(t) witḣ θ(t) = ω, we getṄ c =Ṅ j = 0 along with the additional relations: The above equations are easily solved so that, also considering the normalization condition, we get the islands' occupation: Moreover, λ = √ M so that the ground state energy is given by µ 0 = −Kλ. Thus, we have demonstrated that the fraction of particles belonging to the central island does not depend on the number M of lateral sites.

V. SELECTED APPLICATIONS OF THE THEORY: INTERACTING CASE (U > 0)
The inclusion of repulsive interaction in the model produces in general non-trivial changes to the particles distribution on the network. However, there exist analytically solvable network models which can be used to validate the numerical procedure described in Sec. III A. These models are presented in Sec. V A and V B. Once the iterative procedure has been validated, it is used to study the interacting star-like network model (Sec. V C), which is not analytically solvable.

A. Interacting flower-like network
Hereafter, we study a flower-like network obtained by linking a central island with M lateral sites, so that N = M + 1 represents the total number of islands. Under the assumption of non-vanishing repulsive interaction (i.e., U > 0), the wavefunction of particles constrained to a flower-like network is obtained by solving the following equations: where ψ c is related to the central island, while ψ j (j ∈ {1, ..., M }) are related to the M lateral sites. Equations (61) and (62) can be solved by using the following ansatz: along with the normalization condition N c + M N l = N T . Using the ansatz with the conditionθ = ω in Equations (61) and (62) provides the following relations: with λ = ω/K and u = U/K. The non-linear problem written in Eq. (65)-(66) can be rewritten in terms of the auxiliary variable z = A c /A l as follows: Once the solution of the quartic equation in Eq. (67) is known, physical observables are obtained by using the following relations: In particular, the non-interacting case (see Equations (59) and (60)) is recovered by setting ξ = 0 in Eq. (67) whose acceptable solution under this condition is given by z = √ M . The model can be solved either by direct solution of Eq. (67) or by using the iterative procedure presented in Sec. III A. The outcome of the model is presented in Fig. 5, where a flower-like network with M = 10 lateral islands linked to a central site is studied. In particular, the particles fraction N c /N T allocated on the central site as a function of ξ is presented in left panel of Fig.  5. Repulsive interaction tends to reduce the particles accumulation on a single site and for this reason the occupation of the central island is a decreasing function of ξ. Treating Eq. (67) as described in Appendix C and considering the small-interaction limit (ξ ≪ 1), it is possible to demonstrate that the central island occupation is described by the following approximate relation: as long as the interaction parameter ξ is sufficiently small. Equation (71) provides an excellent approximation of the curve shown in Fig. 5 (left panel) capturing the linear decreasing of the central island occupation for ξ 0.5. Considering the limit M ≫ 1 in Eq. (71), one observes that the competition between connectivity-induced localization and repulsive interaction is controlled by the parameter ξ/ √ M . Moreover, increasing the repulsive interaction also produces a corresponding increasing of the ground state energy µ 0 = −Kλ (see right panel of Fig. 5). We also observe that distinct solution methods show a good agreement and provide a consistent picture of the system's properties.

B. Star-like network with repulsive interaction at the central island
Hereafter, we focus on a star-like network in which the effect of the repulsive interaction between particles manifests only at the central island, while it is negligible elsewhere. This condition can be experimentally realized by considering a central island with reduced volume compared to that of the other network's islands. Thus, the wavefunction of N T Cooper pairs confined to move on a star-like network with p branches and repulsive interaction at the central island obeys the following equations: with n ≥ 1 and j ∈ {1, ..., p}. Using the notation n exp(iθ), we get the following relations: with u = U/K and the boundary condition A (j) 0 = A 0 . Solution of the non-linear problem is obtained by using the ansatz A (j) n = A 0 ρ n with ρ < 1 and n ≥ 1. After the substitution in Equations (74) and (75), we obtain the following relations A 0 can be written in terms of ρ by using the normalization condition. Accordingly, also assuming N ≫ 1, we obtain: while ρ obeys the quartic equation: with ξ = uN T and ρ < 1. A close inspection of Eq.
(79) clearly shows that repulsive interaction effects can be mitigated by the network topology (in graph theory sense). Indeed, by increasing the number of network's branches (i.e. for p ≫ 1) the quantity ξ/(p − 1) 2 becomes negligible and, under this condition, the system shows the tendency to behave as if the interaction at the central island were absent. The above observation corroborates the idea that highly connected islands are affected by topology-induced effective potentials and play the role of attractive centers of the particles' density.
Once the quartic equation has been solved, the occupation of the central island can be computed according to the relation: while the occupation of islands belonging to a generic branch is given by N (j) n = N 0 ρ 2n . The model can be solved either by direct solution of Eq. (79) or by using the iterative procedure presented in Sec. III A. The outcome of the model is presented in Fig.  6, where a three-branches (i.e., p = 3) star-like network with localized interaction is studied. In particular, the particles fraction N 0 /N T allocated on the central site as a function of ξ is presented in left panel of Fig. 6. As already observed before, the repulsive interaction tends to reduce the particles accumulation on a single site and produces a corresponding increasing of the ground state energy µ 0 = −Kλ (see right panel of Fig. 6). Furthermore, the exponential ansatz used to derive Eq. (79) is validated by using the iterative procedure presented in Sec. III A. Using approximate methods described in Appendix C and considering the small-interaction limit (ξ ≪ 1), the central island occupation of a p-branch star-like network with localized interaction is described by the following relation: Equation (81) provides an excellent approximation of the curve shown in Fig. 6 (left panel) for ξ 0.5. When Eq. (81) with p ≫ 1 is considered, one observes that the localization phenomena are controlled by the parameter ξ/ √ p. Interestingly, Equations (71) and (81) take the same form as the limit M ≫ 1 and p ≫ 1, respectively, are considered. The latter observation suggests that, within the considered limit, the occupation of the central island is rather insensitive to the number of islands forming a single branch.

C. Interacting Star-like network
The wavefunction of N T Cooper pairs confined to move on a star-like network with p branches and repulsive interaction obeys the following equations: with n ≥ 1 and j ∈ {1, ..., p}. Within the same notation used in Sec. V B, we get the following relations: with u = U/K and the boundary condition A (j) 0 = A 0 . Solution of the non-linear problem cannot be obtained by using the ansatz A (j) n = A 0 ρ n and thus we have to use the numerical procedure described in Sec. III A. In the following, we focus on the interacting star-like network obtained by linking a central site with three branches (p = 3), each of which made of the connection of 25 islands. Thus, N = 76 represents the total number of islands. The particles distribution of the system with interaction parameter ξ = 1.2 is shown in Figure 7 (left panel). When the aforementioned distribution is compared with that pertaining to the non-interacting case (i.e. ξ = 0), two distinctive features clearly appear. Indeed, interaction produces a non-exponential behavior of the particles distribution accompanied by a lowering of the particles fraction N 0 /N T allocated at the central island. The latter tendency also populates islands belonging to the system's branches. This behavior is affected by the interaction strength. In particular, the particles fraction allocated at the central island is a decreasing function of the interaction parameter ξ (see Figure 7, middle panel), which is an obvious signature of the repulsive interaction between particles. A comparison of the N 0 /N T vs ξ curve with the corresponding shown in Fig. 6 (left panel) evidences a similar behavior as long as the interaction parameter ξ is smaller than ∼ 0.2, while a marked deviation is observed for ξ > 0.2. The latter observation suggests that, as long as ξ 0.1, the N 0 /N T vs ξ curve of an interacting star-like network is well described by Eq. (81), which has been derived considering localized interaction (see Sec. V B). The presence of a finite repulsive interaction between particles is also responsible for an enhancement of the ground state energy of the system compared to the non-interacting case. This statement is easily verified by looking at the right panel of Figure 7 where the dependence of λ on ξ is presented. Indeed, λ is a decreasing function of the interaction strength ξ and thus the ground state energy µ 0 = −Kλ increases due to the repulsive interaction. The localization length as a function of the interaction parameter ξ is reported in the inset of Fig. 7 (right panel). The full-line curve represents the localization length L of a star-like network with interaction localized at the central island (see Sec. V B), while symbols represent the localization length of the interacting star-like network described in the present section. In the latter case, localization length has been extracted by considering the exponential decay of the particles distribution around the central part of the system branches. Localization length increases with ξ evidencing the presence of interaction-driven delocalization phenomena. For ξ > 0.2, the localization length of an interacting starlike network (symbols) grows faster as the interaction increases than the corresponding behavior observed in the case of localized interaction (full line curve). Differences are less relevant for small values of the interaction (ξ < 0.2). From the physics viewpoint, we observe that interacting star-like networks and star-like networks with interaction localized at the central island present a very similar behavior as long as the central island population includes a relevant fraction of the total number of particles (i.e., for small values of ξ). Under this condition, considering an interacting star-like network, islands belonging to the branches of the star are essentially depopulated and particles belonging to this part of the system behave as if the repulsive interaction were absent. Interaction always works at the central island. Intuition clearly indicates that the equivalence between the two models is broken as the interaction parameter ξ increases. Indeed, under this condition, the central island occupation starts to be inhibited favoring the particles migration towards the system's branches. The increased number of particles belonging to the branches makes the repulsive interaction more effective.

VI. DISCUSSION AND CONCLUSIONS
We have presented a generalized Feynman's model to describe the inhomogeneous Cooper pairs distribution in tree-like networks of weakly-coupled superconducting islands. The theory allows the description of the interplay between particles interaction and network connectivity in graph-like artificial structures also evidencing interesting connections with the spectral graph theory. In these systems, the network connectivity defines effective potentials which are able to induce particles localization. Indeed, network's islands with higher connectivity act as localization centers, while the repulsive interaction works in the opposite way. Using a number of representative examples, it has been demonstrated that the particles fraction localized at a given system's site depends on the network's topology and, as expected, it is lowered by the repulsive interaction. In particular, we have demonstrated that the Cooper pairs fraction localized at the central site of a non-interacting flower-like network or at the center of a star-like network cannot exceed fifty percent of the total particles number N T , irrespectively of any other detail. Moreover, the inclusion of repulsive interaction in star-like networks leads to depopulation of the central site and induces non-exponential localization phenomena. Cooper pairs wavefunction is expected to present the same profile exhibited by the superconducting order parameter in view of the close connection between the aforementioned quantities. This expectation is indeed confirmed and we have found that the wavefunction of particles confined on a three-legs star-like network closely follows the superconducting order parameter profile derived in Ref. [35]. These findings suggest interesting analogies between disordered superconducting systems and tree-like networks of weakly-coupled superconducting islands. In both the mentioned systems localization phenomena are accompanied by an enhancement of the critical temperature of the superconducting transition, being the latter property appealing for possible technological applications. Furthermore, since both Josephson arrays and opti- cal lattices are experimental realizations of the Bose-Hubbard model, we expect that atomtronics [53] platforms can be also used to test the localization phenomena reported in this work. distributed according to the relation N (α) n = N 0 ρ 2n , being N 0 the occupation of the central island and with v = V /K > 0. The confining properties of the local potential can be characterized by introducing the localization length L v according to the relation: being the latter solution of the equation exp(−n/L v ) = ρ 2n . A closer inspection to Equations (B3) and (B4) shows that L v is a decreasing function of v. Moreover, particles localization increases when L v decreases. Despite 4th order algebraic equations admit exact solutions (Ferrari-Cardano method), it is here useful searching for approximate roots of Eq. (67) and (79). Approximate roots are expressed via relatively simple mathematical expressions and can be used to understand the interplay between topological confinement effects and repulsive interaction. Both the aforementioned equations can be written in the general form f (s) = 0 with f (s) = s 4 + αs(s 2 − 1) − β (C1) and β > 0. Moreover, the parameter α is proportional to the dimensionless interaction strength ξ. We are interested in finding a real and positive solution. If α is sufficiently small, the exact solution exhibits a small deviation from the solution of the non-interacting problem s 0 = 4 √ β. Let us call δs ≪ s 0 the aforementioned deviation and write the equation f (s 0 + δs) = 0. Disregarding cubic and quartic terms in δs, we obtain the following equation: (6s 2 0 + 3αs 0 )δs 2 + (4s 3 0 + 3αs 3 0 − α)δs + αs 0 (s 2 0 − 1) = 0.
Let us note that δs → 0 for vanishing interaction (i.e., for α → 0). Thus, the approximate solution of the equation f (s) = 0 is given by s 0 + δs.