format: text

text:
# Critical current of a S/F-I-N/S tunnel structure in a parallel magnetic field

## I Introduction

The problem of superconductivity in the presence of spin splitting is currently heavily investigated. In fact, this question was already addressed experimentally and theoretically a while ago for the simplest case of a very thin superconducting film in a parallel magnetic field. For such geometry the quasiparticle energy states in a superconductor are splitted due to the interaction of the field with spin magnetic moments of the quasiparticles, while the orbital effects are negligible. This splitting has been earlier demonstrated by extensive studies of tunnel junctions formed by thin superconducting films (see, e.g., \[1-3\] and references therein).

A more sophisticated case is the proximity effected hybrid systems consisting of a superconductor (S) and a ferromagnet (F). In a ferromagnet, the difference in the energy between spin-up and spin-down bands is an internal property, and that is the physical reason why a singlet Cooper pair, injected from a superconductor into a ferromagnet, acquires a finite momentum. The resulting proximity induced superconducting state of the F layer is qualitatively different from the zero-momentum state: it is spatially inhomogeneous and the order parameter contains nodes where the phase changes by $`\pi `$ (see recent review  and references therein). Particularly, transport properties of the S/F structures have turned out to be quite unusual. Experiments that have been performed on SFS weak links \[5-7\] and SIFS tunnel junctions  directly prove the so-called $`\pi `$-phase superconductivity. Planar tunneling spectroscopy also reveals a $`\pi `$-phase shift in the order parameter, when superconducting correlations coexist with the ferromagnetic order . An inverse phenomenon, the modification of the BCS density of states in mesoscopic S layer under the influence of a ferromagnet is observed in tunnel spectroscopy experiment  as well. For S/F-I-F/S tunnel structures with very thin F layers, it was also predicted  that one can, if there is a parallel orientation of the F layers magnetization, turn the junction into the $`\pi `$-phase state with the critical current inversion, or, if there is an antiparallel orientation of the F layers internal fields, save the 0-phase junction state but enhance the tunnel current. Experimentally, however, the enhancement of the dc Josephson current has not been detected until now. One of the main problems here is the requirement to fix one layer magnetization of the sample, while preserving a possibility to reorient the magnetization of another layer by a weak magnetic field.

In this report we draw attention to another experimental set-up where, as we show below, there is the critical current enhancement, as well as other new effects of the superconductivity in the presence of spin splitting, while the experimental set-up admits of simple variations of an applied field direction. Namely, we consider the dc Josephson effect for tunnel structure where there is spin splitting for two reasons - due to an external magnetic field and due to an internal exchange field of the F layer. The ground state and critical current dependence on the applied field is analyzed for a system of proximity coupled bulk superconductor with thin ferromagnetic and thin normal nonmagnetic metals separated by an insulating barrier (S/F-I-N/S tunnel junction). Based on a microscopic theory of the superconducting proximity effect for S/F and S/N bilayers we demonstrate that the S/F-I-N/S contact with strong enough magnetism of the F layer is in a ground state with $`\pm \pi /2`$ superconducting phase difference on the electrodes. In such a contact low external magnetic field reverses the critical current sign (the $`\pi `$-state of the junction) if the fields (external and internal) have parallel orientation. For the antiparallel fields alignment the critical current behavior is much more complicated: a monotonic increase of the external field causes an enhancement of the critical current, which is then changed by a reverse of current sign.

The article is organized as follows. In Sec. II we outline the model of the junction. In Sec. III the outcomes of the theory of proximity coupled S/F and S/N bilayers have been generalized to include external magnetic field effects. The main results of the paper are presented in Sec.IV. In Sec. V we evaluate the approximations that have been made in the theory and briefly discuss a possible realization of the S/F-I-N/S systems to explore the predicted effects.

## II The model of a junction

The system we are interested in is the S/F-I-N/S tunnel structure in a parallel magnetic field shown in Fig.1. We assume that the first (left) electrode of the junction is formed by a proximity coupled S/F bilayer while the second (right) is formed by a proximity coupled S/N bilayer. The transparency of the insulating layer is small enough to neglect the effect of a tunnel current on the superconducting state of the electrodes. The S/F and S/N interfaces may have arbitrary finite transparency, but it is large compared to the transparency of the junction barrier. The F layer will be treated as a single domain ferromagnet. We also expect the F layer magnetization to be aligned parallel to the interface, so that it does not create a spontaneous magnetic flux penetrating into the S layer. The general relation between the Curie temperature of the F metal and the superconducting transition temperature of the S metal is considered. Transverse dimensions of the junction are supposed to be much less then the Josephson penetration depth, so that (a) all quantities depend only on a single coordinate $`x`$ normal to the interface surface of the materials, and (b) a flux quantum can not be trapped by the junction.

We assume that the applied field penetrates and mostly concentrates into the ”weak” F-I-N part of the junction. We will also neglect the effect of the magnetic field on electron orbits as compared with the effect of the field on electron spins. This assumption is well satisfied for films thinner than their superconducting coherence lengths $`d_{F(N)}<<\xi _{F(N)}`$ . Here $`d_{F,N}`$ is the thickness of the F and N layers, respectively; $`\xi _F\mathrm{~}\sqrt{D_F/2H_e}`$ and $`\xi _N\mathrm{~}\sqrt{D_N/2\pi T_C}`$ are superconducting coherence lengths for F and N metals; $`D_{F,N}`$ is a diffusion coefficient in the F, N metal; $`H_e`$ is the exchange energy. For films of such thicknesses one can also be sure that there is a fairly uniform penetration of a parallel field. (The evaluation of these approximations will be made below, in Sec. V). We consider the ”dirty” limit for all metals, i. e., we assume that $`l_S<<\xi _S`$ and $`l_{F(N)}<d_{F(N)}`$, and that the condition $`H_e<\tau ^1`$ is fulfilled , where $`l_{S,F,N}`$ is the electron mean free path for the S, F, and N layers, respectively; $`\xi _S`$ is the superconducting coherence length in the S metal; and $`\tau `$ is the momentum relaxation time. (We have taken the system of units with $`\mathrm{}=\mu _B=k_B=1`$). We restrict ourselves to the case when the spin-orbit scattering is absent and spin ”up” and spin ”down” electron subbands do not mix with each other. In equations below, we will use the effective coherence length $`\xi \mathrm{~}\sqrt{D_F/2\pi T_C}`$ for the F layer, thus providing the regular crossover to both limits $`T_C>>H_e0`$ and $`H_e>>T_C`$ (see Refs. 14, 15 for details); the relation on film thickness is read as $`d_F<<\xi _F`$ or $`\xi `$ , depending on the relation between $`T_C`$ and $`H_e`$.

## III S/F and S/N sandwiches in the magnetic field

As a first step, we should generalize results of the theory for a proximity coupled S/F (S/N) bilayer of a massive, $`\xi _S<<d_S`$ , superconductor and a thin ferromagnetic (nonmagnetic) metal taking into account the influence of the applied magnetic field $`H`$ . The most direct way to do it is to use the results for a S/F proximity coupled bilayer (Refs. ).

As is well known, the properties of ”dirty” metals can be described by the Usadel equations . Following , we introduce the modified Usadel functions $`\mathrm{\Phi }`$ , used to take into account the normalized confinement on the normal $`G`$ and anomalous $`F`$ Green’s functions. Due to small thickness of the normal (magnetic and nonmagnetic) metal, the proximity effect problem can be reduced to a boundary problem for the S layer and a relation for determining modified Usadel function $`\mathrm{\Phi }`$ of a nonsuperconducting layer (see Refs.  for details). The boundary conditions at the SF (SN) interface are modelled by two parameters: $`\gamma _M`$ and $`\gamma _B`$ . The $`\gamma _M`$ plays the role of an effective pairbreaking parameter near the SF (SN) interface and its value is mainly determined by electron densities of S and F (N) metals. A large value of $`\gamma _M`$ corresponds to high density of quasiparticles in the F (N) metal compared to that in the S near the SF (SN) boundary. In this case the diffusion of quasiparticles into the superconductor leads to a strong suppression of the order parameter in the S region. The scale at which superconductivity reaches the bulk value is $`\xi _S`$ from the SF (SN) interface. In the opposite case, $`\gamma _M<<1`$ , the influence of the F (N) layer on the superconductivity of the S metal is weak, and even vanishes if $`\gamma _M=0`$. The parameter $`\gamma _B`$ $`\mathrm{~}`$ $`R_B`$ , where $`R_B`$ is the SF (SN) boundary resistance, describes the electrical quality of the SF (SN) interface. The case $`\gamma _B<<1`$ corresponds to a vanishing interlayer resistance, i.e., the S and F (N) metals are in good electric contact. The case $`\gamma _B>>1`$ corresponds to low transparency of the SF (SN) barrier. There is also one more parameter, $`H_e`$ , that is the value of the exchange spin splitting in the F layer. In the general case, the problem needs self-consistent numerical calculations of the Usadel equations. However, at some simplifying assumptions, the model admits of analytic solutions that hold all new physics we are interested in.

For F (N) films thinner than their superconducting coherence lengths, $`d_{F(N)}<<\xi _{F(N)}`$ , the critical current of the S/F-I-N/S tunnel contact can be represented through the S layer Green’s functions (see Eq. (3), below). The explicit expressions for these functions can be obtained on the basis of the microscopic theory of proximity effects for S/F and S/N bilayers for a weak, $`\gamma _M<<1`$ , and a strong, $`\gamma _M>>1`$ , proximity effect. We generalize these results assuming the external magnetic field penetrates into the F and N layers. Namely, in the case of a bulk S layer proximity coupled with a thin F layer with a weak proximity effect $`\gamma _M<<1`$ and an arbitrary value of the boundary resistance $`\gamma _B`$ (quite a realistic experimental case) for the function $`\mathrm{\Phi }_{SF}`$ we have the analytical solution (see Eq. (14), Ref. 14 or Eq. (10), Ref. 15):

$$\mathrm{\Phi }_{SF}(\omega )\mathrm{\Phi }_{SF}(\omega ,x=0)=\mathrm{\Delta }_S\left(1\frac{\gamma _M\beta \omega _F}{\gamma _M\beta \omega _F+\omega A_F}\right)\text{ ,}$$
(1)
where $`\mathrm{\Delta }_S`$ is the absolute value of the superconducting order parameter in the bulk of the S layer, $`\omega _F=\omega +i\sigma (H_e+H),`$ $`\omega \omega _n=\pi T(2n+1),n=0,\pm 1,\pm 2,\mathrm{}`$ are Matsubara frequencies , $`\sigma =\pm 1`$ is the spin sign;$`\beta ^2=\left(\omega ^2+\mathrm{\Delta }_S^2\right)^{1/2}/\pi T_C,`$ and $`A_F=\left[1+\gamma _B\omega _F\left(\gamma _B\omega _F+2\omega /\beta ^2\right)/(\pi T_C)^2\right]^{1/2}`$. The sign of the external magnetic field depends on the relative orientation to the internal exchange field. To be definite, we will suppose that $`H_e>0`$ and there will be the case of the parallel external fields orientation if $`H>0`$, while for $`H<0`$ the fields will have antiparallel orientation.

For the S/F bilayer with a strong proximity effect $`\gamma _M>>1`$ and free value of the boundary resistance $`\gamma _B`$ our calculations yield (see Eq. (12) Ref. 15):

$$\mathrm{\Phi }_{SF}(\omega ,x=0)\mathrm{\Phi }_{SF}(\omega )=B(T)\left(\pi T_C+\gamma _B\omega _F\right)/\gamma _M\omega _F$$
(2)
where $`B(T)=2T_C[1(T/T_C)^2][7\zeta (3)]^{1/2}`$ and $`\zeta (3)`$ is the Riemann $`\zeta `$ function.

