text:
# Phenomenological model of propagation of the elastic waves in a fluid-saturated porous solid with non-zero boundary slip velocity

## I Introduction

It has been a common practice in the fluid dynamics to assume that when a fluid flows over an interface with a solid, the fluid molecules adjacent to the solid have zero relative velocity with the respect to solid. So far, this widely used assumption, known as ”no-slip boundary condition”, has been successfully applied to the theoretical modeling of almost all macroscopic experiments. As relevantly noticed by Craig et al. craig , the success of this assumption does not reflect its accuracy, but rather insensitivity of the experiment to a partial-slip boundary condition. It is known that the boundary slip becomes important only when the length scale over which the fluid velocity changes approaches the slip length, that is the distance behind the interface at which the fluid velocity extrapolates to zero, i.e. when the fluid becomes highly confined, e.g., blood flow through capillaries or fluid flow through natural porous rock. Recently, authors of Refs. craig ; zg  presented a convincing experimental evidence of a boundary slip in a Newtonian liquid. They performed direct measurements of the hydrodynamic drainage forces, which show a clear evidence of boundary slip. Also, they found that the boundary slip is a function of the fluid viscosity and the shear rate. These results have important implications for the blood dynamics in capillaries, the permeability of porous media, and lubrication of nano-machines. For example, results of Craig et al. suggest that red blood cells squeeze through narrow capillary walls more easily and induce less shear stress on capillary walls due to the boundary slip. Also, in oil production industry, the residual oil is difficult to produce due to its naturally low mobility. Thus, the enhanced oil recovery operations are used to increase production. It has been experimentally proven and theoretically validated that there is a substantial increase in the net fluid flow through a porous medium if the latter is treated with elastic waves ore1 ; ore2 ; ore3 ; ore4 . We may conjecture that the elastic waves via the pore wall vibration cause boundary slip of the residual oil droplets, which likewise red blood cells, squeeze through pores with less resistance, effectively increasing permeability of the porous medium.

A quantitative theory of propagation of elastic waves in a fluid-saturated porous solid has been formulated in the classic paper by Biot biot . After its appearance, this theory has seen numerous modifications and generalizations. One of the major findings of Biot’s work was that there was a breakdown in Poiseuille flow above a certain characteristic frequency specific to this fluid-saturated porous material. Biot theoretically studied this breakdown by considering the flow of a viscous fluid in a tube with longitudinally oscillating walls under an oscillatory pressure gradient. Biot’s theory can be used to describe interaction of fluid-saturated porous solid with sound for a classic Newtonian fluid assuming no-slip boundary condition at the pore walls holds. However, in the light of recent experimental results of Ref. zg , revision of the classic theory is needed in order to investigate novelties bought about by the boundary slip.

Biot’s theory has been a successful tool for interpreting the experimental data for decades, but there are circumstances when it fails. Gist gist  performed ultrasonic velocity and attenuation measurements in sandstones with a variety of saturating fluids and compared the data with the predictions of the Biot’s theory. Velocity data show systematic deviations from Biot theory as a function of pore fluid viscosity. Ultrasonic attenuation is much larger than the Biot prediction and for several sandstones is nearly constant for a three decade variation in viscosity. This behavior contrasts with synthetic porous media such as sintered glass beads, where Biot theory provides accurate predictions of both velocity and attenuation. The similar claim has been also made by earlier experimental studies e.g. Ref. mochizuki . Even Biot himself acknowledged b62  that the model of a purely elastic solid matrix saturated with a viscous fluid is rather idealistic and there are rare occasions where it gives satisfactory agreement with the experiment. Thus, he generalized theory applying the correspondence principle, and incorporated the solid and fluid attenuation, and dissipation due to non-connected pores. In this work Biot introduced the viscodynamic operator, which is a sum of the viscosity and inertial terms, and it accounts for the dynamic properties of the fluid motions in pores valid for the both low and high frequencies. The theory can model other, more complex, mechanisms such as thermoelastic dissipation, which produce a relaxation spectrum, or local flow mechanisms as the ’squirt’ flow mwk86 , dnn94 . A time domain formulation based on memory variables has been also successfully introduced c93 , cq96 . In order to narrow down the source of discrepancy between Biot theory and experiment, Gist gist  tried to modify natural sandstone by curing a residual saturation of epoxy in the pore space, filling small pores and micro-cracks. This altered rock has a significantly reduced attenuation, demonstrating the dominance of small pores in controlling ultrasonic attenuation. Then he suggested that two non-Biot viscous attenuation mechanisms are needed: local flow in micro-cracks along grain contacts, and attenuation from pore-wall surface roughness. The model model based on these two mechanisms seemed to succeed in improvement of fitting the data. However, this by no means does not exclude other possibilities, and all plausible effects that could explain the discrepancy should be explored. At low frequencies, the fluid flow is of Poiseuille type and the inertia effects are obviously negligible in comparison with the effects of the viscosity. At high frequencies, however, these effects are confined to a thin boundary layer in the vicinity of the pore walls and the inertia forces are dominant. Therefore, we conjecture that in such situations, apart from other yet unknown effects, non-zero boundary slip effect maybe responsible for the deviations between the theory and experiment. This justifies our aim to formulate a model that would account for the boundary slip.

In the section II we formulate theoretical basis of our model and in section III we present our numerical results. In the section IV we conclude with the discussion of the results.

## II The model

In our model we study a Newtonian fluid flowing in a cylindrical tube, which mimics a natural pore, whose walls are oscillating longitudinally and the fluid is subject to an oscillatory pressure gradient. We give analytical solutions of the problem in the frequency domain.

The governing equation of the problem is the linearized momentum equation of an incompressive fluid

$$\rho \frac{\stackrel{}{v}}{t}=p+\mu ^2\stackrel{}{v}.$$
$`(1)`$
Here, $`\stackrel{}{v}`$, $`p`$, $`\rho `$ denote velocity, pressure and mass density of the fluid, while $`\mu `$ is the viscosity coefficient.

Now, let $`u`$ be a velocity of the wall of the tube which oscillates in time as $`e^{i\omega t}`$. The flow of fluid in a cylindrical tube with longitudinally oscillating walls can be described by a singe component of the velocity, namely, its $`z`$-component $`v_z`$ ($`z`$ axis is along the centerline of the tube). We use cylindrical coordinate system $`(r,\varphi ,z)`$ in treatment of the problem. We introduce the relative velocity $`U_1`$ as $`U_1=v_zu`$. Thus, assuming that all physical quantities vary in time as $`e^{i\omega t}`$, we arrive at the following master equation for $`U_1`$

$$^2U_1+\frac{i\omega }{\nu }U_1=\frac{X}{\nu }.$$
$`(2)`$
Here, we have introduced the following notations:

$$\rho X=(p+\rho \frac{u}{t}),$$

which is a sum of the applied pressure gradient and force exerted on the fluid from the oscillating wall of the tube and, $`\nu `$, which is $`\nu =\mu /\rho `$.

The solution of Eq.(2) can be found to be biot

$$U_1(r)=\frac{X}{i\omega }+CJ_0(\beta r),$$

where $`J_0`$ is the Bessel function and $`\beta =\sqrt{i\omega /\nu }`$.

Assuming that the slip velocity is $`U_1(a)=U_s`$ at the wall of the tube, where $`a`$ is its radius, we obtain

$$U_1(r)=\frac{X}{i\omega }\left[1(1+\overline{U}_s)\frac{J_0(\beta r)}{J_0(\beta a)}\right].$$
$`(3)`$
Here,

$$\overline{U}_sU_s\frac{i\omega }{X}=U_s\frac{\nu }{a^2X}(\beta a)^2.$$

Defining the cross-section averaged velocity as

$$\overline{U}_1=\frac{2}{a^2}_0^aU_1(r)r𝑑r,$$

we obtain

$$\overline{U}_1=\frac{Xa^2}{\nu }\frac{1}{(\beta a)^2}\left[1\frac{2(1+\overline{U}_s)J_1(\beta a)}{(\beta a)J_0(\beta a)}\right].$$
$`(4)`$
Following work of Biot biot  we calculate the stress at the wall $`\tau `$,

$$\tau =\mu \left(\frac{U_1(r)}{r}\right)_{r=a}=\frac{\mu \beta X}{i\omega }\left(1+\overline{U}_s\right)\frac{J_1(\beta a)}{J_0(\beta a)}.$$
$`(5)`$
The total friction force is $`2\pi a\tau `$. Following Biot we calculate the ratio of total friction force to the average velocity, i.e.

$$\frac{2\pi a\mu }{\overline{U}_1}=2\pi \mu (\beta a)\left(1+\overline{U}_s\right)\frac{J_1(\beta a)}{J_0(\beta a)}$$

$$\times \left[1\frac{2\left(1+\overline{U}_s\right)J_1(\beta a)}{(\beta a)J_0(\beta a)}\right]^1.$$
$`(6)`$
Simple analysis reveals that (assuming $`\overline{U}_s0`$ as $`\omega 0`$, see discussion below)

$$\underset{\omega 0}{lim}\frac{2\pi a\tau }{\overline{U}_1}=8\pi \mu ,$$

which corresponds to the limiting case of Poiseuille flow. Following Biot biot , we also introduce a function $`F(\kappa )`$ with $`\kappa `$ being frequency parameter, $`\kappa =a\sqrt{\omega /\nu }`$, in the following manner

$$\frac{2\pi a\tau }{\overline{U}_1}=8\pi \mu F(\kappa ),$$

thus,

$$F(\kappa )=\frac{1}{4}\kappa \sqrt{i}\left(1+\overline{U}_s\right)\frac{J_1(\kappa \sqrt{i})}{J_0(\kappa \sqrt{i})}$$

$$\times \left[1\left(1+\overline{U}_s\right)\frac{2J_1(\kappa \sqrt{i})}{\kappa \sqrt{i}J_0(\kappa \sqrt{i})}\right]^1.$$
$`(7)`$
Note, that $`F(\kappa )`$ measures the deviation from Poiseuille flow friction as a function of frequency parameter $`\kappa `$. The Biot’s expression for $`F(\kappa )`$ in the no boundary slip regime can be easily recovered from Eq.(7) by putting $`\overline{U}_s0`$ for all $`\kappa `$’s.

So far, we did not specify $`\overline{U}_s`$, however there are certain physical constraints it should satisfy:

(A) Authors of Ref. craig  demonstrated that the slip length (which, in fact, is proportional to the slip velocity) is a function of the approach rate and they showed why several previous careful measurement of confined liquids have not observed evidence for boundary slip. Under the low approach rates employed in previous measurements slip length is zero and no-slip boundary condition is applicable. Experiments reported in Ref. craig  were performed when half-sphere approached a plane with liquid placed in-between at different approach rates. However, we should clearly realize what term ”approach rate” means in the context of Biot’s theory: Biot investigated fluid flow in the cylindrical tube whose walls are harmonically oscillating in the longitudinal direction as $`x(t)=Ae^{i\omega t}`$, therefore if similar to Ref.craig  experiment would be done for the oscillation tube, the ”approach rate” would be the amplitude of $`\dot{x}(t)`$, i.e. $`i\omega A`$. Thus, when $`\omega 0`$, $`\overline{U}_s`$ should also tend to zero.

(B) At high frequencies the viscous effects are negligible compared to the inertial effects. Thus, fluid at high frequencies behaves as an ideal (non-viscous) pore fluid, which allows for an arbitrary slip. Therefore, when $`\omega \mathrm{}`$, $`\overline{U}_s`$ should tend to zero. When we examine Fig. 6 from Ref.craig  and Fig. 3 (bottom panel) from Ref.zg  closely we should understand that the authors measured ascending part where slip length increases with the increase of the approach (flow) rate. Since, the viscous effects should be negligible compared to the inertial effects for large approach rates, there must be also a descending part – clearly slip length (or the same as slip velocity) cannot grow infinitely as approach rate increases. Viscous effects should give in to the inertial effects and fluid should behave an ideal fluid allowing for an arbitrary slip.

Based upon the above two physical arguments we conclude that $`\overline{U}_s`$ should be a smooth function which tends to zero when $`\omega `$ (or the same as $`\kappa `$) tends to both zero and $`\mathrm{}`$. Therefore, we postulate following phenomenological expression for $`\overline{U}_s`$

$$Re[\overline{U}_s(\kappa )]=Im[\overline{U}_s(\kappa )]=\xi \frac{B\kappa ^4}{(A+\kappa ^2)^5}$$

