text:
# Disease Spreading in Structured Scale-Free Networks

## 1 Introduction

During the last years, there has been a burst of activity in the study of complex networks barabasi02 ; dorogorev . It has been shown that many social and natural systems falou99 ; alexei ; www99 ; strog01 ; montoya02 ; wagner01 ; jeong01 ; spsk ; vazquez  can be represented as a graph where nodes represent individuals or agents and links stand for the physical interactions among them. Surprisingly, many of these networks share some important topological features such as small-world (SW) properties watts98  and scale-free (SF) degree distributions barab99 , where the degree or connectivity $`k`$ of a node is the number of neighbors to which it is linked. Networks displaying scale-free degree distributions are very interesting not only for their relative abundance in Nature but also because of their peculiar statistical properties. In particular, the unbounded fluctuations of the connectivity distribution $`P(k)k^\gamma `$ ($`\gamma 3`$) seem to be a blueprint of all real-world networks studied so far barabasi02 ; dorogorev .

The SW and SF properties have considerable impact on the processes running on top of complex networks. The effects of the complex topological features of SF networks on the dynamics of epidemic disease spreading is perhaps one of the most interesting outcomes in the study of complex networks. Specifically, the understanding of spreading phenomena in these networks can shed light on a large number of practical problems, ranging from computer virus infections newvir ; white  to epidemiology murray ; anders . For instance, it has been recognized only recently pv01a ; pv01b  that random SF networks are completely prone to epidemic spreading allowing the onset of large epidemics whatever the spreading rate of the infection. This radical change with respect to the way a disease is spread in a regular structure is rooted in the diverging connectivity fluctuations of SF networks with $`\gamma 3`$. Similarly, random SF networks have been shown to exhibit extremely robustness to random damages barabasi00 ; newman00 ; havlin01 .

On the other hand, real networks are also characterized by degree correlations that might play a fundamental role in the functional properties of networked systems and on processes running on correlated networks. The study of networks with degree correlations and the extension of previous results obtained for random networks to correlated networks is quite recent berg02 ; marian1 ; marian3 ; n02a ; vw02 ; cond182 . One of the models aimed to account for degree correlations was proposed recently by Klemm and Eguíluz klemm02  (henceforth referred to as deactivation model) where nodes can be deactivated with probability inversely proportional to their connectivity. The model turns out to be more rich than initially thought with a power law connectivity distribution $`P(k)k^\gamma `$ but with $`\gamma `$ laying in a quite wide interval $`2\gamma 4`$ cond183 . Additionally, it has been shown that the topology of the network is essentially linear being the graph a collection of stars of diverse degree connected as a chain by a number of local links, i.e., the graph lacks SW properties. While a threshold may exist structured  due to the peculiar topological properties of the network that makes the epidemic spreading to be dominated by the diffusion of the disease on a linear chain cond183 , the star-like graphs connected as a chain might also lead to new effects in the dynamics of epidemic spreading.

In this paper, we study in detail by large scale numerical simulations two paradigmatic epidemiological models, namely, the Susceptible-Infected-Susceptible (SIS) and the Susceptible-Infected-Removed (SIR) models on top of networks generated using the deactivation model. We found that the existence or not of an epidemic threshold for $`\gamma 3`$ depends on the initial density of infected individuals while for values of $`\gamma >3`$ the epidemic threshold is recovered and it is universal, i.e., it does not depend on the initial density of infected individuals. We also provide analytical arguments in order to explain our numerical findings. Moreover, we show that the linear character of the graph determines the existence of a trivial threshold for the SIR model that can be mapped to a bond percolation problem in one dimension.

## 2 Deactivation model

Let us first summarize the main properties of the deactivation model introduced by Klemm and Eguíluz klemm02 . The growing dynamics of the network is defined in the following way: One starts from a fully connected graph of $`m`$ nodes that are set active. At successive time steps nodes are added one by one following the recipe: $`(i)`$ A new node is connected to all active nodes in the network; $`(ii)`$ one of the active nodes is selected for deactivation with probability

$$p_d(k_i^{\mathrm{in}})=\frac{\left[_{j𝒜}(a+k_j^{\mathrm{in}})^1\right]^1}{a+k_i^{\mathrm{in}}};$$
(1)
where the sum in Eq. (1) runs over the set of active nodes $`𝒜`$, $`a`$ is a model parameter, and $`k_i^{\mathrm{in}}`$ denotes the in-degree of the $`i`$-th node. Finally, $`(iii)`$ the new node just added is set active.

Recently, it has been shown cond183  that most of the topological properties (in particular, the degree distribution) of the network generated using the above rules are very sensitive to the order in which steps $`(ii)`$ and $`(iii)`$ are performed. Thus we shall discriminate in the following two cases cond183 :

* Model A: $`(ii)`$ is performed before $`(iii)`$ .
* Model B: $`(iii)`$ is performed after $`(ii)`$.

The deactivation model is usually run taking $`a=m`$. This makes the deactivation probability inversely proportional to the total connectivity of the nodes $`k=m+k^{\mathrm{in}}`$. The connectivity distribution can be analytically obtained for the limiting cases of lowest cond183  and largest $`m`$ structured  for each model resulting in a power law that reads as,

$$P(k)k^\gamma ,$$
(2)
where the exact value of $`\gamma `$ depends on the model considered and the value of $`m`$ such that,

* model A with $`a=m`$ $`3<\gamma 4`$
* model B with $`a=m`$ $`2\gamma <3`$.

Noticeable and relevant for the physical processes run on top of this model, the dynamics of the deactivation model allows exponents $`\gamma `$ of the connectivity distribution that can lead to an unbounded second moment ($`k^2\mathrm{}`$ for model B) or to finite connectivity fluctuations ($`k^2`$ is finite for model A). More important, however, is the fact that the generated networks lack SW properties. Small-world properties refer to the fact that in many complex networks one can go from one node to any other node of the system visiting a very small number of intermediate nodes. The minimum number of such intermediate nodes for each pair $`(i,j)`$ of nodes in the network is called the minimum path length between $`i`$ and $`j`$, and the diameter of the network is defined as the largest among the shortest paths between any two nodes in the network. For networks displaying SW properties, the graph diameter grows, at most, logarithmically with the system size $`N`$.

The networks generated with the deactivation model consist of a collection of star-like graphs connected as a chain, resembling the topology of a one dimensional lattice. In fact, irrespective of the version considered, both the diameter and the average minimum path distance scale linearly with the network size $`N`$ cond183 . In particular, the mean square displacement of a random walker on the deactivation model scales with time almost as for a one dimensional lattice qualifying for a slightly subdiffusive behavior cond183 . Hence, the existence of an epidemic threshold in the deactivation model is not surprising. However, at the same time, the peculiar characteristics of the graph and, in particular, the presence of stars with very large connectivities introduces new and interesting effects in the dynamics of disease spreading. As we shall discuss in the following section, whether or not a finite prevalence sets in for the SIS model would depend on the initial density of infected individuals.

## 3 The SIS Model

In the standard SIS epidemiological model murray , each node of the network represents an individual and each link is the connection along which the individuals interact and the epidemic can be transmitted. In this model, the individuals can exist in two possible states, namely, susceptible or healthy and infected. Individuals removal due to death or acquired immunization is not allowed and thus they stochastically move through the cycle susceptible $``$ infected $``$ susceptible. The disease transmission is described in an effective way. At each time step, each susceptible node connected to one or more infected nodes gets the infection with probability $`\lambda `$, while the infected nodes are cured and become again susceptible at a rate $`\delta `$ (henceforth set to 1 without lack of generality). In the SIS model whenever the epidemic pervades the system, it gets into an endemic state with a stationary value for the density of infected nodes that acts as the order parameter of the model.

Results on regular structures point out that there is a threshold below which the system does not reach such a stationary state and the infection dies out murray ; pv01b . Moreover, the value of the prevalence at the stationary state and its own existence does not depend on the initial density of infected individuals. One can start with an initial density $`\rho _0=1/N`$ of infected individuals or from, let’s say, half of the lattice infected (in general, $`\rho _0=const`$): The prevalence self-organizes into a stationary state in both cases. It is worth noting that both situations are feasible in practice and thus physically relevant.

On the other hand, the behavior of the SIS model on random SF networks is radically different pv01a ; pv01b  when $`\gamma 3`$. In this case, there is no epidemic threshold in the thermodynamic limit and the networks are completely prone to the spreading of the disease. This drastic change of behavior is due to the diverging connectivity fluctuations of SF networks. Moreover, it has also been shown that for finite-size system an effective threshold is recovered, but its value is significantly overestimated fspv .

In order to study the SIS model in deactivation model networks, we first generate networks using the algorithm described in the preceding section for both settings of the model A and B with $`a=m=3`$. Then we let the system evolve according to the SIS dynamics. Large scale numerical simulations were performed in networks of size up to $`N=10^6`$, averaging over at least 100 different realizations on at least 10 different realizations of the network. Initially, a fraction of nodes $`\rho _0`$ was infected and we let the system relax into the steady state where the prevalence $`\rho `$ attains its stationary value. Two initial conditions for the density of infected individuals were considered. In the first case, we start from a single infected individual, $`\rho _0=1/N`$, and set the observation time $`t_m`$ to be $`10^6`$ time steps, that is, the state of each node is updated up to $`t_m`$ times if there is at least one infected individual. Then, we repeat the process but starting with homogeneous initial conditions. This is achieved by infecting a finite fraction of the network $`\rho _0=const.`$ and monitoring as before whether a stationary state sets in with its corresponding prevalence for the same observation time $`t_m`$. It is worth recalling that the prevalence in the stationary state is computed as the average over all surviving trials in both cases.

The results obtained are depicted in Fig. 1 where the steady state prevalence has been drawn as a function of the spreading rate $`\lambda `$. As expected for model A the epidemic threshold would appear again for both initial conditions since in this case together with the lack of SW properties and the linearity of the network, the second moment $`k^2`$ is finite with $`3<\gamma 4`$. In this case, the phase diagram is the same regardless of the initial density of infected individuals $`\rho _0`$. It is worth noting that whenever an epidemic threshold exists, its value seems to be the same (within statistical errors) for both versions of the deactivation model, depending only on the average connectivity $`k`$. This implies in its turn that the dynamics of the infection spreading is mainly determined by the linear star-like structure of the graphs.

Two radically different behaviors can be noted for the deactivation model B depending on the initial condition. In this case the exponent of the power law decay of the degree distribution is in the interval $`2\gamma <3`$. In particular a best fit analysis for networks of size $`N=10^7`$ yields an exponent $`\gamma =2.65\pm 0.05`$ for $`m=3`$. When the system evolves from a single infected individual, an epidemic threshold $`\lambda _c`$ appears, below which the epidemic cannot pervade the network. In this case the fact that in the thermodynamic limit the connectivity fluctuations are diverging is canceled by the effects introduced by the linear topological nature of the network and the lack of SW properties. As advanced in Ref. cond183 , the SIS dynamics can be reduced in a coarse grained picture to the diffusion of the disease on a linear chain. Thus, the classical picture for regular structures with an epidemic threshold is recovered.

For homogeneous initial conditions, the phase diagram is completely different. In this case, the stars with high connectivity can let the epidemic survive and a final stationary state is reached by a long power-law decay. In order to rule out the presence of finite size effects hiding an abrupt transition we have checked the behavior of the stationary prevalence for several system sizes. Note that larger system sizes are needed in order to correctly depict the prevalence behavior for $`\lambda 1`$, since finite size networks induce finite size corrections to the zero threshold pv01a ; pv01b . The results obtained are plotted in Fig. 2. The straight line is a fit to the form $`\rho (\lambda )\lambda ^\beta `$, with $`\beta 4.5`$, showing that the infection prevalence is assuming a finite stationary value for all values $`\lambda >0`$. This numerical evidence confirms the previous picture obtained for uncorrelated scale-free networks, namely, the absence of an epidemic threshold in the thermodynamic limit. It is interesting to recall that using a dynamical mean field approach for uncorrelated networks with scale-free connectivity distributions given by $`P(k)k^\gamma `$, the same functional form $`\rho (\lambda )\lambda ^\beta `$ is predicted when $`2<\gamma <3`$, although in this case satisfying the relation $`\beta =1/(3\gamma )`$ pv01b . Hence, it seems that correlations, while preserving the general functional form, make the prevalence to decay more faster than it does in the uncorrelated case. This in turn may explain why in finite systems the epidemic threshold is smaller than its counterpart in uncorrelated networks, a fact also observed in random network models marian1 . Below, we provide an analytical argument that help understand the origin of these differences by analyzing the SIS model on an ensemble of disconnected stars.