Qualitatively similar results are obtained for the S/N bilayer. If the influence of the N layer on the superconductivity of the S metal is weak, $`\gamma _M<<1`$ and the value of interface resistance is arbitrary, for $`\mathrm{\Phi }_{SN}(\omega )`$ we have the expression (1) where the substitutions should be made: $`\omega _F\omega _N=\omega +i\sigma H`$ and $`A_FA_N=\left[1+\gamma _B\omega _N\left(\gamma _B\omega _N+2\omega /\beta ^2\right)/(\pi T_C)^2\right]^{1/2}`$. In the opposite case of a strong proximity effect $`\gamma _M>>1`$ , for $`\mathrm{\Phi }_{SN}(\omega )`$ we have the relation (2) where $`\omega _F`$ is changed by $`\omega _N`$ , and, as before, there is an arbitrary value of the SN boundary resistance.

Below we will suppose for simplicity, that the interface parameters involved , $`\gamma _M`$ and $`\gamma _B`$ , have the same values for both SF and SN interfaces. An extension to different SN and SF boundary parameters is straightforward, and will be made towards the end of the Sec. IV. An additional physical approximation is also the assumption that the spin discrimination by the interfaces is unimportant. Generalization for interfaces with different transmission probabilities for quasiparticles is obvious and we do not present it here.

## IV Critical current

The critical current of a S/F-I-N/S tunnel contact can be represented in the form (see, e.g., Ref. )

$$j_C=(eR_N/2\pi T_C)I_C=(T/T_C)Re\underset{\omega >0,\sigma }{}\{G_S\mathrm{\Phi }_S/\omega \}|_F\{G_S\mathrm{\Phi }_S/\omega \}|_N\times $$
(3)
$`\{[1+2\omega _FG_S(\gamma _B/\pi T_C)+\omega _F^2(\gamma _B/\pi T_C)^2]_F\times [1+2\omega _NG_S(\gamma _B/\pi T_C)+\omega _N^2(\gamma _B/\pi T_C)^2]_N\}^{1/2}`$where RN is the resistance of the contact in the normal state; the subscript F (N) labels quantities referring to the left ferromagnetic and right nonmagnetic electrodes. Here we used $`G_S=\omega /(\omega ^2+\mathrm{\Phi }_S\stackrel{~}{\mathrm{\Phi }}_S)^{1/2}`$ , where $`\stackrel{~}{\mathrm{\Phi }}_S(\omega ,H_{eff})=\mathrm{\Phi }_S^{}(\omega ,H_{eff})`$ and $`H_{eff}`$ is $`(H+H_e)`$ or $`H`$ .

### A Strong proximity effect

Let us first present results for the case of a strong proximity effect, $`\gamma _M>>1`$ , and a vanishing interface resistance, $`\gamma _B=0`$. If there is a strong suppression of the order parameter near S/F (S/N) boundary then, using the expressions for $`\mathrm{\Phi }_{SF}`$ and $`\mathrm{\Phi }_{SN}`$ , we can recast the Eq. (3) in the form:

$$j_C(H,H_e)2(T/T_C)B_M^2(T)\underset{\omega }{}\frac{\omega ^2H(H_e+H)}{\omega ^2\{[(\omega ^2H(H_e+H)]^2+\omega ^2(H_e+2H)^2\}}$$
(4)
where $`B_M(T)=B(T)\pi T_C/\gamma _M`$ , and in proceeding these relations, we have taken into account that the value of $`\mathrm{\Phi }_{SF(SN)}`$ is small, $`\mathrm{\Phi }_{SF(SN)}\mathrm{~}\gamma _M^1<<1`$. Performing the substitution $`(H_e+H)H_{eR}`$ and $`H\pm H_{eL}`$ in Eq. (4) , one can obtain the results for the (S/F)L-I-(F/S)R tunnel structure for (anti)parallel orientation of the $`H_{eR}`$ and $`H_{eL}`$ fields (see Eq. (15) in Ref. 15 and Eq. (5) in Ref. 12). If $`H_e=0`$ and $`H=0`$ one can also restore the result for S/N-I-N/S junction (see, e.g., Eq. (31) of Ref. 17 and Eq. (46) of Ref. 18). In the general case, as it is seen from Eq. (4), in the geometry with the parallel fields by increasing the external magnetic field one can inverse the supercurrent sign of the S/F-I-N/S tunnel structure. In other words, at large enough applied field, $`H(H_e+H)>>\omega ^2\mathrm{~}(\pi T_C)^2`$ , the critical current changes its sign in comparison with that when $`H0`$ . In particular, from Eq. (4) we obtain in the limit of large field values:

$`j_C(H,H_e)\text{ }\mathrm{~}B^2(T)/\gamma _M^2{\displaystyle \underset{\omega }{}}1/\{\omega ^2H(H_e+H)\}\text{ }\mathrm{~}\text{ }(T_CT)^2/\{\gamma _M^2H(H_e+H)\},`$i. e. , a crossover of the S/F-I-N/S junction from the 0-phase state to the $`\pi `$-phase state takes place.

The most interesting results have been obtained for the geometry with the opposite field orientation: $`H<0`$. Simple substitution into Eq. (4) gives for $`H=H_e`$ that $`j_C(H_e,H_e)=j_C(0,H_e)`$ , while if $`H=H_e/2`$ we directly obtain that $`j_C(H_e/2,H_e)>j_C(0,H_e),`$ $`j_C(H_e,H_e)`$ . Comparing, we also find that $`j_C(H_e/2,H_e)j_C(H_e/2,H_e)`$ and $`j_C(H_e,H_e)j_C(H_e,H_e)`$ . As it follows from Eq. (4), for large enough values of the external field, the critical current again inverses its sign and the crossover of the S/F-I-N/S junction from the 0-phase state to the $`\pi `$-phase state takes place. However, while increasing the field in the interval $`0<|H|<H_e`$ , a nonmonotonic current versus field dependence is observed. The main features of this dependence are shown on Fig. 2, where the Josephson current amplitude of the S/F-I-N/S junction versus the applied magnetic field is given for the case of a strong proximity effect, $`\gamma _M=10`$ and $`\gamma _B=0`$ , and different exchange field intensity. (Here and below, dashed curves are used for the results of the numerical calculations for the parallel fields orientation and full curves are used for the results in the opposite geometry). It can be seen, that the dc current greatly depends on relative orientation of the fields.

To reveal the underlying physics behind the critical current behavior on field, let us consider a ground state and a phase shift across the contact at zero external field. As follows from Eq. (2), the modified Usadel function $`\mathrm{\Phi }_{SF}(\omega )`$ can be written in the form (the case $`\gamma _B=0`$)

$`\mathrm{\Phi }_{SF}(\omega )=B(T)(\pi T_C/\gamma _M){\displaystyle \frac{\omega iH_e}{\omega ^2+H_e^2}}`$Taking into account that the typical value $`\omega `$ ~$`\pi T_C`$ , one can see that in the limit $`H_e>>\pi T_C`$ the unified correlation function of the S/F bilayer acquires an additional $`\pm \pi /2`$ phase shift in comparison with the similar function for the S/N bilayer. When the external field is applied, $`H0`$ and $`\gamma _M>>1`$ , $`H_e>>\pi T_C`$ , we may evaluate the critical current as

$$j_C\text{ }\mathrm{~}\text{ }Re\underset{\omega }{}F_{SF}(\omega )F_{SN}(\omega )\text{ ~}Re\underset{\omega }{}\mathrm{\Phi }_{SF}(\omega )\mathrm{\Phi }_{SN}(\omega )/\omega ^2\text{ ~}$$
(5)
$`\text{~}Re{\displaystyle \underset{\omega }{}}{\displaystyle \frac{\mathrm{exp}\{\pm i[\pi /2+\mathrm{arctan}(H/\omega )]\}}{\omega ^2(\omega ^2+H^2)^{1/2}}}`$i.e., for $`H0`$ , the value of the supercurrent may be low enough - the exchange field blocks a contact. So, the S/F-I-N/S contact with strong enough magnetism of the F layer can have a ground state with the superconducting phase difference on the electrodes about $`\pi /2`$. Note, that this mechanism of the $`\pi /2`$ junction realization should not be confused with the $`\pi /2`$-phase state induced by current fluctuations in the junction plane, the situation that has been treated earlier in Ref.  for SFS structures with a thick F layer. We postpone a more detailed description of the $`\pi /2`$ state of S/F-I-F/S and S/F-I-N/S hybrid systems to a forthcoming publication , but here discuss the effects of the external magnetic field. As one can see from expression (5), if the contact is in magnetic field, the additional phase shift of the Usadel function $`\mathrm{\Phi }_{SN}(\omega )`$ has been induced. Depending on the relative orientation, the applied field can prompt the total phase shift across the contact to a $`\pi `$-state of the junction or to a $`0`$-phase state. For the latter case one can expect an essential enhancement of the critical current amplitude. Such type of behavior we can see on Figs. 2-7.

As seen on Fig. 2, a low external magnetic field reverses the critical current sign if the fields (external and internal) have parallel orientation. For antiparallel fields alignment the critical current behavior is much more complicated. In some interval of the applied field the enhancement of the dc Josephson current takes place in comparison with the case of $`H=0`$, while the 0-phase state of the junction is held. By increasing the magnetic field one can turn the current though a few extremum and then adopt the junction to the $`\pi `$-phase state superconductivity with the opposite current direction.

The results above are obtained for the case of a vanishing interface resistance, $`\gamma _B=0`$. After simple but cumbersome algebra, one can also obtain expressions for $`\gamma _M>>1`$ , and $`\gamma _B0`$ . We will not extract here the formulae but illustrate our calculations on Fig. 3, where the dc current dependence versus the applied magnetic field is shown for the S/F-I-N/S junction with a strong influence of the F (N) layer on superconducting layers , $`\gamma _M=10`$ , and different interface transparency: $`\gamma _B=1÷10`$. A critical current amplitude of the junction depends greatly on bilayers parameters, and on the interface transparency $`\gamma _B`$ , in particular, decreasing if the transparency decreases, as it is seen on Fig. 3.

### B Weak proximity effect

Using the expression (1), after simple transformations one can obtain for the dc Josephson current in the case of vanishing interface resistance, $`\gamma _B=0`$ , and a weak proximity effect, $`\gamma _M<<1`$:

$$j_C(H,H_e)2(T/T_C)Re\underset{\omega }{}\frac{\mathrm{\Delta }_S^2}{\mathrm{\Omega }^2}\times \{1+2i\sigma \omega \gamma _M\beta (H_e+2H)\mathrm{\Omega }^2$$
(6)
$`\gamma _M^2\beta ^2[H^2+(H_e+H)^2+4\omega ^2\mathrm{\Omega }^2H(H_e+H)]\mathrm{\Omega }^2\}^{1/2}\text{ , }`$where $`\mathrm{\Omega }^2=\mathrm{\Delta }_S^2+\omega ^2`$ . In proceeding these relations, we have taken into account that $`\gamma _M\beta `$ is small, $`\gamma _M\beta <<1`$. Simple analysis shows that $`j_C(0,H_e)=j_C(H_e,H_e)`$ as it should be from the symmetry of the fields geometry. If $`\gamma _M\beta H_e<<1`$ , one can also obtain that

$`j_C(H_e/2,H_e)2(T/T_C){\displaystyle \underset{\omega }{}}{\displaystyle \frac{\mathrm{\Delta }_S^2}{\mathrm{\Omega }^2}}\{1(\gamma _M\beta H_e)^2{\displaystyle \frac{\mathrm{\Delta }_S^2\omega ^2}{2\mathrm{\Omega }^2}}\}^{1/2}>j_C(0,H_e)`$The dependence of the Josephson current amplitude of S/F-I-N/S junction on the applied magnetic field for a weak proximity effect, $`\gamma _M=0.1`$ , and different exchange field values of the F layer is shown on Fig. 4. As in the case of strong proximity effect, the current greatly depends on the relative orientation of the fields. For parallel geometry, a ground state with a phase shift across the contact about $`\pi /2`$ is found at low field, while with the field increasing the ground state is changed to the $`\pi `$ phase superconductivity. For the antiparallel orientation the current behavior is essentially nonmonotonic and the increasing of the field turns the current through a few extrema until the junction crossovers into the $`\pi `$-phase state with opposite current direction. Shown on Fig. 5. are representative dc current dependences versus the applied magnetic field for the junction with $`\gamma _M=0.1`$ , $`H_e=1.5\pi T_C`$ , and different interface transparency. As is expected, the critical current amplitude decreases with increasing the interface resistance.

As one can see on Figs. 2-5, for antiparallel geometry the $`\pi `$ phase superconductivity takes place when the external field value is about the internal one, $`HH_e`$ . This result, far from being surprising, is related to the quasiparticle spin splitting in the effective magnetic fields in bilayers. Namely, the S/F-I-N/S structure in antiparallel magnetic field $`H=H_e`$ has the same spin splitting as S/N-I-F/S one in zero magnetic field. In other words, for $`H>H_e`$ we should have the critical current versus field dependence as those we have for $`H>0`$ . That is what we see on Figs. 2-5.

### C Asymmetric bilayer interface parameters

To illustrate the main features of the Josephson effect for the S/F-I-N/S tunnel structure caused by an applied magnetic field, we took, for simplicity, the same interface parameters for both SF and SN interfaces. As the nonsuperconducting metals, - ferromagnetic and nonmagnetic ones, - should be different, this approximation may be unrealizable experimentally. A further step consists in extending the theory to a real system, where the boundary conditions for S/F and S/N bilayers can be different.

Generalization of the results to SF and SN interfaces with different transmission probabilities for quasiparticles is obvious. We present the results for two main cases: (a) the structure with a strong proximity effect for the SN boundary and a weak proximity effect for the SF boundary, and vice versa: (b) a strong proximity effect for the SF boundary and a weak proximity effect for the SN boundary. Let us first discuss the analytical results for vanishing interface resistance $`\gamma _B=0`$ . We obtain

$$j_C(H,H_e)\text{ }\mathrm{~}\text{ }(T/T_C)B_{MN}(T)\underset{\omega }{}\frac{\mathrm{\Delta }_S}{\mathrm{\Omega }(\omega ^2+H^2)}\{1\gamma _{MF}\beta H(H_e+H)\mathrm{\Omega }^2\}$$
(7)
if $`\gamma _{MN}>>1,`$ and $`\gamma _{MF}<<1`$ ; here $`B_{MN}(T)=B(T)\pi T_C/\gamma _{MN}`$ . As before, in the geometry with the parallel fields, the increase of external magnetic field inverses the supercurrent direction of the S/F-I-N/S tunnel structure at large enough applied field, $`\gamma _{MF}H(H_e+H)>>\omega ^2\mathrm{~}(\pi T_C)^2`$ . In the geometry with the opposite field orientation, simple substitution into Eq. (7) gives that $`j_C(H_e,H_e)=j_C(0,H_e)`$ , while at $`H=H_e/2`$ we again obtain that $`j_C(H_e/2,H_e)>j_C(0,H_e),`$ $`j_C(H_e,H_e)`$ . At low values of the field, $`|H|<H_e`$ , the nonmonotonic current versus field dependence can be observed; for large enough values of the external field, $`|H|>H_e`$ , the critical current inverses its sign and the crossover of the S/F-I-N/S junction from the 0-phase state to the $`\pi `$-phase state takes place.

Qualitatively similar results are obtained for a junction with $`\gamma _{MN}<<1`$ , $`\gamma _{MF}>>1`$ . If $`\gamma _B=0`$ our analytical calculations give

$$j_C(H,H_e)\text{ }\mathrm{~}\text{ }(T/T_C)B_{MF}(T)\underset{\omega }{}\frac{\mathrm{\Delta }_S}{\mathrm{\Omega }[\omega ^2+(H+H_e)^2]}\{1\gamma _{MN}\beta H(H_e+H)\mathrm{\Omega }^2\}$$
(8)
where $`B_{MF}(T)=B(T)\pi T_C/\gamma _{MF}`$ .

More general cases, $`\gamma _B0`$ , are illustrated by numerical calculations, shown on Figs. 6 and 7. In Fig.6 the dependence of the Josephson current amplitude of the junction on the applied magnetic field is given for a strong proximity effect for the SN boundary, $`\gamma _{MN}>>1`$ , a weak proximity effect for the SF boundary, $`\gamma _{MF}<<1`$, and different interface transparency $`\gamma _B=0÷2`$ . (Here and for Fig. 7 the numerical calculations have been made for simplicity for SF and SN interfaces with the same resistivity, i. e., $`\gamma _{BF}=\gamma _{BN}=\gamma _B`$ ). Opposite situation is illustrated by the results on Fig. 7. As is seen, the main features of the current-field behavior are robust enough and are held for the S/F-I-N/S junction with asymmetric bilayer parameters, too. However, one can also see, that the current amplitude depends on the conditions on the SF and SN interfaces.

In an experiment, the exchange and applied fields may be noncollinear. When these is an angle $`\theta `$ between the exchange and applied fields directions, the expression for the critical current can be written in the form

$`j_C(\theta )=j_C^P\mathrm{cos}^2(\theta /2)+j_C^A\mathrm{sin}^2(\theta /2)`$where $`j_C^{P(A)}`$ is the current for parallel (antiparallel) fields alignment. The dependence of the critical current on the angle between the fields is illustrated by Fig. 8 for S/F-I-N/S structure with $`\gamma _M=0.1`$, $`\gamma _B=5`$ for the SF and SN boundaries and$`H_e=\pi T_C`$.

## V Conclusion

The question of superconductivity in the presence of spin splitting is currently heavily investigated, both theoretically and experimentally. However, in spite of a larger number of theoretical predictions, most of them have not been detected experimentally until now. There are some set-up problems, and one of them is the requirement to fix one layer magnetization, while preserving the possibility to reorient the magnetization of another layer by a weak magnetic field. In this report we have analyzed a superconducting tunnel structure where there is spin splitting for two reasons: spin splitting caused by an external magnetic field and by an internal exchange field. The advantage of such structure is the possibility to get a convenient governing and adjustable parameter - an external magnetic field. We demonstrate that in the contacts under consideration the critical current can reverse its sign (the $`\pi `$-state of the junction) for the parallel field orientation. For the antiparallel field alignment the critical current behavior is much more complicated: for $`|`$$`H|<H_e`$ a monotonic increase of the external field causes an enhancement of the critical current; further field’s increase $`|`$$`H|>H_e`$ , reverses the current direction. We reveal the physics behind such critical current versus external field behavior. Namely, we show that in the case of a strong enough magnetism of the F layer a superconducting phase difference on the junction electrodes is about $`\pi /2`$. If such a contact is in magnetic field, an additional phase shift on the electrodes has been induced, and, depending on the relative orientation, the applied field prompts the junction to the $`\pi `$ \- or to the $`0`$ \- phase state. From the fundamental point of view, these results provide new effects of the superconductivity in the presence of a spin splitting in hybrid superconductor-ferromagnet-normal metal structures and may lead to further understanding of the proximity effects in ferromagnetic materials. The feature important for practical applications is that the superconducting properties of S/F-I-N/S junctions can be simply varied by changing the external magnetic field direction.

Some simplifying assumptions have been made, such as the model admitted of analytical solutions. We consider the dirty limit. However, when the thickness of the F and N layers is much less then its coherence length, $`d_{F(N)}<<\xi _{F(N)}`$ , the results are, in most cases, applicable for clean F and N metals, too. Indeed, if scattering centers are on average separated by a distance smaller than the coherence length, then the dirty limit may be applied. The physical presence of interfaces with atomic roughness ensures that there are such scattering processes for very thin layers.

Let us estimate the orbital effects. The magnetic field H induces screening currents and leads to some suppression of the condensate function. For ferromagnetic and nonmagnetic layers in the dirty limit the depairing rate due to the Meissner current is determined by the energy $`Dp_S^2`$ , where $`p_SHd_{F(N)}/\varphi _O`$ is the condensate momentum (here $`Hd_{F(N)}`$ is the magnitude of the vector potential into the F (N) layer and $`\varphi _O`$ is the flux quantum). So, the orbital effects can be neglected if

$`Dp_S^2\mathrm{~}(D/\xi _{F(N)}^2)(H/H_{C2})^2(d_{F(N)}/\xi _{F(N)})^2\mathrm{~}\mathrm{\Delta }_S(H/H_{C2})^2(d_{F(N)}/\xi _{F(N)})^2<<H,H_e\text{ , }`$where $`H_{C2}`$ ~$`\varphi _O/\xi _{F(N)}^2`$ is the upper critical field of F (N) layer. We see, that the nanoscale nonsuperconducting layers thickness, $`d_{F(N)}<<\xi _{F(N)}`$ , ensures that there is a fairly uniform penetration of an applied field and that orbital depairing effects are minimized. However, for nanoscale hybrid structures a strong electric field arising near metal-metal boundaries may be a source of spin-flip processes . Decreasing the average (ferro)magnetic field experienced by the Cooper pair, the spin-flip scattering is much more destructive than the orbital effect .

It was also assumed that the applied field is concentrated into the ”weak” F-I-N part of the junction. In a real situation, an external magnetic field can partly penetrate into the S layers and cause some suppression of the order parameter $`\mathrm{\Delta }_S`$ . This issue was already addressed a while ago and well understood (see, e.g., de Gennes , Chap.VII). While, in the general case, it is necessary to take into account self-consistently corrections to the pair potential of the S layer, this fact, however, does not qualitatively influence the results of the paper, if the magnetic flux that penetrates into the contact is much less that a flux quantum. The condition is fulfilled when the transverse dimensions of the junction are much less then the Josephson penetration depth.

Let us briefly comment on an experimental set-up, as well. Several recent works (see, e. g., \[5-7,25\]) have found that dilute ferromagnetic $`Cu_xNi_{1x}`$ alloys (with x $`<`$0.6) are most adopted as materials for the F layers, since their weak ferromagnetism is less devastating to superconductivity. Magnetic hysteresis measurements show for $`Cu_xNi_{1x}`$ films with x $``$ 0.4$`÷`$0.5 the coercive field of ~100 Oe (see, e.g.,  and references therein ). One can also pin the magnetization of the $`CuNi`$ layer by adjacent $`FeMn`$ layer via exchange bias, and in such a way enlarge the interval of an applied magnetic field, where the magnetization remains fixed. On the other hand, various niobium oxides, one of them NbO, possess the properties of a normal nonmagnetic metal. Such an N layer is usually formed during the fabrication process (see, e. g.,  and references therein). So, conventional, e. g., magnetron sputtering technique can be used to grow the tunnel superconductor-ferromagnet-insulator-nonmagnetic metal-superconductor structures.

The authors are grateful to M. A. Belogolovskii for reading the paper and useful discussions.

Fig. 1. The S/F-I-N/S system in a parallel magnetic field.

Fig. 2. Critical current of the S/F-I-N/S tunnel junction vs external magnetic field at high interface transparency, $`\gamma _B=0`$ , and strong proximity effect, $`\gamma _M`$ = 10, for different values of the exchange field of the F layer: $`H_e/\pi T_C=`$ 0.5, 0.7, 0.8 and 1.0 (curves 1, 2, 3 and 4, respectively). Here and on Figs. 3-7, solid curves illustrate the case of antiparallel orientation of the internal F layer exchange field and the external magnetic field; dashed curves illustrate the case of parallel orientation of these fields. The numerical results where obtained for $`T=0.1T_C`$ .

Fig. 3. Critical current of the S/F-I-N/S tunnel junction vs external magnetic field at strong proximity effect, $`\gamma _M=10`$, $`H_e=\pi T_C`$ and different values of SF (SN) boundary transparency: $`\gamma _B=`$ 1, 2, 7, and 10 (curves 1, 2, 3 and 4, respectively).

Fig. 4. Critical current vs external magnetic field for the S/F-I-N/S tunnel junction with $`\gamma _B=5`$ and weak proximity effect $`\gamma _M`$ = 0.1 for different values of the exchange field of the F layer: $`H_e/\pi T_C=`$ 1.0, 1.5, and 2.0 (curves 1, 2, and 3, respectively).

Fig. 5. Critical current of an S/F-I-N/S tunnel junction vs exchange energy and weak proximity effect $`\gamma _M=0.1`$, $`H_e=1.5\pi T_C`$ and different values of SF (SN) boundary transparency: $`\gamma _B=`$ 3, 4, and 5 (curves 1, 2, and 3, respectively).

Fig. 6. Critical current for nonsymmetric S/F-I-N/S tunnel structure with a strong proximity effect for the SN boundary, $`\gamma _{MN}=10`$ , and a weak proximity effect for the SF boundary, $`\gamma _{MF}=0.1`$ , $`H_e=\pi T_C`$ and different values of SF (SN) interface transparency: $`\gamma _B=`$ 0.5, 1, and 2 (curves 1, 2 and 3, respectively).

Fig. 7. Critical current of nonsymmetric S/F-I-N/S tunnel structure with a strong proximity effect for the SF boundary, $`\gamma _{MF}=10`$ , and a weak proximity effect for the SN boundary, $`\gamma _{MN}=0.1`$ , $`H_e=\pi T_C`$ and high values of SF (SN) interface transparency: $`\gamma _B=`$ 0.2, 0.5, 0.7 (curves 1, 2 and 3, respectively).

Fig. 8 The dependence of the critical current on the angle $`\theta `$ between the fields for S/F-I-N/S structure with $`\gamma _M=0.1`$, $`\gamma _B=5`$for SF and SN boundaries and$`H_e=\pi T_C`$ ; $`\theta `$ $`=0^o,30^o,60^o,90^o,120^o,150^o,180^o`$ (curves 1 – 7, respectively).


text:
# Virtual Time Horizon Control via Communication Network Design11footnote 1Invited review paper for Volume on Computational Complexity and Statistical Physics, Santa Fe Institute Studies in the Sciences of Complexity series, Oxford University Press, 2003.

## 1 Introduction to the scalability problem of massively parallel discrete-event simulations

The description and understanding of complex systems dynamics is in most cases impossible via analytic methods. The density of problems that are rigorously solvable with analytic tools is vanishingly small in the set of all problems. The only way one can obtain system level understanding of such problems is through direct simulation. The class of complex systems considered here are made of a large number of interacting individual elements with a finite number of attributes, or local state variables, each assuming a countable number (typically finite) of values. The dynamics of the local state variables are discrete events occurring in continuous time. Between two consecutive updates, the local variables stay unchanged. Another important assumption we make is that the interactions in the underlying system to be simulated have finite range. Examples of such systems include: magnetic systems (spin states and spin flip dynamics), surface growth via molecular beam epitaxy (height of the surface at a given point, measured from a growth level, molecular deposition and diffusion dynamics); epidemiology (health of an individual, health state change due to infection, or recovery); financial markets (wealth state, buy/sell dynamics), wireless communications, or queueing systems (number of jobs, job arrival dynamics).

Often, as the case studied here, the dynamics of such systems is inherently stochastic and asynchronous. The simulation of such systems is rather non-trivial and in most cases the complexity of the problem requires simulations on distributed architectures, defining the field of Parallel Discrete-Event Simulations (PDES) . Conceptually, the computational task is divided among $`N`$ processing elements (PE-s), where each processor evolves the dynamics of the allocated piece. Due to the interactions among the individual elements of the simulated system (spins, atoms, packets, calls, etc.) the PE-s must coordinate with a subset of other PE-s during the simulation. For example, the state of a spin can only be updated if the state of the neighbors is known. However, some neighbors might belong to the computational domain of another PE, thus message passing will be required in order to preserve causality. In the PDES schemes we analyze, update attempts are self-initiated  and are independent of the configuration of the underlying system . Although these properties simplify the analysis of the corresponding PDES schemes, they can be highly efficient  and are readily applicable to a large number problems in science and engineering. Further, the performance and the scalability of these PDES schemes become independent of the specific underlying system i.e., we learn the generic behavior of these complex computational schemes.

The update dynamics, together with the information sharing among PE-s, make the parallel discrete-event simulation process a complex dynamical system itself. In fact, it perfectly fits the type of complex systems we are considering here: the individual elements are the PE-s, and their states (local simulated time) evolve according to update events which are dependent on the states of the neighboring PE-s.

With the number and size of parallel computers on the rise, the problem of designing efficient parallel algorithms, or update schemes becomes increasingly important. In passing, we can mention a few examples of large parallel computers: the 9472-node ASCII Red at Sandia, the 12288-node QCDSP Teraflop Machine at Brookhaven, the 8192-node IBM ASCII White with 12.3 Teraflops. IBM is currently building the Blue Gene/L with 200 Teraflops, with over $`10^6`$ nodes. As a matter of fact the largest supercomputer ever built is by Nature itself: the brain, which does an immense parallel computing task to sustain the individual. In particular the human brain has $`10^{11}`$ PE-s (neurons) each with an average of $`10^4`$ synaptic connections, creating a bundle on the order of $`10^{15}`$ “wires” jammed into a volume of approximately $`1400\text{cm}^3`$.

The design of efficient parallel update schemes is a rather challenging problem, due to the fact that the dynamics of the simulation scheme itself is a complex system, whose properties are hard to deduce using classical methods of algorithm analysis. Here we present a less conventional approach to the analysis of efficiency and scalability for the class of massively parallel conservative PDE-s schemes, by exactly mapping the parallel computational process itself onto a non-equilibrium surface growth model . This allows us to translate the questions about efficiency and scalability into questions formulated in terms of certain topological properties of this non-equilibrium surface. Then, using methods from statistical mechanics, developed some time ago to study the dynamics of such surfaces (in a completely different context), we solve the scalability problem of the computational PDES-s scheme . Similar connections between computational schemes and complex systems behavior have recently been made  for rollback-based PDES-s algorithms  and self-organized criticality .

Since one is interested in the dynamics of the underlying complex system, the PDES scheme must simulate the ‘physical time’ variable of the complex system. When the simulations are done on a single processor machine, a single (global) time stream is sufficient to “label” or time-stamp the updates of the local configurations, regardless whether the dynamics of the underlying system is synchronous or asynchronous. When simulating asynchronous dynamics on distributed architectures, however, each PE generates its own physical, or virtual time, which is the physical time variable of the particular computational domain handled by that PE. Due to the varying complexity of the computation at different PE-s, at a given wall-clock instant the simulated virtual times of the PE-s can differ, a phenomenon called “time horizon roughening”. We denote the simulated, or virtual time at PE $`i`$ measured at wall-clock time $`t`$, by $`\tau _i(t)`$. For non-interacting subsystems the wall-clock time $`t`$ is directly proportional to the (discrete) number of parallel steps simultaneously performed on each PE, also called the number of Monte-Carlo steps (MCS) in dynamic Monte Carlo simulations. Without altering the meaning, $`t`$ from now on will be used to denote the number of discrete steps performed in the parallel simulation. The set of virtual times $`\{\tau _i(t)\}_{i=1}^N`$ forms the virtual time horizon of the PDES-s scheme after $`t`$ parallel updates.

In conservative PDES-s schemes , a PE will only perform its next update if it can obtain the correct information to evolve the local configuration (local state) of the underlying physical system it simulates, without violating causality. Otherwise, it idles. Specifically, when the underlying system has nearest-neighbor interactions, each PE must check with its “neighboring” PEs (mimicking the interaction topology of the underlying system) to see if those are progressed at least up to the point in virtual time where the PE itself did . Based on the fundamental notion of discrete-event systems that the state of a local state variable remains unchanged between two successive update attempts, the above rule guarantees the causality of the simulated dynamics . A simple example illustrating this is shown in Figure 1.

One can consider, for example, a magnetic system as the underlying physical system, where the spins are arranged in the sites of a one-dimensional lattice, and that a single spin is handled by a single PE (for more realistic and efficient implementations see ). In Fig. 1, which shows the distribution of the virtual simulated times at a given wall-clock instant $`t`$, the only PE that can update from the set $`\{i1,i,i+1\}`$ is in site $`i`$, since the state of the neighboring spins at sites $`i\pm 1`$ are already known. However, PE-s $`i\pm 1`$ cannot update their spin states at wall clock instant $`t`$, because, the state of the neighboring spin $`i`$ at their simulated times (at $`\tau _{i1}`$ and $`\tau _{i+1}`$) is not known yet. In other words PE $`i`$ can only update at wall-clock instant $`t`$ if $`\tau _i(t)\mathrm{min}\{\tau _{i1}(t),\tau _{i1}(t)\}`$, i.e., its virtual time is a local minimum among the virtual times of its neighboring PE-s. It is easy to see that the same conclusion holds for arbitrary PE topologies. Let the topology for the communication among the processing elements be symbolized by a graph $`G(V,E)`$, where $`V`$ is the vertex set of $`N`$ nodes and $`E`$ is the edge set of $`G`$. Given a node $`iV(G)`$, we denote by $`S_i^{(1)}`$ the set of first order neighbors on $`G`$ of $`i`$. Then, node (PE) $`i`$ can update its state, in the conservative PDES-s scheme iff:

$$\tau _i^{(G)}(t)\underset{jS_i^{(1)}}{\mathrm{min}}\{\tau _j^{(G)}(t)\},i=1,..,N.$$
(1)
In the following, the set of (active) nodes which obey condition (1) at time $`t`$, will be denoted by $`𝒜(t)`$. Now we are in the position to formulate the scalability problem of PDES-schemes for systems with asynchronous dynamics. For the PDES-s scheme to be fully scalable, the following two criteria must be met: (i) the virtual time horizon must progress on average at a nonzero rate, and (ii) the typical spread of the time horizon should be finite, as the number of PE-s $`N`$ goes to infinity. When the first criterion is ensured for large enough times $`t`$, the simulation is said to be computationally scalable, and it just means that when increasing the size of the computation to infinity, while keeping the average computational domain/load on a single PE the same, the simulation will progress at a nonzero rate. However, as we will show below, increasing the system size, the spread in the time horizon can diverge, severely hindering frequent data collection about the state of the simulated system. Specifically, when one requires to take a measurement of some physical property of the simulated system at (virtual, or simulated) physical time $`\tau `$, we have to wait (in wall-clock time) until all the virtual simulated times at all the PE-s pass through the value of $`\tau `$. Thus, in order to collect system-wide measurements from the simulation, we incur a waiting time proportional to the spread, or width of the fluctuating time horizon. For PDES-s schemes for which the spread diverges with system size, the waiting time for the measurements will also diverge, and the scheme is not measurement scalable. When condition (ii) is fulfilled for large enough times $`t`$, we say that the PDES-s scheme is measurement scalable.

The scalability criteria above can simply be formalized in terms of the properties of the virtual time horizon, $`\{\tau _i^{(G)}(t)\}_{i=1}^N`$. The average of the time horizon after $`t`$ parallel steps is:

$$\overline{\tau }^{(G)}(t)=\frac{1}{N}\underset{j=1}{\overset{N}{}}\tau _j^{(G)}(t).$$
(2)
At a given (wall-clock) time $`t`$ the only PE-s that can make progress, i.e., are not in idle, are those which have simulated virtual times obeying condition (1). Thus, the rate of progress of the time horizon average becomes:

$$\overline{\tau }^{(G)}(t+1)\overline{\tau }^{(G)}(t)=\frac{1}{N}\underset{l𝒜(t)}{}\left[\tau _l^{(G)}(t+1)\tau _l^{(G)}(t)\right].$$
(3)
The difference in the square brackets on the r.h.s of (3) is the physical time elapsed between two consecutive events in the physical domain simulated by the $`l`$th PE, and it is determined by the physical process responsible for the stochastic dynamics of the simulated complex system. If we replace the time intervals in square brackets in (3) with their (obviously finite) average value $`\mathrm{\Delta }`$, we obtain that the average progress rate of the time horizon, or average utilization $`u^{(G)}(t)=\overline{\tau }^{(G)}(t+1)\overline{\tau }^{(G)}(t)`$ is proportional to the number of non-idling, or active PE-s. The average $``$ is taken over the stochastic event dynamics, which is assumed to be the same at all sites. For many cases, the $`\mathrm{\Delta }`$ factor is independent on $`N`$ (due to the finite range of the interaction in the complex system), so the computational efficiency, or average utilization of the simulation can simply be identified with the average density of the active PE-s:

$$u^{(G)}(t)=\frac{|𝒜^{(G)}(t)|}{N}$$
(4)
where $`|𝒜^{(G)}(t)|`$ denotes the number of elements of the set $`𝒜^{(G)}(t)`$. Thus, the PDES-s scheme is computationally scalable, if there exists a constant $`c>0`$, such that:

$$u^{(G)}(\mathrm{})=\underset{\genfrac{}{}{0pt}{}{t\mathrm{}}{N\mathrm{}}}{lim}\frac{|𝒜^{(G)}(t)|}{N}>c.$$
(5)
The measurement scalability of the PDES-s scheme, is characterized by the spread of the virtual time horizon. Instead of dealing with the actual spread (difference between the maximum and minimum values) we shall consider the average “width” of the time horizon defined as:

$$[w^{(G)}]^2(t)=\frac{1}{N}\underset{j=1}{\overset{N}{}}\left[\tau _j^{(G)}(t)\overline{\tau }^{(G)}(t)\right]^2.$$
(6)
A PDES-s scheme is measurement scalable, if there exists a constant $`M>0`$, such:

$$[w^{(G)}]^2(\mathrm{})=\underset{\genfrac{}{}{0pt}{}{t\mathrm{}}{N\mathrm{}}}{lim}\frac{1}{N}\underset{j=1}{\overset{N}{}}\left[\tau _j^{(G)}(t)\overline{\tau }^{(G)}(t)\right]^2<M.$$
(7)
In reality, the number of PEs $`N`$ or the simulation time $`t`$ can never be taken to infinity, so for practical purposes, the scalability is deduced from the scaling behavior of the quantities for long times and for large number of PEs. The setup presented above is perfectly suitable to establish a mapping between non-equilibrium surface growth models  and conservative PDES-s schemes. We discuss extensively this mapping in the next section.

The paper is organized as follows: in Section 2 we discuss the scalability of the computational phase and the failure of scalability of the measurement phase of the basic conservative scheme on regular topologies. Then we show how a simple modification of the communication topology (from a regular lattice to a small-world structure) leads to a fully scalable PDES-s scheme. In section 3 we study the scalability problem on scale-free network topologies. Section 4 is devoted to conclusions.

## 2 Scalability of conservative PDES-s schemes on regular and small-world topologies

In many large complex systems the stochastic event dynamics can be characterized by a Poisson distributed stream. For example, in an Ising magnet with single spin-flip Glauber dynamics  the spin-flip attempts are Poisson distributed events, or in wireless cellular communications the call arrivals also obey Poisson statistics , etc. In the following we restrict ourselves to such Poisson distributed stochastic processes for event dynamics, however numerical simulations show, that our conclusions for scalability hold for a large class of other stochastic distributions, as well. The evolution of the virtual time horizon incorporating condition (1) for Poisson asynchrony is given by the equation:

$$\tau _i^{(G)}(t+1)=\tau _i^{(G)}(t)+\eta _i(t)\underset{jS_i^{(1)}}{}\theta (\tau _j^{(G)}(t)\tau _i^{(G)}(t))$$
(8)
Here $`\theta (x)`$ is the Heaviside step function, and $`\eta _i(t)`$ the Poisson distributed virtual time increment at PE $`i`$, and time $`t`$ is. These increments are drawn at random, independently of $`i`$ and $`t`$, and of the existing time horizon.

### 2.1 The basic conservative scheme on regular topologies

Next, we consider the basic conservative scheme, which is defined on regular, square lattice communication topologies, in $`d`$ dimensions, so that $`N=L^d`$. For brevity, we drop from the superscript $`(G)`$ in the notation for $`\tau _i(t)`$. In particular, we first illustrate our analysis on the simplest regular topology, that of a regular one-dimensional lattice, with periodic boundary conditions ($`G`$ is a ring). Later we discuss the general, $`d`$-dimensional case. The evolution equation on the ring becomes simply:

$$\tau _i(t+1)=\tau _i(t)+\eta _i(t)\theta (\tau _{i1}(t)\tau _i(t))\theta (\tau _{i+1}(t)\tau _i(t)).$$
(9)
with the boundary conditions $`\tau _{N+1}=\tau _N=\tau _0`$. The total number of active sites/PE-s is thus given by $`|𝒜(t)|=_{i=1}^N\theta (\tau _{i1}(t)\tau _i(t))\theta (\tau _{i+1}(t)\tau _i(t))`$ so the average utilization (4) becomes:

$$u(L,t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\theta (\tau _{i1}(t)\tau _i(t))\theta (\tau _{i+1}(t)\tau _i(t))$$
(10)
The average $``$ is performed over the random variables $`\{\eta _i(t^{})\}_{\genfrac{}{}{0pt}{}{i=1,..,L}{t^{}=1,..,t}}`$ which have an exponential distribution, $`Prob\{x<\eta x+\delta x\}=_x^{x+\delta x}𝑑ye^y`$. In spite of the simple appearance of the dynamics (9), and the exponential (or Poisson) stochastic dynamics at nodes, calculating the average utilization (10) is very difficult. Even a rigorous proof for the existence of the lower bound (5) using direct methods is still an open problem.

Here we present a different approach, by mapping first the problem to a non-equilibrium surface grown via molecular beam epitaxy (where atoms or molecules are deposited from vapors, or beams onto the surface) model. The quantities brought in analogy are: the $`i`$-th PE is the site $`i`$ in the substrate; the number of parallel updates $`t`$ is the number of deposited monolayers; $`\tau _i(t)`$ is the height $`h_i(t)`$ at site $`i`$ and time $`t`$; a virtual time increment of $`\eta _i(t)`$ at PE $`i`$ in the $`t`$-th step, corresponds to a material “rod” of length $`\eta _i(t)`$ deposited onto the surface, see Figure (2). The length of the rod is a Poisson distributed random variable. During the $`t`$-th update, the rods are deposited only in local minima of the surface. The utilization of the PDES-s scheme corresponds to the the density of local minima of the growing surface. Even though the length of the rods are random independent variables, the fact that they can only be deposited in local minima will generate lateral correlations into the surface fluctuations, and makes the problem hard to compute exactly.

The rods are deposited onto the surface in a parallel update scheme: after all local minima are updated (deposited onto) the the time $`t`$ is incremented by unity. We will coin the surface growth analog of our basic conservative PDES-s scheme , the ‘MPEU’ model (Massively Parallel Exponential Update Model).

Both the utilization (density of minima) and the width of the time horizon are quantities characterizing the fluctuations of the growing surface. The type of fluctuations can be classified into universality classes, each class having distinct statistical properties. Studying the PDES-s scheme as a surface growth model, we can describe its fluctuations and identify the surface growth universality class it belongs to. In order to do that we first introduce the slope variables, $`\varphi _i=\tau _i\tau _{i1}`$. Provided $`\tau _i(t)`$ is a local minimum, a rod deposited of length $`\eta _i`$ corresponds to taking an amount of $`\eta _i`$ from $`\varphi _{i+1}`$ and adding it to $`\varphi _i`$, since $`\varphi _i(t+1)=\tau _i(t)\tau _{i1}(t)+\eta _i(t)`$ and $`\varphi _{i+1}(t+1)=\tau _{i+1}(t)\tau _i(t)\eta _i(t)`$, i.e., in the surface of slopes, $`\{\varphi _i\}_{i=1}^L`$ the dynamics is biased surface diffusion, given by the equation:

$$\varphi _i(t+1)\varphi _i(t)=\eta _i(t)\theta (\varphi _i(t))\theta (\varphi _{i+1}(t))\eta _{i1}(t)\theta (\varphi _{i1}(t))\theta (\varphi _i(t))$$
(11)
with the constraint $`_{i=1}^L\varphi _i=0`$ generated by the periodic boundary conditions in the $`\tau `$ variables. In terms of the local slope variables, the expression for the average density of minima, or average utilization becomes: $`u(L,t)=\frac{1}{L}_{i=1}^L\theta (\varphi _i(t))\theta (\varphi _{i+1}(t))`$. Translational invariance implies (no node is statistically special) $`u(L,t)=\theta (\varphi _i(t))\theta (\varphi _{i+1}(t))`$ for any $`i=1,2,..,L`$. From (9) it follows that $`\tau _i(t+1)\tau _i(t)=u(L,t)`$, thus, the average rate of propagation of the MPEU surface is identical to the average utilization of the PDES-s scheme. It is also easy to see that in the slope language it is identical to the average current in the ring. Next we perform a naive coarse graining by using the representation $`\theta (\varphi )=lim_{\kappa 0}\frac{1}{2}\left[1+\mathrm{tanh}(\varphi /\kappa )\right]`$, and keeping only the terms up to the order $`(\varphi /\kappa )`$. Performing this coarse graining, one obtains:

$$\varphi _i(t+1)\varphi _i(t)=\frac{1}{4\kappa }\varphi _{i+1}(t)2\varphi _i(t)+\varphi _{i1}(t)\frac{1}{4\kappa ^2}\varphi _i(\varphi _{i+1}\varphi _{i1}).$$
(12)
Strictly speaking all of the $`(\varphi /\kappa )^n`$, $`n=1,2,\mathrm{}`$ terms are divergent, but taking the proper continuum limit and introducing an appropriately scaled bias, the only relevant terms can be shown to be those appearing in Eq. (12) . In the continuum limit, one thus, obtains for the coarse-grained field:

$$\frac{}{t}\widehat{\varphi }=\frac{^2}{x^2}\widehat{\varphi }\lambda \frac{}{x}\widehat{\varphi }^2$$
(13)
where $`\lambda `$ is a parameter related to the coarse-graining procedure. The above nonlinear partial differential equation (13) is known as the nonlinear biased diffusion, or the Burgers equation . Returning to the coarse-grained equivalent of the height, or virtual times, $`\widehat{\tau }`$, via $`\widehat{\varphi }=\widehat{\tau }/x`$, we obtain the Kardar-Parisi-Zhang (KPZ) equation :

$$\frac{\widehat{\tau }}{t}=\frac{^2\widehat{\tau }}{x^2}\lambda \left(\frac{\widehat{\tau }}{x}\right)^2.$$
(14)
To capture the fluctuations, one typically adds a delta-correlated noise term $`\xi (x,t)`$, to the right hand side, conserved for Eq. (13) (i.e., $`\xi 𝑑x=0`$), and non-conserved for Eq. (14). It is important to note that we obtained the KPZ equation as a result of a coarse-graining procedure. While this results in the “loss” of some of the microscopic details for the original growth model on the lattice, Eq. (14) with noise added captures the long-wavelength behavior of the MPEU model. Thus, we claim that virtual time horizon for the basic conservative PDES scheme exhibit kinetic roughening and it belongs to the KPZ universality class. Identifying the universality class of a model is one of the main objectives and used extensively in surface science to classify fluctuation statistics. Our procedure above indicates that the long-wavelength statistics of the fluctuations of the time horizon for the basic conservative PDES-s scheme is in fact captured by the nonlinear KPZ equation. In one dimension a steady-state for the surface fluctuations is reached (in the long time limit) for any finite system-size and it is governed by the Edwards-Wilkinson (EW) Hamiltonian $`_{EW}𝑑x\left(\frac{\widehat{\tau }}{x}\right)^2`$ (see, e.g., ). The corresponding surface is a simple random-walk surface, where the slopes are independent random variables in the steady state. This means that out of the four local configurations of slopes around a point (down-up, down-down, up-up, up-down) only one contributes in average to a minimum (down-up), and since they are all equally likely, we conclude that $`u_{EW}(L\mathrm{},t\mathrm{})=1/4=0.25`$. Our numerical simulations for the MPEU model, see Fig.3a) indicate a value of $`u(L\mathrm{},t\mathrm{})=0.24641\pm (7\times 10^6)`$ a value close, but definitively not the same as for the simple random walk surface. The reason for the obvious difference is that the coarse-grained version and the original microscopic model are not identical over the whole spectrum of wavelengths of the fluctuations.

The coarse-graining procedure preserves the statistics of the long-wavelength modes, but it looses some information on the short-wavelength ones. In particular the density of minima is heavily influenced by the short wavelengths (by how “fuzzy” the interface is). However, the density of minima cannot vanish in the thermodynamic limit: a zero density of local minima would imply that it is zero on all length-scales which would contradict the fact that it belongs to the EW universality class. The fact that the steady-state of the MPEU model belongs to the EW universality class guarantees that the local slopes are short-range correlated \[Fig.3c)\], and that the finite-size corrections for the density of local minima (average propagation rate of the surface) follows a universal scaling form :

$$u(L,\mathrm{})u(\mathrm{},\mathrm{})+\frac{\text{const.}}{L^{2(1\alpha )}}.$$
(15)
Here $`\alpha `$ is the roughness exponent (equals to $`1/2`$ for the EW universality class), characterizing the macroscopic surface-height fluctuations, as described in detail in the next paragraph. Figure 3 confirms this scaling behavior. Further, calculating the spread in the average utilization the steady state as function of system size, $`\sigma _L^2=u^2(L,\mathrm{})u(L,\mathrm{})^2`$, we obtain $`\sigma _LL^{1/2}`$. These findings suggest that the utilization is a self-averaging macroscopic quantity: its full distribution $`P_L(u)`$ for large $`L`$ is a Gaussian \[Fig. 3b)\].