This function has a maximum at $`\kappa _{}=\sqrt{2A/3}`$. By fixing $`B`$ at $`5^5A^3/(4\times 3^3)`$ we force it to achieve its maximum equal to unity at $`\kappa =\kappa _{}`$ (when $`\xi =1`$). We introduced $`\xi `$ (usually $`0\xi 1`$) as a sort of ”weight” of the boundary slip effect into the solution. Putting, $`\xi =0`$ this way would easily allow us to recover no-slip case, while the increase of $`\xi `$ we would be able to trace effect of non-zero slip onto the all physical quantities. We plot this function in Fig. 1 for $`A=25`$ and $`\xi =1`$.

It is worthwhile to mention that, at first glance, it seems that non-zero boundary slip, which appears at ”high approach rates”, should be attributed to the non-linear effects. However, thorough interpretation of the experimental results of Ref. craig , in the context of the oscillatory tube (see points (A) and (B) above) allows us to conclude that the non-zero boundary slip can be also incorporated into the linear Biot’s theory.

Next, we plot both $`F_r(\kappa )=Re[F(\kappa )]`$ and $`F_i(\kappa )=Im[F(\kappa )]`$ in Fig. 2 for the three cases: when there is no boundary slip ($`\xi =0`$), and $`\xi =0.05,\mathrm{\hspace{0.33em}0.1}`$. We gather from the plot that the $`\xi =0`$ is identical to Fig. 4 from Ref.biot  as it should be. However, we also notice a noticeable difference from the classic case when $`\xi `$ is non-zero in the intermediate frequencies domain. Of course, according to our definition of the phenomenological form of $`\overline{U}_s(\kappa )`$, even for non-zero $`\xi `$ (non-zero boundary slip), when $`\kappa 0`$ and $`\kappa \mathrm{}`$, asymptotically $`F(\kappa )`$ behaves as classic Biot’s solution biot , i.e.

$$\underset{\kappa 0}{lim}F_r(\kappa )=1,$$
$`(8)`$
and

$$\underset{\kappa 0}{lim}F_i(\kappa )=0,$$
$`(9)`$
$$\underset{\kappa \mathrm{}}{lim}F(\kappa )=\frac{\kappa }{4}\sqrt{i}=\frac{\kappa }{4}\left(\frac{1+i}{\sqrt{2}}\right).$$
$`(10)`$
Since in our phenomenological model we allow for the deviations in the intermediate frequency domain, it is easy to foresee that these will have an impact on all of the predictions of the Biot’s theory precisely in that frequency range. Namely, all observable quantities predicted by the Biot’s theory, such as phase velocities and attenuation coefficients of the both rotational and dilatational waves will be affected by the introduction of boundary slip into the model in the intermediate frequency range.

Biot biot  showed that the general equations which govern propagation of rotational and dilatational high-frequency waves in a fluid-saturated porous medium are the same as in the low-frequency range provided the viscosity is replaced by its effective value as a function of frequency. In practice, it means replacing the resistance coefficient $`b`$ by $`bF(\kappa )`$.

The equations describing dynamics of the rotational waves are biot

$$\frac{^2}{t^2}(\rho _{11}\stackrel{}{\omega }+\rho _{12}\stackrel{}{\mathrm{\Omega }})+bF(\kappa )\frac{}{t}(\stackrel{}{\omega }\stackrel{}{\mathrm{\Omega }})=N^2\stackrel{}{\omega },$$
$`(11)`$
$$\frac{^2}{t^2}(\rho _{12}\stackrel{}{\omega }+\rho _{22}\stackrel{}{\mathrm{\Omega }})bF(\kappa )\frac{}{t}(\stackrel{}{\omega }\stackrel{}{\mathrm{\Omega }})=0,$$
$`(12)`$
where, $`\rho _{11},\rho _{12}`$ and $`\rho _{22}`$ are mass density parameters for the solid and fluid and their inertia coupling; $`\stackrel{}{\omega }=\mathrm{curl}\stackrel{}{u}`$ and $`\stackrel{}{\mathrm{\Omega }}=\mathrm{curl}\stackrel{}{U}`$ describe rotations of solid and fluid with $`\stackrel{}{u}`$ and $`\stackrel{}{U}`$ being their displacement vectors, while the rigidity of the solid is represented by the modulus $`N`$. Substitution of a plane rotational wave of the form

$$\omega =C_1e^{i(lx+\chi t)},\mathrm{\Omega }=C_2e^{i(lx+\chi t)},$$
$`(13)`$
into Eqs.(11) and (12) allows us to obtain a characteristic equation

$$\frac{Nl^2}{\rho a^2}=E_riE_i,$$
$`(14)`$
where $`l`$ is wavenumber, $`\chi =2\pi f`$ is wave cyclic frequency, $`\rho =\rho _{11}+2\rho _{12}+\rho _{22}`$ is the mass density of the bulk material and $`a`$ is a pore radius.

The real and imaginary parts of Eq.(14) can be written as

$$E_r=\frac{(\gamma _{11}\gamma _{22}\gamma _{12}^2)(\gamma _{22}+ϵ_2)+\gamma _{22}ϵ_2+ϵ_1^2+ϵ_2^2}{(\gamma _{22}+ϵ_2)^2+ϵ_1^2},$$
$`(15)`$
and

$$E_i=\frac{ϵ_1(\gamma _{12}+\gamma _{22})^2}{(\gamma _{22}+ϵ_2)^2+ϵ_1^2},$$
$`(16)`$
where $`\gamma _{ij}=\rho _{ij}/\rho `$, $`ϵ_1=(\gamma _{12}+\gamma _{22})(f_c/f)F_r(\kappa )=(\gamma _{12}+\gamma _{22})(f_c/f)F_r(\delta \sqrt{f/f_c})`$, $`ϵ_2=(\gamma _{12}+\gamma _{22})(f_c/f)F_i(\kappa )=(\gamma _{12}+\gamma _{22})(f_c/f)F_i(\delta \sqrt{f/f_c})`$. The function $`F(\kappa )`$ was written here more conveniently as a function of frequency $`f`$ , i.e. $`F(\kappa )=F(\delta \sqrt{f/f_c})`$ biot , where $`\delta `$ is a factor dependent on pore geometry. For the hollow cylinder-like pores, $`\delta =\sqrt{8}`$ biot  and we use this value throughout the paper. $`f_c`$ is the critical frequency above which the Poiseuille flow breaks down, and it equals $`b/(2\pi \rho _2)`$. Here $`\rho _2`$ denotes the product of porosity and fluid mass density.

In order to obtain phase velocity and attenuation coefficient of the rotational waves, we put $`l=Re[l]+iIm[l]`$. Thus, the phase velocity is then $`v_r=\chi /|Re[l]|`$. Introducing a reference velocity as $`V_r=\sqrt{N/\rho }`$, we obtain the dimensionless phase velocity as

$$\frac{v_r}{V_r}=\frac{\sqrt{2}}{\left[\sqrt{E_i^2+E_r^2}+E_r\right]^{1/2}}.$$
$`(17)`$
To obtain the attenuation coefficient of the rotational waves, we introduce a reference length, $`L_r`$, defined as $`L_r=V_r/(2\pi f_c)`$. The length $`x_a`$ represents the distance over which the rotational wave amplitude is attenuated by a factor of $`1/e`$. Therefore we can construct the dimensionless attenuation coefficient as $`L_r/x_a`$,

$$\frac{L_r}{x_a}=\frac{f}{f_c}\frac{\left[\sqrt{E_i^2+E_r^2}E_r\right]^{1/2}}{\sqrt{2}}.$$
$`(18)`$
The equations describing dynamics of the dilatational waves are biot

$$^2(Pe+Qϵ)=\frac{^2}{t^2}(\rho _{11}e+\rho _{12}ϵ)+bF(\kappa )\frac{}{t}(eϵ),$$
$`(19)`$
$$^2(Qe+Rϵ)=\frac{^2}{t^2}(\rho _{12}e+\rho _{22}ϵ)bF(\kappa )\frac{}{t}(eϵ),$$
$`(20)`$
where, $`P,Q`$ and $`R`$ are the elastic coefficients, $`e=\mathrm{div}\stackrel{}{u}`$ and $`ϵ=\mathrm{div}\stackrel{}{U}`$ are the divergence of solid and fluid displacements. Again, substitution of a plane dilatational wave of the form

$$e=C_1e^{i(lx+\chi t)},ϵ=C_2e^{i(lx+\chi t)},$$
$`(21)`$
into Eqs.(19) and (20) allows us to obtain a characteristic equation

$$(zz_1)(zz_2)+iM(z1)=0,$$
$`(22)`$
where $`z=l^2V_c^2/\chi ^2`$, $`V_c^2=(P+R+2Q)/\rho `$ represents the velocity of a dilatational wave when the relative motion between fluid and solid is absent, $`z_{1,2}=V_c^2/V_{1,2}^2`$ with $`V_{1,2}`$ being the velocities of the purely elastic waves with subscripts 1,2 referring to the two roots of Eq.(22), and finally $`M=(ϵ_1+iϵ_2)/(\sigma _{11}\sigma _{22}\sigma _{12}^2)`$ with $`\sigma _{11}=P/(P+R+2Q)`$, $`\sigma _{22}=R/(P+R+2Q)`$ and $`\sigma _{12}=Q/(P+R+2Q)`$.

Eq.(22) has two complex roots $`z_I`$ and $`z_{II}`$. Phase velocities of the two kinds of dilatational waves can be defined as

$$\frac{v_I}{V_c}=\frac{1}{Re[\sqrt{z_I}]},\frac{v_{II}}{V_c}=\frac{1}{Re[\sqrt{z_{II}}]},$$
$`(23)`$
while the corresponding attenuation coefficients can be also introduced as

$$\frac{L_c}{x_I}=Im[\sqrt{z_I}]\frac{f}{f_c},\frac{L_c}{x_{II}}=Im[\sqrt{z_{II}}]\frac{f}{f_c}.$$
$`(24)`$
## III Numerical results

In order to investigate the novelties brought about into classical Biot’s theory of propagation of elastic waves in porous medium biot  by the inclusion of boundary slip, we have studied the full parameter space of the problem.

In all forthcoming results, we calculate phase velocities and attenuation coefficients for the case 1 from Table I taken from Ref. biot , which is $`\sigma _{11}=0.610`$, $`\sigma _{22}=0.305`$, $`\sigma _{12}=0.043`$, $`\gamma _{11}=0.500`$, $`\gamma _{22}=0.500`$, $`\gamma _{12}=0`$, $`z_1=0.812`$, and $`z_2=1.674`$.

We calculated normalized phase velocity of the plane rotational waves, $`v_r/V_r`$, and the attenuation coefficient $`L_r/x_a`$ using our more general expression for $`F(\kappa )`$ (which takes into account non-zero boundary slip) given by Eq.(7).

In Fig. 3 we plot phase velocity $`v_r/V_r`$ as a function of frequency for the three cases: the solid curve corresponds to $`\xi =0`$ (no boundary slip), while long-dashed and short-dashed curves correspond to $`\xi =1.0`$ and $`\xi =1.5`$ respectively (we have used these large values in order to emphasize the effect of the boundary slip). Note that the $`\xi =0`$ case perfectly reproduces the curve 1 in Fig. 5 from Ref.biot . For $`\xi =1.0`$ and $`1.5`$ we notice a deviation from the classic behavior in the form of an decrease of phase velocity in the domain of intermediate frequencies.

Fig. 4 shows the attenuation coefficient $`L_r/x_a`$ of the rotational wave as a function of frequency for the three values of $`\xi `$: the solid curve corresponds to $`\xi =0`$ (no boundary slip), while long-dashed and short-dashed curves represent the cases $`\xi =1.0`$ and $`\xi =1.5`$ respectively. Note that $`\xi =0`$ case coincides with curve 1 in Fig. 6 from Ref.biot . We observe that, in the intermediate frequency range, with increase of boundary slip our model yields substantially higher values of the attenuation coefficient than in the classic no-slip case, which is in qualitative agreement with the experimental data gist , mochizuki .

We also calculated normalized phase velocities of the plane dilatational waves, $`v_I/V_c`$ and $`v_{II}/V_c`$, and the attenuation coefficients $`L_c/x_I`$ and $`L_c/x_{II}`$ using our more general expression for $`F(\kappa )`$ given by Eq. (7).