### 3.1 SIS model on an ensemble of disconnected stars

Let us study a graph made up of an ensemble of disconnected stars. To be more precise, we consider a graph where all nodes with degree $`k>1`$ are connected to nodes with degree $`k=1`$. If $`P(k)`$ ($`k1`$) is the degree distribution then $`P(1)=_{k>1}kP(k)`$. Since the stars are isolated the stationary prevalence is given by

$$\rho =\underset{k>1}{}P(k)n(k)v(k),$$
(3)
where $`n(k)`$ is the stationary number of infected nodes on a star of degree $`k`$ and $`v_k`$ is the probability that at $`t=0`$ there was at least one infected node in a star with $`k`$ leaves.

Now, to compute $`n(k)`$ we focus on the SIS dynamics on a single star. The spreading takes place from the central node to the leaves, then from the leaves to the central node and so on. For stars with a large number of leaves ($`k1`$) we can approximate the average number of infected nodes by the number of infected leaves. Hence, the number of infected nodes at step $`t+1`$ is

$$n(k,t+1)=\left[1(1\lambda )^{n(k,t)}\right]\lambda k.$$
(4)
The first factor in the right hand side is the probability that the central node gets infected, receiving the disease from at least one leaf. The second factor is the number of infected leaves given that the center was infected. This linear map has always the trivial solution $`n(k,t)=0`$. To investigate its stability we assume $`n(k,t)1`$ resulting

$$n(k,t+1)=\lambda k\mathrm{ln}\frac{1}{1\lambda }n(k,t).$$
(5)
Thus, there is a critical star size

$$k_c(\lambda )=\left(\lambda \mathrm{ln}\frac{1}{1\lambda }\right)^1,$$
(6)
such that for $`k<k_c`$ the prefactor in the right hand side of Eq. (5) is smaller than 1 and, therefore, $`n(k)=0`$. On the contrary, for $`k>k_c`$ there is an exponential growth of $`n(k,t)`$ indicating that the solution $`n(k,t)=0`$ is not stable, i.e. $`n(k)>0`$. Moreover, the critical degree depends on $`\lambda `$, with the limiting cases $`k_c(0)=\mathrm{}`$ and $`k_c(1)=0`$.

Going back to Eq. (3) we obtain that

$$\rho =\underset{k>k_c}{}P(k)n(k)v(k).$$
(7)
Let us now distinguish between two different initial conditions considered before. In one case we start from only one infected node. The probability that one node of a star of $`k`$ leaves is infected is

$$v^{(1)}(k)=\frac{k+1}{N}.$$
(8)
The substitution of this expression in Eq. (7) yields

$$\rho ^{(1)}=\underset{k>k_c}{}P(k)n(k)\frac{(k+1)}{N}\frac{k^2}{N}.$$
(9)
Now, for $`2<\gamma <3`$ the second moment grows at most as $`k^2=𝒪(N^{3\gamma })`$. Hence, in the thermodynamic limit the prevalence is asymptotically zero. However, a different result is obtained if we infect a finite fraction $`\rho _0`$ of the nodes. In this case the probability that at least one node of a star with $`k`$ leaves is infected is given by

$$v^{(2)}(k)=1(1\rho _0)^{k+1}\rho _0.$$
(10)
From this expression and Eq. (7) it follows that

$$\rho ^{(2)}\rho _0\underset{k>k_c}{}P(k)n(k).$$
(11)
Thus, there is a finite prevalence in the stationary state.

These two opposite scenarios are the same observed in the numerical simulations of the preceding section. When there is a finite fraction of infected nodes at $`t=0`$, a finite prevalence is obtained. On the contrary, when the initial infection is concentrated in only one node, the stationary state reflects the topological nature of the graph. For the ensemble of disconnected stars we obtain a stationary state with no infected node, corresponding to the SIS model in a zero-dimensional system. On the other hand, for the deactivation model, we get a phase transition at a finite infection rate as expected for a one-dimensional system marro99 .

## 4 SIR model

The SIR model assumes that individuals can exist in three classes: susceptible, infected and removed. The main difference with the SIS model is that once an individual gets infected it is removed and, therefore, it can not catch the infection again. If we start from a single infected node the size of the outbreak of the disease is given by the number of nodes that can be reached assuming that each link is occupied with a probability $`\lambda `$. Thus, the SIR model is equivalent to a bond percolation problem with bond occupation probability $`\lambda `$ newman02a . Moreover, the size of the outbreak is just the size of the giant component.

The study of the SIR model in random SF networks confirmed the epidemiological picture obtained for the SIS model in complex networks with power-law connectivity distributions lloydsir ; moreno02 . For instance, it has been shown that the effective epidemic threshold is inversely proportional to the connectivity fluctuations $`k^2`$ and hence it is vanishing in the thermodynamic limit $`N\mathrm{}`$ for $`\gamma 3`$. The high heterogeneity of SF networks also causes that the relative incidence of an outbreak strongly depends on the connectivity of the first infected nodes moreno02 . This dependency on the initial seed should not be confused with the dependency on the initial conditions found previously since for the SIR model the differences in the relative incidence of an epidemic outbreak is an intrinsic effect of the large heterogeneity in the connectivity distribution of SF networks and is not related to the presence of any topological constrain.

We have made numerical simulations of the bond percolation problem on top of the structured networks generated using the two variants of the deactivation model, A and B. Once the graphs are generated each link is removed with a probability $`1\lambda `$. Then, the size of the giant component of the resulting graph is computed. All the results reported below were obtained taking an average over 10 graph realizations of the deactivation model and 10 realizations of the link removal procedure. In Fig. 3 we plot the size of the giant component $`S`$ as a function of $`\lambda `$ for different graph sizes $`N`$ for models A and B, with $`m=3`$. Notice that in both cases the qualitative picture is the same.

Between the limiting cases $`S(0)=0`$ and $`S(1)=1`$, there is an intermediate range of $`\lambda `$ where $`S`$ goes from a value close to 0 to a value close to 1. The width of this interval decreases appreciable with increasing $`N`$. Moreover, the point at which the transition takes place systematically shifts to larger values of $`\lambda `$ approaching 1. Since the deactivation model has essentially a one dimensional topology we expect that in the large $`N`$ limit $`S=0`$ for any $`\lambda <1`$. This hypothesis can not be confirmed by direct observation of Fig. 3 but it can be checked using a finite size scaling analysis. When we plot $`S`$ as a function $`1\lambda `$ in a log-linear scale we observe equidistant shifts of the curves after increasing $`N`$ from $`10^4`$ to $`10^5`$ and then to $`10^6`$. This scenario corresponds to a scaling of the form $`S=f[(1\lambda )N^\alpha ]`$, where $`f(x)`$ is a scaling function that is independent of $`N`$ and $`\alpha `$ is a scaling exponent.

In Fig. 4 we plot $`S`$ versus $`(1\lambda )N^\alpha `$ in a log-linear scale using the value of $`\alpha `$ that gives the best data collapse. This was achieved for $`\alpha =0.15(1)`$ and $`\alpha =0.14(1)`$ for models A and B, respectively. Within the numerical error, the exponent $`\alpha `$ is the same for both models. Moreover the scaling function $`f(x)`$ is also unique, up to a constant factor. These results point out that with regard to the bond percolation problem, or the SIR model, the differences in the degree distributions of model A and B are qualitatively irrelevant.

## 5 Discussion and conclusions

We have study the dynamics of infectious diseases in structured scale-free networks generated using the deactivation model introduced by Klemm and Eguíluz. By means of large scale numerical simulations, we have shown that the existence or not of an epidemic threshold for the SIS model depends on the initial density of infected individuals when the connectivity fluctuations of the network are unbounded. This effect is completely new and reflects the peculiar topological nature of the networks generated with the deactivation model algorithm. Guided by the analytical solution to the SIS model in an ensemble of disconnected stars, we may summarize the SIS dynamics as follows.

The linear topological nature of the networks provoke the dynamics of the epidemic spreading to be almost purely diffusive, like in a one-dimensional lattice. However, for model B the small value of $`\gamma `$, that leads to the divergence of the connectivity fluctuations, allows the existence of stars with very large connectivities such that for any value of $`\lambda `$ the infection can always get trapped in these stars. Hence, starting from homogeneous initial conditions, the probability of hitting a star in which the epidemic survives is not zero. This is not certainly the case if one starts at a single infected node. In this case, there is no time for the infection to spread through the network. It would be of further interest to study this diffusion process in more details, for example, by looking at the distribution of waiting times for each node in the cicle infected $``$ susceptible $``$ infected. On the other hand, when the connectivity fluctuations are finite (model A), the probability of finding a star with degree $`k>k_c(\lambda )`$ is zero below $`\lambda _c`$ and thus the epidemic threshold is recovered for all possible initial conditions.

Finally, the study of the SIR model confirms that the dominant factor in this case is the chain structure of the network with star-like nodes connected locally. The threshold for this epidemiological model coincides, as it should, with the critical point of the corresponding one-dimensional bond percolation problem that gives an epidemic threshold $`\lambda _c=1`$.

In summary, we have provided evidences that in this kind of networks the dominant factor determining the behavior of both the SIS and the SIR models is the linear topology of the network. The existence of high degree nodes introduces new effects such as the dependency on the initial conditions for the deactivation model B with $`\gamma 3`$. As a final remark, we should say that the results for the SIS and the SIR models here reported cannot be directly extended to all correlated networks. The peculiar topological features of the networks generated with the deactivation model make them unique within the class of correlated networks. Some recent works have begun to marian3 ; cond182 ; cond450  address the SIS and the SIR models in correlated networks establishing the general conditions for the existence or not of an epidemic threshold. The extensive numerical exploration of these models in such networks is still to do.

###### Acknowledgements.

We thank A. Vespignani and R. Pastor-Satorras for helpful comments and discussions. This work has been partially supported by the European Commission - Fet Open project COSIN IST-2001-33555.


format: text

text:
# The Environmental Dependence of Star Formation in the Nearby Universe

## 1. Introduction

Studies of the effects of the large scale environment around a galaxy on the properties of its star formation are essential for understanding how galaxies evolve. It is well known since long time that galaxies with large star formation rates avoid dense environments, and many recent works have confirmed that star formation is lower for galaxies in clusters and groups than for those in the field (e.g., Ellingson et al. 2001, Martínez et al. 2002). Several different process have been proposed to explain why and how star formation is affected by the environment, like mergers, ram pressure, and tidal interactions. While all of them may take place in some form or another, the main process that affects star formation in galaxies is not clear yet.

In the current era of large redshift surveys, enormous amounts of homogeneous data on nearby and distant galaxies are becoming available and may be used to study statistically many aspects of star formation in galaxies. Here we use the first 100k Data Release of the 2dF Galaxy Redshift Survey (2dFGRS; Colless et al. 2001) to examine some properties of star-forming galaxies in the nearby universe.

## 2. A Nearby Galaxy Sample

The sample analyzed here was extracted from the 2dFGRS. We adopt a volume limited sample because, in this case, the selection function is in principle uniform and the variations in galaxy number density are due to galaxy clustering only. We have considered galaxies with velocities between 600 km s-1 and 15000 km s-1 brighter than the extinction corrected magnitude $`b_J=18.50`$ (corresponding to $`M_{b_J}^{lim}=17.38+5\mathrm{log}h`$) and located within the SGP and NGP strips covered by the survey. Galaxy distances were computed from the redshifts, neglecting peculiar velocities and $`k`$-corrections.

The redshift completeness and its magnitude dependence vary over the area of the survey but, nevertheless, they can be represented by a mask that gives the selection function at each position in the survey area (Colless et al. 2001, Norberg et al. 2001). We have used numerical simulations to study how the survey geometry and its current incompleteness affect our estimates of the galaxy density (Mateus & Sodré, in preparation). These simulations allowed us to estimate the mean local correction factor, $`𝒞`$, that should be applied to the observed density to correct it of local incompleteness and border effects ($`𝒞1`$). To minimize these corrections and, at the same time, to have a sample as large as possible, we have considered only regions for which $`𝒞2.5`$.

## 3. Determination of the Local Galaxy Density