In the following we show numerical results supporting our claim that the MPEU model belongs to the KPZ universality class. One of the fundamental characteristic quantities strongly influenced by the long-wavelength modes is the average width of the height fluctuations, as given in Eq. (6). As the surface grows due to deposition, after an initial transient the growth of the width will grow as a power law $`w^2(L,t)t^{2\beta }`$ along with the lateral surface correlations $`\xi _{||}(L,t)t^{1/z}`$, until the correlations reach the system size ($`\xi _{||}=L`$) at a crossover time $`t_\times `$ (see, e.g., Ref. ). After the crossover time $`t_\times `$ (for any finite system $`L`$) the surface fluctuations are governed by a steady-state distribution and the width scales as

$$w^2(L,\mathrm{})L^{2\alpha }.$$
(16)
The exponent $`\beta `$ is called the growth exponent, $`\alpha `$ is the called roughness exponent and $`z`$ is called the dynamic exponent in the surface growth literature . It is easy to show that the three exponents are not all independent, and $`\alpha =z\beta `$ holds . Also, these scaling forms allow one to collapse all the different curves for the width onto a single function in the scaling regime, expressing the dynamic scaling property of the width: $`w^2(L,t)=L^{2\alpha }f(t/L^z)`$ ($`f`$ is easy to read of after comparing it to the scaling behavior). For the KPZ interface, the analytically obtained exact values for the exponents are: $`\beta =1/3`$, $`\alpha =1/2`$ and $`z=3/2`$. Fig. 4 shows the numerically measured scaling properties for the width of the MPEU model. the numerically obtained values for the exponents $`\beta =0.326\pm 0.005`$, $`\alpha =0.49\pm 0.01`$ for large system sizes ($`L=10^5`$) confirm the KPZ behavior, including the dynamical scaling property (inset). Another confirmation for the EW universality class in the steady state, comes from measuring the full width distribution $`P(w^2)`$. For systems belonging to the EW universality class and having the same type of boundary conditions imposed, the width distribution has a universal scaling form $`P(w^2)=\frac{1}{w^2}\mathrm{\Phi }\left(\frac{w^2}{w^2}\right)`$ with