In Figs. 5 and 6 are similar to Fig. 3 except that now we plot phase velocities $`v_I/V_c`$ and $`v_{II}/V_c`$ as a function of frequency for $`\xi =0,\mathrm{\hspace{0.17em}0.05},\mathrm{\hspace{0.17em}0.1}`$. Note that the solid curves on the both graphs reproduce curves 1 in Figs. 11 and 12 from Ref. biot . We gather from the graph that introduction of the non-zero boundary slip leads to decrease of $`v_I/V_c`$, while $`v_{II}/V_c`$ increases in the in the domain of intermediate frequencies as boundary slip, $`\xi `$, increases.

In Figs. 7 and 8 we plot the attenuation coefficients $`L_c/x_I`$ and $`L_c/x_{II}`$ in a similar manner as in Fig. 4, but now for $`\xi =0,\mathrm{\hspace{0.17em}0.05},\mathrm{\hspace{0.17em}0.1}`$. Again we observe that the $`\xi =0`$ case reproduces curves 1 in Figs. 13 and 14 from Ref.biot . We gather from Figs. 7 and 8 that introduction of the non-zero boundary slip yields higher than in the no boundary slip case values of the attenuation coefficients for both types of the dilatational waves in the domain of intermediate frequencies, which, again, is in qualitative agreement with the experimental results gist , mochizuki .

It is worthwhile to note that for no-zero $`\xi `$ the asymptotic behavior of the elastic waves in the $`f/f_c\mathrm{}`$ as well as $`f/f_c0`$ limit is identical to the classic behavior established by Biot biot . This is a consequence of the assumptions of our phenomenological model, which is based on robust physical arguments.

## IV discussion

In this paper we have studied the effect of introduction of the boundary slip in the theory of propagation of elastic waves in a fluid-saturated porous solid originally formulated by Biot biot . Biot’s theory does not account for boundary slip effect, however, the boundary slip becomes important when the length scale over which the fluid velocity changes approaches the slip length, i.e. when the fluid is highly confined, for instance, fluid flow through porous rock or blood vessel capillaries. In the light of recent convincing experimental evidence of a boundary slip in a Newtonian liquid zg , it is necessary to take into account this effect into the Biot’s theory where appropriate. We have studied the effect of introduction of boundary slip upon the function $`F(\kappa )`$ that measures the deviation from Poiseuille flow friction as a function of frequency parameter $`\kappa `$. This function crucially enters Biot’s equations which describe dynamics of fluid-saturated porous solid. Therefore, a revision of Biot’s theory was needed in order to incorporate boundary slip effect into the all measurable predictions of this theory such as phase velocities and attenuation coefficients of the both rotational and dilatational waves. We have performed such analysis, and in summary, we found that the introduction of the non-zero boundary slip into the Biot’s theory of propagation of the elastic waves in a fluid-saturated porous solid results in

* the decrease, as compared to the no-slip limiting case, of the phase velocities of both rotational waves ($`v_r/V_r`$) and dilatational wave of the first kind ($`v_I/V_c`$) in the domain of intermediate frequencies. On contrary, the phase velocity of the dilatational wave of the second kind ($`v_{II}/V_c`$) experiences an increase as compared to the no-slip limiting case in the domain of intermediate frequencies.
* in the domain of intermediate frequencies the attenuation coefficients of both the rotational ($`L_r/x_a`$) and dilatational waves ($`L_c/x_I`$ and $`L_c/x_{II}`$) are increased as compared to the no-slip limiting case as the boundary slip increases, which is in qualitative agreement with the experimental data gist , mochizuki .
* behavior of all physical quantities which describe the elastic waves in the asymptotic limits of both small and large frequencies is not affected by the introduction of the non-zero boundary slip. The deviation occurs only in the domain of intermediate frequencies, as prescribed by our phenomenological model.

The investigation of properties of elastic waves is important for a number of applications. The knowledge of phase velocities and attenuation coefficients of elastic waves is necessary, for example, to guide the oil-field exploration applications, acoustic stimulation of oil producing fields (in order to increase the amount of recovered residual oil), and the acoustic clean up of contaminated aquifers ore1 ; ore2 ; ore3 ; ore4 . Therefore, our results would be useful for various applications in oil production as well as in ecology.

From the recent experimental results of Ref.zg  we gathered that there are physical situations were the no-slip boundary condition becomes invalid. We have formulated a phenomenological model of elastic waves in the fluid-saturated porous medium based on Biot’s linear theory and certain physically justified assumptions on the variation of boundary slip velocity with frequency, $`\overline{U}_s(\kappa )`$. Since, there are no experimental measurements of $`\overline{U}_s(\kappa )`$, for a cylindrical tube, on which ”Biot-like” theory relies, there is a certain freedom of choice, which could be used to obtain a better fit of experimental data with the theory in cases where classic Biot’s theory fails to do so. If fact, our model predicts higher than the Biot’s theory values of attenuation coefficients of the both rotational and dilatational waves in the intermediate frequency domain, which is in qualitative agreement with the experimental data gist , mochizuki . Therefore, the introduction of the boundary slip yields three-fold benefits:

* Better agreement of theory with an experimental data, since, the parametric space of the model is larger (includes effects of boundary slip).
* Possibility to identify types of porous medium and physical situations where boundary slip is important.
* Constrain model parameters that are related to the boundary slip.

We would like to close this paper with the following two remarks: First, the reported slip length varies from tens of nano-meters craig  to a few microns zg , depending on a type of a liquid and differences in the experimental set up. The latter experimental work clearly demonstrates the applicability of our model to the elastic wave phenomena in a usual porous media found in Nature. Second, in this work we postulated a phenomenological expression for the boundary slip velocity based on some physically sound arguments (in fact, similar, successful, approach has been used by the authors of Ref.tn01  in the context of solar physics). However, perhaps it may be possible to derive an expression based on more general ”first principles”, which should be a subject of a future work.

Figure Captions

Fig. 1: Behavior of $`Re[\overline{U}_s(\kappa )]=Im[\overline{U}_s(\kappa )]`$ as function of frequency parameter, $`\kappa `$, for $`\xi =1`$, $`A=25`$.

Fig. 2: Behavior of $`F_r(\kappa )`$ (thick curves) $`F_i(\kappa )`$ (thin curves) as function of frequency parameter, $`\kappa `$, according to Eq.(7). Solid curves correspond to the case when $`\xi =0`$ (no boundary slip), while long-dashed and short-dashed curves correspond to $`\xi =0.05`$ and $`\xi =0.1`$ respectively.

Fig. 3: Behavior of dimensionless, normalized phase velocity of the rotational wave, $`v_r/V_r`$, as a function of frequency. Solid curve corresponds to the case when $`\xi =0`$ (no boundary slip), while long-dashed and short-dashed curves correspond to $`\xi =1.0`$ and $`\xi =1.5`$ respectively.

Fig. 4: Behavior of dimensionless, normalized attenuation coefficient of the rotational wave, $`L_r/x_a`$, as a function of frequency. Solid curve corresponds to the case when $`\xi =0`$ (no boundary slip), while long-dashed and short-dashed curves correspond to $`\xi =1.0`$ and $`\xi =1.5`$ respectively.

Fig. 5: Behavior of dimensionless, normalized phase velocity of the dilatational wave, $`v_I/V_c`$, as a function of frequency. Solid curve corresponds to the case when $`\xi =0`$ (no boundary slip), while long-dashed and short-dashed curves correspond to $`\xi =0.05`$ and $`\xi =0.1`$ respectively.

Fig. 6: The same as in Fig. 5 except for the curves are for $`v_{II}/V_c`$.

Fig. 7: Behavior of dimensionless, normalized attenuation coefficient of the dilatational wave, $`L_c/x_I`$, as a function of frequency. Solid curve corresponds to the case when $`\xi =0`$ (no boundary slip), while long-dashed and short-dashed curves correspond to $`\xi =0.05`$ and $`\xi =0.1`$ respectively.

Fig. 8: The same as in Fig. 7 except for the curves are for $`L_c/x_{II}`$.


format: text

text:
# Towards quantum well hot hole lasers

## I Introduction

Hot hole lasers OQE23-1991 ; Hovenier-MPSSW-1997apl  emit in the THz (far-infrared) with an unusually broad gain spectrum, allowing amplification and generation of laser pulses on a picosecond time scale. The THz band has important potential applications in (e.g.) medical imaging and office communications. Bulk hot-hole lasers have been realised in p-Ge, producing gains of $``$0.25cm-1 around 4THz. Of the III-V materials, both GaAs and InSb are suitable for hot hole lasers; although their performance is not as good as in Ge Kinsler-W-2000jap-pre . Investigation of these is the most useful for industrial applications because of the existing ability to grow high quality III-V structures.

We present a discussion of likely modes of operation of quantum well hot-hole lasers in GaAs, together with predictions from Monte Carlo simulations for one design. The simulations use an infinite-well $`k.p`$ bandstructure, and include optical phonons, acoustic phonons, and piezoelectric phonons MonteCarlo ; Kane-1966ss ; Ridley-QPS . As yet they do not include ionised impurity scattering, but this is a small effect and should not affect the character of the results significantly.

## II Quantum Wells

Bulk hot-hole lasers can be described using the heavy and light hole valence bands: an electric field accelerates the heavy holes in a streaming motion to high energies $`E>E_{LO}`$; from where (ideally) they scatter into light hole cyclotron orbits formed by the magnetic field; and then emit a photon and return to low energy in the heavy hole band to repeat the cycle. In a quantum well each valence band breaks up into a set of subbands: heavy hole subbands HH1, HH2, HH3, …; and similarly for light holes LH1, LH2, etc. The non parabolicity of the bulk bandstructure leads to a variety of possible quantum well bandstructures, depending on the orientation of the crystal axes in the well material. Figure 1 shows the lower subbands of two simple cases schematically. Note that the  well is not symmetric in $`x`$ and $`y`$; and that for the  well the HH dispersions are no longer even approximately parabolic, with the HH2 having a noticeable local maximum at the origin.

The  well has two heavy-hole subbands (HH1, HH2) that need to be considered. For wells over about $`50`$Å wide, the two lowest energy, and hence most heavily populated subbands will be HH1 and HH2. This means that in contrast to the bulk case we do not need a magnetic field to confine the light holes in cyclotron orbits. The LH1 is well above the HH2 for a 100Å well, and so has minimal effect on hot-hole laser operation. Along one direction ($`y`$) both HH subbands are very flattened at the base, but with a LH-like curvature for larger $`k_y`$; whereas along the other ($`x`$) HH2 sits between HH1 and LH1. This means that the heavy-hole distributions will be roughly rectangular in outline; and an increasing electric field will shift the distributions along its direction. In the $`x`$ direction the two HH subbands (nearly) cross at a point about 15meV above the bottom of HH1, allowing for fast inter-subband scattering. Holes that get further up HH1 to an energy equal to $`E_{LO}`$ above the bottom of HH2 will quickly emit an LO phonon and return to populate either HH2 or HH1. Note that HH2 is flat to larger $`k_y`$ values than HH1, making it relatively easy to get inversion where HH2 is flat but HH1 is not.

Fig. 2(a) shows the difference in the distribution functions between HH2 and HH1, for an electric field of 250V/cm along the $`x`$ direction. We see regions of inversion for $`k_x0`$ and $`k_y5.00\times 10^8`$m-1. There is a small displacement in the direction of the field, and the amount of inversion decreases with increasing field strength. If the electric field is applied along the $`y`$ direction, we see only one inversion peak, because the other is wiped out as the HH1 distribution is shifted underneath it; and a field applied along $`x=y`$ has a similar, but not so marked, effect. Figure 2(b) shows the optical gain due to the inversion obtained and shown in (a). We see optical gain occurring over the range of energies between the small splitting at the HH1-HH2 anti-crossing ($`0`$meV) and the separation of the subband minima ($`15`$meV).