We characterize the environment of a galaxy by the three-dimensional galaxy number density in its neighborhood. This density is calculated from the distance $`r`$ to its $`5`$th nearest neighbor in the sample and corrected with the local value of $`𝒞`$. Since we are neglecting peculiar velocities, galaxies in dense environments, where the velocity dispersion is high, will probably have their local densities underestimated. However, only one of the cataloged galaxy clusters in the survey region (De Propis et al. 2002), S0006, is inside our selected volume. Hence, our sample is more representative of the field than of the cluster galaxy population.

The local densities (of galaxies brighter than $`M_{b_J}^{lim}`$) were computed for all the galaxies in the sample above. It can be shown that our results are not affected, qualitatively, by the density corrections.

## 4. Spectral Classes

Spectra obtained in modern redshift surveys like the 2dFGRS have much more information than that needed to obtain just redshifts, and can be used to provide informations on the galaxy type or star formation activity. The use of equivalent widths (EWs) of certain spectral features are particularly useful in the study of survey data because they are not affected by reddening and are relatively insensitive to the instrumental resolution.

Following Balogh et al. (1999), we use the EWs of the \[O II\]$`\lambda `$3727 doublet and of $`H\delta `$ to classify galaxies in distinct spectral classes. The EW of \[O II\] measures the activity of star formation, whereas the EW($`H\delta `$) is sensitive to the age of star formation bursts. EWs in the rest-frame of the objects were computed and 146 of them that presented bad quality spectra or were classified as AGN (i.e, with EW($`H\alpha `$) $`>`$ 10 Å and EW(\[N II\]) $`>`$ 0.55 EW($`H\alpha `$); Lewis et al. 2002) were removed from the analysis. The remaining sample, with $`N_{ALL}=1234`$ galaxies, will be analyzed hereafter.

We consider star-forming galaxies those with EW(\[O II\]) larger than 5 Å, roughly the detection limit of this spectral feature in the survey spectra. Notice that we adopt here positive and negative values for the EW of emission and absorption lines, respectively. The spectral classes defined in the plane EW($`H\delta `$) vs. EW(\[O II\]) and discussed in this work are:

* non star-forming galaxies (P: EW(\[O II\]) $``$ 5 Å)
passive galaxies, without evidence of current significant star formation (in general E or S0);
* starburst galaxies (SB: EW(\[O II\]) $`>`$ 5Å, EW($`H\delta `$) $`>`$ 0)
galaxies where a large fraction of the light comes from a starburst that started less than $``$200 Myr ago; these galaxies are called SSB (short starburst) by Balogh et al. (1999);
* ordinary star-forming galaxies (SF: EW(\[O II\])$`>`$5Å, EW($`H\delta `$)$`<`$0)
includes most normal spirals and irregulars, that have been forming stars for several hundred million years.

With these definitions, the total number of galaxies with star formation is $`N_{TSF}=635`$, of which $`N_{SB}=66`$ and $`N_{SF}=569`$, and of those without evidence of ongoing star formation is $`N_P=599`$. Considering only galaxies with star formation, SBs correspond to 10.4% of the total.

## 5. Type vs. Local Density

Our main results are summarized in Figure 1. In all plots the data of each class was divided in six density bins with approximately the same number of galaxies per bin. Fig. 1a shows that the fraction of star-forming galaxies decreases steadily with increasing local density. At the same time, the fraction of non star-forming galaxies increases with increasing density. These trends may be seen as a consequence of the well-known density - morphology relation proposed by Dressler (1980). The fraction of SBs, however, is essentially independent of local density. Since the fraction of star-forming galaxies decrease with increasing density, in dense environments SBs become relatively more frequent among star-forming galaxies.

The relation between the mean value of EW(H$`\alpha `$) and local density, shown in Fig. 1b, also provides interesting information about the behavior of star formation in distinct environments. Passive galaxies do not have significant H$`\alpha `$ emission, as expected. For ordinary star-forming galaxies the EW(H$`\alpha `$) is essentially constant with density. For the SBs, however, the EW(H$`\alpha `$) presents a clear decrease at $`\rho 0.06h^3`$ Mpc-3. At lower densities these galaxies have mean EW(H$`\alpha `$) larger than those of SF galaxies. For higher densities, however, it drops from $``$ 30 Å to $``$ 10 Å, below the mean value of SFs, $``$ 20 Å.

## 6. Discussion

In the hierarchical scenario, as galaxy clustering evolves, the density around a galaxy tends to increase, in all environments. Higher density probably means more interactions. Our results give support to this broad view because star formation seems to be affected by the environment in all range of densities covered by our sample, as evidenced by the decrease in the fraction of star-forming galaxies with increasing local density shown in Fig. 1. Lewis et al. (2002), analyzing a different sample extracted from the 2dFGRS, find a density threshold, below which the star formation is insensitive to the galaxy environment. This is in variance with our results, which suggest that the processes that affect the star formation activity in a galaxy do act everywhere, although they may be different in distinct environments.

Most of star-forming galaxies have their star formation activity inhibited or even interrupted as density increases, but the fraction of SBs seems to be roughly independent of local density, contrarily to the expectation that denser environments favor the triggering of a starburst. However, since the fraction of star-forming galaxies decreases in dense regions, the fraction of SB galaxies relative to all star-forming galaxies is actually enhanced. This behavior have been noticed already in a sample of galaxies in the Shapley supercluster (Cuevas 2000). From a study of 8 Abell clusters, Moss & Whittle (2000) have also concluded that the fraction of SBs among spirals increases from regions of lower to higher local galaxy surface density. They noticed that in most starbursts the H$`\alpha `$ emission is circumnuclear and the galaxy morphology is disturbed, indicating the action of tidal interactions; these galaxies are the local counterparts of the Butcher-Oemler objects. Here we have shown that the fraction of SBs is essentially the same in all environments. This favors mechanisms like tidal interactions for triggering starbursts because they act everywhere, and not only in high-density environments, like ram-pressure stripping or high-speed interactions between galaxies (‘harassment’).

We have also found that the EW(H$`\alpha `$) of starbursts suffers a strong reduction in high density regions. This reduction occurs at a characteristic intergalactic distance ($`\rho ^{1/3}`$) of $`2.5h^1`$ Mpc. This is a quite large separation, when compared with the distances between the members of our Local Group.

We have inspected visually the sample of starburst galaxies using Digital Sky Survey images. A clear trend with morphology is present: at low densities most of the starbursts are apparently normal spirals or galaxies with an apparently compact morphology, whereas at high densities the starburst class contains a large fraction of tidally distorted spirals. In the lowest density bin only 1 out of 11 galaxies is clearly a tidally distorted spiral; in the highest density bin, however, 6 out of 11 are galaxies of this kind. This suggests that tidal interactions in dense environments may indeed favor the triggering of starbursts in spiral galaxies.

### Acknowledgments.

We thank FAPESP, CNPq, PRONEX, and CCINT/ USP for supporting this work.

## References

Balogh, M.L., Morris, S.L., Yee, H.K.C., Carlberg, R.G., Ellingson, E. 1999, ApJ, 527, 54

Colless, M., et al.(The 2dFGRS Team) 2001, MNRAS, 328, 1039

Cuevas, H. 2000, PhD thesis, IAG-USP

Dressler, A. 1980, ApJ, 236, 351

Dressler, A., & Gunn, J. E. 1983, ApJ, 270, 7

Ellingson E., Lin H., Yee H. K. C., Carlberg R. G., 2001, ApJ, 547, 609

Lewis, I. J. et al (The 2dFGRS Team) 2002, MNRAS, 334, 673

Madgwick

Martínez, H. J., Zandivarez, A., Domínguez, M., Merchán, M. E.,, & Lambas, D. G., 2002, MNRAS, 333, L31

Moss, C. & Whittle, M. 2000, MNRAS, 317, 667

Norberg, P., et al. (The 2dFGRS Team) 2001, MNRAS, 328, 64


text:
# Superconductivity with hard-core repulsion: BCS-Bose crossover and s-/d-wave competition

## I  Introduction

Superconductivity is possible when carriers attract. The attraction in metals is normally due to the phonon coupling. The strong coulomb repulsion is normally well screened in metals and under some conditions it can be overcome by the attraction leading to superconductivity. A locally attractive (non-retarded) model, such as the BCS one, can thus give a good qualitative description of the superconducting state since the correlation length $`\xi `$ is much larger than both the screening length and the interparticle density.

This is not the case when the carriers are non conventional. A relevant example is the single-band Kondo lattice with a RKKY interaction. For large Kondo coupling the ground state of this system is formed by local singlets: it is thus insulating at half-filling. At small doping the carriers are bachelor spins: they feel a hard core repulsion, since they cannot sit on the same site at the same time. The RKKY antiferromagnetic exchange interaction gives an attraction among holes in the singlet channel, and thus it could lead to a superconducting state Andrei ; nozieres . But in this case the hard-core repulsion cannot be neglected in the description of the $`s`$-wave superconductivity. More details on the Kondo lattice model are discussed in Reference nozieres, .

In this paper we consider a simple model where the role of the hard-core repulsion can be analyzed in a transparent way. For simplicity and relevance to physical systems we consider a two dimensional square lattice with hopping and attractive interaction among nearest neighbor sites. Different range for the interaction and the hopping produces interesting results, that will be analyzed. Similar models have been studied by different authors. Ohkawa and Fukuyama proposed an ansatz to solve the mean-field equations for superconductivity with local repulsion in Kondo systems Fukuyama . Later Micnas, Ranninger, and Robaszkiewicz Ranninger  and Bastide, Lacroix, and Rosa Simoes BastideLacroix , applied that technique to calculate the superconducting critical temperature. Similar equations in a different context were considered by Aligia and coworkers.Aligia

We will reconsider the problem at zero temperature and investigate the competition among different symmetries of the order parameter as one spans the crossover from BCS superconductivity to Bose condensation Legget ; NSR ; ranninger ; RanderiaPRL ; Haussmann ; PistolesiStrinati ; zwerger ; Levin ; Hertog ; andrenacci . In particular we will clarify an apparent disagreement with the prediction of Randeria Randeria2D  on the presence of bound states and superconducting instabilities in 2D systems. We will also discuss the unusual nature of the superconducting state: the main peculiarity is the presence of a radial line of nodes in $`k`$ space for the gap function $`\mathrm{\Delta }(𝐤)`$.

The plan of the paper is the following. In Section II we define in detail the model we are considering. In Section III we solve the two body problem and we discuss the relevance to the Bose limit. In Section IV we set up the mean field equations and find the resulting phase diagram for the ground state. The competition between $`s`$\- and $`d`$-wave superconductivity is analyzed in Section V. In Section VI the effect of a mismatch between the hopping and the interaction range is considered by introducing a hopping $`t^{}`$ among next-neighbor sites.

## II The model

We consider fermions hopping on a 2D square lattice of spacing $`a`$ with nearest-neighbor sites hopping matrix elements $`t`$. (In some cases we will consider also next-neighbor hopping $`t^{}`$.) The fermions interact with a local repulsive interaction parameterized by $`U`$, that we will set to infinity. They furthermore feel an attraction $`V`$, which we assume to act only between nearest neighbors (in a real RKKY example that would not be the case) with opposite spins. The relevant parameters at zero temperature are only $`V/t`$ and the filling of the band $`0<n<2`$. The spin does not play an essential role. Since we are interested in the singlet state, we can disregard the parallel-spin part of the interaction, the mean field equations would not be modified by that. The model Hamiltonian is thus:

$`H`$ $`=`$ $`{\displaystyle \underset{𝐤\sigma }{}}(t_𝐤\mu )c_{𝐤\sigma }^{}c_{𝐤\sigma }^{}+U{\displaystyle \underset{i}{}}(n_i{\displaystyle \frac{1}{2}})(n_i{\displaystyle \frac{1}{2}})`$ (1)
$`V{\displaystyle \underset{i,\delta }{}}(n_i{\displaystyle \frac{1}{2}})(n_{i+\delta }{\displaystyle \frac{1}{2}})`$
where the vector $`\delta `$ spans the four nearest neighbor sites of the square lattice, $`\sigma `$ is the spin projection in a chosen direction and $`t_𝐤=2t[\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)]`$ is the single-particle dispersion relation.

The Hamiltonian in the form of Eq. (1) is particle-hole symmetric: for half-filling $`\mu =0`$. For convenience it can also be written in a non p-h symmetric, but shorter form:

$$H=\underset{𝐤\sigma }{}(t_𝐤\mu _o)c_{𝐤\sigma }^{}c_{𝐤\sigma }^{}+U\underset{i}{}n_in_iV\underset{i,\delta }{}n_in_{i+\delta },$$
(2)
where $`\mu _o=\mu +(U4V)/2`$. The two forms are equivalent but from a constant. The interaction term can be written in Fourier space as usual:

$$\underset{𝐤_1+𝐤_2=𝐤_3+𝐤_4}{}V(𝐤_2𝐤_3)c_{𝐤_1}^{}c_{𝐤_2}^{}c_{𝐤_3}^{}c_{𝐤_4}^{}$$
(3)
with $`V(𝐤)=UVw(𝐤)`$ and $`w(𝐤)=_\delta e^{i𝐤\delta }=2[\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)]`$.

In the following we investigate the superconducting order in the ground state of the Hamiltonian (1). In order to treat the strong repulsion correctly we consider first the two body problem, that is relevant in the dilute limit.

## III Two-body problem

If the attractive interaction is strong enough to bind a pair we expect to reach the Bose condensed limit for a low density of carriers. It is thus crucial to know whether a bound state is present or not. In our case the hard-core repulsion will fight against an $`s`$-wave bound state, while it will not affect the $`d`$-wave state. In the dilute limit the choice between $`s`$\- and $`d`$-wave is determined simply by the lowest 2-body bound state, if present.

We thus begin our investigation by solving the two body problem: we find the threshold for the appearance of bound states. The two-body problem is exactly solvable by diagonalization of a 5 by 5 matrix, corresponding to the central and 4 neighbor sites involved in the equations of motions. This can be done directly in real space, but we will proceed in momentum space in order to use the same approach for the mean field equations in the next Section. In $`k`$-space the simplification comes from the separability of the interaction potential. As a matter of fact, $`V_{\mathrm{𝐤𝐤}^{}}=V(𝐤𝐤^{})`$, defined above can be written in a separable form as follows:

$$V_{\mathrm{𝐤𝐤}^{}}=\underset{\alpha \beta }{}w_\alpha (𝐤)𝒱_{\alpha \beta }w_\beta (𝐤^{})$$
(4)
where

$$\{\begin{array}{ccc}\hfill w_0(𝐤)& =& 1\hfill \\ \hfill w_1(𝐤)& =& [\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)]/\sqrt{2}\hfill \\ \hfill w_2(𝐤)& =& [\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)]/\sqrt{2}\hfill \\ \hfill w_3(𝐤)& =& \mathrm{sin}(k_xa)\hfill \\ \hfill w_4(𝐤)& =& \mathrm{sin}(k_ya)\hfill \end{array},$$
(5)
$`𝒱_{\alpha \beta }=v_\alpha \delta _{\alpha \beta }`$, with $`v_0=U`$ and $`v_\alpha =2V`$ for $`\alpha 0`$. We have chosen a symmetrized and anti-symmetrized combination of cosines in order to exploit the symmetry of the lattice.

The Green’s function for the two-body problem on the basis $`|\mathrm{𝐪𝐐}=c_{𝐪+𝐐/2}^{}c_{𝐪+𝐐/2}^{}|0`$ is diagonal in $`𝐐`$. We consider only the case $`𝐐=0`$ Qdep . The corresponding equation of motion is thus that of a particle in an external potential:

$`(\omega 2t_𝐤)𝒢_{\mathrm{𝐤𝐩}}(\omega ){\displaystyle \underset{𝐪}{}}V_{𝐤,𝐪}𝒢_{\mathrm{𝐪𝐩}}(\omega )=\delta _{\mathrm{𝐤𝐩}}.`$ (6)
Since $`V`$ is separable the above integral equation reduces to an algebraic one for 5 parameters. The substitution of the separable form of $`V`$ gives:

$$(\omega 2t_𝐤)𝒢_{\mathrm{𝐤𝐩}}(\omega )\underset{\alpha }{}w_\alpha (𝐤)A_{\alpha 𝐩}(\omega )=\delta _{\mathrm{𝐤𝐩}},$$
(7)
where

$$A_{\alpha 𝐩}(\omega )=v_\alpha \underset{𝐪}{}w_\alpha (𝐪)𝒢_{\mathrm{𝐪𝐩}}(\omega ).$$
(8)
The Green’s function is thus:

$$𝒢_{\mathrm{𝐤𝐩}}(\omega )=𝒢_𝐤^o(\omega )\left[\delta _{\mathrm{𝐤𝐩}}+\underset{\alpha }{}w_\alpha (𝐤)A_{\alpha 𝐩}(\omega )\right],$$
(9)
where $`𝒢_𝐤^o(\omega )=(\omega 2t_𝐤)^1`$. Substituting (9) into (8) we obtain the following closure equation for $`A_{\alpha 𝐩}`$:

$$A_{\alpha 𝐩}(\omega )=v_\alpha w_\alpha (𝐩)𝒢_𝐩^o(\omega )+v_\alpha \underset{\beta }{}B_{\alpha \beta }(\omega )A_{\beta 𝐩}(\omega ),$$
(10)
where

$$B_{\alpha \beta }(\omega )=\underset{𝐪}{}w_\alpha (𝐪)w_\beta (𝐪)𝒢_𝐪^o(\omega ).$$
(11)
The set of linear equations (10) for $`A`$ can be solved by inversion of the 5x5 matrix, $`1𝒱B(\omega )`$

$$A_{\alpha 𝐩}(\omega )=\underset{\beta }{}𝒯_{\alpha \beta }(\omega )w_\beta (𝐩)𝒢_𝐩^o(\omega )$$
(12)
where we have defined $`𝒯(\omega )=[1𝒱B(\omega )]^1𝒱`$. The final expression for the Green’s function is:

$$𝒢_{\mathrm{𝐤𝐩}}(\omega )=𝒢_𝐤^o(\omega )\left[\delta _{\mathrm{𝐤𝐩}}+𝒢_𝐩^o(\omega )\underset{\alpha \beta }{}𝒯_{\alpha \beta }(\omega )w_\alpha (𝐤)w_\beta (𝐩)\right].$$
(13)
Note that $`𝒯`$ is the usual T-matrix of scattering theory:

$$T_{\mathrm{𝐤𝐩}}=\underset{\alpha \beta }{}𝒯_{\alpha \beta }w_\alpha (𝐤)w_\beta (𝐩).$$
(14)
The problem can now be solved for any value of $`U`$, but since we are interested in the $`U=\mathrm{}`$ limit it is convenient to write $`𝒯`$ as:

$$𝒯=(𝒱^1B)^1$$
(15)
where $`(𝒱^1)_{00}=1/U=0`$. In this way we eliminate $`U`$ from the outset. One can readily verify that for $`U=\mathrm{}`$ the resulting local Green function $`𝒢(𝐱_i=0)`$ actually vanishes. As it should, the hard-core repulsion prevents two fermions with opposite spin projection to sit on the same site.

Let us come back to the explicit evaluation of $`𝒯`$ in two dimensions. The $`B`$ matrix is in block form, $`𝒢_𝐪^o`$ is in fact even under both $`k_xk_x`$ and $`k_yk_y`$. $`w_2`$ is odd under $`k_xk_y`$, and $`w_3`$ and $`w_4`$ are odd under $`k_xk_x`$ and $`k_yk_y`$ respectively. These symmetries leave few non vanishing elements and split the matrix in block form, corresponding to the irreducible representations of the symmetry group of the lattice:

$$B=\left(\begin{array}{ccccc}B_{00}& B_{01}& 0& 0& 0\\ B_{01}& B_{11}& 0& 0& 0\\ 0& 0& B_{22}& 0& 0\\ 0& 0& 0& B_{33}& 0\\ 0& 0& 0& 0& B_{33}\end{array}\right).$$
(16)
The bound states are given by the solution of the equation $`det(𝒱^1B)=0`$. We obtain thus the following equations:

$`{\displaystyle \frac{1}{2V}}`$ $`=`$ $`{\displaystyle \frac{B_{01}^2}{B_{00}}}B_{11},`$ (17)
$`{\displaystyle \frac{1}{2V}}`$ $`=`$ $`B_{22}={\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{[\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)]^2}{\omega 2t_𝐤}},`$ (18)
$`{\displaystyle \frac{1}{2V}}`$ $`=`$ $`B_{33}={\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\mathrm{sin}^2(k_xa)}{\omega 2t_𝐤}}.`$ (19)
Eqs. (17), (18), and (19) refer to $`s`$-, $`d`$-, and $`p`$-wave bound states, respectively. An additional $`s`$-wave bound state has been eliminated by setting $`U=\mathrm{}`$. Note that $`U`$ does not affect the $`p`$\- and $`d`$-solutions, since the wavefunction vanishes at the origin. In contrast, the $`s`$-wave bound state is strongly affected by $`U`$. As a matter of fact, $`B_{00}(\omega )`$, $`B_{01}(\omega )`$, and $`B_{11}(\omega )`$ all diverge logarithmically for $`\omega 2D=8t`$, with $`2D`$ the single particle band width. Using Eqs. (20) and (21) below one can verify that these divergences disappears in Eq. (17), due to mutual cancellation: the $`s`$-wave bound state is present only for $`V`$ larger than a threshold $`V_t^s`$, while for a purely attractive potential binding always occurs in two-dimensions.

It is not difficult to prove the following exact relations among the different $`B`$’s by summing and subtracting $`\omega `$ and cosine terms in the numerator:

$`\omega B_{00}(\omega )+4\sqrt{2}tB_{01}(\omega )=1,`$ (20)
$`\omega B_{01}(\omega )+4\sqrt{2}tB_{11}(\omega )=0.`$ (21)
Substituting Eqs. (20) and (21) into Eq. (17) the $`s`$-wave bound state equation becomes:

$$\frac{2t}{V}=\frac{\omega }{8t}+\frac{1}{8tB_{00}(\omega _B)}s\mathrm{wave},$$
(22)
where explicitly $`B_{00}(\omega )=_𝐤(\omega 2t_𝐤)^1`$. At the threshold $`V_t^s`$ the bound state energy lies at the band bottom $`\omega =8t`$ where $`B_{00}`$ diverges: hence the threshold $`V_t^s`$ for a $`s`$-wave bound state is simply $`V=2t`$. The $`s`$-wave solution appears non analytically at the threshold: $`\omega _B(V)2D4t\mathrm{exp}\{\pi V_t^s/(VV_t^s)\}`$. For values of $`V>V_t^s`$ Eq. (22) can be solved numerically. We plot the result in Figure 1 together with the $`d`$-wave bound state energy obtained by numerical solution of Eq. (18). As it can be seen from the picture the $`d`$-wave threshold $`V_c^d7.35t`$ is higher than the $`s`$-wave one, (despite the fact that the latter is affected by the hardcore). This reflects the cost in kinetic energy (curvature of the wave functions) of changing the sign as we go around $`𝐱=0`$.

We thus conclude that for $`V>V_t^s`$ and for low density the ground state of the system is a Bose condensate of fermions pairs in a $`s`$-wave bound state. The Bose limit is safely achieved only if the bound state energy, $`|\omega _B+2D|`$, is much larger than the interparticle distance energy scale, i.e. $`ϵ_Ftn`$. Since $`\omega _B+2D`$ vanishes exponentially near $`V_c`$ the Bose region in the $`n`$-$`V`$ diagram rapidly shrinks as $`VV_c`$ (cf. Fig. 2 in the following).

This picture for dilute pairs is quite clear and it will not be modified by the introduction of a small $`t^{}`$ (see following discussion). The opposite BCS limit is on the other hand more subtle.

## IV Mean field equations for superconductivity

In order to study the onset of superconductivity for small coupling and large density, we consider the usual BCS mean field theory.

As in the two-body problem, the integral equation can be reduced to a simple algebraic equation. The $`s`$-wave solution is strongly modified by the hard-core interaction, and the order parameter has a $`𝐤`$ dependence that comes from the mixing of the $`w_0`$ and $`w_1`$ terms, as in the two body problem. This mixing is crucial in order to ensure the vanishing of the superconducting wavefunction at the origin.

The mean-field theory equations describing the superconducting state of the Hamiltonian (1) are obtained by the usual decoupling:

$$H_{MF}=\underset{𝐤\sigma }{}\xi _𝐤c_{𝐤\sigma }^{}c_{𝐤\sigma }^{}+\underset{k}{}\left[\mathrm{\Delta }_𝐤c_𝐤^{}c_𝐤^{}+cc\right].$$
(23)
The minimization procedure gives the familiar equations for the ground state:

$`n`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}\left[1{\displaystyle \frac{\xi _𝐤}{E_𝐤}}\right],`$ (24)
$`\mathrm{\Delta }_𝐤`$ $`=`$ $`{\displaystyle \underset{𝐤^{}}{}}V(𝐤𝐤^{}){\displaystyle \frac{\mathrm{\Delta }_𝐤^{}}{2E_𝐤^{}}},`$ (25)
where $`\xi _𝐤=t_𝐤\mu _o+V(0)n/2`$, $`E_𝐤=\sqrt{\xi _𝐤^2+|\mathrm{\Delta }_𝐤|^2}`$, and for consistency we retained the Hartree term. Putting together the Hartree shift and the original shift of the chemical potential we have $`\xi _𝐤=t_𝐤\mu ^{},`$ with $`\mu ^{}=\mu +\frac{1}{2}(U4V)(n1)`$. (There is no Fock exchange term because of our interaction between opposite spins only.)

We concentrate now on the gap equation (25). We exploit again the separability of the interaction to write:

$$\mathrm{\Delta }_𝐤=\underset{\alpha }{}v_\alpha w_\alpha (𝐤)\underset{𝐤^{}}{}w_\alpha (𝐤^{})\frac{\mathrm{\Delta }_𝐤^{}}{2E_𝐤^{}}=\underset{\alpha }{}A_\alpha w_\alpha (𝐤),$$
(26)
where the new coefficients $`A`$ satisfy:

$$A_\alpha =v_\alpha \underset{\beta }{}W_{\alpha \beta }[A]A_\beta ,$$
(27)
and we have defined the matrix $`W`$, analogous of the matrix $`B`$ for the two-body problem:

$$W_{\alpha \beta }[A]=\underset{𝐤}{}\frac{w_\alpha (𝐤)w_\beta (𝐤)}{2E_𝐤}.$$
(28)
The $`A`$ dependence comes through $`\mathrm{\Delta }`$ in $`E`$. In matrix form Eq. (27) reads simply $`(1+𝒱W)A=0`$. As before it is convenient to write the equation for $`A`$ in terms of $`𝒱^1`$, in order to perform easily the limit $`U\mathrm{}`$. Eq. (27) thus becomes:

$$(𝒱^1+W[A])A=0.$$
(29)
In order to find $`\mathrm{\Delta }_𝐤`$ we need to solve the nonlinear equation (29), analogous of the linear equation (10) for the two-body problem. Let $`\lambda _i`$ and $`v_\alpha ^i`$ be the eigenvalues and eigenvectors of $`𝒱^1+W`$. We project $`A_\alpha =_ic_iv_\alpha ^i`$ on the $`v_\alpha ^i`$, so that

$$\underset{i}{}\lambda _i[c]c_iv^i=0$$
(30)
$`c_i`$ is non vanishing only if $`\lambda _i=0`$. The corresponding coefficient $`c=c_i`$ is fixed by the equation $`\lambda _i[c]=0`$. The form of the solution for $`A`$ is simply $`A_\alpha =cv_\alpha ^i`$, and for $`\mathrm{\Delta }_𝐤=c_\alpha v_\alpha ^iw_\alpha (𝐤)`$.

Due to the symmetry of the crystal, the matrix $`W`$ is already in the form (30) for $`\alpha =`$ 2, 3, and 4. The conditions $`\lambda _i=0`$ in these cases give the two equations

$`{\displaystyle \frac{1}{2V}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{[\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)]^2}{2E_𝐤}},`$ (31)
$`{\displaystyle \frac{1}{2V}}`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\mathrm{sin}^2(k_xa)}{2E_𝐤}},`$ (32)
that correspond to the $`d`$\- and $`p`$-wave solutions. In the “$`s`$” subspace $`\{w_0,w_1\}`$ we are instead left with a 2 by 2 matrix:

$$𝒱^1+W[A]=\left(\begin{array}{cc}W_{00}& W_{01}\\ W_{01}& 1/2V+W_{11}\end{array}\right).$$
(33)
The condition of vanishing of the determinant gives:

$$\frac{1}{2V}=\frac{W_{01}^2}{W_{00}}W_{11},$$
(34)
When a solution is found the form of the gap function will be given by the corresponding eigenvector of (33): $`v=(W_{01},W_{00})`$, thus

$$\mathrm{\Delta }_𝐤=\mathrm{\Delta }\left[1r\frac{\mathrm{cos}k_xa+\mathrm{cos}k_ya}{\sqrt{2}}\right]$$
(35)
where

$$r=W_{00}/W_{10}.$$
(36)
The mixing of the two $`w`$’s is crucial in order to take into account correctly the hard-core repulsion. This was already clear from the two-body problem, where the exact solution gives a wave function that is a superposition of $`w_0`$ and $`w_1`$. At the MF level we could, in principle, search for solutions with only one of these two components different from zero: such a procedure is wrong. The simplest way to see that is to evaluate the anomalous correlation function $`(𝐱𝐲)=c_{}(𝐱)c_{}(𝐲)`$. This function is the equivalent of the wavefunction for the two-body problem: it gives the probability of finding two fermions on the same site, thus for infinite $`U`$ it is bound to be zero at the origin. We can verify that our solution correctly gives $`(0)=0`$:

$`(0)={\displaystyle \underset{𝐤}{}}_𝐤={\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\mathrm{\Delta }_𝐤}{2E_𝐤}}`$ (37)
$`\left[W_{01}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{1}{E_𝐤}}W_{00}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{w_1(𝐤)}{E_𝐤}}\right]=0`$
It is clear that if, out of the two order parameters, only one is present, $`(0)`$ would be always non-vanishing. For $`p`$\- and $`d`$-wave, the symmetry alone is sufficient to guarantee the vanishing of the wavefunction.

The condition (37) has important consequences on the superconducting order. It is clear that in order to fulfill (37) $`\mathrm{\Delta }_𝐤`$ must change sign in the Brillouin zone. In particular, since $`w_1(𝐤)=t_𝐤/(2\sqrt{2}t)`$, $`\mathrm{\Delta }_𝐤`$ is constant on surfaces of fixed energy:

$$\mathrm{\Delta }_𝐤=\mathrm{\Delta }(\xi _𝐤)=\mathrm{\Delta }[1+r(\xi _𝐤+\mu ^{})/(2\sqrt{2}t)].$$
(38)
Thus the line of nodes corresponds to a given energy $`\xi _𝐤`$ instead of given direction, as it would be in $`p`$\- or $`d`$-wave superconductivity. The energy for which $`\mathrm{\Delta }`$ vanishes is $`\xi _N=2\sqrt{2}t/r\mu ^{}`$. Other consequences of (37) will be discussed in more details in Section IV.4, after the solution of the self-consistent equations.

Before proceeding, we find a simpler expression for the relevant $`W`$’s. Adding and subtracting $`\mu ^{}`$ in the numerator of $`W_{10}`$ we have:

$$W_{10}=\frac{\mu ^{}}{2\sqrt{2}t}W_{00}\frac{1}{2\sqrt{2}t}\underset{𝐤}{}\frac{\xi _𝐤}{2E_𝐤}.$$
(39)
We use Eq. (24) fixing $`n`$ to eliminate the second term. We obtain:

$$2\mu ^{}W_{00}+4\sqrt{2}tW_{10}=n1,$$
(40)
that is a generalization of Eq. (20) for the two-body problem. Similarly for $`W_{11}`$ we have:

$$W_{11}=\frac{\mu ^2}{8t^2}W_{00}+\frac{\mu ^{}(1n)}{8t^2}+b$$
(41)
where

$$b=\frac{1}{16t^2}\underset{𝐤}{}\frac{\xi _𝐤^2}{E_𝐤}=\frac{\mu ^{}}{16t^2}+\frac{1}{16t^2}\underset{𝐤}{}\frac{\xi _𝐤}{E_𝐤}(\xi _𝐤E_𝐤)$$
(42)
In conclusion we obtain the following set of equations:

$`{\displaystyle \frac{1}{V}}`$ $`=`$ $`2b(\mu ^{},\mathrm{\Delta },r){\displaystyle \frac{(n1)^2}{16t^2W_{00}(\mu ^{},\mathrm{\Delta },r)}}`$ (43)
$`{\displaystyle \frac{1}{r}}`$ $`=`$ $`{\displaystyle \frac{\mu ^{}}{2\sqrt{2}t}}+{\displaystyle \frac{n1}{4\sqrt{2}tW_{00}(\mu ^{},\mathrm{\Delta },r)}}`$ (44)
This form is more convenient for the following discussion since $`b`$ remains finite when $`\mathrm{\Delta }0`$.

Equations (43) and (44), together with Eq. (24), form the complete set of equations for the three unknowns $`\{\mu ^{},\mathrm{\Delta },r\}`$ (the gap is defined in Eq. (35)). This will be our starting point to discuss the physics of the $`s`$-wave superconducting state.

### IV.1 Threshold line for $`s`$-wave superconductivity

We have seen that in the two-body problem a pair is bound only if $`V`$ exceeds a threshold potential $`V_t^s`$. For small densities it is clear that for $`V>V_t^s`$ the system is superconducting and in the Bose limit, i.e. with a coherence length much smaller that the average interparticle distance. As shown rigorously by Randeria et al. Randeria2D , in the low density limit and in 2D, the presence of a bound state in the two-body problem is a necessary and sufficient condition for $`s`$-wave superconducting instability to take place. We thus expect that for $`V<V_t^s`$ and $`n0`$ the system is not superconducting. A separation line, starting at $`V_t^s`$ for $`n0`$, as function of $`n`$, separates the $`s`$-wave superconducting region and the normal behavior. The equation for this line of second-order transition can be found by setting $`\mathrm{\Delta }=0`$ in Eq. (43) $`n`$ being fixed by the chemical potential $`\mu ^{}`$ ($`4t<\mu ^{}<4tD`$). We obtain:

$$\frac{1}{V}=\frac{\mu ^{}}{8t^2}+\frac{1}{8t^2}\underset{𝐤}{}\left[|\xi _𝐤|\xi _𝐤\right].$$
(45)
When $`\mathrm{\Delta }=0`$, $`\mu ^{}+D=ϵ_F^o`$ is the Fermi energy of an ideal gas measured from the band edge. For small doping $`ϵ_F^o=2\pi tn`$, and since for small $`n`$ the integral in (45) is quadratic in $`ϵ_F`$ we can neglect it and find a simple analytic expression for the threshold curve:

$$\frac{V_t^s}{t}=2+\pi n+O(n^2).$$
(46)
The complete curve is shown in Fig. 2. The critical line increases smoothly from $`V/t=2`$ for $`n=0`$ to $`V/t4.93`$ at half-filling ($`n=1`$). This can be easily understood by the fact that the particles feel more and more the presence of the hard-core repulsion.

### IV.2 Near the threshold: BCS limit

Near the transition line the system behaves like a BCS superconductor for weak coupling. In particular, the superconducting gap is much smaller than $`ϵ_F^o`$, the Fermi energy of the free gas, and $`\mu ^{}`$ remains $`D+ϵ_F^o`$. In order to prove this fact we assume that these properties hold, and we evaluate Eq. (43). $`W_{00}`$ has the usual $`\mathrm{ln}\mathrm{\Delta }`$ divergent behavior and it can be evaluated in leading order by assuming a constant density of states and constant $`\mathrm{\Delta }_𝐤`$ near the Fermi surface ($`t_𝐤=\mu ^{}`$). It is convenient to write the equation in terms of $`VV_t^s(n)`$, by subtracting the critical line equation. In this way we obtain:

$$\mathrm{\Delta }(\mu ^{})2\sqrt{D^2\mu ^2}\mathrm{exp}\left\{\frac{(1n)^2VV_t(n)}{16t^2[VV_t(n)]\rho (\mu ^{})}\right\},$$
(47)
where $`\rho (\mu ^{})`$ is the density of states per spin evaluated at $`\mu ^{}`$. Near the threshold $`\mathrm{\Delta }/ϵ_F^o1`$ and thus the superconductor is of BCS type. Such a regime persists as long as $`\mathrm{\Delta }(\mu ^{})/(\mu ^{}+D)`$ is small. The dashed line in Fig. 2 indicates the values of $`n`$ and $`V`$ for which $`\mathrm{\Delta }(\mu ^{})/(\mu ^{}+D)=0.1`$. The “BCS region” is restricted between that line and the critical line $`V_t^s`$ – a narrow region compared to the remainder of the phase diagram. The behavior of the chemical potential is shown in Fig. 4.