$$\mathrm{\Phi }(x)=\frac{\pi ^2}{3}\underset{n=1}{\overset{\mathrm{}}{}}(1)^{n1}n^2e^{\frac{\pi ^2}{6}n^2x},$$
(17)
where the above scaled distribution is for the case of periodic boundary conditions. Figure 4b) is a confirmation that MPEU indeed belongs to the steady state of the EW class, implying that the average utilization (density of local minima) approaches a non-zero, finite value in the thermodynamic limit (5) as reflected by Eq. (15). Therefore, the basic conservative scheme is computationally scalable. \[For an in-depth and systematic analytical calculation of the density of minima (utilization) for a number of surface growth models (including the EW class) see Ref. .\] The measurement phase of the basic conservative scheme, however, is not scalable, as indicated by the power-law divergence of the width in the long-time large $`L`$ limit, Eq. (16). For higher-dimensional topologies, using universality arguments, the conclusion remains the same: the basic conservative PDES-s is computationally scalable, but the measurement phase may not be, depending on the upper critical dimension  of the surface (see Ref. ).

### 2.2 The conservative scheme on small world networks

From the previous section it follows that the average width of the fluctuations scales in the steady state as $`w^2(L,t`$$`=`$$`\mathrm{})L^{2\alpha }=L`$ , i.e., it grows linearly with the system size. This means that the basic conservative PDES-s scheme is not measurement scalable. Standard methods to control the width of the virtual time horizon in a PDES scheme utilize some kind of a windowing technique . That is, the height of the local simulated time at any PE cannot progress beyond an appropriately chosen and regularly updated “cap”, measured from the global minimum of the time horizon . Thus, a PDES scheme with a moving window relies on frequent global synchronizations or communications, which, depending on the architecture, can get costly for large number of PEs. Here we show how to modify the original conservative scheme such that the scheme is also measurement scalable without global “intervention” .