The  well HH2 subband has a distinct local maximum at $`k=0`$; and also LH1 is rather close to HH2, and should not be neglected. Our code does not yet allow for these more complicated dispersions, but we can see that a lasing cycle might be as follows: A HH1 of moderate positive $`k_x`$ is accelerated by an electric field $`+F_x`$, where it might scatter into HH2 near the anti-crossing. Then it will scatter back onto the inverted part of HH2, either by acoustic phonon scattering or LO phonon emission. Here it will be accelerated backwards past $`k=0`$ to the point of inflection ($`k_{LHi}`$) where the HH2 effective mass is infinite; and subsequently emit a photon and drop down to HH1, from where it will eventually scatter to moderate positive $`k_x`$ and will repeat the cycle. Note that HH1 also has a local maxima, but it is smaller and has less effect; and its point of inflection where holes can collect is offset from the HH2 one, and so should not interfere with population inversion too much. In contrast to the  well, we expect the emission spectrum of this  well to be peaked and centered at $`E_{HH2}(k_{LHi})E_{HH1}(k_{LHi})`$, where the HH2’s will accumulate.

Quantum well hot hole lasers also benefit from modulation doping. One significant source of unhelpful scattering is that due to holes scattering off either other holes or impurities. We can halve this contribution to the scattering simply by moving impurities used to add the holes to the device to outside the active region. Also, we can vary the well widths to vary the inter-subband separations to greater or less than $`E_{LO}`$ to enhance or reduce optical phonon scattering as required; or even use the well’s crystal orientation to adjust the densities of states (DOS) of the subbands. The  well has HH1, HH2 with flat $`E(k)`$ in the $`k_y`$ direction, leading to an enhanced density of states, hence making those regions into relatively preferred destinations for scatterings. This differs markedly from the DOS effects in bulk material Kinsler-W-2000jap-pre .

## III Conclusions

We have shown the potential for quantum well hot-hole lasers by an investigation of the possible bandstructure in combination with computer simulation. Although we have not yet explored the full parameter space of possible designs, our first attempt (the 100Å  quantum well) produced simulations which showed inversion over a range of electric field strengths and directions, with the optimum being for a 250V/cm field in the in-plane $`x`$ direction. Further, the 100Å well has a bandstructure which promises a good inversion. We aim to continue to test different designs and field combinations, including the addition of magnetic fields, in order to determine the most practical designs for experimental investigation.

Acknowledgements: This work is funded by the European Commission via the program for Training and Mobility of Researchers.


text:
# The transition between hole-pairs and four-hole clusters in four-leg 𝑡-𝐽 ladders

## I Introduction

Strongly correlated electrons confined to move on ladders formed from coupled chains have been an active topic of research in recent years. These examples of lightly doped spin liquids can be very efficiently analyzed using the density matrix renormalization group (DMRG) method, which allows long ladders to be accurately simulated. Ladders with an even number of legs have spin liquid groundstates at half-filling and the spin liquid character remains upon doping due to the binding of holes to form pairs at low concentrations. The underlying physics of the ladder systems has close similarities to that of a doped resonant valence bond (RVB) phase which is widely believed to be realized in the high-$`T_c`$ superconductors. The numerically accurate DMRG simulations of the strongly correlated $`t`$-$`J`$ model give a clear picture of the behavior of hole-pairs but only in qualitative terms as the density varies.

In this paper our aim is to analyze the DMRG results and to use them to determine the effective interactions between hole-pairs which govern their behavior. Since the spin sector of even leg ladders remains gapped, there are only charge degrees of freedom at low energies and these can be represented in terms of suitable hard core boson models. These in turn can be analyzed much more easily and completely. In an earlier paper we described such a representation or mapping of the hole-pairs on a two-leg ladder to hard core bosons on a single chain. This enabled us to determine the effective repulsive interactions between the hard core bosons by fitting to the hole density distribution determined by DMRG for the $`t`$-$`J`$ model. In the present paper we extend this analysis to four-leg ladders. The extra transverse degrees of freedom allow the formation of larger clusters of two hole-pairs (i. e. four hole clusters) when the hole density exceeds a critical concentration $`\delta _c\frac{1}{8}`$. A cluster represents not a bound state of the hole-pairs but a finite energy resonance which becomes populated when the chemical potential reaches a certain threshold.

The outline of the paper is as follows. First we treat the low concentration region where the density profiles obtained from DMRG simulations show that hole-pairs simply repel each other. In this region we extend the mapping onto the model of hard core bosons that we used earlier for the two-leg ladder and obtain a parameterization of the repulsive interactions and Luttinger liquid exponent $`K_\rho `$. We compare the evolution of $`K_\rho `$ with hole density in the two and four leg ladders. In both cases $`K_\rho 1`$ as the hole concentration $`\delta 0`$ indicating predominantly superconducting correlations. To treat the formation of the four hole clusters it is necessary to introduce additional transverse degrees of freedom. This is done by mapping the hole-pairs on a four-leg ladder onto a hard core boson model on a two-leg ladder. Choosing a potential with a repulsive tail but with a finite energy resonance on a single rung allows us to reproduce the DMRG results for the four-leg $`t`$-$`J`$ model, both for the kink in the chemical potential and for the change in the density profile at the critical density. The paper ends with a concluding section.

## II Hole Pairs and Clusters

We consider the $`t`$-$`J`$ ladder Hamiltonian

$`H`$ $`=`$ $`t{\displaystyle \underset{𝐢,𝐣,\sigma }{}}𝒫(c_{𝐢\sigma }^{}c_{𝐣,\sigma }+c_{𝐣,\sigma }^{}c_{𝐢,\sigma })𝒫`$ (2)
$`+J{\displaystyle \underset{𝐢,𝐣}{}}(𝐒_𝐢𝐒_𝐣{\displaystyle \frac{1}{4}}n_𝐢n_𝐣)`$
for a four-leg ladder (4LL) configuration as depicted in Fig. 1.

Here $`\sigma =(,)`$ denotes the spin index and $`𝐢,𝐣`$ the summation over nearest-neighbor sites with $`𝐢=(x,y)`$ as site index. The operators $`c_{𝐢,\sigma }^{}`$ and $`c_{𝐢,\sigma }`$ create or destroy electrons with spin $`\sigma `$ at position $`𝐢`$ respectively. The corresponding density operator is $`n_{𝐢,\sigma }=c_{𝐢,\sigma }^{}c_{𝐢,\sigma }`$ and $`n_𝐢=n_{𝐢,}+n_{𝐢,}`$. The projection operator $`𝒫_𝐢(1n_{𝐢,}n_{𝐢,})`$ prohibits double occupancy of a site and $`𝐒_𝐢`$ denotes the spin at site $`𝐢`$. We use $`J=0.35t`$ throughout. If not noted otherwise, we use the DMRG method to compute ground state properties which implies that open boundary conditions on the ends of the ladder are used.

We start by examining the internal structure of a hole-pair (HP) on a $`t`$-$`J`$ 4LL, by studying the hole-hole correlation (HHC) function. Since the DMRG computations are performed with open boundaries, special care must be taken in the measurement of the HHC-function, which we define as

$$g_{y,y^{}}(x,x^{})=\frac{n_{x,y}^hn_{x^{},y^{}}^h}{n_{x,y}^hn_{x^{},y^{}}^h}n^h.$$
(3)
Here $`n_{x,y}^h=1n_{x,y}`$ denotes the hole density operator acting on position $`(x,y)`$ and $`n^h`$ is the average hole density. Introducing relative and center of mass coordinates, $`r=(0,1,\mathrm{})`$ and $`R_0=(1,3/2,2,\mathrm{})`$ in the long direction, we define with

$$g(r,R_0)=\underset{y,y^{}}{}g_{y,y^{}}(R_0\frac{r}{2},R_0+\frac{r}{2})$$
(4)
the HHC between different rungs. To obtain a good approximation for the HHC function $`g(r)`$ on an infinite ladder, we measure $`g(r,R_0)`$ for an ($`L\times 4`$) ladder in the middle of the system, i. e. for $`R_0L/2`$. In this way we have measured $`g(r)`$ on a ($`32\times 4`$)-ladder with open boundaries. The result is plotted in Fig. 2. It shows clearly that in the ground state the two holes are bound. The correlation function decreases exponentially as can be seen from the logarithmic plot in the inset. The correlation length is $`1.358`$ and hence larger than in the two-leg ladder (2LL) case where we have obtained $`1.184`$ . This suggests that the binding energy of a HP is reduced as the effective dimension is enhanced. With $`E(N)`$ as the ground state energy with $`N`$ holes doped into the ladder we have computed the binding energy defined as

$$E_b=2E(1)E(0)E(2)$$
(5)
for two and 4LL using systems with $`40`$ and $`24`$ rungs respectively. As expected, we have obtained a lower binding energy $`E_b=0.1062t`$ for the 4LL than for the 2LL case, $`E_b=0.1487t`$.

If the HP behave as hard core bosons (HCB) at low hole doping one should find as many maxima in the hole density profile of a $`t`$-$`J`$ 4LL with open boundaries as the number of HP. In fact, we observe this only below a critical doping of $`\delta _c\frac{1}{8}`$ in contrast to the results obtained for two-leg ladders, where this equality holds to much higher dopings. Figures 3 and 4 shows the hole density profiles for the groundstate of a $`(16\times 4)`$ ladder for $`2`$ up to $`12`$ holes. The plot with $`6`$ holes doped into the ladder corresponds to a hole density per site $`\delta <\frac{1}{8}`$ and shows three maxima as expected. For $`8`$ holes doped into the ladder ($`\delta =\frac{1}{8}`$) we could expect to find four maxima, but we observe only three, suggesting that larger clusters of four holes are forming in some way as will be discussed later.

## III Low density of hole-pairs

In this section we examine the $`t`$-$`J`$ 4LL at low hole doping $`\delta `$, i. e. $`\delta <\delta _c`$.

In order to set up an effective bosonic model we first check the symmetry of the ground state at low hole dopings using a Lanczos method to get the lowest lying energy states with exact diagonalisation for a $`6\times 4`$ $`t`$-$`J`$ ladder with periodic boundary conditions along the legs. For $`2`$ and $`4`$ holes doped into the half filled ladder the ground state has even parity under reflection along and perpendicular to the legs.

Since the spin part of the ground state wavefunction is a singlet and thus antisymmetric and the spin excitations are gapped, the charge degrees of freedom can be described by an even parity wavefunction where HP are considered as effective (hard core) bosons. At low hole dopings, $`\delta <\delta _c`$, only one HP is found on a rung. Hence, we can map HP in the $`(L\times 4)`$ $`t`$-$`J`$ ladder with periodic boundary conditions in the long direction on HCB on a single closed chain as shown in Fig. 5. This effective model is completely analogous to the model for 2LL introduced in Ref. .

Since holes on a ladder can pair with more weight on adjacent rungs, our effective model incorporates the possibility that the “center of mass” of a HP can lie on a rung or between two rungs as shown in Fig. 5. Note that for even (odd) distance $`r`$ along the legs between two holes, the center of mass lies on a rung (between two rungs). The HHC function $`g(r)`$ is connected to the probability of finding the center of mass of a pair on a rung, $`w_{\mathrm{int}}`$, or between two rungs, $`w_{\mathrm{half}}`$, by

$`w_{\mathrm{int}}`$ $`=`$ $`{\displaystyle \underset{r_{\mathrm{even}}}{}}g(r)`$ (6)
$`w_{\mathrm{half}}`$ $`=`$ $`{\displaystyle \underset{r_{\mathrm{odd}}}{}}g(r).`$ (7)
The same occupation probabilities can be obtained with the one boson Hamiltonian

$`H_\mathrm{B}`$ $`=`$ $`t^{}{\displaystyle \underset{R=\frac{1}{2},1,\mathrm{}}{\overset{L}{}}}(B_R^{}B_{R+\frac{1}{2}}+B_{R+\frac{1}{2}}^{}B_R)+ϵ{\displaystyle \underset{R=\frac{1}{2},\frac{3}{2},\mathrm{}}{\overset{L}{}}}N_R`$ (8)
for a boson $`B_R^{}`$ which moves on a closed chain with length $`L`$, i. e. $`2L`$ sites, under the action of an alternating on-site potential which is $`0`$ or $`ϵ`$ on integer or half integer sites respectively. Here $`N_R=B_R^{}B_R`$ and $`R=(1/2,1,\mathrm{},L)`$.

Figure 5 shows the mapping of HP from the ladder to effective bosons on a single chain. Note that the center of mass coordinate $`R`$ of the pair determines the position of the effective boson. Here the ratio of the probabilities $`w_{\mathrm{half}}/w_{\mathrm{int}}`$ to find the boson on a site with half integer or integer $`R`$ depends only on $`ϵ/t^{}`$ and can easily be obtained as