### IV.3 Bose limit and intermediate regime

As discussed above, for $`V>V_t^s`$ and small density we inevitably reach the Bose limit. This is confirmed by the mean field equations. Indeed for $`n0`$ the quantity $`\mathrm{\Delta }/|\mu ^{}+D|`$ vanishes and $`\mu ^{}<D`$. One can thus verify that the Eq. (43) reduces to Eq. (22) for the bound state energy of the two-body problem with $`\mu ^{}=\omega _B/2`$. At the same time $`\mathrm{\Delta }`$ is fixed by Eq. (24) that gives $`\mathrm{\Delta }t\sqrt{n}`$. This behavior of $`\mathrm{\Delta }`$ can be seen in Fig. 3 for small $`n`$. For intermediate couplings and densities the numerical solution crosses over between the BCS and Bose limits described above. In Fig. 5 we plot the lines of constant $`\mathrm{\Delta }(\mu ^{})/(\mu ^{}+D)`$. We have seen that this ratio characterizes both the BCS (small positive value) and the Bose limit (small negative value). The phase diagram for the crossover from BCS to Bose condensation in the purely attractive case was given in Ref. andrenacci,  for both $`s`$-wave (local attraction) and $`d`$-wave (nearest neighbour attraction). A comparison of the results shows that the effect of the repulsion is to eliminate a large portion of the $`s`$-wave BCS region, which becomes normal. The $`d`$-wave case will be discussed in Section V.

All these curves are plotted in the range of densities (0,1): it should be stressed that the have no real meaning near half filling since the mean field approach misses the insulating state. We expect that the approximation remains valid as far as $`z_kv`$, the weight of the quasiparticle peak, remains of the order of 1. We also remember that in the corresponding Kondo problem the limit $`n=1`$, corresponds to no particles at all.

### IV.4 Nodes of $`\mathrm{\Delta }(\xi )`$

We have seen above that $`\mathrm{\Delta }`$ depends only on $`\xi `$ and that it vanishes for $`\xi =\xi _N`$. As a consequence we can write $`\mathrm{\Delta }(\xi )`$ as follows:

$$\mathrm{\Delta }(\xi )=\frac{\mathrm{\Delta }r}{2\sqrt{2}t}[\xi \xi _N]$$
(48)
where, substituting Eq. (44) in the expression for $`\xi _N`$ we have:

$$\xi _N=\frac{1n}{2W_{00}(\mu ^{},\mathrm{\Delta },r)}.$$
(49)
For small value of $`\mathrm{\Delta }`$, $`\xi _Nt/\mathrm{ln}(\mathrm{\Delta }/t)`$, thus it vanishes logarithmically in $`\mathrm{\Delta }`$. This means that approaching the threshold line $`V_t^s`$ the lines of nodes approach the Fermi surface. It turns out that keeping it at a logarithmically small distance (in $`\mathrm{\Delta }`$) suffices to fulfill the hard-core condition (37).

The presence of a line of nodes near the Fermi surface modifies the usual shape of the occupation number distribution $`n_𝐤`$. Specifically, it introduces a particle-hole asymmetry for $`|\xi |\xi _N`$. For large values of the energy the $`n`$ distribution is anomalous. Usually, far from the Fermi surface, $`n_𝐤`$ is either 0 or 1, depending on the sign of $`\xi `$. In contrast our $`\mathrm{\Delta }(\xi )`$ is $`\xi `$ when $`|\xi |`$ is large with respect to both $`\mathrm{\Delta }`$ and $`\xi _N`$, leading to

$$n(\xi )\frac{1}{2}\left[1\mathrm{sgn}(\xi )\frac{1}{\sqrt{1+\mathrm{\Delta }^2r^2/8t^2}}\right].$$
(50)
This behavior can be recognized in Figure 6 where the form of $`n(\xi )`$ is shown for $`V/t=4.5`$ and for two different values of $`n`$. Inspection of that figure shows that $`n(\xi )`$ remains symmetrical around the chemical potential only for $`|\xi |\mathrm{\Delta }(\mu ^{})`$. The asymmetry for $`|\xi |>\xi _N`$ is clearly apparent in the inset.

The aforementioned asymmetry explains the anomalous behavior of the chemical potential in the superconducting state seen in Fig. 4. The chemical potential is normally shifted down with respect to the free gas value! The BCS distribution function is particle-hole symmetric, and states are “missing” below the lower band edge: $`\mu `$ must go down in order to retain the same density. In contrast here the lower spectral density near $`\xi =\xi _N`$ leads to an increase of $`\mu `$ in a narrow range where $`\mathrm{\Delta }`$ and $`\xi _N`$ are comparable, as shown in Figures 4 and 7.

## V $`d`$-wave superconductivity

We consider now the possibility of $`d`$-wave superconductivity. By symmetry, the hard-core part of the interaction is not seen by the two-body wave function, thus the scattering states should always feel an attraction at finite energy. This is enough to induce BCS superconductivity at any finite density and for any value of the coupling constant $`V`$. In contrast the Bose limit can be reached only if a bound state is present in the two body problem, $`V>V_t^d=7.35t`$. Remember that in this limit $`s`$-wave always wins, we thus restrict to $`V<V_t^d`$.

In the opposite limit, for small values of the coupling $`VV_t^s`$, the situation is again clear: we expect a standard BCS $`d`$-wave superconductors, since there is no $`s`$-wave solution there. (The presence of a small hopping at nearest neighbors is discussed in the next section and it will not change the final outcome.)

The situation is much more interesting at intermediate couplings, $`V_t^s<V<V_t^d`$. In this region there is competition between the two possible superconducting orders: the lowest in energy will be the actual ground state. We analyze qualitatively this competition as a function of the density $`n`$ for a fixed value of the coupling $`V`$. For small density ($`n0`$) the $`d`$-wave order parameter has a BCS-like form with an effective interaction $`V_{eff}=Vw_2^2_{FS}`$ where $`w_2^2_{FS}`$ is the average over the Fermi surface of the angular factor. Since $`w_2(𝐤)(k_x^2k_y^2)`$, $`V_{eff}`$, which is of order $`w_2^2`$, is reduced by a factor $`n^2`$. The gap $`\mathrm{\Delta }_𝐤`$ is exponentially small and the resulting energy gain vanishes extremely rapidly with $`n`$ since $`\delta E\rho (\mu )\mathrm{\Delta }_𝐤^2_{FS}/2`$.

The $`s`$-wave solution, instead, is in the Bose limit, thus its energy gain is simply the number of bosons multiplied by the bound state energy: $`\delta E=n\omega _B/2`$. For $`n0`$ the $`s`$-wave is always preferred to the $`d`$-wave.

The opposite holds near the threshold for $`s`$-superconductivity: there $`d`$-wave still lowers energy, and is thus preferred. We expect a first order transition in between. The question is “where”? Will the BCS $`s`$-region be completely sweeped out by the $`d`$-wave BCS state?

In order to answer this question we need an explicit calculation of the ground state energy gain for the superconducting state:

$$\delta E=2\underset{𝐤}{}t_𝐤n_𝐤2\underset{𝐤<𝐤_F}{}t_𝐤\underset{𝐤}{}\mathrm{\Delta }_𝐤_𝐤.$$
(51)
The last term in (51) is $`\mathrm{\Delta }_d^2/(2V)`$ for the $`d`$-wave symmetry, and $`\mathrm{\Delta }^2r^2/(2V)`$ for $`s`$-wave. $`\mathrm{\Delta }_d`$ is obtained by solving Eq. (18) and Eq. (24) with $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_dw_2(𝐤)`$. $`\mathrm{\Delta }`$ and $`r`$ are the quantities introduced in (35) and discussed in the previous Section. Numerical results are shown in Fig. 8. One can clearly see the crossing of the two energy values indicating the first order transition. For small density, as expected, the $`s`$-wave state is always lower in energy (left side of the diagram). Note however that the value of the energy for which the transition takes place depends strongly on the value of $`V`$ considered. We can distinguish two cases:

(i) $`V`$ slightly above $`V_t^s(n=0)`$ the $`s`$-wave superconductivity spans the whole crossover from the Bose limit (for $`n0`$) to the BCS one. The first order transition to the $`d`$-phase occurs deep in the $`s`$-BCS region. Note that this happens at tiny values of $`\delta E`$ (for instance for $`V=2.5t`$, $`\delta E_c/t<10^{10}`$). Between the $`s`$\- and the $`d`$-wave superconducting phases there exists a region that for all practical purposes is normal. Increasing the density at fixed $`V`$ one thus find that superconductivity is reentrant.

(ii) At larger values of $`V`$, the $`s`$\- and $`d`$-state energies cross for a value of the energy gain reasonably large (of the same order of the maximum value achieved along the $`s`$-wave evolution). Thus a nearly monotonic behavior of $`\delta E`$ results (cf. the line $`V=6t`$ ).

In conclusion we found that the $`d`$-wave competes with the $`s`$-wave only for large values of $`V`$. When $`V`$ is near the threshold a quasi-normal region between the $`s`$-wave and the $`d`$-wave can survive, and the BCS-Bose crossover in $`s`$-wave can take place fully.

## VI Next neighbor hopping and $`d`$-wave $`s`$-wave competition

In this Section we consider the effect of the inclusion of a nearest-neighbors hopping integral $`t^{}`$. The effect of this term on the superconducting critical temperature has been investigated in Ref. Ranninger  and in more details in Ref. BastideLacroix . Here we consider ground state properties and in agreement with Refs. Ranninger ; BastideLacroix  we find that the introduction of $`t^{}`$ strongly modifies the $`s`$-wave phase diagram. As a matter of fact, the threshold line $`V_t^s`$ is washed out by an arbitrary small $`t^{}`$. A transparent proof of that has been given in Ref. BastideLacroix . We reproduce here their arguments for zero temperature. Equation (34) for $`\mathrm{\Delta }`$ can be rewritten as follows:

$$\frac{1}{2VW_{00}}=w_1^2w_1^2=(w_1w_1)^2,$$
(52)
where we have defined:

$$g(𝐤)\underset{𝐤}{}\frac{g(𝐤)}{E_𝐤}/\underset{𝐤}{}\frac{1}{E_𝐤}.$$
(53)
When $`\mathrm{\Delta }`$ is small, the weighted average over the Brillouin zone is peaked on the Fermi surface. If $`w_1(𝐤)`$ has the same $`𝐤`$ dependence of the kinetic energy, the righthand side of (52) vanishes at leading order in $`\mathrm{\Delta }`$, and the logarithmic divergence in $`W_{00}`$ will be ineffective: superconductivity exists for $`V`$ larger than a threshold $`V_t`$. In contrast, when the spectrum $`t_𝐤`$ is such that $`w_1(𝐤)`$ varies on the Fermi surface, the right-hand side of (52) is finite for $`\mathrm{\Delta }0`$ and superconductivity is possible however small $`V`$.

This is discussed extensively in BastideLacroix . Here we want to point out two interrelated facts: (i) the relation of superconductivity with the presence of a bound state in the two-body problem. (ii) The unusual nature of the superconducting state created by a small $`t^{}`$ in the region of the phase diagram where the system would be normal but for $`t^{}`$.

Randeria et al. Randeria2D  proved, within a two-dimensional continuum model, that the existence of an $`s`$-wave bound state for the two-body problem is a necessary and sufficient condition for the existence of $`s`$-wave superconductivity at low density. In the model at hand the presence of a small $`t^{}0`$ makes superconductivity always possible while it maintains a threshold for the two-body bound state. That threshold stems from the fact that $`w_1`$ has no fluctuations near the band edge – hence no $`\mathrm{log}\omega _B`$ term in the pair propagator. This result seems to contradict the prediction of Ref. Randeria2D .

The solution of this inconsistency comes from the $`𝐤`$ dependence of $`\mathrm{\Delta }`$ at the Fermi surface. In fact, even if the order parameter is of $`s`$-wave symmetry, it actually changes sign 8 times on the Fermi surface and it has a zero average on it. This can be verified by solving the set of Eqs. (31) and (24) in power series of $`t^{}`$. One finds the following expression for $`\mathrm{\Delta }_𝐤`$:

$$\mathrm{\Delta }_𝐤=\mathrm{\Delta }\frac{t^{}}{\mu ^{}}\left[v(𝐤)v_{FS}\right]+O(t_{}^{}{}_{}{}^{2}/\mu ^{})$$
(54)
where $`t_𝐤=2t\sqrt{2}w_1(𝐤)+t^{}v(𝐤)`$. More generally, it is clear that for any perturbation $`v(𝐤)`$ to the nearest neighbor hopping, $`\mathrm{\Delta }_𝐤`$ will change sign on the Fermi surface, and since it has $`s`$-wave symmetry this must happen at least 8 times.

Due to the unusual behavior of $`\mathrm{\Delta }_𝐤`$ the results of Randeria et al. do not apply to this $`s`$-wave state which is, in practice, close to an angular momentum $`l=4`$ state (which is isotropic as regards cubic symmetry.)

Hence $`s`$-wave superconductivity survives the hard core for arbitrary small values of the attraction. But since the superconducting state is close to $`l=4`$, the resulting value of $`\mathrm{\Delta }`$ (or $`T_c`$) is extremely small, even compared to the $`d`$-wave solution. We checked explicitly this fact for our model. We considered values of $`t^{}<t/2`$, since at this value the Fermi surface becomes double sheeted. We calculated the effective $`s`$-wave coupling

$$V_{eff}^s/V=w_1^2w_1^2\left[v^2v^2\right]t^2/(8t^2)$$
(55)
as a function of the Fermi energy and of the coupling $`t^{}`$. We compared $`V_{eff}^s`$ with the effective $`d`$-wave coupling: $`V_{eff}^d/V=w_2^2`$. In all cases we found that $`V_{eff}^d>V_{eff}^s`$.

We thus conclude that the introduction of a $`t^{}<t/2`$, or of other similar perturbations, does not change qualitatively the above discussion, since $`d`$-wave superconductivity always hides this anomalous $`s`$-wave state.

## VII Conclusions

We have studied the appearance of superconductivity on a 2D lattice in presence of a hard-core repulsion and of a nearest neighbor attraction. We constructed the mean field solution and compared it to the two-body problem in vacuum. Even if the procedure is restricted to mean field, it correctly forbids the double occupancy of the same site of the superconducting wavefunction. The main results are the following: (i) the $`s`$-wave solution is suppressed for small values of the coupling at any density. The introduction of additional small hopping integrals cannot change this, since the $`d`$-wave solution is always preferred. (ii) the $`d`$-wave solution is always possible, but the actual value of the energy gain becomes extremely small at low density. Thus for small $`V`$ and $`n`$ the critical temperatures are tiny: in practice at any temperature the system is normal. (iii) For $`V`$ larger than the threshold for the $`s`$ two body-bound state the system exhibits a crossover from the Bose-Einstein condensation of fermions pairs to a BCS behavior as the density is increased. A first order transition to the $`d`$-wave state occurs. The situation is summarized in the phase diagram of Fig. 9.

When only nearest neighbor attraction is present the competition between $`d`$-wave and $`s`$-wave at $`T_c`$ has been considered recently by Wallington and Annett Wallington . They found that for any coupling the $`s`$-wave has the lowest critical temperature, but for a small region near half-filling, where the van Hove singularity of density of states weighted with the angular factors stabilizes the $`d`$-wave phase. The presence of the hard-core repulsion allows $`d`$-wave superconductivity to appear on a much larger portion of the phase diagram. The correct treatment of the repulsion is thus crucial to study the competition between $`s`$\- and $`d`$-wave.

The technical procedure developed can be applied to the Kondo-lattice problem sketched in the introduction. However, one has to keep in mind that in that problem superconductivity and Kondo hybridization should be considered on the same footing, since each of the two order parameter reacts on the other. We will not try thus to apply the above result directly to the Kondo lattice model: we leave the detailed solution of the superconductivity-Kondo hybridization in presence of the hard-core repulsion for future work PistolesiNozieres .


text:
# Packing-limited growth of irregular objects

## I Introduction

In a previous work dodds2002pa , we examined packings formed by spheres growing in $`d`$ dimensions, immediately stopping upon contact with another sphere. We termed this to be packing-limited growth (PLG). Here, we address the question of what happens when we consider irregular (i.e., non-spherical) objects of similar shape growing in $`d=2`$ dimensions. We focus in particular on collisions of growing rectangles. By combining rectangles, we are able to form and examine a range of shapes whose growth patterns vary broadly in radial anisotropy. Whereas in dodds2002pa  we found universality classes depending only on dimension $`d`$, we find here that marked non-universal behavior arises when shapes are varied.

Our problem finds several motivations. First, this is a relatively unexplored kind of packing. Packings are generally static or randomly generated mandelbrot1977 : physical mechanisms are only seldomly connected with the creation of packings brilliantov94 ; andrienko94 . Furthermore, packings are typically of monodisperse objects of the same form or taken from a small set of forms (huillet1998, ; conway1999, ). From a physical point of view, two-dimensional packings of growing objects may be of use in modeling or understanding certain biological and ecological patterns. We consider the geometric approach presented here as a preliminary step towards describing how shape alters the size distributions of populations.

In section II we describe the various objects we construct from rectangles. Although not all of these objects have obvious physical parallels, they present a range of typological cases from which direct applications and comparisons may be sought. In section III, we provide a scaling analysis that empirically appears to be exact for $`d4`$ in the case of hyperspheres dodds2002pa  (details of the calculations for this section are given in Appendix A). We report the results of our numerical investigations in section IV along with a discussion of the failings of the scaling theory. We conclude the paper in section V and outline our algorithm for packing rectangles in Appendix B.

## II Irregular objects considered

The building block shape of the objects we consider is the rectangle. We denote the aspect ratio as $`a`$, defining all rectangles such that $`a1`$. We record rectangle size $`r`$ as half the length of the long side, the dimensions then being $`2r`$ and $`2r/a`$. Examples of each object are provided in the packings of Figures 1 and 2.

For $`a=1`$, we have squares and as $`a\mathrm{}`$, rectangles effectively become line segments and the packing becomes one of fitting $`d=1`$ objects into a $`d=2`$ volume. Rectangles afford a basic example of anisotropic growth since the long side grows at a rate slower than the short side by a factor of $`1/a`$. The growing square is distinguished from the rectangle since the former expands uniformly perpendicular to its edges. However, in comparison to disks, the growth rate of edge points relative to the center of both squares and rectangles is non-uniform.

The simplest combination of two rectangles is a cross (Figures 2A and 2B) and a natural generalization is the $`2n`$-spoke object (with $`n=1`$ being a rectangle and $`n=2`$ being a cross, see Figure 2C). Each $`2n`$-spoke object is then part of a family of shapes indexed by $`a`$.

The last object we consider is an eight-pointed star as shown in Figure 2D. This object is formed by two squares overlaid at angle of $`\pi /4`$ to each other.

All objects are packed using the approach of Manna manna92a : objects are added sequentially and allowed to instantly grow so as to just reach the existing packing structure. We have observed and argued dodds2002pa  that the value of the size distribution exponent and all other related exponents is independent of the growth dynamic.

In the context of packing, spheres make for straightforward calculations since the contact point between any two colliding spheres occurs along the line through their centers. Rectangles constitute a relatively simple generalization of spheres from a numerical point of view, hence our use of them here.

## III Scaling theory

In dodds2002pa , we argued that a simple scaling assumption may be made regarding the form of $`P(r;n)`$, the distribution of sphere sizes after $`n`$ objects have been packed. We take $`P(r;n)`$ to be described by a power law for radii above a cutoff value $`r_c`$ and uniform below:

$$P(r;n)=\{\begin{array}{cc}\frac{\alpha 1}{\alpha }r_c^1\hfill & \text{for }r<r_c\hfill \\ \frac{\alpha 1}{\alpha }r_c^1\left(\frac{r}{r_c}\right)^\alpha \hfill & \text{for }rr_c.\hfill \end{array}$$
(1)
The tail of the distribution is fixed and the distribution fills in ($`r_c`$ decreases with $`n`$) with the $`(n+1)`$th object being chosen from within the pore space. The size distribution $`P(r;n)`$ is connected to the probability of inserting an object of size $`r`$ after $`n`$ objects have been deposited, $`P_{\mathrm{ins}}(r;n)`$. Since the probability of adding a sphere of vanishing radii must be proportional to the total surface area of the existing spheres $`S(n)`$, we were able to estimate $`P_{\mathrm{ins}}(r;n)`$ as dodds2002pa

$$P_{\mathrm{ins}}(r;n)=\frac{S(n)}{\mathrm{\Phi }(n)},0rr_c,$$
(2)
where $`\mathrm{\Phi }(n)`$ is the pore space volume. Using equation (1) to calculate $`S(n)`$ and $`\mathrm{\Phi }(n)`$ and requiring that $`P_{\mathrm{ins}}(r;n)`$ be normalized, we obtained the estimate $`\alpha =11/4`$ for disks.

More generally in $`d=2`$ dimensions, if we now write area as $`A=\rho r^2`$ and perimeter as $`L=\sigma r`$, we find the scaling exponent to be

$$\widehat{\alpha }=\frac{9\sigma +4\rho }{3\sigma +2\rho }.$$
(3)
Setting $`\sigma =2\pi `$ and $`\rho =\pi `$, we recover $`\alpha =11/4`$ for disks. Note that the probability of adding a non-spherical object is still proportional to $`S(n)`$ in the limit of $`r0`$, regardless of the object’s shape. For example, it is irrelevant that an added square of side length $`2r`$ may be oriented such that its center is between $`r`$ and $`\sqrt{2}r`$ away from its point of contact with the existing packing. Also note that as $`a\mathrm{}`$, the theoretical prediction tends to $`\widehat{\alpha }=3`$, in agreement with a previous determination of an upper bound for polydisperse packings aste1996 . Details of the calculations finding $`\widehat{\alpha }`$ for the objects we consider here are to be found in Appendix A.

All other exponents depend on $`\alpha `$ via simple scaling relations dodds2002pa . In particular, the pore volume $`\mathrm{\Phi }(n)`$ decays as $`n^\beta `$ with the connection between $`\beta `$ and $`\alpha `$ being

$$\alpha =1+\frac{2}{1+\beta }.$$
(4)
We use this equation to calculate theoretical estimates of $`\beta `$, denoting them by $`\widehat{\beta }`$. Since direct measurement of $`\beta `$ is a significantly more robust exercise than determining $`\alpha `$ from $`P(r)`$, we employ equation (4) in estimating $`\alpha `$ in the following section.

## IV Numerical results

For each shape, we generate statistics for single packings with $`10^5`$ objects. Some example distributions taken from rectangle packings for $`a=1`$, 10, and 100 are shown in Figure 3. The distribution $`N(r)`$ is the unnormalized frequency distribution corresponding to $`P(r)`$. In the case of plain rectangles, we see that as $`a`$ increases, both the observed $`\alpha `$ and the theoretical estimate $`\widehat{\alpha }`$ increase while $`\beta `$ and $`\widehat{\beta }`$ accordingly decrease. Furthermore, for all the shapes considered we find that $`\alpha <\widehat{\alpha }`$ for all finite $`a`$.

We note that measuring power law exponents is not a trivial procedure and here we take some care to ensure the validity of our results. While simple regression is the basic tool of analysis, the presence of, for example, crossovers and finite size cutoffs can substantially degrade the level of precision attainable. The method of measuring exponents we use here is based on examining a smoothed version of the derivative of the distribution as viewed in double logarithmic space (similar approaches are to be found in voigt97 ; dodds2001pa ). For example, for $`\mathrm{\Phi }(n)`$, we perform regression analysis on $`\mathrm{log}_{10}\mathrm{\Phi }(n)`$ versus $`\mathrm{log}_{10}n`$ over a sliding, variable width window of values of $`\mathrm{log}_{10}n`$. Writing the upper and lower limits of this window as $`\mathrm{log}_{10}n_1`$ and $`\mathrm{log}_{10}n_2`$, we have $`w=\mathrm{log}_{10}n_2/n_1`$ being the width. In general, depending on the number of orders of magnitude spanned by the data, we would preferably choose $`w`$ in the range $`0.5w3`$. Here, we fix $`w=1`$ and calculate the “local” exponent $`\beta (n_1)`$ for each window. We find in all cases that $`\beta (n_1)`$ tends towards a constant value. This indicates the scaling law is robust and further allows us to estimate $`\beta `$ along with an error based on the fluctuations observed for $`\beta (n_1)`$. Measured values of $`\alpha `$ and $`\beta `$ are recorded in Table 1.

As we have noted above, our scaling theory approach appears to be exact for $`d4`$ dimensions in the case of hyperspheres and an overestimate of the true value of $`\alpha `$ for $`d<4`$. We therefore do not expect the scaling theory to be correct for non-spherical objects in $`d=2`$ dimensions. Indeed, in all cases, we observe the theoretical estimate $`\widehat{\alpha }`$ is an overestimate of the measured $`\alpha `$. This direction of error makes sense in light of the scaling assumption made by equation (2). The actual form of $`P_{\mathrm{ins}}(r;n)`$ is not precisely uniform but rather rolls over less steeply than a step function around $`r=r_c`$. Given the manipulations that leads to (3), it can be argued that $`\widehat{\alpha }\alpha `$ must hold dodds2002pa .

In Figure 4, we show how $`\alpha `$ varies as a function of both aspect ratio $`a`$ of the constituent rectangles and the object’s particular shape. The strongest influence is evidently aspect ratio with $`\alpha `$ varying from $`2.564`$ (disks, $`a=1`$) up to $`3`$ (rectangles, $`n`$-spokes, $`a=\mathrm{}`$). We observe a secondary effect due to the details of the shape. Square and the eight-pointed star packings have a value of $`\alpha `$ increased above that of disks. Cross packings have a higher value of $`\alpha `$ than do rectangles with the same aspect ratio.

## V Concluding remarks

We have extended a model describing the interaction of growing disks dodds2002pa  to the problem of growing irregular and anisotropic objects. In this model, the exponent characterizing the size distribution of objects is found to be independent of growth dynamics. However, the main result presented here is that the exponent is highly shape-dependent, adopting a continuous range of values $`2.564\alpha 3`$. Ultimately, understanding how geometry impacts on the structure of plant communities will require imposing a notion of packing-limited growth onto a reasonable set of dynamics. The results here demonstrate that in so doing we must keep in mind that shape matters.

###### Acknowledgements.

PSD acknowledges the support of the Columbia Earth Institute.

## Appendix A Shape parameters for mean field calculations

In this section, we derive the formulas for the area and perimeter of irregular objects in $`d=2`$. For general rectangles, we have area given by

$$A_{\text{rect}}=4r^2/a,$$
(5)
and perimeter as

$$L_{\text{rect}}=4(1+1/a)r.$$
(6)
Using equation (3), we therefore have the mean field estimate of $`\alpha `$ for rectangles as

$$\widehat{\alpha }_{\text{rect}}=3\frac{2}{5+3a},$$
(7)
where writing the result in this fashion makes plain the limiting value of $`\widehat{\alpha }_{\text{rect}}=3`$ for $`a\mathrm{}`$. For the $`2n`$-spoke objects described above, area grows as

$$A_{n\text{-spoke}}=2n/a^2(2a\mathrm{cot}\pi /2n)r^2,$$
(8)
and perimeter is given by

$$L_{n\text{-spoke}}=4n/a(a+1\mathrm{cot}\pi /2n)r,$$
(9)
leading to the estimate

$$\widehat{\alpha }_{n\text{-spoke}}=\frac{13+9a(9+2/a)\mathrm{cot}\pi /2n}{5+3a(3+1/a)\mathrm{cot}\pi /2n}.$$
(10)
Note that for $`n=1`$, we recover the result for rectangles in equation (7). Finally, for the eight-pointed star we have

$$A_{\text{8-point}}=2\left[2+(1\mathrm{tan}\pi /8)^2\right]r^2,$$
(11)
and

$$L_{\text{8-point}}=16(1\mathrm{tan}\pi /8)r,$$
(12)
which yields

$$\widehat{\alpha }_{\text{8-point}}=\frac{4240\mathrm{tan}\pi /8+2\mathrm{tan}^2\pi /8}{1514\mathrm{tan}\pi /8+\mathrm{tan}^2\pi /8}=11/8,$$
(13)
the same as for disks and squares.

## Appendix B Rectangle collisions

We describe an arbitrary rectangle in the $`x`$-$`y`$ plane as follows. The sides of the rectangle are $`2r`$ and $`2r/a`$ where $`r>0`$ is half the length of the long side of the rectangle and $`a1`$ is take the aspect ratio. The rectangle is centered at $`(x_0,y_0)`$ and rotated at an angle $`\theta `$. We take $`\theta `$ as the angle between the direction of the positive $`x`$-axis and the short axis of the rectangle (i.e., parallel to the side with length $`2r/a`$) so that $`0\theta <\pi `$.

A rectangle described by $`(r,a)`$ and $`\theta =0`$ and $`(x_0,y_0)=(0,0)`$ satisfies the equation

$$\mathrm{max}\left|\left[\begin{array}{c}ax\\ y\end{array}\right]\right|=\mathrm{max}\left|\left[\begin{array}{cc}a& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]\right|=r.$$
(14)
This is a remapping of a rectangle into a square of side length $`2r`$. For an arbitrary rectangle, we can map it onto a basic square by recentering it at the origin, removing the rotation, and undoing the dilation by the aspect ratio $`a`$:

$$\mathrm{max}\left|\left[\begin{array}{cc}a& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right]\left[\begin{array}{c}xx_0\\ yy_0\end{array}\right]\right|=r.$$
(15)
Therefore, to determine whether or not an arbitrary point $`(x,y)`$ lies on or within a given rectangle, we need to check whether or not

$$\mathrm{max}\left|\left[\begin{array}{cc}a& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right]\left[\begin{array}{c}xx_0\\ yy_0\end{array}\right]\right|r.$$
(16)
In packing growing rectangles, the above is used to check if a newly seeded rectangle has been placed in pore space and not within an existing rectangle.

For a rectangle that passes this test, the next calculation is to determine how large it may grow preserving its aspect ratio and orientation so that it just reaches an existing rectangle.

To do this, we consider one arbitrary rectangle that does not cover the origin and find the size of largest rectangle centered at the origin with the short side along the $`x`$-axis such that the rectangles just touch. We will then generalize to any configuration by appropriate rotations.

First, the “growing” rectangle could hit the already existing one at any of the latter’s corner points. The four corner points are given by

$`\left[\begin{array}{c}x_0\\ y_0\end{array}\right]+r\left[\begin{array}{c}\pm 1/a\mathrm{cos}\theta \mathrm{sin}\theta \\ \pm 1/a\mathrm{sin}\theta +\mathrm{cos}\theta \end{array}\right],`$ (21)
$`\left[\begin{array}{c}x_0\\ y_0\end{array}\right]+r\left[\begin{array}{c}a\mathrm{cos}\theta +\mathrm{sin}\theta \\ a\mathrm{sin}\theta \mathrm{cos}\theta \end{array}\right].`$ (26)
The other possible collisions are between the corners of the added rectangle and the sides of the existing rectangle. There may be 0, 1, 2, 3 or 4 such interceptions (1 or 3 if the rectangles touch at corners). To calculate these points, we parametrize each side of the existing rectangle and find the interceptions with $`y=\pm ax`$. For an arbitrary line segment described by

$$\begin{array}{c}x=b_1+b_2t\text{and}y=c_1+c_2t,\end{array}$$
(27)
with $`1t1`$, the intersection with $`y=\pm ax`$ occurs when

$$\begin{array}{c}t=\frac{c_1\pm ab_1}{c_2ab_2}.\end{array}$$
(28)
Parametrizing the sides of an arbitrary rectangle gives the following:

$`\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}x_0\\ y_0\end{array}\right]+r\left[\begin{array}{c}\mathrm{sin}\theta \\ \mathrm{cos}\theta \end{array}\right]+{\displaystyle \frac{rt}{a}}\left[\begin{array}{c}\mathrm{cos}\theta \\ \mathrm{sin}\theta \end{array}\right],`$ (37)
$`\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}x_0\\ y_0\end{array}\right]+r\left[\begin{array}{c}\mathrm{sin}\theta \\ \mathrm{cos}\theta \end{array}\right]+{\displaystyle \frac{rt}{a}}\left[\begin{array}{c}\mathrm{cos}\theta \\ \mathrm{sin}\theta \end{array}\right],`$ (46)
$`\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}x_0\\ y_0\end{array}\right]+{\displaystyle \frac{r}{a}}\left[\begin{array}{c}\mathrm{cos}\theta \\ \mathrm{sin}\theta \end{array}\right]+rt\left[\begin{array}{c}\mathrm{sin}\theta \\ \mathrm{cos}\theta \end{array}\right],`$ (55)
$`\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}x_0\\ y_0\end{array}\right]+{\displaystyle \frac{r}{a}}\left[\begin{array}{c}\mathrm{cos}\theta \\ \mathrm{sin}\theta \end{array}\right]+rt\left[\begin{array}{c}\mathrm{sin}\theta \\ \mathrm{cos}\theta \end{array}\right].`$ (64)
The first two equations describe the short sides and the latter two describe the long ones.

Using equation (28) in equations (37) to (64), we respectively have eight possible solutions (2 for each side):

$`t={\displaystyle \frac{(y_0+r\mathrm{cos}\theta )\pm a(x_0r\mathrm{sin}\theta )}{r/a\mathrm{sin}\theta r\mathrm{cos}\theta }},`$ (65)
$`t={\displaystyle \frac{(y_0r\mathrm{cos}\theta )\pm a(x_0+r\mathrm{sin}\theta )}{r/a\mathrm{sin}\theta r\mathrm{cos}\theta }},`$ (66)
$`t={\displaystyle \frac{1/a(y_0+r/a\mathrm{sin}\theta )\pm (x_0+r/a\mathrm{cos}\theta )}{r/a\mathrm{cos}\theta r\mathrm{sin}\theta }},`$ (67)
$`t={\displaystyle \frac{1/a(y_0r/a\mathrm{sin}\theta )\pm (x_0r/a\mathrm{cos}\theta )}{r/a\mathrm{cos}\theta r\mathrm{sin}\theta }}.`$ (68)
Each of these has to be tested to see if $`1t1`$. If so, then upon substituting the values of $`t`$ determined in in equations (65) through (68) into equations (37) through (64), we have

$`x={\displaystyle \frac{(x_0\mathrm{sin}\theta y_0\mathrm{cos}\theta r)}{\mathrm{sin}\theta a\mathrm{cos}\theta }},`$ (69)
$`x={\displaystyle \frac{(x_0\mathrm{sin}\theta y_0\mathrm{cos}\theta +r)}{\mathrm{sin}\theta a\mathrm{cos}\theta }},`$ (70)
$`x={\displaystyle \frac{(x_0\mathrm{cos}\theta +y_0\mathrm{sin}\theta +r/a)}{\mathrm{cos}\theta \pm a\mathrm{sin}\theta }},`$ (71)
$`x={\displaystyle \frac{(x_0\mathrm{cos}\theta +y_0\mathrm{sin}\theta r/a)}{\mathrm{cos}\theta \pm a\mathrm{sin}\theta }},`$ (72)
where $`y=\pm ax`$ in all cases.

We are now able to write down how to determine the size of the largest rectangle that may fit in given it is centered at the origin with an angle $`\varphi `$ and there is one other rectangle already present at $`(x_0,y_0)`$ oriented at an angle $`\theta `$.

1. Determine whether or not the new rectangle lies within the old one using equation (16).
2. If it is not enclosed, rotate the coordinate system by an angle $`\varphi `$ so the central rectangle sits square with the axes. All that needs to be done is to move the existing rectangle from $`(x_0,y_0)`$ to $`(\mathrm{cos}\varphi x_0+\mathrm{sin}\varphi y_0,\mathrm{sin}\varphi x_0+\mathrm{cos}\varphi y_0)`$, and change its angle of rotation to $`\theta ^{}=\theta \varphi (mod\pi )`$.
3. Find the position of the four corner points using equation (21) with $`\theta ^{}`$ and $`(x_0^{},y_0^{})`$.
4. Determine which, if any, intersection points with the diagonals of the new rectangle exist by calculating $`t`$ in equations (65) through (68).
5. Each value of $`t`$ that satisfies $`1t1`$, is a valid intersection point. The positions of all such points are then taken from equations (69) through (72) along with $`y=\pm ax`$.
6. For each valid $`(x,y)`$ pair, calculate $`\mathrm{max}(|ax|,|y|)`$ which gives us the half-width of the shortest side of the rectangle centered at the origin that passes through the point $`(x,y)`$. Take the minimum over all such values to find the half-width of the shortest side of the actual largest rectangle that may be inserted without overlapping.