The divergence of the width of the surface fluctuations is closely related to the fact that the lateral surface correlations also grow with the system size. In particular, for the one-dimensional EW surface in the steady state, for large $`L`$ (and fixed $`l`$)

$$\widehat{\tau }_i\widehat{\tau }_{i+l}\xi _{||}(L,\mathrm{})|l|,$$
(18)
where $`\widehat{\tau }_i`$ are the coarse-grained height fluctuations measured from the mean and $`\xi _{||}(L,\mathrm{})L`$. Thus, $`w^2(L,\mathrm{})=\widehat{\tau }_i^2\xi _{||}(L,\mathrm{})L`$. The “height-height” correlations can be characterized by introducing the structure factor for the heights:

$$S^{(\tau )}(k)=\frac{1}{L}\stackrel{~}{\tau }_k\stackrel{~}{\tau }_k$$
(19)
where $`k=2\pi n/L`$, $`n=0,1,2,..,L1`$ is the wave-vector, and $`\stackrel{~}{\tau }_k=_{j=0}^{L1}e^{ikj}(\tau _j\overline{\tau })`$ is the discrete spatial Fourier transform of the fluctuations of virtual time horizon. Then

$$\widehat{\tau }_i\widehat{\tau }_{i+l}=\frac{1}{L}\underset{k}{}e^{ikl}S^{(\tau )}(k)$$
(20)
and

$$w^2(L,\mathrm{})=\frac{1}{L}\underset{k}{}S^{(\tau )}(k).$$
(21)
Since the universality class for the time horizon evolution is EW, it follows that the expected behavior for the steady-state structure factor for small wave-numbers is

$$S^{(\tau )}(k)\frac{1}{k^2}$$
(22)
(see, e.g., Eq. (11) in Ref. ).

Indeed, this is also confirmed by our direct simulation results, shown Fig. 5a). This form of the structure factor implies that there are no lengthscales other than the lattice constant and the systems size, and thus, the correlation length and the width diverge in the thermodynamic limit, as also can be seen by directly evaluating Eq. (21).

To de-correlate the surface fluctuations, we modify the communication topology in the following way : for every node $`i`$, at the onset of the simulation, we introduce one extra quenched random communication link $`r(i)`$. Together with the existing regular topology, these extra communication links will form a small-world graph . Note that in our specific construction of the small-world network, each node has exactly one random connection and $`r(r(i))`$$`=`$$`i`$, so that there are exactly $`L/2`$ random links distributed. The updating on PE $`i`$ will obey the following probabilistically chosen condition:

$`\tau _i\{\begin{array}{c}\mathrm{min}\{\tau _{i1},\tau _{i+1},\tau _{r(i)}\}\text{with probability}p\hfill \\ \mathrm{min}\{\tau _{i1},\tau _{i+1}\}\text{with probability}\mathrm{\hspace{0.33em}\hspace{0.33em}1}p\hfill \end{array}.`$ (25)
The PE actually performs the update (generate the virtual time of the next update, or deposit the rod at $`i`$ in the MPEU surface) if condition (25) is fulfilled. This means that for sites that would normally be updated within the basic conservative scheme, i.e., $`\tau _i\mathrm{min}\{\tau _{i1},\tau _{i+1}\}`$ the PE will make an extra check for the condition $`\tau _i\tau _{r(i)}`$ with probability $`p`$. The parameter $`p`$ allows us to tune the scalability properties of the corresponding PDES scheme on the quenched small-world network continuously from the pure basic conservative scheme ($`p=0`$) to the “fully” small-world conservative scheme ($`p=1`$).

These occasional extra checkings through the quenched random links are not necessary for the faithfulness of the simulation. It is merely used to synchronize the PEs in such a way that the fluctuations of the time horizon remain bounded in the infinite system-size limit. Most importantly, as the width is reduced from “infinity” (or some large number proportional to $`L`$ for finite number of PEs) to a finite controlled value, the utilization still remains bounded away from zero.

To support this statement, we first use the same coarse-graining procedure used to derive the KPZ equations as the continuum counterpart of the MPEU model. For the small-world topology we obtain

$$\frac{\widehat{\tau }}{t}=\gamma (p)\widehat{\tau }+\frac{^2\widehat{\tau }}{x^2}\lambda \left(\frac{\widehat{\tau }}{x}\right)^2+\xi (x,t).$$
(26)
with $`\gamma (p)=0`$ for $`p=0`$, and $`\gamma (p)>0`$ for $`0<p1`$. This implies that the extra checking along the random links introduces a strong relaxation (first term on the rhs of (26)) for the long-wavelength modes of the surface fluctuations, resulting in a finite width. A more transparent picture is gained if we look at the steady-state structure factor (19). Restricting our attention to the linear terms in Eq. (26) we obtain

$$S^{(\tau )}(k)\frac{1}{\gamma +k^2}.$$
(27)
In this approximation, the lateral correlation length $`\xi _{||}`$ scales as $`1/\sqrt{\gamma }`$, and remains finite (and system-size independent) in the thermodynamic limit for all $`p>0`$ (i.e., for an arbitrary small probability to utilize the random links). Figure 5b) shows the structure factor for the small-world network with $`p=0.1`$, confirming the prediction of Eq. (27) for small wave numbers. Consequently, the height-height correlations decay exponentially

$$\widehat{\tau }_i\widehat{\tau }_{i+l}\xi _{||}e^{|l|/\xi _{||}},$$
(28)
and the width remains finite, $`w^2(L,\mathrm{})\xi _{||}`$, where $`\xi _{||}`$ is independent of the systems size for all $`p>0`$. Further, for the structure factors of the local slopes (the Fourier transform of the slope-slope correlations) one obtains

$$S^{(\varphi )}(k)=\frac{1}{L}\stackrel{~}{\varphi }_k\stackrel{~}{\varphi }_k=k^2S^{(\tau )}(k)\frac{k^2}{\gamma +k^2}=1\frac{\gamma }{\gamma +k^2}.$$
(29)
Both terms above yield short-range correlations (delta function for the first term and exponential decay for the second one), thus, the slopes remain short-range correlated, resulting in a non-zero density of local minima. Figure 6 shows two snapshots of the virtual time horizons for the basic conservative scheme $`p=0`$, and the small-world scheme with $`p=0.1`$. Figure 7a) shows the scaling of the steady state width with the system size for various $`p`$ values and Fig.7b) shows the scaling of the average, steady state utilization with the system size for the same set of $`p`$ values. Notice that when increasing $`p`$ (from $`p`$$`=`$$`0`$ to $`p`$$`=`$$`0.01`$), the width instantaneously drops from a linear divergence to a saturated value, while at the same time, the utilization hardly changes. In fact, an infinitesimally small $`p`$ will make the width bounded, but at an only infinitesimal expense to the utilization. For example, for a hypothetically infinite system, taking $`p`$$`=`$$`0.01`$, the width is reduced from infinity to about $`40`$, while the utilization from $`0.2464`$ only to about $`0.246`$; for $`p`$$`=`$$`0.1`$, the width is further reduced to about $`5`$, while the utilization only to $`0.242`$. By further increasing $`p`$, the width further reduces, and at $`p`$$`=`$$`1`$ it is about $`1.46`$, whereas the utilization decreases to $`0.141`$, still clearly bounded away from zero in the thermodynamic limit.

## 3 Scalability of the Conservative PDES-s scheme on scale-free network topologies