$$\frac{w_{\mathrm{half}}}{w_{\mathrm{int}}}=\frac{ϵ^2+(4t^{})^2ϵ\sqrt{ϵ^2+(4t^{})^2}}{ϵ^2+(4t^{})^2+ϵ\sqrt{ϵ^2+(4t^{})^2}}.$$
(9)
From this equation and Eq. (6) we obtain for $`ϵ`$

$$ϵ=4t^{}\mathrm{sinh}\left(\frac{1}{2}\mathrm{ln}\left|\frac{_{r_{\mathrm{even}}}g(r)}{_{r_{\mathrm{odd}}}g(r)}\right|\right).$$
(10)
With $`g(r)=g(r,R_0)`$ obtained as explained in Sec. II for $`R_0L/2`$ from a ($`32\times 4`$) $`t`$-$`J`$ ladder we get $`w_{\mathrm{half}}/w_{\mathrm{int}}=1.18`$ and $`ϵ=0.33t^{}`$. This result shows that the HP is mainly centered between two rungs.

Once the boson density $`N_R`$ for the model (8) has been computed, we obtain the hole rung density $`n_x^h=_{y=1}^4n_{x,y}^h`$ on rung $`x`$ by the convolution

$`n_x^h`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{_{r_{\mathrm{even}}}g(r)(N_{x\frac{r}{2}}+N_{x+\frac{r}{2}})}{_{r_{\mathrm{even}}}g(r)}}`$ (12)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \frac{_{r_{\mathrm{odd}}}g(r)(N_{x\frac{r}{2}}+N_{x+\frac{r}{2}})}{_{r_{\mathrm{odd}}}g(r)}}.`$
Next we generalize the one boson model to a finite density including the interactions between the HCB and also the effect of open boundary conditions and write the effective Hamiltonian for the HCB as

$`H_{\mathrm{eff}}`$ $`=`$ $`H_\mathrm{B}+V_{\mathrm{int}}+V_\mathrm{b}`$ (13)
The potential $`V_{\mathrm{int}}`$ gives the interaction between HCB, i. e. HP in the $`t`$-$`J`$ model. Since the simulations are for finite systems, we have to take into account the interaction of the HP with the boundaries. The potential $`V_\mathrm{b}`$ has been introduced to describe this effect.

### A Computing $`V_\mathrm{b}`$ and $`V_{\mathrm{int}}`$

Our procedure is to compute the HCB density $`N_R`$, convolute it with $`g(r)`$ according to Eq. (12) and then to compare it with the hole rung density $`n_x^h`$ of the corresponding $`t`$-$`J`$ system. In this way we obtain $`V_\mathrm{b}`$ and $`V_{\mathrm{int}}`$ by fitting the density profile of the effective model to that of the $`t`$-$`J`$ model.

However there is a problem in the density profiles of the $`t`$-$`J`$ ladders. The open boundary conditions induce density oscillations at the boundaries which cannot be neglected as in the case of the 2LL. Due to the higher number of degrees of freedom the HP in the 4LL gets more distorted, when it approaches the open ends of the system. To circumvent this problem we use smooth bounday conditions (SBC) as proposed in Ref.  for both, the $`t`$-$`J`$ and the effective model. The SBC introduce smoothly decreasing energy parameters into the Hamiltonian as the open ends of the ladder are approached. We use SBC only in the long direction which extend into the $`t`$-$`J`$ ladder until the fifth rung from the open ends. We choose the same smoothing function as used in Ref.  for the Hubbard model. Details are explained in Appendix A.

We obtain $`V_\mathrm{b}`$ by considering one HP in the $`t`$-$`J`$ model and choose an exponentially decreasing form for the potential term $`V_\mathrm{b}`$

$`V_\mathrm{b}`$ $`=`$ $`v_b{\displaystyle \underset{R}{}}N_R\left(e^{\frac{R1}{\xi _b}}+e^{\frac{LR}{\xi _b}}\right)`$ (14)
with the two parameters $`v_\mathrm{b}`$ and $`\xi _\mathrm{b}`$.

Figure 7 shows the results of the fit and the optimal parameter values are displayed in Table I. Note that the rather large value of $`v_\mathrm{b}`$ is a consequence of the smoothing function, which also decreases the energy parameters of $`V_\mathrm{b}`$, i. e. $`v_b`$, near the boundaries.

To obtain the interaction potential $`V_{\mathrm{int}}`$ we proceed in the same way as for $`V_\mathrm{b}`$ and choose a hard core form

$`V_{\mathrm{int}}`$ $`=`$ $`{\displaystyle \underset{R}{}}{\displaystyle \underset{R^{}>R}{}}v_{\mathrm{int}}(|RR^{}|)N_RN_R^{}`$ (15)
$`v_{\mathrm{int}}(\mathrm{\Delta }R)`$ $`=`$ $`\{\begin{array}{cc}\mathrm{}\hfill & \mathrm{\Delta }R<\mathrm{\Delta }R_{\mathrm{min}}\hfill \\ v_1\hfill & \mathrm{\Delta }R=\mathrm{\Delta }R_{\mathrm{min}}\hfill \\ ve^{(RR_{\mathrm{min}}\frac{1}{2})/\xi }\hfill & \mathrm{\Delta }R>\mathrm{\Delta }R_{\mathrm{min}}\hfill \end{array}.`$ (19)
We consider a ($`24\times 4`$) $`t`$-$`J`$ ladder with four holes and the corresponding effective model with two HCB and use three fit parameters for the interaction potential and the additional parameter $`\mathrm{\Delta }R_{\mathrm{min}}`$, which broadens the hard core of the bosons, to model the interaction between the pairs.

Using the fits to the density profiles as shown in Fig. 7 the parameter values quoted in Table II are obtained. Figure 9 shows the result together with that from Ref.  for the two leg ladder.

We have tested the results by comparing the density profiles for various numbers of HP and ladder lengths. We find good agreement between the density profiles obtained from the $`t`$-$`J`$ model and the effective model for hole dopings $`\delta 0.08`$, as can be seen from Fig. 8.

### B Luttinger liquid parameter $`K_\rho `$

Having determined the effective HCB model that describes the low energy properties of the HP in the $`t`$-$`J`$ model, we can express the charge density and the superconducting correlation functions in terms of Luttinger liquid parameter $`K_\rho `$ as

$`N_RN_0\overline{N}^2`$ $``$ $`\mathrm{const}\times R^2+\mathrm{const}\times \mathrm{cos}(\pi (N_b/L)R)R^{2K_\rho }`$ (20)
$`B_R^{}B_0`$ $``$ $`\mathrm{const}\times R^{1/2K_\rho },`$ (21)
with $`N_b`$ as the number of HCB on the chain and $`\overline{N}=N_b/2L`$ the mean HCB density per site. These relations show that the superconducting correlations $`B_r^{}B_0`$ are dominant if $`K_\rho >\frac{1}{2}`$. For HCB in one dimension, $`K_\rho `$ can be obtained from the relations

$`K_\rho `$ $`=`$ $`\pi v_\mathrm{c}N\left({\displaystyle \frac{^2E_0}{N_\mathrm{b}^2}}\right)^1`$ (22)
$`v_\mathrm{c}K_\rho `$ $`=`$ $`{\displaystyle \frac{\pi }{N}}{\displaystyle \frac{^2E_0(\mathrm{\Phi })}{\mathrm{\Phi }^2}}|_{\mathrm{\Phi }=0}.`$ (23)
Here $`E_0`$ denotes the ground state energy for a closed ring of length $`L`$, i. e. $`N=2L`$ sites, with $`N_\mathrm{b}`$ HCB and $`E_0(\mathrm{\Phi })`$ is the ground state energy of the system penetrated by a magnetic flux $`\mathrm{\Phi }`$ which modifies the hopping by the usual Peierls phase factor, $`t^{}t^{}\mathrm{exp}(\pm i\mathrm{\Phi }/N)`$. From these two equations the charge velocity $`v_\mathrm{c}`$ can be eliminated. We used exact diagonalization for HCB chains with lengths between $`32`$ and $`220`$ and with $`N_\mathrm{b}=2`$.

The Luttinger liquid parameter $`K_\rho `$ for the interaction potential given by the parameters in Table II approaches the universal value $`K_\rho =1`$ at half filling. The results are shown in Fig. 10 where $`K_\rho `$ is plotted as a function of the hole doping per rung in the corresponding $`t`$-$`J`$ ladder, $`\delta _r=4\delta =2N_b/N`$. For $`N_\mathrm{b}/N0`$, corresponding to a very dilute HCB gas, we have $`K_\rho =1+O(N_\mathrm{b}/N)`$ consistent with Ref. . Up to $`\delta 0.06`$ ($`\delta _r0.23`$) the superconducting correlations are dominant, since $`K_\rho >\frac{1}{2}`$.

We have tested numerically the influence of the various parameters in Eq. 15 on $`K_\rho `$ and $`v_c`$. A longer ranged interaction leads to a lower $`K_\rho `$ and a larger $`v_c`$ leaving the product $`K_\rho v_c`$ nearly unchanged. As Fig. 9 shows, the interaction between HP on 4LL is longer ranged than on 2LL. Accordingly, the Luttinger liquid parameter for the 4LL is lower in the investigated density range. As can be seen from Figure 10 the long range charge density wave correlations tend to overcome the superconducting correlations already at lower hole dopings with increasing dimensionality, i. e. when going from 2LL to 4LL.

## IV The formation of four hole clusters

In this section we examine the four-hole clusters (FHC) which form on a $`t`$-$`J`$ 4LL at hole dopings higher than a critical doping $`\delta _c\frac{1}{8}`$.

First we clarify the question whether a change in the reflection symmetry occurs for the ground state wave function with increasing doping. We consider the real space off-diagonal singlet pairing expectation value (SPE) defined as

$`s_{𝐢,𝐣}^h`$ $`=`$ $`h2|c_{𝐢,}^{}c_{𝐣,}^{}c_{𝐢,}^{}c_{𝐣,}^{}|h`$ (24)
for nearest neighbor sites $`𝐢,𝐣`$. Here $`|h`$ denotes the ground state for $`h`$ holes doped into the ladder. We have calculated the SPE on a $`12\times 4`$ $`t`$-$`J`$ ladder for $`h`$ up to 12, i. e. a hole doping $`\delta =\frac{1}{4}`$. Figure 11 shows the results for $`h=4`$. The SPE show a characteristic sign-pattern which would change with the reflection symmetry over the doping, i. e. when going from $`hh+2`$. We observe no change as $`h`$ is further increased and conclude that the ground state symmetry remains unchanged in the investigated density range. This result suggests that HP enter the same band at dopings below and above $`\delta _c`$. Figure 11 reveals that at low hole doping HP enter in states where they reside predominantly on the outer 2LL since the magnitudes of the SPE between next neighbors on these ladders are largest. But this does not imply that two holes sit on the same rung. Instead the HHC function $`g_{y,y^{}}(x,x^{})`$ shows that two holes can most probably be found in a diagonal configuration on a $`2\times 2`$ plaquette on one of the two coupled 2LL.

To set up an effective model which can describe the formation of the FHC, we again model the charge degrees of freedom of HP by effective (hard core) bosons moving in a spin liquid but introduce additional transverse degrees of freedom by considering the $`t`$-$`J`$ 4LL as two coupled 2LL. This allows us to map HP – represented by their center of mass – from the $`(L\times 4)`$ $`t`$-$`J`$ ladder on HCB on a $`((L1)\times 2)`$ ladder as depicted in Fig. 12. In this representation a HP is centered on a $`2\times 2`$ plaquette on one of the two coupled 2LL. The Hamiltonian is given by