The Internet is a spontaneously grown collection of connected computers. The number of (only) webservers by February 2003 reached over 35 million . The number of PC-s in use (Internet users) surpassed 660 million in 2002, and it is projected to surpass one billion by 2007 .

The idea for using it as a giant supercomputer is rather natural: many computers are in an idle state, running at best some kind of screen-saver software, and the “waisted” computational time is simply immense. Projects such as SETI@home  or the GRID consortium  are targeting to harness the power lost in screen-savers.

Most of the problems solved currently with distributed computation on the Internet are “embarrassingly parallel”  , i.e., the computed tasks have little or no connection to each other similar to starting the same run with a number of different random seeds, and at the end collecting the data to perform statistical averages. However, before more large-scale, complex problems can be solved in real time on the Internet a number of challenges have to be solved, such as the task allocation problem which is rather complex by itself .

Here we ask the following question: Given that task allocation is resolved and the PE communication topology on the internet is a scale-free network, what are the scalability properties of a PDES-s scheme on such networks? Here we present numerical results, for the PDES-s update scheme, as measured on a model of scale-free networks, namely the Barabási-Albert model . This network is created through the stochastic process of preferential attachment: to the existing network at time $`t`$ of $`N`$ nodes, attaches the $`N+1`$-th node with $`m`$ links (“stubs”) at time $`t+1`$, such that each stub attaches to a node with probability proportional to the existing degree (at $`t`$) of the node. Here we will only present the $`m=1`$ case, when the network is a scale-free tree. Once we reach a given number of nodes in the network, we stop the process and use the random network instance to run the MPEU model on top of it, using the evolution equation (8) for the time horizon. While in case of regular topologies, the degree of a node is constant, e.g. for $`d`$-dimensional “square” lattices, $`P^{(L^d)}(k)=2d\delta _{k,2d}`$ , for the BA network, it is a power law in the asymptotic ($`N\mathrm{}`$) limit : $`P^{BA}(k)2m^2k^3`$.

The condition (1) for a site to be updated, i.e., that its virtual time is a local minimum, is a local property, and thus we expect that the utilization itself be correlated with local structural properties of the graph, such as the degree distribution. To get a more detailed picture, we define two more quantities:
1) connectivity-utilization:

$$u_k(N,t)=\frac{|𝒜_k(t)|}{N}$$
(30)
which is the fraction of active nodes of degree $`k`$, and
2) the relative connectivity-utilization:

$$r_k(N,t)=\frac{|𝒜_k(t)|}{N_k}$$
(31)
which is the fraction of active nodes of degree $`k`$ within the set of all nodes of degree $`k`$. From the abobe definitions we find the following ovious relations: $`_ku_k(N,t)=u(N,t)`$ and $`_kr_k(N,t)N_k/N=_kr_k(N,t)P^{BA}(k)=_ku_k(N,t)=u(N,t)=r_k(N,t)_{network}`$ at all times. Figure 8 shows the steady state ($`t\mathrm{}`$, in the MPEU model on a fixed BA network of $`N`$ nodes) values of the average utilization as function of the network size $`N`$. The inset in Fig. 8 is analogous to Fig. 3(a) which showed the same quantity on a ring. Notice that strictly speaking, the PDES-s scheme is computationally non-scalable. An empirical fit suggest that $`u^{}(N)=u(N,t`$$`=`$$`\mathrm{})\left[\mathrm{ln}\left(aN^b\right)\right]^1`$ with $`a3.322`$ and $`b=0.902`$, i.e., the computation is only logarithmically (or marginally) non-scalable. For a system of $`N`$$`=`$$`10^3`$ nodes we have found a steady state utilization (for the worst case scenario) of $`u^{}(10^3)=0.1328`$ ($`13.3\%`$ efficiency), while for a system of a million nodes, $`N`$$`=`$$`10^6`$, the utilization dropped only to $`u^{}(10^6)=0.073`$ ($`7.3\%`$ efficiency), by less than half of its value! For practical purposes the PDES-s scheme can be considered computationally scalable, and this type of non-scalability we will call logarithmic (or marginal) non-scalability.

Figure 9 shows the scaling of the width of the fluctuations for the time horizon as function of time, and the scaling of its value in the steady-state as function of system size (inset). Notice, that while the steady state width diverges to infinity, it only does so logarithmically, $`w^2(N,t`$$`=`$$`\mathrm{})\left[\mathrm{ln}\left(cN^d\right)\right]`$ with $`c1.25`$ and $`d=0.401`$. Some specific values: $`w^2(10^3,t`$$`=`$$`\mathrm{})3.01`$, $`w^2(10^5,t`$$`=`$$`\mathrm{})4.78`$. This means that the measurement phase of the PDES-s scheme on a scale-free network is non-scalable either, however, it is so only logarithmically, and for practical purposes the scheme can be considered scalable. Overall, the PDES-s update scheme has logarithmic (or marginal) non-scalability on scale-free networks. If one examines the connectivity-utilization and relative connectivity-utilization in the steady state, as shown in Fig. 10, one finds that with good approximation $`u_k^{}(N)k^3`$, and $`r_k^{}(N)=const.u^{}(N)`$ for $`kk_\times `$ and $`r_k^{}(N)k^3`$ for $`k>k_\times `$, with $`k_\times 1/u^{}(N)=\mathrm{ln}\left(aN^b\right)\mathrm{ln}N`$, being a crossover degree.

## 4 Conclusions

We studied the fundamental scalability problem of conservative PDES schemes where events are self-initiated and have identical distribution on each PE. First, we considered the scalability of the basic conservative scheme for systems with short-range interactions on regular lattices. By exploiting a mapping  between the progress of the simulation and kinetic roughening in non-equilibrium surfaces, we found that while the average progress rate of the PEs $`u(\mathrm{},\mathrm{})`$ is a finite non-zero value, the spread of the progress of the PEs about the mean $`w^2(\mathrm{},\mathrm{})`$ diverges. The former property makes the measurement phase of the algorithm non-scalable. In order to make the measurement part of the simulation scalable as well, we introduced  quenched random connections between PEs (exactly one for each) so that the resulting random links on the top of the regular short-range connections formed a small-world-like connection topology. Invoking the same conservative protocol used at an arbitrarily small rate through the random links was sufficient to achieve full scalability: the PEs progress at a non-zero, near uniform rate without requiring global synchronization . The above construction of a fully scalable algorithm for simulating large systems with asynchronous dynamics and short-range interactions is an example for the enormous “computational power and synchronizability”  that can be achieved by small-world couplings. The suppression of critical fluctuations of the virtual time horizon is also closely related to the emergence of mean-field-like phase transitions and phase ordering in non-frustrated interacting systems . In particular, the fluctuations exhibited by the virtual time horizon with small-world synchronization should exhibit very similar characteristics to the fluctuations of the order parameter in the XY-model placed on a small-world network .

Second, we studied the scalability properties for a causally constrained PDES scheme hosted by a network of computers where the network is scale-free following a “preferential attachment” construction . Here the PEs simply have to satisfy the general criterion Eq. (1) in order to advance their local time. Despite some nodes in the network having abnormally large connectivity (as a result of the scale-free nature of the degree distribution), we found that the computational phase of the algorithm is only marginally non-scalable. The utilization exhibited slow logarithmic decay as a function of the number of PEs. At the same time, the width of the time horizon diverged logarithmically slowly, rendering the measurement phase of the simulations marginally non-scalable as well. An intriguing question to pursue is how the logarithmic divergence of the surface fluctuations observed here can be related to the collective behavior (in particular, the finite-size effects of the magnetic susceptibility) of Ising ferromagnets on scale-free networks  with the same degree distribution.

## Acknowledgements

Discussions with P.A. Rikvold, B.D. Lubachevsky, Z. Rácz and G. Istrate are gratefully acknowledged. Z.T. was supported by the DOE under contract W-7405-ENG-36. This research is supported in part by NSF through Grant No. DMR-0113049 and the Research Corporation through Grant No. RI0761.


text:
# On the physical significance of 𝑞-deformation to isovector pairing interactions in nuclei

## Abstract

The quantum deformation concept is applied to a study of isovector pairing correlations in nuclei of the mass $`40A100`$ region. While the non-deformed ($`q1`$) limit of the theory provides a reasonable global estimate for strength parameters of the pairing interaction, the results show that the $`q`$-deformation plays a significant role in understanding higher-order effects in the interaction.

1. Introduction — In addition to purely mathematical examinations of quantum algebraic concepts (see e.g. ), studies focused on quantum (or $`q`$-) deformation  in various fields of physics has been the focus of considerable attention in recent years. Studies of interest include applications in string/brane theory, conformal field theory, statistical mechanics, and metal clusters . A feature of any quantum theory is that in the $`q1`$ limit one recovers the “classical” (non-deformed) results.
The earliest applications of $`q`$-deformation in nuclear physics were related to a description of rotational bands in axially deformed nuclei  using the second order Casimir operator of $`SU_q(2)`$. Although optimum values for the deformation parameter did not differ much from the “classical” $`q1`$ limit, an overall improved fit to the experimental excitation energies was achieved with $`q1`$. Nevertheless, such applications can contradict the physical interpretation of the generators of the algebra when the generators are associated with fundamental symmetries of the system . In particular, when one $`q`$-deforms $`SO(3)SU(2)`$, rotational invariant of the system is compromised. But the situation can be very different if one considers the algebra of a many-body interaction. In such scenarios, the $`q`$-deformation accounts for non-linear contributions of higher-order interactions  without affecting physical observables associated with a subset of the generators of the algebra that remain unchanged by the deformation. This approach also allows one to construct $`q`$-deformed nuclear Hamiltonians with exact solutions. By considering a $`q`$-deformed generalization of some nuclear structure models, the role of the deformation can be explored by comparing the “classical” and $`q`$-deformed results with the experimental data.
It is well-known that effective interactions in nuclei are dominated by pairing and quadrupole terms. The former accounts for the formation of fermion pairs that give rise to a superconducting/pairing gap in nuclear spectra, and the latter is responsible for the strong enhanced electric quadrupole transitions in collective rotational bands. Indeed, within the framework of the harmonic oscillator shell-model with degenerate single-particle energies, these two limiting cases have a very clear algebraic structure in the sense that the spectra exhibit a dynamical symmetry. In the pairing limit (appropriate for near closed-shell nuclei), the Kerman-Klein quasi-spin $`SU(2)`$ group  together with its dual, the unitary symplectic group $`Sp(2\mathrm{\Omega })`$ , where $`2\mathrm{\Omega }`$ is equal to the degeneracy of the shell, allows one to introduce the seniority quantum number that can be used to classify the spectra. On the other hand, in the quadrupole limit the symplectic group $`Sp(6,)`$  governs a shape-determined dynamics.
The inclusion of a proton-neutron isovector (isospin $`T=1`$) pairing interaction enlarges the group structure from the Kerman-Klein $`SU(2)Sp(2)`$ to $`SO(5)Sp(4)`$ . And a further extension to the $`SO(8)`$ group  allows one to also take the competing isoscalar ($`T=0`$) proton-neutron pairing mode into account . In this letter we focus on isovector pairing correlations in nuclei and, in addition, include a proton-neutron isoscalar term in the interaction that is diagonal in an isospin basis. The simple two-body $`Sp(4)`$ pairing model gives a very reasonable estimate for the nuclear interaction strength, which we assume to be constant for all nuclei within a major shell. The new feature reported on in this article is an extension of the theory to include non-linear deviations from the pairing solution as realized through a $`q`$-deformation of the underlying $`Sp(4)`$ symmetry group.
Our results show that the $`q`$-deformation parameter is decoupled from the interaction strength. In addition, since the $`q1`$ results are uniformly superior to those of the non-deformed limit, the results also suggest that the deformation has physical significance over-and-above the simple pairing gap concept, extending to the very nature of the nuclear interaction itself. In short, our results suggest that the $`q`$-deformation has physical significance beyond what can be achieved by simply tweaking the parameters of a two-body interaction. The results also underscore the need for additional studies to achieve a more comprehensive understanding of $`q`$-deformation in nuclear physics.

2. The $`q`$-deformed pairing model — The $`𝔰𝔭_q(4)`$ deformed algebra  is constructed in terms of $`q`$-deformed creation and annihilation operators, $`\alpha _{jm\sigma }^{}`$ and $`\alpha _{jm\sigma }`$, each of which creates and annihilates a nucleon with isospin $`\sigma `$ ($`+\frac{1}{2}`$ for proton, $`\frac{1}{2}`$ for neutron) in a single-particle state of total angular momentum $`j`$ (half-integer) with projection $`m`$ along the $`z`$-axis. The $`q`$-deformed fermion operators are defined through their anticommutation relations, $`\{\alpha _{jm\sigma },\alpha _{jn\sigma }^{}\}_{\pm 1}=q^{\pm \frac{N_{2\sigma }}{2\mathrm{\Omega }}}\delta _{mn},\{\alpha _{jm\sigma },\alpha _{kn\tau }^{}\}=0(jk,\sigma \tau )`$, and $`\{\alpha _{jm\sigma }^{},\alpha _{kn\tau }^{}\}=\{\alpha _{jm\sigma },\alpha _{kn\tau }\}=0`$ , where the $`q`$-anticommutator is $`\{A,B\}_\kappa =AB+q^\kappa BA`$ and $`N_{2\sigma }=\left(_{jm}\alpha _{jm\sigma }^{}\alpha _{jm\sigma }\right)_{q=1}`$. In a model with degenerate single-particle levels, the dimension of the space for given $`\sigma `$ is $`2\mathrm{\Omega }=_j(2j+1)`$. The non-deformed (“classical”) version of the theory is obtained in the limit $`q1`$. The generators of $`Sp_q(4)`$ are constructed as a bilinear product of fermion operators coupled to total angular momentum and parity $`J^\pi =0^+`$,

$`T_\pm `$ $`={\displaystyle \frac{1}{\sqrt{2\mathrm{\Omega }}}}{\displaystyle \underset{jm}{}}\alpha _{jm,\pm {\scriptscriptstyle \frac{1}{2}}}^{}\alpha _{jm,{\scriptscriptstyle \frac{1}{2}}},`$
$`A_\mu ^{}`$ $`={\displaystyle \frac{1}{\sqrt{2\mathrm{\Omega }(1+\delta _{\sigma \tau })}}}{\displaystyle \underset{jm}{}}(1)^{jm}\alpha _{jm,\sigma }^{}\alpha _{j,m,\tau }^{},`$ (1)
$`A_\mu `$ $`=(A_\mu ^{})^{},\mu =\sigma +\tau `$
in addition to the proton (neutron) number operators $`N_{\pm 1}`$, which remain undeformed. These generators of the Cartan subalgebra can also be realized in terms of the isospin projection operator $`T_0=\frac{1}{2}(N_{+1}N_1)`$ and the total nucleon number operator $`N=N_{+1}+N_1`$. The $`𝔰𝔭_q(4)`$ algebra contains four distinct $`𝔰𝔲_q(2)`$ $`q`$-deformed subalgebras (see table 1). The commutation relations between the generators are symmetric with respect to the exchange $`qq^1`$ .

As for the microscopic “classical” approach , the most general Hamiltonian of a system with $`Sp_q(4)`$ dynamical symmetry and conserved proton and neutron particle numbers can be expressed in terms of the generators (1)

$`H_q`$ $`=\epsilon _qNF_q(A_{+1}^{}A_{+1}+A_1^{}A_1)G_qA_0^{}A_0`$
$`2C_q([\frac{N}{2}\mathrm{\Omega }]^2[\mathrm{\Omega }]^2)D_q[T_0]^2\frac{E_q}{2}\left(\{T_+,T_{}\}\left[\frac{N}{2\mathrm{\Omega }}\right]\right),`$ (2)
where $`[X]^2=\mathrm{\Omega }\left[\frac{X}{2\mathrm{\Omega }}\right]\left(\left[X+1\right]_{\frac{1}{2\mathrm{\Omega }}}+\left[X1\right]_{\frac{1}{2\mathrm{\Omega }}}\right)\stackrel{q1}{}X^2`$ and by definition $`[X]_\kappa =\frac{q^{\kappa X}q^{\kappa X}}{q^\kappa q^\kappa }\stackrel{q1}{}X`$. In principle, the parameters $`\gamma _q=\{\epsilon _q,F_q,G_q,C_q,D_q,E_q\}`$ can be different from their non-deformed counterparts $`\gamma =\{\epsilon ,F,G,C,D,E\}`$. The model describes the motion of $`N_+`$ valence protons and $`N_{}`$ valence neutrons in the mean-field of a doubly-magic nuclear core. The basis states are constructed by the action of the pair-creation operators $`A_{0,\pm 1}^{}`$, on the vacuum. The quantum numbers used to specify the basis set $`|n_1,n_0,n_1)`$ count the number of $`pp`$, $`pn`$ and $`nn`$ pairs, respectively, and are associated with the corresponding limits of $`𝔰𝔭_q(4)`$ (see table 1). The basis vectors model $`0^+`$ states with dominant isovector pair correlations. In the mass $`40A100`$ region they can be used to describe $`0^+`$ ground states of even-$`A`$ nuclei and the lowest isobaric analog $`0^+`$ state in the odd-odd nuclei with a $`J0`$ ground state.
From a “classical” perspective, the deformation introduces higher-order, many-body terms into a theory that starts with only one-body and two-body interactions, the latter including an isovector pairing interaction (parameters $`F,G`$) together with a proton-neutron isoscalar force diagonal in an isospin basis. The way in which the higher-order effects enter into the theory is governed by the $`[X]`$ form. Since $`\varkappa `$ and $`q`$ are related to one another, $`q=e^\varkappa `$, everything is tied to the deformation with $`[X]=\frac{\mathrm{sinh}(\varkappa X)}{\mathrm{sinh}(\varkappa )}=X(1+\varkappa ^2\frac{X^21}{6}+\varkappa ^4\frac{3X^410X^2+7}{360}+\mathrm{})X`$ in the $`\varkappa 0`$ ($`q1`$) limit. The deformation is applied to the same region of nuclei where the “classical” model has already proven to provide for a reasonable description of the $`0^+`$ states under consideration. Thus the $`q`$-deformation does not remedy the non-deformed model but complementarily can improve it.

3. Novel properties of $`q`$-deformation — To explore the physics of $`q`$-deformation, we fit the eigenvalues of the deformed Hamiltonian to the relevant experimental $`0^+`$ state energies for groups of 36 and 100 nuclei. This was carried out in two steps. First we determined a set of $`\gamma `$ parameters for the non-deformed limit ($`q=1`$) of the theory that yielded a best overall fit to the data. Then another global fit was made with the $`\gamma _q`$ and $`q`$ parameters allowed to vary. The $`\gamma _q`$ set that was found differed very little from that of the non-deformed case. In short, varying the deformation parameter affected the pairing strengths very little. This observation underscores the fact that the deformation represents something fundamentally different, a feature that cannot be “mocked up” by allowing the strengths of the non-deformed interaction to absorb its effect.
Once the general $`\gamma `$ parameters were determined, the deviation of the predicted $`0^+`$ state energy from the corresponding experimental number, $`\left|H_qE_{\mathrm{exp}}\right|^2`$, was minimized with respect to $`q`$ for each nucleus. This procedure yielded either two symmetric solutions for $`\varkappa `$ ($`H_q=E_{\mathrm{exp}}`$), that is one physical solution $`|\varkappa |`$, or one value of $`q`$ at the minimum of the $`q`$-deformed energy $`H_q`$ (see Fig. 1). In the second case the minimum occurs at the “classical” energy ($`\varkappa =0`$) and its difference from the experimental value can be attributed to the presence of other types of interactions that are not in the model.

The higher-order terms, which correspond to many-body interactions, can be recognized through the expansion of the eigenvalues of the $`q`$-deformed Hamiltonian (2). As an example, in the limit of identical pairing interactions, ı.e. $`𝔰𝔲_q^\pm (2)`$ subalgebra (table 1), and in the $`pn`$ pairing limit, ı.e. $`𝔰𝔲_q^0(2)`$ limit, the ratio between the $`q`$-deformed $`E_q^\pm `$ ($`E_q^{pn}`$) and “classical” $`E^\pm `$ ($`E^{pn}`$) energies is

$`R^{pp(nn)}=\frac{E_q^\pm }{E^\pm }=`$ $`1+\frac{\varkappa ^2}{6}\left(\frac{8n_{\pm 1}^25+4\mathrm{\Omega }^2}{8\mathrm{\Omega }^2}+\left(\frac{E^\pm }{Fn_{\pm 1}}\right)^2\right)+\mathrm{},`$ (3)
$`R^{pn}=\frac{E_q^{pn}}{E^{pn}}=`$ $`1+\frac{\varkappa ^2}{6}\left(\frac{n_0^214\mathrm{\Omega }^2}{4\mathrm{\Omega }^2}+\left(\frac{E^{pn}}{Gn_0}\right)^2\right)+\mathrm{},`$ (4)
where the expansions include higher-order terms that may not be negligible and the like-particles limit (3) can be compared with the earlier studies . While the quadratic coefficient in $`R^{pp(nn)}`$ is positive, the one in $`R^{pn}`$ is negative. This leads to a decrease of the binding energy of the $`pn`$ pairs as $`|\varkappa |`$ increases from zero. As the deformation parameter increases from the “classical” limit, the like-particle pairing is strengthen, yielding a larger pairing gap.

4. Analysis of the role of the $`q`$-deformation — The analysis yields solutions for the deformation parameter $`|\varkappa |`$ that fall on a smooth curve that tracks with the energy of the lowest $`2^+`$ states (see Fig. 2). These energies are largest near closed shells where the pairing effect is essential for determining the low-lying spectrum and decrease with increasing collectivity and shape deformation. Similar properties are suggested for the $`q`$-deformation and this result, even though qualitative, gives some insight into the understanding of the nature of the $`q`$-deformation. The observed smooth behavior of the deformation parameter reveals its functional dependence on the model quantum numbers.
The many-body nature of the interaction is most important around closed shells and the regions with $`N_+N_{}`$. For these nuclei the $`q`$-parameter has significant values and the experimental energies can be reproduced exactly. An interesting point is that $`q`$ tends to peak for even-even nuclei when $`N_+=N_{}`$ where strong pairing correlations are expected (see Fig. 2).
“Classical” values of the $`q`$-deformation parameter ($`q1`$) are found in nuclei with only one or two particle/hole pairs from a closed shell. This is an expected result since the number of particles is insufficient to sample the effect of higher-order terms in a deformed interaction. For these nuclei the non-deformed limit gives a good description. Around mid-shell ($`N2\mathrm{\Omega }`$) the deformation adds little improvement to the theory with the experimental values remaining close to the “classical” limit. This suggests that for these nuclei the pairing interaction is not sufficient for their description. The results imply that even though the $`q`$-parameter gives additional freedom for all the nuclei, it only improves the model around regions of dominant pairing correlations.

A $`q`$-deformed extension of the $`Sp(4)`$ model, which is the underlying symmetry for describing isovector ($`T=1`$) pairing correlations in atomic nuclei, has been investigated. When compared to experimental data the theory shows a smooth functional dependence of the deformation parameter $`q`$ on the proton and neutron numbers, which resembles the behavior of the lowest $`2^+`$ state energies. Since a $`q`$-deformation of $`Sp_q(4)`$ introduces higher-order, non-linear terms in the $`np,`$ $`pp`$ and $`nn`$ pairing interactions into the nuclear Hamiltonian, the outcome suggests the presence and importance of higher-order pairing correlations in nuclei, especially for nuclei just beyond closed shells and with $`N_+N_{}`$. The results also show that the deformation is decoupled from parameters that are used to characterize the two-body interaction itself, which means the latter can be assigned best-fit global values for the model space under consideration without compromising overall quality of the theory. Moreover, the specific features of the nuclear structure can be investigated through the use of a local $`q`$ value that varies smoothly with nuclear mass number. In summary, the concept of quantum deformation has been linked to the smooth behavior of physical phenomena in atomic nuclei.

This work was supported by the US National Science Foundation, Grant Numbers 9970769 and 0140300.