$`H`$ $`=`$ $`t_l{\displaystyle \underset{i=1}{\overset{L2}{}}}{\displaystyle \underset{j=1}{\overset{2}{}}}(B_{i,j}^{}B_{i+1,j}+B_{i+1,j}^{}B_{i,j})`$ (27)
$`t_r{\displaystyle \underset{i=1}{\overset{L1}{}}}(B_{i,1}^{}B_{i,2}+B_{i,2}^{}B_{i,1})`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j,i^{},j^{}}{}}v(|ii^{}|,|jj^{}|)N_{i,j}N_{i^{},j^{}}`$
where $`B_{i,j}^{}`$ and $`B_{i,j}`$ create and destroy a HCB on site $`(i,j)`$ respectively and $`N_{i,j}`$ is the corresponding density operator. The matrix elements $`t_l`$ and $`t_r`$ stand for nearest neighbor hopping along the legs and rungs of the 2LL. The last term describes the interaction between HP and is determined by the parameter $`v(r_x,r_y)`$ with $`r_x`$ and $`r_y`$ as the relative separation of two HCB along and perpendicular to the legs respectively.

As hopping matrix element along the legs, $`t_l`$, we use the value obtained in Ref.  for 2LL, $`t_l=0.303t`$. The energy difference of the bonding and antibonding band is given by $`2t_r`$. We have computed this energy difference at wavevector $`k_x=0`$ with exact diagonalisation for two holes on a $`6\times 4`$ leg $`t`$-$`J`$ ladder using periodic boundary conditions in the long direction and obtained $`t_r=0.092t`$ so that $`t_l/t_r3`$.

We are able to reproduce qualitatively the formation of FHC but the formation depends crucially on the form of the interaction between the HP. We know from Ref.  that on a 2LL HP repel each other and hence choose a repulsive interaction for $`r_x>0`$ which tends to zero when $`r_x`$ increases. Further we expect this repulsion to be less for hole pairs on different legs of the effective 2LL, i. e. $`v(r_x,1)<v(r_x,0)`$. At short distance we can expect a certain energy lowering due to a lower cost in magnetic energy in the configuration with four holes on neighboring sites. Thus we set the constraint that $`v(r_x,1)`$ has a local minimum at $`r_x=0`$.

The solution of the two particle problem for the effective HCB model with this type of interaction shows that in the ground state the two HCB are unbound. But there are also resonant cluster states which can be occupied at higher energies. Hence, it can be expected that with increasing HCB density beyond a certain value of the chemical potential a resonant state is occupied. This occurs when the expense in kinetic energy due to the loss of some degrees of freedom is compensated by a smaller interaction energy. Above the critical doping a HP can cluster with another on the same rung but then it tunnels from one suitable HP to another HP thereby gaining additional kinetic energy.

A set of parameter values for $`v(r_x,r_y)`$ is given in Table III. With these we calculated the HCB density per rung, $`N_i=_{j=1}^2N_{i,j}`$, for $`3`$ up to $`6`$ HCB on a $`15\times 2`$ ladder shown in Fig. 13.

These results correspond to the hole density profiles shown in Fig. 3 and Fig. 4 for the $`16\times 4`$ $`t`$-$`J`$ ladder for $`3`$ up to $`6`$ HP. As can be seen, the formation of FHC is reproduced qualitatively. The number of maxima in the density profiles does not change, if a fourth HCB (HP in the $`t`$-$`J`$ model) is added and for $`6`$ HCB (HP) we observe three well separated maxima, each containing two HCB (HP). The density profiles for $`5`$ HCB and HP look a bit different. Here the two additional HP cluster with the others but, due to the repulsive interaction between them, they are pushed towards the open ends of the sample in x-direction. We observe the same effect in the effective model, where the two outer maxima have higher weights than the maxima in the center.

We tested the effective model also for longer systems. Figure 14 shows the rung density profiles for $`5`$ up to $`8`$ HP on a $`(24\times 4)`$ $`t`$-$`J`$ ladder. The FHC start to form when going from $`5`$ to $`6`$ HP doped into the ladder, which corresponds again to a critical doping $`\delta \frac{1}{8}`$. By filling the ladder further with HP the number of maximum can also decrease. For $`7`$ HP we observe only four maxima. The area attached to each maxima suggests that two FHC form near the open ends by occupying seven rungs each, whereas the remaining three HP share the inner ten rungs, where one pair is loosely bound with the other two. Finally, for $`8`$ HP we observe $`4`$ FHC which occupy $`6`$ rungs each. Also this “contraction” of HP can be reproduced with the effective model using the interaction parameters from Table III. Figure 15 displays the corresponding rung density profiles for $`5`$ up to $`8`$ HCB on a $`23\times 2`$ ladder. They show the same evolution with HCB number as the hole density profiles in Fig. 14 with HP number. We cannot expect the effective model to reproduce the behavior of the $`t`$-$`J`$ model in all details, but in view of the simplicity of the model, the results are satisfactory. After scanning a wide parameter range for $`v(r_x,r_y)`$ we came to the conclusion that the formation of FHC can occur when the interaction leads to resonant states.

Since we observe a qualitative change in the density profiles, we have computed the chemical potential $`\mu (h)`$ of the HP as a function of hole doping

$`\mu (h)`$ $`=`$ $`E(h+1)E(h1)`$ (28)
for both models in order to check whether it shows a singularity at the critical doping. Here, $`E(h)`$ denotes the ground state energy for a system with $`h`$ holes or $`h/2`$ HCB in the $`t`$-$`J`$ and effective model respectively.

Figure 16 displays the chemical potential for HCB in the bosonic model for a $`23\times 2`$ ladder as a function of the hole doping $`\delta `$ in the corresponding $`t`$-$`J`$ model. The hole doping is related to the HCB number $`N_b`$ by $`\delta =N_b/2L`$. The chemical potential shows a deviation from the quadratic form observed at low densities right at the position where the formation of the FHC sets in, i. e. when $`N_b`$ is increased from $`5`$ to $`6`$. The exact position of this kink can be tuned by choosing appropriate interaction parameter values $`v(r_x,r_y)`$. The difference in the chemical potential between the critical and zero doping, $`\mathrm{\Delta }\mu 0.42t`$, is to a good approximation given by the difference between the energies of the first resonant state and the ground state, $`\mathrm{\Delta }E=0.43t`$, obtained by solving the two particle problem. The critical doping can be shifted to lower values by making the interaction longer ranged.

Figure 17 shows the chemical potential for hole pairs obtained from the corresponding $`24\times 4`$ $`t`$-$`J`$ ladder. The dashed dotted line shows a fit to a second order polynomial for $`\delta <\delta _c`$. As expected, the chemical potential $`\mu (h)`$ deviates from this polynomial for $`\delta >\delta _c`$.

Finally, we note that while the effective model reproduces the evolution of the system from HP to FHC nicely, it does not reproduce the density profiles quantitatively in the FHC density regime as can be seen, for example, by comparing the upper curves in Fig. 14 and 15. Therefore, further refinement of the effective model would be required to obtain the accuracy necessary to calculate parameters such as $`K_\rho `$ from the effective model in this density regime.

## V Conclusions

In this work our previous study of the low energy properties of the hole doped $`t`$-$`J`$ ladders with two legs is extended to the case of four legs. In both cases density matrix renormalization group (DMRG) results show that the holes bind in pairs at low densities and a finite spin gap is preserved. To analyze the DMRG results in more detail we introduce a hard core boson model on a single chain to describe the low energy degrees of freedom. The effective interactions between the hard core bosons which represent hole-pairs, are determined by fitting the density profiles to those obtained by DMRG methods for the $`t`$-$`J`$ model. This effective repulsive interaction is longer ranged for the four-leg ladder and this in turn reduces the value of the Luttinger liquid parameter $`K_\rho `$. As a consequence the region of predominantly superconducting correlations is reduced in the wider ladder. Whereas the hole-pairs in the two-leg ladder simply repel each other, in the wider four-leg ladder a modification of the hole density profiles appears beyond a critical hole doping which can be simply interpreted as the formation of four hole clusters. To reproduce this behavior in the hard core boson model, it is necessary to introduce an extra transverse degree of freedom by replacing the single chain with a two-leg ladder. It is also necessary to modify the interaction potential to incorporate a four hole cluster as a finite energy resonance. This is achieved by replacing the monotonically decreasing repulsive interaction on the single chain by one with a potential minimum at short range. The longer range repulsive tail of the interaction potential causes the chemical potential to add a hole-pair to rise as the density is increased. When the chemical potential exceeds a critical value approximately equal to the resonance energy, the four hole clusters are formed leading to a kink in the chemical potential as a function of hole density.

## VI Acknowledgements

SRW acknowledges the support of the NSF under grant DMR98-70930. MT was supported by the Swiss National Science Foundation. The DMRG calculations have been performed on the SGI Cray SV1 of ETH Zurich.

## A Smooth boundary conditions

The SBC introduce smoothly decreasing energy parameters into the Hamiltonian as the edges of the lattice are approached. The result of this operation is that, instead of having a sharp and rigid boundary, the boundary extends itself into the system and its exact size is not fully determinable. In this way we talk of the bulk of the system as the region where the energy parameters are constant, and of the boundary as the region over which the parameters are smoothly turned off.

For the $`t`$-$`J`$ ladders we use SBC only in the long direction and denote the width of the boundary with $`m`$.

By applying SBC to the Hamiltonians (2) and (8) we replace the energy parameters $`t`$, $`J`$, $`t^{}`$, $`ϵ`$, $`v_\mathrm{b}`$, and $`v_{\mathrm{int}}`$ by the site dependent parameters $`t_{𝐢,𝐣}`$, $`J_{𝐢,𝐣}`$, $`t_R^{}`$, $`ϵ_R`$, $`v_{\mathrm{b},R}`$, and $`v_{\mathrm{int},R,R^{}}`$ respectively. All these parameters are scaled according to the smoothing function $`f`$ defined as

$$f(d)=\{\begin{array}{cc}0\hfill & d=0\hfill \\ \frac{1}{2}\left[1+\mathrm{tanh}\frac{d\frac{1}{2}}{d(1d)}\right]\hfill & 0<d<1\hfill \\ 1\hfill & d>=1.\hfill \end{array}$$
(A1)
and plotted in Fig.18. Here $`d`$ denotes the distance to the nearest outermost rung of the system in units of $`m`$. For the $`t`$-$`J`$ model we use

$`d_{𝐢,𝐢^{}}`$ $`=`$ $`d_{x,y,x^{},y^{}}=\mathrm{min}({\displaystyle \frac{x+x^{}}{2}}1,L{\displaystyle \frac{x+x^{}}{2}})/m`$ (A2)
with $`L`$ as the number of sites in the long direction of the $`t`$-$`J`$ ladder and scale the energy parameters in such a way that $`t_{𝐢,𝐣}/t=f(d_{𝐢,𝐣})`$ and $`J_{𝐢,𝐣}/J=f(d_{𝐢,𝐣})`$, where $`t`$ and $`J`$ are the bulk values.

For the one dimensional chain with length $`N`$ of the effective model we define

$`d_{R,R^{}}`$ $`=`$ $`\mathrm{min}({\displaystyle \frac{R+R^{}}{2}}1,N{\displaystyle \frac{R+R^{}}{2}})/m`$ (A3)
$`d_R`$ $`=`$ $`\mathrm{min}(R1,NR)/m`$ (A4)
and scale the energy parameters in the way that $`t_R^{}/t^{}=f(d_{R+\frac{1}{4}})`$, $`ϵ_R/ϵ=f(d_R)`$, $`v_{\mathrm{b},R}/v_\mathrm{b}=f(d_R)`$, and $`v_{\mathrm{int},R,R^{}}/v_{\mathrm{int}}=f(d_{R,R^{}})`$. Here again, $`t^{}`$, $`ϵ`$, $`v_\mathrm{b}`$, and $`v_{\mathrm{int}}`$ are the bulk values.

To determine an appropriate value for the width $`m`$ of the boundaries we have considered the curvature of the rung hole density profile of a $`(24\times 4)`$ $`t`$-$`J`$ ladder. For $`m<4`$ the curvature is far from being a smooth function of the rung coordinate but for $`m=4`$ it is. So we have used $`m=4`$ for the computations with SBC.


text:
# Scale-free topology of e-mail networks

(February 11, 2002)

## Abstract

We study the topology of e-mail networks with e-mail addresses as nodes and e-mails as links using data from server log files. The resulting network exhibits a scale-free link distribution and pronounced small-world behavior, as observed in other social networks. These observations imply that the spreading of e-mail viruses is greatly facilitated in real e-mail networks compared to random architectures.

Complex networks as the World Wide Web or social networks often do not have an engineered architecture but instead are self-organized by the actions of a large number of individuals. From these local interactions non-trivial global phenomena can emerge as, for example, small-world properties Watts and Strogatz (1998) or a scale-free distribution of the degree Barabási and Albert (1999). These global properties have considerable implications on the behavior of the network under error or attack Albert et al. (2000), as well as on the spreading of information or epidemics Pastor-Satorras and Vespignani (2001a). Here we report that networks composed of persons connected by exchanged e-mails show both the characteristics of small-world networks and scale-free networks.

The nodes of an e-mail network correspond to e-mail addresses which are connected by a link if an e-mail has been exchanged between them. The network studied here is constructed from log files of the e-mail server at Kiel University, recording the source and destination of every e-mail from or to a student account over a period of 112 days 111Three e-mail addresses has been excluded as artifacts, since they reached a large quantity of the students only because all students used the same server (e.g. by service e-mails from the university computer center). Therefore their large degree was caused solely by sampling log files from only one e-mail server.. The resulting network consists of $`N=`$ 59,912 nodes (including 5,165 student accounts) with a mean degree of $`<k>`$ = 2.88 and contains several separated clusters with less than 150 nodes and one giant component of 56,969 nodes (mean degree $`<k_{\text{large}}>`$ = 2.96). The degree distribution $`n(k)`$, i.e. the distribution of the number $`k`$ of a node’s next neighbors, obeys a power law

$$n(k)k^{1.81},$$
(1)
with exponential cut-off (Fig. 1).

Most of the scaling exponents reported so far for the degree distributions of computer and social networks lie in the range of -2.0 to -3.4 Albert and Barabási (2002). One exception is the social network of co-authorships in high energy physics, for which Newman found an exceptionally small scaling exponent of -1.2 Newman (2001). Similar to our work are studies of networks of phone calls made during one day. These phone-call networks show scale-free behavior of the degree distribution as well, with an exponent of -2.1 Abello et al. (1999); Aiello et al. (2000).

Let us briefly discuss how our result on e-mail networks may be influenced by the measurement process. The sampling of the network has been restricted to one distinct e-mail server. Therefore, only the degrees of accounts at this server are known exactly. Here, these internal accounts correspond to e-mail addresses of local students, whereas the external nodes are given by all other e-mail addresses. We resolve the degree distribution of internal accounts only (Fig. 2), and find that it can be approximated by a power-law $`n_{\text{int}}(k)k^{1.32}`$ as well (mean degree $`<k_{\text{int}}>=`$ 14.86). Since the degrees of external nodes typically are underestimated, this exponent is smaller than for the whole network. For the same reason, there are fewer nodes with small degree in the distribution of students’ degrees. Note that the cut-off of both distributions is about the same. Therefore, external sources addressing almost all internal nodes (e.g. advertisement or spam) do not bias the degree statistics. Thus, it can be concluded that the e-mail network exhibits scale-free behavior.

Furthermore, the e-mail network shows the properties of a “small world” Watts and Strogatz (1998), i.e. a high probability that two neighbors of one node are connected themselves (clustering) and a small average length $`\mathrm{}`$ of the shortest path between two nodes. The clustering is measured by the clustering coefficient $`C`$ of a network which is defined in the following way: The clustering coefficient $`C_\nu `$ of a node $`\nu `$ is given by the ratio of existing links $`E_\nu `$ between its $`k_\nu `$ first neighbors to the potential number of such ties $`\frac{1}{2}k_\nu (k_\nu 1)`$. By averaging $`C_\nu `$ over all nodes one arrives at the clustering coefficient $`C`$ of the network

$$C=C_\nu _\nu =\frac{2E_\nu }{k_\nu (k_\nu 1)}_\nu .$$
(2)
To verify the small-world properties, experimental results are compared to the respective values derived for random networks, as well as for networks with links assigned randomly but obeying the same degree distribution. We will call a network a “random network” if the probability $`p`$ that an arbitrary pair of nodes is connected is constant $`p=k/(N1)`$, leading to a Poissonian degree distribution. In this case, the clustering coefficient is just this probability $`C_{\text{rand}}=p`$. Additionally, we deduced an estimate $`C^{}`$ for an upper bound of the clustering coefficient of a network with identical degree distribution, but randomly assigned links. Hence, $`C^{}`$ gives an upper bound of the clustering that is expected from the degree distribution alone. Employing the generating function approach for networks with arbitrary degree distributions Newman et al. (2001) and assuming that fluctuations of the mean degree of a node’s neighborhood can be ignored, results in:

$$C^{}=\frac{1}{kN}\left(\frac{k^2}{k}1\right)^2.$$
(3)
This equation exactly applies in the case of the Poissonian degree distribution of a random network ($`C_{\text{rand}}=C^{}`$). Next we compare the experimental average shortest path length $`\mathrm{}`$ with the respective value $`\mathrm{}^{}`$ of a network with the same degree distribution and randomly assigned links. With the generating function approach we obtain Davidsen et al. (2002):

$$\mathrm{}^{}\frac{\mathrm{log}\left(\frac{N}{k}\right)}{\mathrm{log}\left(\frac{k^2k}{k}\right)}+1.$$
(4)
For the Poissonian distribution of a random network, Eqn. (4) simplifies to:

$$\mathrm{}_{\text{rand}}\frac{\mathrm{log}N}{\mathrm{log}k}.$$
(5)
Note that only $`N`$, $`k`$ and $`k^2`$ are used for the estimates $`\mathrm{}^{}`$ and $`\mathrm{}_{\text{rand}}`$ as in the derivation of Eqn. (3).

In the following, the experimental values are compared to the above estimates. The experimental clustering coefficient $`C=`$ 0.156 is about one magnitude larger than one would expect solely from the degree distribution ($`C^{}`$ = 0.0187). For a random network, the respective clustering coefficient is $`C_{\text{rand}}=4.8210^5`$.

The mean shortest path length in the giant component was determined to $`\mathrm{}`$ = 4.95 $`\pm `$ 0.03 with the Dijkstra algorithm Rosen (2000). It is larger than the mean shortest path length in a network with the same degree distribution 222To be specific, $`\mathrm{}`$ has to be calculated using the degree distribution of the giant component. Since this distribution does not differ significantly from the degree distribution of the whole network, employing the latter yields the same result within numerical accuracy used here. ($`\mathrm{}^{}`$ = 3.43) since more links are consumed for forming local clusters. It is still smaller than the path length of a random network $`\mathrm{}_{\text{rand}}`$ = 10.10 (where each pair of nodes is connected with a constant probability leading to the same mean degree) because of the highly connected ”hubs” present in a scale-free network.

To further investigate the emergence of the scale-free degree distribution, we study the e-mail network as a directed graph, where an e-mail corresponds to a directed link pointing from the sender to the receiver. Although the e-mail network has to be treated as an undirected graph in the context of virus spreading (see below), it seems reasonable that the sending and receiving of e-mails are governed by different processes. Again, the analysis is done for the distributions of all nodes and of internal nodes only, where for the latter, the in- and out-degree can be determined exactly. The distribution of the in-degree $`i`$, i.e. a node’s number of different nodes it has received e-mails from, are very similar for all nodes and for internal nodes, respectively (Fig. 3). They can both be approximated by a power law $`n(i)i^{1.49}`$ over about two orders of magnitude. Deviations of the two distributions for in-degrees $`i<6`$ are due to the underestimation of the degree of external nodes. One explanation for an in-degree exponent of about -1.5 is the assumption of stochastic multiplicative growth as in the model of Huberman and Adamic Huberman and Adamic (1999); Adamic and Huberman (2000). They proposed that the number of links a node received at a time step is a random fraction of the number of links it already has received. The treatment of the out-degree is more difficult. For the whole network, the distribution of out-degree $`j`$, i.e. a node’s number of links to other nodes, shows pronounced scale-free behavior $`n(j)j^{2.03}`$ (Fig. 4). However, the corresponding distribution for internal nodes is broad but does not show scale-free behavior over a sufficient range. This may be caused by the limited size of the sample but may also point to the systematic error caused by students possibly using different (external) accounts for sending e-mails. The out-degree scaling exponent of the whole network lies in a quite common range for communication and social networks, as, e.g., the movie actors’ network or the phone call network Albert and Barabási (2002), where the principle of preferential attachment can be used for modeling Barabási and Albert (1999). It applies to the assumption that the probability $`p_j`$ that a link originates in the set of nodes with out-degree $`j`$ is proportional to the number of links already starting in this set $`f[j]`$:

$$p_jjf[j].$$
(6)
This corresponds to Simon’s general model for such copy and growth processes Simon (1955); Bornholdt and Ebel (2001). Let us briefly apply this model to the e-mail network. From our data we estimated the ratio of the growth rate of nodes to the growth rate of links to $`\alpha =`$0.597 nodes per links which is sufficient to calculate the scaling exponent $`\gamma `$ Bornholdt and Ebel (2001):

$$\gamma =1+\frac{1}{1\alpha }.$$
(7)
Thus, the preferential linking model leads to a steep exponent of -3.48 not in accordance with observation. On the other hand, a model based only on transitive linking Davidsen et al. (2002), i.e. on the assumption that two nodes are more likely to be linked if they have a common neighbor, can in principle reproduce the small-world properties and a broad degree distribution but leads to a too high clustering and does not yield a power-law degree distribution for this particular network. From this we conjecture, that including both preferential and transitive linking may consistently model the e-mail network.

What are the implications of the above results for the spreading of e-mail viruses? The occurrence of e-mail viruses has become a well-known phenomenon in today’s communication experience. An e-mail virus or e-mail worm is a program attached to an e-mail which, when opened by the recipient, causes the recipient’s e-mail program to remail numerous infected e-mails to e-mail addresses found in the address book or in stored e-mails. Hence, for the propagation of e-mail viruses the network is undirected. This is different for chain e-mails, where each recipient is asked to forward the chain e-mail to other addresses. E-mail viruses can cause serious damage to computer networks by destroying data at infected computers or by overloading e-mail servers and other infrastructure. In May 2000, for instance, the e-mail worm “I love you” infected more than 500,000 individual systems worldwide CERT Coordination Center, Carnegie Mellon University  and obstructed 21 % of the computer workplaces in Germany Bundesministerium des Inneren (2001) (Federal Ministry of the Interior of Germany).

In scale-free networks, the threshold for the propagation rate above which an infection of the network spreads and becomes persistent is very much lower than in other disordered networks, or even vanishes Pastor-Satorras and Vespignani (2001b). This means that the self-organized structure of the e-mail network facilitates the spreading of computer viruses, as well as of any other information. In addition, the e-mail network is quite robust in case of “failures” of random nodes if, for instance, some participant does not answer e-mails for a while or uses anti-virus software. However, it is sensitive to the loss of highly connected participants because of the power-law degree statistics Albert et al. (2000). Hence uniformly applied immunization of nodes is less likely to eradicate infections until almost all participants are immunized, whereas targeting prevention efforts at the highly connected sites significantly suppresses epidemic outbreaks and prevalence Pastor-Satorras and Vespignani (2001c); Newman (2002).

These observations suggest helpful and advantageous applications, but also point to the inherent dangers of e-mail networks. The security of e-mail communication can be improved by identifying highly connected hub addresses and monitoring them for viruses more strictly, e.g. in corporate e-mail networks to prevent the damaging and costly spreading of e-mail viruses. In a different application, making use of the high clustering, commercial e-mail providers can identify communities of users more easily Kleinberg and Lawrence (2001) and focus marketing more efficiently. In general, communication by e-mail can be interfered with as well as utilized more extensively due to the non-trivial topological features of the e-mail network that we found here. Exploring the web of e-mails does not only extend our knowledge of social and communication networks but it also shows how vulnerable and exploitable these systems can be.

In conclusion, we have shown that an e-mail network, where nodes are given by e-mail addresses and links by exchanged messages, exhibits both small-world properties and scale-free behavior. The e-mail network is studied in terms of undirected, as well as directed networks. Spreading of e-mail viruses is considered, based on the appropriate viewpoint of an undirected graph. The scale-free nature of the e-mail network strongly eases persistence and propagation of e-mail viruses but also points to effective countermeasures.

###### Acknowledgements.

We thank A.-L. Barabási, J. Davidsen, S. N. Dorogovtsev, and A. Vespignani for useful discussions and comments. H. E. acknowledges support by the Studienstiftung des deutschen Volkes (German National Merit Foundation).


text:
# 1 Introduction

## 1 Introduction

It has long been believed that quantum gravity effects become important at distances comparable to the Planck length, $`\mathrm{}_P`$, and a number of arguments have been presented to support this idea . To show this sort of thing, one mainly needs to show that gravity effects become important at some fixed length scale depending only on the constants $`c`$, $`G`$ and $`\mathrm{}`$. Dimensional analysis does the rest, since $`\mathrm{}_P`$ is the only quantity with dimensions of length that one can construct from these constants. However, in situations where a second length scale becomes relevant, one cannot use dimensional analysis to settle all controversies. For example, Ng and van Dam  have recently argued that quantum gravity effects cause surprisingly large uncertainties in the measurement of a large distance $`\mathrm{}`$, namely

$$\mathrm{\Delta }\mathrm{}\mathrm{}^{1/3}\mathrm{}_P^{2/3}$$

where the symbol $``$ means that we are ignoring a constant factor of order unity. Amelino-Camelia  has gone even further, arguing that

$$\mathrm{\Delta }\mathrm{}\mathrm{}^{1/2}\mathrm{}_P^{1/2}.$$

Uncertainties on this scale are on the brink of being experimentally detectable, lending extra interest to the issue. However, in what follows, we reanalyze the Ng–van Dam thought experiment and show that by modifying its design we can dramatically reduce the uncertainty of distance measurements. Our modified thought experiment gives

$$\mathrm{\Delta }\mathrm{}\mathrm{}_P.$$

## 2 Ng–van Dam Thought Experiment

The elements of the Ng–van Dam thought experiment are straightforward, and the aim is to show that through a simple application of the uncertainty principle, together with limits imposed by general relativity, we are led to the conclusion that a fundamental distance uncertainty arises that may be far larger than the Planck scale.

The argument proceeds as follows. First consider two nearby objects in free fall approximately at rest relative to one another: an observer consisting of a clock and light emitter, and a mirror. If the observer wants to know the distance to the mirror, he may simply emit a burst of light, wait a time $`t`$ for the light to return, and conclude that the mirror is a distance $`\mathrm{}=ct/2`$ away.

Now we are interested in the uncertainty of this measurement. Following an argument due to Wigner  we treat the clock as a free quantum mechanical particle and impose the uncertainty condition $`\mathrm{\Delta }q\mathrm{\Delta }p\mathrm{}`$. Writing $`\mathrm{\Delta }p=m\mathrm{\Delta }v`$ where $`m`$ is the mass of the clock, we thus obtain the following bound on the uncertainty of the clock’s position at time $`t`$:

$$\begin{array}{ccc}\mathrm{\Delta }q(t)& =& \mathrm{\Delta }(q+tv)\hfill \\ & =& \sqrt{(\mathrm{\Delta }q)^2+(t\mathrm{\Delta }v)^2}\hfill \\ & & \mathrm{\Delta }q+\frac{\mathrm{}t}{m\mathrm{\Delta }q}.\hfill \end{array}$$
(1)
To minimize the position uncertainty at time $`t`$, we find that the optimal position uncertainty at time zero should be $`\mathrm{\Delta }q=\sqrt{\mathrm{}t/m}`$. Plugging this back into equation (1), we find the minimum uncertainty at time $`t`$ to be

$$\mathrm{\Delta }q(t)\sqrt{\mathrm{}t/m}.$$
(2)
It is also convenient to write this in terms of the distance to be measured and the Compton wavelength of the clock, $`\mathrm{}_C=\mathrm{}/mc`$:

$$\mathrm{\Delta }q(t)\mathrm{}^{1/2}\mathrm{}_C^{1/2}.$$
(3)
This uncertainty in the position of the clock contributes to the uncertainty in $`\mathrm{}`$, the distance between the clock and mirror. We can ignore the uncertainty in the position of the mirror, which behaves similarly, and obtain this lower bound on $`\mathrm{\Delta }\mathrm{}`$:

$$\mathrm{\Delta }\mathrm{}\mathrm{}^{1/2}\mathrm{}_C^{1/2}.$$
(4)
So far we have only considered the effects of quantum mechanics and the speed of light, with no mention of the effects of general relativity. Next, Ng and van Dam consider the details of the clock itself. They take the clock to consist of two parallel mirrors a distance $`d`$ apart, and consider a tick of the clock to be the time $`2d/c`$ that it takes light to travel back and forth between them. Since we now have the length scale $`d`$ and the mass scale $`m`$ of the clock, we can now begin to consider general relativity effects. In particular, Ng and van Dam assert that the size of the clock, $`d`$, must be larger than its Schwarzschild radius $`\mathrm{}_S=Gm/c^2`$. If the tick of the clock is a lower bound on the accuracy of its time measurements, this requirement implies that

$$\mathrm{\Delta }\mathrm{}\mathrm{}_S.$$
(5)
Finally, squaring the uncertainty from equation (4) and multiplying the result by equation (5), we obtain $`(\mathrm{\Delta }\mathrm{})^3\mathrm{}\mathrm{}_C\mathrm{}_S`$. Note that $`\mathrm{}_C\mathrm{}_S`$ is equal to $`\mathrm{}_P^2`$, the Planck length squared. Thus the primary result obtained from the Ng–van Dam thought experiment is that the minimum uncertainty in this kind of measurement satisfies

$$\mathrm{\Delta }\mathrm{}\mathrm{}^{1/3}\mathrm{}_P^{2/3},$$
(6)
a bound depending only on the distance $`\mathrm{}`$ to be measured and the Planck length.

Ng and van Dam have suggested that this uncertainty is inherent in any distance measurement, and that it may be apparent in the latest generation of interferometers designed for detecting gravitational waves. As explained in Section 3, we believe that this conclusion is unwarranted, and that changes in the details of the clock used in the thought experiment will have an effect on the final uncertainty result, and thus lead to a breakdown of the Ng–van Dam uncertainty as a fundamental property of distance measurement.

Before turning to this we should comment a bit on the Wigner clock thought experiment. Our derivation of equation (2) is rigorous if in the second step we assume that the fluctuations in $`q`$ and $`p`$ are uncorrelated, or more precisely, that the expectation value $`pq+qp`$ equals $`2pq`$. This is true if, for example, the wavefunction for the clock’s center of mass is initially Gaussian. If we allowed the wavefunction to be arbitrary, we could arrange for $`\mathrm{\Delta }q(t)`$ to be arbitarily small at some chosen time $`t`$, simply by taking a wavefunction close to a delta function at time $`t`$ and evolving it backwards to get the wavefunction at time zero. However, for our purposes this counts as ‘cheating’, since we can only reduce $`\mathrm{\Delta }q(t)`$ arbitrarily this way if we already know the time $`t`$ that is to be measured.

## 3 A Modified Thought Experiment

As was pointed out by Adler, Nemenman, Ovenduin, and Santiago , the result of the above thought experiment is only valid when the clock evolves according to the free Schrödinger equation, and this need not be the case. Ng and van Dam  have, in turn, asserted that if the clock is bound through some potential, then it must be bound to something, and that this something can then be considered as a part of the clock. We will explore this line of reasoning through a specific example, and argue that attaching the clock to an external object does indeed affect the final uncertainty result.

For simplicity, let us take the clock to be attached to one end of a rod in free fall. We assume this rod has equilibrium length $`x`$, elastic constant $`Y`$, and mass $`m`$. We take the mass of the clock to be negligible compared to $`m`$. It is thus the rod’s mass, rather than that of the clock, which limits the spreading of the clock’s wavefunction with the passage of time. By making $`m`$ large we can make this spread as slow as we like. However, special relativity puts a bound on the elastic constant $`Y`$, since a very rigid medium would allow sound waves to propagate with speed greater than that of light. This means that the ends of the rod will undergo zero-point fluctuations which also contribute to the uncertainty of any position we measure using the clock.

To calculate these, recall that the speed of sound along the rod will be $`\sqrt{Yx/m}`$. Since this must be less than $`c`$, we have

$$Ymc^2/x.$$
(7)
Now, a one-dimensional rod with one end attached to an immovable wall can perform oscillations about its equilibrium length $`x`$ exactly as a harmonic oscillator with spring constant $`k=Y/x`$ and oscillator mass $`m_{osc}=m/3`$. Thus, a rod in space will have a mode of oscillation that allows each end to perform harmonic oscillations about its respective equilibrium position with equivalent spring constant $`k=2Y/x`$ and oscillator mass $`m_{osc}=m/6`$. Treating the problem quantum-mechanically, the ground state of such an oscillator has position uncertainty

$$\mathrm{\Delta }x=\sqrt{\frac{\mathrm{}}{2(km_{osc})^{1/2}}}=\sqrt{\frac{\mathrm{}}{2}\left(\frac{3x}{mY}\right)^{1/2}}$$
(8)
Using equation (7) and dropping factors of order unity, we obtain

$$\mathrm{\Delta }x\sqrt{\mathrm{}x/mc}$$
(9)
There are other states where $`\mathrm{\Delta }x`$ is smaller at one moment, but not for the entire duration of the experiment, at least if the rod’s period of oscillation is short compared to the time $`t`$ which the clock is to measure — and if it were not, the rod would fail to slow the spread of the clock’s wavefunction. We may thus take this value of $`\mathrm{\Delta }x`$ as an approximate lower bound on the uncertainty of the clock’s position with respect to rod’s center of mass. In terms of the Compton wavelength of the rod, it follows that

$$\mathrm{\Delta }xx^{1/2}\mathrm{}_C^{1/2}.$$
(10)
The rod’s center of mass will also have an uncertainty in its position, given by equation (3). Both this and $`\mathrm{\Delta }x`$ will contribute to the uncertainty of any distance we measure by sending a burst of light from the clock and measuring the time it takes for the light to return. Assuming these uncertainties are uncorrelated, we obtain

$$\mathrm{\Delta }\mathrm{}(\mathrm{}^{1/2}+x^{1/2})\mathrm{}_C^{1/2}.$$
(11)
Next we must rethink the requirements imposed by general relativity on the clock–rod system. The requirement made in the original Ng–van Dam thought experiment was that the clock must be larger than its own Schwarzschild radius. While this is still true, it is not important now, since we no longer need the clock to be heavy to prevent its wave function from spreading with time. Instead, the important relativity consideration is that the rod be longer than its Schwarzschild radius:

$$x\mathrm{}_S$$
(12)
Although we do not have spherical symmetry in this example, the Schwarzschild radius remains useful as a quick and dirty estimate of the length scale at which the rod would collapse to form a black hole. Indeed, the ‘hoop conjecture’  says roughly that any system compressed within its Schwarzschild radius must give rise to a singularity.

Plugging this lower bound on $`x`$ into equation (10), we obtain

$$\mathrm{\Delta }\mathrm{}(\mathrm{}^{1/2}+\mathrm{}_S^{1/2})\mathrm{}_C^{1/2}.$$
(13)
We can minimize this uncertainty by making the rod very heavy, so that $`\mathrm{}_C0`$, leaving us with

$$\mathrm{\Delta }\mathrm{}\mathrm{}_P.$$
(14)
It is interesting to compare this argument to that given by Ng and van Dam. First, unlike their argument, ours does not assume any limitation on the clock’s accuracy solely due to its size. Second, while their result was obtained by multiplying two independent lower bounds on $`\mathrm{\Delta }\mathrm{}`$, one from quantum mechanics and the other from general relativity, ours arises from an interplay between competing effects. On the one hand, we wish to make the rod as heavy as possible to minimize the quantum-mechanical spreading of its center of mass. To prevent it from becoming a black hole, we must also make it very long. On the other hand, as it becomes longer, the zero-point fluctuations of its ends increase, due to the relativistic limitations on its rigidity. We achieve the best result by making the rod just a bit longer than its own Schwarzschild radius.

We should also mention and rebut a possible objection to our argument. One could try to eliminate the effect of zero-point fluctuations on the position of the clock by attaching the clock to the midpoint of the rod rather than the end. This would indeed eliminate the effect of the vibrational mode we have been considering. However, there are other modes in which the midpoint of the rod oscillates. If we take these into account, some of the numbers in our formulas change, but we again find that $`\mathrm{\Delta }\mathrm{}`$ is greater than $`\mathrm{}_P`$ times some constant of order unity.

## 4 Conclusions

It appears that by attaching the clock used in the Ng–van Dam thought experiment to a heavy object, the surprisingly large uncertainties arising in distance measurements can be avoided. It is thus our opinion that these large uncertainties cannot be a fundamental feature of the union of quantum mechanics and general relativity.

Following Ng and van Dam, most recent papers on this topic suggest laser interferometry as a possible way to actually detect the uncertainties in measurements of distance caused by quantum gravity effects . While this is tempting, it is worth noting that interferometers do not measure distances so much as periodic changes in distance. Indeed, the design of an interferometer is very different from the Ng–van Dam thought experiment or our modified version, so any attempt to detect quantum gravity effects at LIGO will require a separate analysis.