text:
# A simple model for a transverse dune field

## 1 Introduction

There are over hundred different types of dunes that have been classified by geomorphologists according to their shape and size . Dunes can be found in deserts, in coastal areas, on the sea-bottom , and even on Mars . Although no general theoretical model has been developed up to now to explain this wide variety of dune shapes, it is well known that the dune morphology is mainly determined by the amount of available sand and by the behavior of the wind over the year. The simplest type of dunes occurs when there is not very much sand, and when the direction of the blowing wind doesn’t change with time. These dunes have a crescent-like shape and are called “barchans”. Since the pioneering work of Bagnold in 1941 , several field measurements led to a good understanding of barchan dunes. In particular, it has been observed that these dunes move without changing their shape with a velocity which is inversely proportional to their height . Fields with hundreds of barchans are found for instance in Peru , Namibia  and Marroco . In a recent work , a two-dimensional model to describe the collective motion of barchan dunes has been presented, which was able to reproduce the same kind of spatial correlations between the dunes as in real dune fields.

As the sand availability becomes larger, parallel dunes perpendicular to the wind direction, which are called “transverse dunes”, arise for a wind blowing steadily from the same direction. Ensembles of transverse dunes, with crest to crest distances ranging from a few meters to over 3 km, are found to cover about 40% of all terrestrial sand seas, being common in the Northern Hemisphere , and also dominating the sand landscape on Mars . The last few years have seen an increasing interest in the investigation of the geomorphological aspects and the evolution of transverse dunes . In spite of this effort, very little is known about these dunes, and several phenomena are still not well understood, as for example the effect of coalescence between dunes with different velocities .

This work reports on a simple one-dimensional model for the evolution of transverse dunes in a field, which reproduces the observation that these dunes have generally the same height and velocity of propagation. This is independent on the initial configuration of the field and on coalescence effects between close dunes with different velocities. In our model each dune is considered to have zero width and infinite length, but a variable height, which determines the associated propagation velocity and that also controls the interaction between neighbouring dunes. The paper is organized as follows. In section 2, the model is presented; the results for both cases with and without coalescence are discussed in section 3; conclusions are made in section 4.

## 2 The Model

We consider a set of $`N`$ dunes distributed along the $`x`$ axis, each one with a variable height $`h_i`$, but zero width. The dunes are allowed to move in the $`x`$ direction with a velocity $`v_idx_i/dt`$ given by:

$$v_i=\frac{a}{h_i},$$
(1)
where $`x_i`$ is the position of the $`i`$th dune, $`i=1,2,\mathrm{},N`$, and $`a`$ is a constant phenomenological parameter that contains information about the wind speed, the grain size, etc.

The behavior of the sand flux on a dune field is crucial to the evolution of the field dynamics, as interactions between dunes strongly depend on the sand transport among them. For a dune that migrates shape invariantly through sand avalanches on its slip face, the sand flux $`\varphi `$ over the dune behaves as $`d\varphi /dh=`$constant, where $`h(x)`$ is the height profile of the dune. For the model studied in this paper, we consider the height $`h_i`$ of the dune at the position of the crest, and set for the $`i`$th dune:

$$\varphi _i=bh_i,$$
(2)
where $`b`$ is our second phenomenological parameter, again depending on wind speed, grain size, etc. Note that equation (2) means that the rate of sand transport over the crest is larger for the bigger dunes than for the smaller ones. As a consequence, exchange of sand is expected between neighbouring dunes with different heights, as a bigger dune tends to loose sand more rapidly than a smaller one.

To model the inter-dune interaction, we first observe that for a two-dimensional slice of a dune, the mass of sand in the dune, or equivalently, its area $`A`$, increases with the height of the dune as $`Ah^2`$. On the other hand, the amount of sand transported depends on the difference between the sand fluxes on the crests of the neighbouring dunes, i.e. if $`dA_i`$ is the volume change for the $`i`$th dune during a time interval $`dt`$, then $`dA_i(dh_i)^2(\varphi _i\varphi _{i1})`$. Thus, we write:

$$\frac{dh_i}{dt}=c\times \sqrt{|\varphi _i\varphi _{i1}|}\times \text{sign}(\varphi _i\varphi _{i1}),$$
(3)
where $`c`$ is the third parameter of our model that depends on the shape of the dune, and the last factor on the right side ensures that more sand is transported from the larger dune to the smaller one. As in a real transverse dune field, we set the sand influx to point in the positive wind ($`x`$) direction and assume it to be constant. We call the sand influx $`\varphi _0`$, which can be seen as an effective dune fixed at the position $`x_0=0`$, with a constant height $`h_0=\varphi _0/b`$. Equations (1), (2) and (3) are repeatedly iterated, which gives the time evolution of the dune field.

## 3 Results and discussion

We studied the evolution of a set of $`N`$ dunes with initial positions $`x_i(0)`$ and heights $`h_i(0)`$, for different configurations with $`N=10`$, as shown in Figure 1. In the first group of simulations (I), the dunes are equally spaced along the $`x`$ axis at $`t=0`$, and have the same heights $`h_i(0)=H,i=1,2,\mathrm{},N`$ (Figure 1(a)). In the second group (II), we take random values for the heights $`h_i(0)`$ chosen from a homogeneous distribution between 0 and 1, but still the dune spacing is regular (Figure 1(b)). And in the third group (III), both the values for the heights and the positions of the dunes along the $`x`$ axis are chosen randomly (Figure 1(c)).

The trajectory of the dunes along the $`x`$ axis and the changes in their heights are shown as a function of the time in Figures 2(a) and 2(b) respectively, for particular realizations corresponding to each one of the initial states I, II and III. The values of the phenomenological parameters are $`a=300\text{m}^2/\text{year}`$, $`b=1500\text{year}^1`$ and $`c=0.05`$, and the sand influx is chosen to be $`\varphi _0=2.25\times 10^4`$m/year $`4.3`$cm/min. In the simulations, each iteration corresponds to $`0.01`$ year. As can be seen from these figures, after a certain number of iterations $`tT_{\mathrm{}}`$, all dunes reach the same height $`h_i(T_{\mathrm{}})=\varphi _0/b=15`$m, and consequently the same velocity $`v_i(T_{\mathrm{}})=a/h(T_{\mathrm{}})=(a\times b)/\varphi _0=20`$m/year, independently on the initial configuration.

In the discussion above, the distances between the dunes were large enough that their relative order was the same during the whole evolution of the field. However, if a smaller, faster migrating dune is very closely behind a larger, slower one, they can reach the same position after some finite time. Dune coalescence has been studied for instance in Ref. , in which a critical distance $`d_{\mathrm{min}}`$ between two neighbouring dunes was considered, below which they collide and become one single dune. In our one-dimensional model, we set $`d_{\mathrm{min}}=0`$. For simplicity, dune coalescence is implemented here through the emergence of a single dune, with height $`h_\mathrm{c}=\sqrt{h_i^2+h_{i+1}^2}`$. Consequently, the number of dunes $`N`$ in the field decreases by unity, whenever two dunes coalesce.

The main plot of Figure 3 shows the variable $`x_i(t)`$ for a set of $`N(0)=20`$ dunes with random values of $`x_i(0)`$ and $`h_i(0)`$ for $`t=0`$. The values of $`\varphi _0`$ and of the parameters $`a`$, $`b`$ and $`c`$ are the same as in Figures 2(a) and 2(b), but now their distances in the initial state are much smaller. The number of dunes as a function of time for the realization in this figure is shown in the upper-left-hand-corner inset. Note that this number decreases monotonically until the last coalescence between dunes has occured, and the final number of dunes, $`N(\mathrm{})=3`$ has been reached. As we can see from the bottom-right-hand-corner inset, coalescence of dunes doesn’t modify the result that the dunes in the field do all have the same height in the final state.

For a given distribution of dunes in the field at $`t=0`$, the number $`N(\mathrm{})`$ of dunes in the final state tends to be larger for smaller values of the parameter $`a`$, the other parameters $`b`$, $`c`$ and $`\varphi _0`$ being fixed. This is because $`a`$ controls the velocity of the dunes, and coalescence is more frequent in a field with faster dunes. In the same way, the closer the dunes are in the initial state, the more they are going to coalesce on average. As an example, we show in Figure 4 the behavior of $`N(\mathrm{})`$ as a function of the length $`L`$ of the system, i.e. the initial distance between the first dune and the last one in the field, for a set having a fixed number of $`N(0)=50`$ dunes, using the same values of the phenomenological parameters as in Figure 1. As we can see in this figure, the smaller (larger) the value of $`L`$, the more (less) the dunes tend to coalesce. The statistical error bars are obtained from different realizations of dunes with random heights and positions at time $`t=0`$. The full line shown in Figure 4 corresponds to the fit:

$$N(\mathrm{})=N(0)[1q\mathrm{exp}(L/L^{})],$$
(4)
where $`N(0)=50`$, $`q1`$, and $`L^{}1.34\times 10^3`$m.

In Figure 5, we show for three different values of $`L`$ the corresponding curves $`N_L(t)`$, averaged over five realizations of $`(x_i(0),h_i(0))`$, as a function of time $`t`$. The values of $`L`$ shown are $`L=100`$m (circles), $`L=500`$m (triangles) and $`L=1000`$m (squares). The dashed line corresponds to a straight line in a linear-log fit to the curve corresponding to $`L=100`$m. We can observe that $`N_L(t)`$ decays approximately logaritmically with time, independently on $`L`$. After a variable time interval $`T_cT_{\mathrm{}}`$, coalescence stops and an asymptotic value $`N(\mathrm{})`$ is reached.

As all dunes attain the same height $`h_0\varphi _0`$ for $`t=T_{\mathrm{}}`$, the time to reach the final state depends on the strength of the sand flux in the field, $`\varphi _0`$. A plot of $`T_{\mathrm{}}`$ as a function of $`\varphi _0`$ for the initial configuration I of Figure 1a is presented in Figure 6. As we can see from this figure, $`T_{\mathrm{}}`$ scales after a transient interval of variation in $`\varphi _0bh_i(0)`$:

$$T_{\mathrm{}}\sqrt{|\varphi _0\varphi |},$$
(5)
i.e. the closer the value of $`\varphi _0`$ to $`\varphi =bh_i(0)=bH`$, the more rapidly the system will reach the final state. From equation (3), we see that the phenomenological parameter $`c`$ controls the rate of sand transport between neighbouring dunes. We have found that the number of iterations to reach the final state decreases as $`1/c`$, which in fact is independent on the initial height $`H`$ of the dunes.

## 4 Conclusions

We presented a simple one-dimensional model for a transverse dune field which reproduces the observation that the dunes in a field attain approximately the same heights and velocities. The fact that our results do not depend on the initial configuration of the field and on the details of coalescence effects between dunes is of particular interest, since the genesis of dune fields is still not understood by geographers and geologists. Of course, serious limitations emerge from considering dunes as points in a line. Further work is required to compare the results presented here with real dune field observations.

We acknowledge V. Schwaemmle for his contribution in the initial stage of this work. We acknowledge A. O. Sousa for useful comments. E. J. R. Parteli acknowledges support from CAPES - Brasília/Brazil.


format: text

text:
# Nonlinear Bessel beams

## 1 Introduction

Diffractive spreading of waves is a classical phenomenon in wave dynamics and an inherent feature of beam propagation. Much attention has been devoted to the possibility of counteracting the dispersive spreading by focusing effects due to medium nonlinearities e.g. the Kerr effect, cf. . However, it has also been pointed out, , that non-diffracting beams are possible also in linear media. In particular, the Helmholz equation that governs the linear diffractive dynamics of a wave beam allows classes of diffraction-free solutions. In addition to the plane wave solutions, the two-dimensional counterparts, the cylindrically symmetric Bessel solutions, also propagate with preserved form, while also allowing for a concentrated beam profile. The drawback from an application point of view is the fact that these beams have infinite energy, and consequently cannot be realized physically. Various ways to circumvent this problem have been suggested, the most obvious being to truncate the Bessel beam at some radius e.g. by a Gaussian truncation, forming the so called Bessel-Gauss beams, . While such a truncation clearly reintroduces diffraction, the beam broadening could be made small if the propagation length is kept smaller than the corresponding diffraction length of the Bessel-Gauss beam. In particular, since the Bessel beam diffracts sequentially, starting with the outer lobes, cf. , the central part of the beam remains intact for a certain distance of propagation.

Recently, there has been growing interest in nonlinear effects in connection with Bessel and Bessel-Gauss beams, . Of special interest for the present investigation is the attention given to media with an intensity dependent refractive index, i.e., Kerr media, see e.g. . The work carried out in  considers the limit of weak nonlinearity, which makes it possible to use a perturbation approach involving an expansion around the lowest order linear (stationary) Bessel solution for solving the evolution equation, being the nonlinear Schrödinger equation.

In the present paper, we investigate in more detail the nonlinear generalisation of the linear diffraction-less Bessel beam solutions as well as the nonlinear dynamics of Bessel-Gauss beams. Stationary solutions in two dimensions are determined by the Bessel equation modified by a nonlinear term, i.e., the radially symmetric nonlinear Schrödinger equation. The modified Bessel solutions, the “nonlinear Bessel beams”, are studied using approximate analytical and numerical methods. The results show that the nonlinearity primarily affects the central high intensity parts of the beam profile, which become radially compressed or decompressed depending on whether the nonlinearity is focusing or defocusing, respectively. The beam profile for large radii remains a Bessel function with a phase shift being the only remaining effect of the nonlinearity. However, for the defocusing nonlinearity an amplitude threshold exists, above which no solutions decaying to zero exist.

The dynamical properties of Gaussian-truncated Bessel beams in the presence of a Kerr nonlinearity are also studied. An exact analytical solution was previously found for the linear dynamics of the Bessel-Gauss beams, . Based on the virial theorem, which gives an exact analytical description of the variation of the beam width in the root-mean-square (RMS) sense, important information about the effect of the nonlinearity on the beam dynamics can be obtained. In particular, it is found that a focusing nonlinearity tends to cause an evolution stage where the central parts of the Bessel beams are initially compressed. Depending on the strength of the nonlinearity, different scenarios are possible, e.g. the subsequent evolution may involve an essentially diffraction-dominated behaviour, but for increasing nonlinearity, two forms of collapse may appear. Either a part of the beam collapses, while the RMS width of the beam still increases, or above a certain threshold, the RMS width goes to zero in a finite distance. Numerical simulations of the dynamics illustrate the different scenarios.

## 2 The nonlinear Schrödinger equation

The propagation of an optical wave in a nonlinear Kerr medium is described by the nonlinear Schrödinger equation. This implies that the slowly varying wave envelope, $`\psi (z,r)`$, of a cylindrically symmetric beam satisfies the following equation

$$\mathrm{i}\frac{\psi }{z}=\frac{1}{2k_0}\left(\frac{^2\psi }{r^2}+\frac{1}{r}\frac{\psi }{r}\right)+\kappa |\psi |^2\psi ,$$
(1)
where $`z`$ is the distance of propagation, $`k_0`$ is the wave number in vacuum, and $`\kappa `$ is the nonlinear parameter. Additional physical effects, e.g. attenuation and gain, can be modelled by using complex coefficients in Eq. (1). The obtained complex equation, which in one dimension has analytical soliton solutions, , is the cylindrical generalisation of the Pereira-Stenflo equation. It has been investigated using a variational approach, , but further work is needed to fully describe the complex case. In the present work, the coefficients are assumed to be real. It is convenient to introduce the normalisation $`\stackrel{~}{r}=r/a_0`$, where $`a_0`$ is a characteristic width of the beam, and $`\stackrel{~}{z}=z/L_D`$, where $`L_D2k_0a_0^2`$ is the Rayleigh length. Eq. (1) then takes the form

$$\mathrm{i}\frac{\psi }{\stackrel{~}{z}}=\frac{^2\psi }{\stackrel{~}{r}^2}+\frac{1}{\stackrel{~}{r}}\frac{\psi }{\stackrel{~}{r}}+\stackrel{~}{\kappa }|\psi |^2\psi ,$$
(2)
where $`\stackrel{~}{\kappa }=L_D\kappa `$. For simplicity we will suppress the tilde in the subsequent expressions.

We begin the analysis by looking for stationary solutions of Eq. (2). For this purpose we write $`\psi =\psi (z,r)=A(r)\mathrm{}^{\mathrm{i}\delta z}`$, which leads to the eigenvalue equation

$$\frac{d^2A}{dr^2}+\frac{1}{r}\frac{dA}{dr}+\delta A+\kappa A^3=0.$$
(3)
This equation is to be solved subject to the boundary conditions that the solution should be finite when $`r=0`$ and go to zero as $`r\mathrm{}`$. The lowest order solution in the physical situation when the nonlinearity balances the diffraction, i.e., the focusing case with $`\kappa >0`$, corresponds to the so called Townes soliton , which has essentially the same properties and $`sech`$-shaped form as the one-dimensional soliton solution, cf. . In particular, this solution only exists for negative eigenvalues, which are uniquely related to the maximum amplitude.

The Bessel beams are the solutions of the linear Schrödinger equation ($`\kappa =0`$) and are given by $`A(r)=A_0J_0(\sqrt{\delta }r)`$. Clearly, well-behaved solutions exist only for positive eigenvalues $`\delta `$. In contrast to the nonlinear case, the linear eigenvalue problem has a continuous set of (positive) eigenvalues, which are independent of the amplitude of the beam profile. The first task of the present analysis is to analyse the nonlinear Bessel beams, being the solutions of Eq. (3) for positive $`\delta `$ and $`\kappa 0`$. We note that by introducing

$$\overline{r}=\sqrt{|\delta |}r,\overline{A}=\sqrt{|\kappa /\delta |}A,$$
(4)
only the signs of $`\delta `$ and $`\kappa `$ remain in Eq. (3). Thus, without loss of generality, it will be assumed that $`\delta =1`$ and $`\kappa =\pm 1`$.

## 3 Nonlinear Bessel beams

In order to analyse the properties of the nonlinear Bessel beams analytically it is instructive to start by examining the central part of the pulse, which is determined by $`\sqrt{\delta }r1`$. Since the initial derivative of $`A(r)`$ must be zero, we have the approximation $`A(r)=A_0+O(r^2)`$, which implies that to second order in $`r`$, Eq. (3) can be approximated by the linear equation

$$\frac{d^2A}{dr^2}+\frac{1}{r}\frac{dA}{dr}+(\delta +\kappa A_0^2)A=0,$$
(5)
with the corresponding solution

$$A=A_0J_0(\sqrt{\delta +\kappa A_0^2}r).$$
(6)
This solution is valid for small $`r`$ only, but nevertheless gives important information about the nonlinear modifications of the Bessel beam. In the focusing case ($`\kappa >0`$), the main lobe tends to be compressed and we expect that the amplitude of the nonlinear Bessel beam will oscillate more rapidly than the linear Bessel beam. On the other hand, in the defocusing case ($`\kappa <0`$), the main lobe should become wider and the nonlinear solution should oscillate slower than the linear one. In particular, if the amplitude is chosen to fulfil

$$\delta +\kappa A_0^2<0,$$
(7)
the expression under the square root will be negative. This corresponds to the modified Bessel function, which is growing with $`r`$. Thus the presence of a defocusing nonlinearity can qualitatively change the behaviour of the solution.

It is clear that if the second derivative of $`A`$ is positive initially, it will remain positive. This can be seen by rewriting the equation as

$$\frac{d^2A}{dr^2}+\frac{1}{r}\frac{dA}{dr}+\delta _{\mathrm{eff}}A=0,$$
(8)
where $`\delta _{\mathrm{eff}}\delta +\kappa A^2`$. If $`\delta _{\mathrm{eff}}`$ is negative at $`r=0`$, the solution will be growing for small $`r`$ and $`\delta _{\mathrm{eff}}`$ will be further decreased. This increases the derivative of the solution, and implies that a sufficiently strong defocusing nonlinear term will give rise to a monotonically growing solution, which is not compatible with the condition at infinity.

A more accurate description of the main lobe of the solution can be obtained using variational analysis and Ritz optimisation, cf. . When using variational analysis, it is important to find a good set of trial functions that gives tractable calculations while maintaining sufficient accuracy. A trial function that should approximate the main lobe of the nonlinear Bessel beam reasonably well is

$$A_T=A_0J_0\left(j_0\frac{r}{r_0}\right),$$
(9)
where $`j_0`$ is the first zero of the Bessel function. This also has the advantage that the exact linear result is recovered in the limit $`\kappa 0`$. In the variational procedure, we assume $`r_0`$ to be given and consider $`A_0`$ as a free parameter. The Lagrangian corresponding to Eq. (3) is

$$L=_0^{r_0}L[A_T]r,$$
(10)
where

$$L[A]=\frac{r}{2}\left(\frac{dA}{dr}\right)^2\frac{r\delta A^2}{2}\frac{r\kappa A^4}{4}.$$
(11)
Using the ansatz (9) we obtain

$$=\frac{1}{4}\left[j_0^2A_0^2J_1^2(j_0)r_0^2\delta A_0^2J_1^2(j_0)C_1r_0^2\kappa A_0^4\right],$$
(12)
where

$$C_1=_0^1J_0^4(j_0x)xx7.62\times 10^2.$$
(13)
The variation with respect to $`A_0`$, i.e., $`/A_0=0`$ yields

$$A_0^2=\frac{(j_0^2r_0^2\delta )J_1^2(j_0)}{2C_1r_0^2\kappa }r_0^2=\frac{j_0^2J_1^2(j_0)}{\delta J_1^2(j_0)+2C_1\kappa A_0^2}.$$
(14)
By setting $`\kappa =0`$, the linear result $`\sqrt{\delta }r_0=j_0`$ is recovered. If $`\kappa `$ is positive the value of $`r_0^2`$ is decreased, which corresponds to compression of the main lobe. A negative $`\kappa `$ gives the opposite effect. This result confirms the previously obtained picture of the effects of the nonlinearity on the main lobe of the linear Bessel solution.

The result of the variational analysis is compared with numerical solutions in Figs. 1 and 2. Different $`r_0`$-values have been used, and they are easily identified in the figures, since the variational approximation is zero when $`r=r_0`$. For clarity the plotted curves have been normalised with respect to their amplitudes at $`r=0`$. It is seen that the Bessel ansatz represents a rather good approximation in the focusing case, but that the presence of the nonlinearity changes the shape for small $`r`$, making it more peaked than the Bessel profile. In the defocusing case, the nonlinear solution instead has a flatter form than the linear Bessel function, Fig. 2. It is also seen that for increasing $`r_0`$, or equivalently increasing $`A_0`$, the approximation deteriorates. This is due to the fact that in this case there exists a threshold value for the amplitude in order to have well-behaved solutions, cf. (7). The critical value for the amplitude is $`A_0=1`$. The variational result also predicts this behaviour, although, as is inferred from Eq. (14), the critical value is found to be slightly different

$$A_0>\sqrt{\delta J_1^2(j_0)/(2C_1\kappa )}1.33.$$
(15)
Clearly, as $`A(0)`$ approaches the threshold value, which is unity, we expect the accuracy of the variational approximation to deteriorate.

We now turn to an investigation of the overall behaviour of the nonlinear Bessel beam profiles. The fact that the main influence of the nonlinear term is a rescaling of the radial coordinate makes it reasonable to look for an approximate solution of the form

$$A(r)A_0J_0(f(r)),$$
(16)
where $`f(r)`$ is a function of $`r`$, $`\delta `$, $`\kappa `$, and $`A_0`$. It is difficult to determine $`f`$ using analytical methods, but by noticing that the linear solution can be written as

$$A=A_0J_0\left(_0^r\sqrt{\delta }r^{}\right),$$
(17)
and by comparing with Eq. (5) it seems reasonable that a good approximation should be obtained by the implicit expression

$$A=A_0J_0\left(_0^r\sqrt{\delta +\kappa A^2(r^{})}r^{}\right).$$
(18)
Although this is, in fact, an integral equation for $`A(r)`$, it nevertheless provides a very simple formula for finding $`A`$ numerically. The corresponding approximate solution is compared to the numerical solution of the full equation in Fig. 3. When the amplitude is low the two curves are identical, since the ansatz then reduces to the Bessel function. In the case of a focusing nonlinearity, there is good agreement between the two approaches, but it is also seen that a phase shift appears between the curves for increasing $`A_0`$. Quite good agreement is seen also in the defocusing case. In particular, the initial flattening is well modelled. The phase shift is now of the opposite sign.

This approximate solution implies that the argument of the Bessel function increases approximately as $`_0^r\sqrt{\delta +\kappa A^2}r^{}`$, which is a nonlinear generalisation of the linear case. Thus the main effect of a focusing nonlinearity is to increase the curvature of the peaks by increasing the growth rate of the argument, making the solution radially compressed. In the defocusing case the curvature is decreased, which is most clearly seen in the main lobe.

Finally, Figs. 4 and 5 further illustrate the nonlinear deformations of the linear diffraction-less Bessel solutions. The numerically obtained curves clearly show the features discussed above; the radial compression of the central lobe in the focusing case and the radial expansion in the defocusing case. The expansion effect in the latter case rapidly increases as the amplitude approaches the critical value $`A_0=1`$, above which no stationary solutions are possible. The phase shifting effect of the nonlinearity on the Bessel-like oscillations is also seen, the shift changing sign with the sign of the nonlinearity.

## 4 Analysis of nonlinear Bessel-Gauss beams

The linear diffraction properties of Bessel-Gauss beams have been analysed and solved analytically, see . In the present section we will analyse the nonlinear dynamics of beams, which initially have a profile in the form of a Bessel function truncated by a Gaussian. Since a general solution of this problem cannot be given, we will use the virial theorem to obtain analytical information and numerical simulations for determining the evolution of the beam profile.

The virial theorem, see e.g. , provides exact and explicit information about the dynamic variation of the width of the beam, $`\sigma `$, defined in the RMS sense as $`\sigma ^2r^2`$, where

$$f(r)\frac{_0^{\mathrm{}}f(r)|\psi (z,r)|^2rr}{_0^{\mathrm{}}|\psi (z,r)|^2rr}.$$
(19)
The virial theorem asserts that

$$\frac{d^2\sigma ^2}{dz^2}=8\frac{H}{I}=\text{constant},$$
(20)
where $`I`$ and $`H`$ are invariants of the two-dimensional nonlinear Schrödinger equation, Eq. (1), and are defined as follows

$`I`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}|\psi (z,r)|^2rr,`$ (21)
$`H`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}\left[\left|{\displaystyle \frac{\psi (z,r)}{r}}\right|^2{\displaystyle \frac{\kappa }{2}}|\psi (z,r)|^4\right]rr.`$ (22)
This means the invariants correspond to the (integrated) beam intensity and the Hamiltonian. Thus, the virial theorem implies that $`\sigma ^2`$ must be a second order polynomial in $`z`$, with coefficients determined by the initial beam profile, $`\psi (0,r)`$. For initial phase functions that do not depend on $`r`$, the linear term in $`z`$ vanishes, and the beam width varies as

$$\sigma ^2(z)=\sigma ^2(0)\left(1+𝑠𝑖𝑔𝑛(H)\frac{z^2}{L_0^2}\right)$$
(23)
where $`L_0`$ is a characteristic length given by

$$L_0^2=\left|\frac{4H}{\sigma ^2(0)I}\right|=\left|\frac{4_0^{\mathrm{}}\left[\left|\frac{\psi (z,r)}{r}\right|^2\frac{\kappa }{2}|\psi (z,r)|^4\right]rr}{_0^{\mathrm{}}r^2|\psi (0,r)|^2rr}\right|.$$
(24)
Clearly this approach cannot be used for analysing the linear, or the nonlinearly modified, stationary Bessel beam solutions of the nonlinear Schrödinger equation since all integrals involved in the virial theorem are infinite. However, for a physical beam, with finite integral content, the virial theorem is useful. In general it is seen that with weak nonlinear focusing effects, the Hamiltonian is positive and the RMS width will increase quadratically with a characteristic diffraction length given by $`L_0`$. When the amplitude of the beam increases, the Hamiltonian decreases and eventually changes sign. This implies that the RMS width goes to zero after a finite length equal to $`L_0`$—the well known nonlinear collapse phenomenon, where $`L_0`$ now plays the role of collapse length.

For Bessel-Gauss beams, , the initial profile is of the form

$$\psi (0,r)=A_0J_0\left(\frac{r}{r_0}\right)\mathrm{exp}\left(\frac{r^2}{2\rho _0^2}\right).$$
(25)
Inserting this into Eq. (24) we obtain the following expression for the characteristic length $`L_0`$

$$\left(\frac{r_0}{2L_0}\right)^2=\frac{S_1(\mu )\mathrm{\Lambda }S_2(\mu )}{S_3(\mu )}$$
(26)
where $`\mu =r_0^2/\rho _0^2`$, $`\mathrm{\Lambda }=\kappa A_0^2r_0^2`$, and the integrals $`S_n`$, $`n=1,2,3`$, are given by

$`S_1(\mu )`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}[\mu xJ_0(x)+J_1(x)]^2\mathrm{exp}(\mu x^2)xx,`$ (27)
$`S_2(\mu )`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}J_0^4(x)\mathrm{exp}(2\mu x^2)xx,`$ (28)
$`S_3(\mu )`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}x^2J_0^2(x)\mathrm{exp}(\mu x^2)xx.`$ (29)
Since we are primarily interested in the case when the Gauss function truncates the outer parts of the Bessel function, we have $`\mu =r_0^2/\rho _0^21`$. In this limit, the asymptotic values of the integrals, $`S_n`$, $`n=1,2,3`$, are obtained analytically as

$`S_1(\mu )`$ $``$ $`{\displaystyle \frac{1}{2\sqrt{\pi }}}{\displaystyle \frac{1}{\sqrt{\mu }}},`$ (30)
$`S_2(\mu )`$ $``$ $`D_1{\displaystyle \frac{3}{4\pi ^2}}\mathrm{ln}\mu ,`$ (31)
$`S_3(\mu )`$ $``$ $`{\displaystyle \frac{1}{4\sqrt{\pi }}}{\displaystyle \frac{1}{\mu ^{3/2}}},`$ (32)
Since $`S_2(\mu )`$ goes rather slowly towards infinity as $`\mu `$ becomes small, it is necessary to determine the constant $`D_1`$ in order to have good accuracy for finite $`\mu `$. Using numerical evaluation of the integral we find $`D_10.202`$. This implies that the characteristic length can be approximated as

$$\left(\frac{r_0}{2L_0}\right)^22\mu \left[1\mathrm{\Lambda }\sqrt{\mu }\left(D_2\frac{3\mathrm{ln}\mu }{2\pi ^{3/2}}\right)\right],$$
(33)
with $`D_20.715`$. In the linear case, the characteristic length is seen to scale simply as $`L_0\rho _0`$, i.e., the diffraction is determined solely by the truncation radius and as $`\rho _0\mathrm{}`$, the non-diffracting Bessel beam is recovered. For increasing values of the nonlinearity parameter, $`\mathrm{\Lambda }`$, but for a fixed truncation radius, the value of $`L_0`$ increases and for a certain critical value of $`\mathrm{\Lambda }`$, the nonlinearity balances the diffraction to give a Bessel-Gauss beam, which is diffraction-less in the RMS sense.

## 5 Dynamics of Bessel-Gauss beams

When a truncated linear Bessel-Gauss beam propagates, the diffraction initially affects only the outermost part of the pulse, where the truncation has changed it from the Bessel shape. The central parts are initially diffraction balanced and remains so until the “diffraction front” propagating inwards from the outer parts eventually reach the inner lobes and also these parts start to diffract outwards. On the other hand, the nonlinear effect is strongest at the centre of the beam, where the intensity is highest, and with a focusing nonlinear term, the main lobe will start to compress. In fact, it will start to compress irrespective of the degree of nonlinearity since the linear diffraction is already balanced. If the nonlinear effect is weak the compression will eventually stop, and diffraction will become the dominating effect. This evolution is illustrated in Fig. 6, where an FDTD simulation using $`A_0=1`$ and $`\mu =0.01`$ is shown. Although the central parts initially compress, the virial theorem predicts beam broadening in the RMS sense. Clearly this is no contradiction since the broadening of the outer parts more than compensate the compression of the centre. In order to further illustrate how the main lobe is compressed, the intensity at $`r=0`$ has been plotted for different initial amplitudes as a function of propagation distance in Fig. 7. The curves have been normalised with respect to their amplitudes at $`r=0`$ and $`\mu =0.01`$. We emphasise the oscillating behaviour for the highest amplitudes. If the nonlinear effect is sufficiently strong, the virial theorem predicts a collapse of the beam to zero RMS width. According to Eq. (33), the amplitude threshold for this type of behaviour is $`A_0>A_c2.26`$ for $`\mu =0.01`$. It is well known that the nonlinear evolution of two-dimensional beams may lead to a break up of the beam into a diffracting background profile with a monotonously compressing filament, that collapses in a finite distance of propagation. Thus, the width of the filament goes to zero and the intensity becomes infinite whereas the beam width in the RMS sense still increases. The simulations for the present case of Bessel-Gauss beams show that when the amplitude is increased further above $`A_0=1`$, the small second peak of Fig. 7 will start to dominate and eventually the simulations indicate a collapse of the second peak, although the RMS width still tends to infinity. In fact, even if the RMS width remains constant, the beam should still be able to undergo partial collapse.

Much effort has been devoted to the study of two-dimensional collapse phenomena induced by the Kerr nonlinearity, see e.g.  and references therein. In particular, it has been found that the virial theorem poses a sufficient but not necessary condition for the occurrence of a singularity where the amplitude becomes infinite. Thus the appearance of a partial collapse singularity below the threshold for a global collapse, as predicted by the virial theorem, is in accordance with earlier results.

## 6 Conclusions

Based on the nonlinear Schrödinger equation in cylindrical geometry we have studied the modification of the diffraction-less linear Bessel beams caused by the nonlinear Kerr effect. The stationary as well as the dynamic properties of the solutions to this equation have been analysed and both analytical and numerical techniques have been used. The investigation shows that the nonlinearity primarily affects the main and inner lobes of the Bessel beams. In the case of the stationary solutions, the central region of the nonlinear Bessel beam tends to become radially more narrow or more extended depending on whether the nonlinearity is focusing or defocusing, respectively. Asymptotically the solutions are still of the same oscillating form as the diffraction-less Bessel beams, the only remaining feature of the nonlinearity being a phase shift as compared to the linear case. However, in the case of the defocusing nonlinearity, there is a finite amplitude threshold for well-behaved solutions to exist. Above this limit the nonlinear diffraction effect becomes larger than the linear effect, which counteracts the diffraction, and no solutions are possible which vanish at infinity.

The properties of Gaussian-truncated Bessel beams have also been studied in the presence of the Kerr nonlinearity. It has been shown, using the virial theorem, that a non-diffracting situation in the RMS sense is possible to obtain by balancing nonlinear focusing and linear diffraction. However, this situation does not correspond to a stationary case of the beam profile. Significant redistribution of the beam occurs and using numerical simulations, the dynamic interaction between linear diffraction and nonlinear focusing has been analysed for varying degrees of nonlinearity. It has been found that, in particular the central parts of the beam may become significantly distorted and may even partially collapse even though the beam width, defined in the RMS sense, remains constant or even increases. This result is in agreement with the classical picture of dynamic self focusing of two-dimensional beams in nonlinear Kerr media.


text:
# Quantifying the Bull’s Eye Effect

## 1 Introduction

Distortions in redshift maps of galaxy positions have been known for decades. Important work was done on linear redshift distortions during the mid 1980’s (Kaiser, 1987), and more recently work has been done toward utilizing the effect to measure the linear redshift distortion parameter $`\beta \mathrm{\Omega }_{0}^{}{}_{}{}^{0.6}/b`$ (Hamilton, 1998; Peacock et al., 2001). This measure can yield an estimate of the mass density, $`\mathrm{\Omega }_0`$, if one has independent knowledge of the bias parameter, $`b`$. In addition, Hui, Kofman & Shandarin (2000) have explored the probability distribution function of density in relation to redshift distortions.

It has been noted that one consequence of redshift distortions, known as the Bull’s Eye effect (Praton & Schneider, 1994; Praton et al., 1997), depends upon the value of the mass density and may be independent of $`b`$ (Melott et al., 1998). As pointed out by Bridle et al. (2003), $`\mathrm{\Omega }_0`$ is a critical parameter in testing our cosmological models but is currently not known with high certainty. Measuring the strength of the Bull’s Eye effect may give an alternative estimate of this important cosmological parameter without requiring knowledge of the value of the bias parameter. The only data this approach requires is redshifts and so we plan to utilize data from galaxy redshift surveys such as 2dF (Colless et al., 2001) and Sloan Digital Sky Survey (SDSS) (Stoughton et al., 2002).

Redshift distortions occur because galaxies have local (peculiar) velocities which either add to or subtract from their line of sight redshifts due to the general Hubble expansion. The Bull’s Eye effect is composed of two effects of redshift distortion: small scale stretching (known as “Fingers of God”) and large scale compression, both along the line of sight. Stretching is caused by the small scale random velocities of galaxies near the center of collapsed clusters. Here, some galaxies are moving toward the observer (decreasing their redshift) and others are moving away (increasing their redshift), so the net effect is to spread the distribution out in the line of sight direction. Compression is the result of large scale coherent motion such as infall of galaxies toward the center of a cluster or wall. Here, galaxies on the near side of the cluster are moving away, while those on the far side are moving toward the observer, so the net effect is to compress the distribution along the line of sight. For more on both of these effects, see the thorough review of redshift distortions by Hamilton (1998), as well as Hui, Kofman & Shandarin (2000).

The combination of these two components results in a visual effect in which the observer sees some structures pointing directly at him, and others that seem to encircle him; this is the origin of the name “Bull’s Eye.” Figure 1 gives a schematic depiction of the two components of the effect (based upon a similar figure in Hamilton (1998)). Figure 2 shows how the effect appears in redshift space maps of galaxy positions. These plots were made using a random subset of points from a two dimensional slice of the full three dimensional simulations (see Section 2 for details on the simulations). The number of points was chosen to give approximately the same density as that of $`L_{}`$ galaxies in the Universe. Each plot contains approximately the same number of points. The biased plot appears more sparse due to greater clustering. Smoothing reduces this difference between biased and unbiased by removing points in low density areas (voids), as may be seen in Fig. 3.

The strength of the compression part of the Bull’s Eye effect depends on the value of $`\mathrm{\Omega }_0`$. In fact, for our purposes here the stretching part is noise. See Fig. 2 for an example of the effect with differing values of $`\mathrm{\Omega }_0`$. The connection between the strength of the effect and the value of $`\mathrm{\Omega }_0`$ can be understood by using the Zel’dovich approximation (Zel’dovich, 1970). This approximation relates the final position, $`\stackrel{}{r}`$, of a particle to its initial position, $`\stackrel{}{q}`$:

$$\stackrel{}{r}=a(t)[\stackrel{}{q}D(t)_\stackrel{}{q}\mathrm{\Phi }(\stackrel{}{q})]$$
(1)
where $`D(t)`$ is the growing mode solution for density perturbations, $`a(t)`$ is the cosmic scale factor, and $`\mathrm{\Phi }`$ is the primordial gravitational potential. One can relate $`D`$ to $`\mathrm{\Omega }_0`$ through

$$f=\frac{d\mathrm{log}D}{d\mathrm{log}a}$$
(2)
where $`f`$ can be approximated as $`\mathrm{\Omega }_0^{0.6}`$ (Peebles, 1980).

The velocity $`\stackrel{}{v}`$ of an element can be found by differentiating Eq. 1, which yields:

$$\stackrel{}{v}=H\stackrel{}{r}a(t)\dot{D}_\stackrel{}{q}\mathrm{\Phi }(\stackrel{}{q})$$
(3)
This velocity can be translated into an effective distance in redshift space

$$d_z=v/H=r_zfa(t)D(t)\frac{}{z}\mathrm{\Phi }(\stackrel{}{q})=aq_z(1+f)a(t)D(t)\frac{}{z}\mathrm{\Phi }(\stackrel{}{q})$$
(4)
Where $`z`$ denotes the redshift (line of sight) direction. Hence the displacement exactly along the line of sight is multiplied by a factor $`(1+f)`$. Since $`f\mathrm{\Omega }^{0.6}`$, it is apparent that the compression (along the line of sight) depends upon $`\mathrm{\Omega }_0`$.

Our ultimate goal is to refine a method which can determine numerically the strength of the redshift distortion associated with a particular value of $`\mathrm{\Omega }_0`$. This method will then be applied to galaxy redshift survey data in order to determine the value of $`\mathrm{\Omega }_0`$ for our universe. Here we use N-body simulations rather than surveys in order to develop methods which probe the effect and can then be modified and applied to surveys.

Throughout the paper we will use $`\mathrm{\Omega }_0`$ to denote the mass density in baryons plus cold dark matter at a redshift of zero. The effect of a cosmological constant or “dark energy” component is negligible for this effect. Also note that here “high” $`\mathrm{\Omega }_0`$ indicates the critical value, while “low” $`\mathrm{\Omega }_00.3`$.

We present two methods and their results in Section 2 and discuss the results in Section 3.

## 2 Methods and Results

We have developed two methods of quantifying the Bull’s Eye effect using N-body simulations. Each method is used to show two things: first, that we can measure the effect (for a given model); second, that we can distinguish between different values of $`\mathrm{\Omega }_0`$ using that measure of the effect. It should be emphasized that the purpose of this study is to work out methods which can probe the effect. Therefore, we have used simulations, wherein there is less uncertainty and parameters may be easily varied in order to test our methods for a variety of possible cases. Future work will utilize knowledge gained through this study to measure the effect using actual redshift survey data.

The simulations used here are numerical models for the gravitational dynamics of collisionless particles in an expanding background. All of the simulations used are done with a particle-mesh (PM) code with $`256^3`$ particles in an equal number of grid points (Melott, 1986; Melott et al., 1988). The physical size of the simulation box is $`512Mpc`$ on a side (for the value of $`h`$ given below). The simulations were normalized to an amplitude $`\sigma _8=0.93`$ at redshift moment $`z=0`$, and assume a Harrison-Zel’dovich primordial power spectrum. The typical correlation length is about three cells ($`6Mpc`$). All models were interpreted with an assumed Hubble constant $`h=H_0/100kms^1Mpc^1=\frac{2}{3}`$. When power spectra are parameterized in $`Mpc`$, the shape is dependent upon $`\mathrm{\Gamma }h`$, where $`\mathrm{\Gamma }=\mathrm{\Omega }_0h`$ is a formal shape parameter in the power spectrum formula of Bardeen, Bond, Kaiser, & Szalay (1986). Here, we decouple $`\mathrm{\Gamma }`$ and $`\mathrm{\Omega }_0`$ since it is known with reasonable certainty that the shape of the power spectrum is not too far from $`\mathrm{\Omega }_0h0.2`$ (Percival et al., 2002). Our intent in doing so is to assess the value of the Bull’s Eye effect in measuring $`\mathrm{\Omega }_0`$.

We do not wish to test our method only against currently favored cosmologies. We ran a simulation with $`\mathrm{\Omega }_0=1`$ and $`\mathrm{\Gamma }h=0.15`$; such models have been called $`\tau `$CDM in the past. This is our “high $`\mathrm{\Omega }_0`$” model. We also ran a simulation with $`\mathrm{\Omega }_0=0.34`$, $`\mathrm{\Gamma }h=0.15`$, which is most consistent with a variety of findings at this time. This is our “low $`\mathrm{\Omega }_0`$” model. There are, in addition, a variety of alternative models with the cosmological constant, $`\mathrm{\Lambda }`$, all of which have very small and totally linear effects on large-scale velocities, especially for $`\mathrm{\Omega }_00.3`$ (Hamilton, 2001). Therefore, we omit this in favor of wider exploration of parameter shifts with large effects.

For each simulation, redshift space versions were constructed, where the redshift distortion is included along one of the three orthogonal directions. In addition, we have versions where bias is included. There are many ways to produce a biased numerical simulation. In this case, particles were tagged if, in the initial conditions, they lay in regions of density enhanced 1$`\sigma `$ above the mean on a scale of R = 1 $`Mpc`$ (Cloud in Cell smoothing (Hockney & Eastwood, 1981)). Using only these particles to define the final state defines the biased density field. The bias factor $`b`$ = 1.8 for this procedure, where $`b`$ is defined as the ratio of $`\sigma _8`$ in the biased density field to that in the underlying simulation from which it was constructed. While the value of the bias parameter in our universe is not known, our $`b`$ is comfortably larger than any indicated by current data so that we can put a safe upper limit on the extent to which this is likely to be a problem for our method.

For our analysis, we used three realizations of each simulation (differing only in their initial random number seed). Each of these realizations has three redshift space versions (with the redshift distortion along one of the three orthogonal directions) and one real space version. The analysis is performed using two dimensional slices of the three dimensional simulation density field, in Cartesian coordinates. A slice is one grid cell thick and is 256 cells on a side. From each redshift space simulation box we take 256 slices containing the redshift distortion. Specifically, for redshift along the x and y directions, we take slices at constant z (x-y planes); for redshift along the z direction we take slices at constant x (y-z planes). For the real space boxes we take 256 slices along each direction. Therefore, we have (3 realizations)x(3 redshift directions)x(256) redshift slices and (3 realizations)x(3 directions)x(256) real space slices available.

From this set of slices we compile an ensemble to use for the analysis. This ensemble is obtained by keeping only every fourth slice from each simulation box in order to reduce correlation between slices due to proximity. This yields an ensemble of 3x3x(256/4)=576 slices (for a given model), which, when taken together, comprises slightly less volume than that of the 2dF survey (Colless et al., 2001).

Slices of the density field have been used since the surveys we will be applying our methods to (such as 2dF and SDSS) include thin slice regions. This also provides a large number of realizations to our ensemble. Cartesian coordinates have been used in order to facilitate development of a method which can then be adapted to apply to polar data such as that from redshift surveys. A Cartesian approach can, in principle, be used directly for a deep enough survey where a plane parallel approximation may be used and may lead to a statistic which is usable in polar coordinates. Again we emphasize that our purpose here is to work out a method in a simpler situation which may then lead to a method suitable for use with real data.

Before any analysis is done, the density slices are smoothed by convolution with a Gaussian $`e^{r^2/2\lambda ^2}`$. This smoothing serves two purposes. First, it ensures that the density does not depend upon grid size in the simulations, since the grid is smaller than the typical smoothing length. Second, smoothing helps to compensate for greater noise in real data as compared to simulations, due to a lower density of galaxies in real data. Typical galaxy separation in the universe is $`3Mpc/h`$ ($`4.5Mpc`$ for $`h=2/3`$) while in our simulations the separation is $`2Mpc`$. The vast majority of galaxies are clustered, as evidenced by the fact that most galaxies have a nearest neighbor closer than the mean separation (Ryden & Turner, 1984; Linder, et al., 1996). Therefore, smoothing emulates the real (clustered) galaxy distribution better.

The smoothing length, $`\lambda `$, is set to be a fraction of the correlation length, $`r_0`$, of the slice. The correlation length is used to set the smoothing length because it is stable and easy to measure. All the data presented here are for a fraction of $`1/2`$. We note that the presence of bias changes the correlation length and so will imply an effect on the smoothing length. Therefore, investigation of variability of the smoothing length is accounted for by examining the biased simulations.

### 2.1 Path length statistics

Our first approach involves using contour crossings of the 2D slices of the 3D density field. Contour plots of a typical slice are shown in Fig. 3. Real and redshift space (with redshift in the vertical direction) are shown for comparison, with high $`\mathrm{\Omega }_0`$, biased high $`\mathrm{\Omega }_0`$ and low $`\mathrm{\Omega }_0`$. In these plots the contours are made for a fixed value of the filling factor, $`\varphi `$, which is defined as the fraction of the total number of sites which have densities above a given threshold. That is, the actual density level used varies between slices and $`\mathrm{\Omega }_0`$ values for a specific value of $`\varphi `$. Using filling factor instead of density threshold reduces our sensitivity to bias (Melott et al., 1988). Biasing increases the amplitude of fluctuations. Therefore, using a fixed fraction of the density above a threshold tends to normalize biased to unbiased density fields (as long as biased density increases monotonically with unbiased density). This can be seen in Fig. 3 by comparing the biased and unbiased contours which appear similar here. Note, however, that biasing also increases the correlation length and therefore (since we use a fraction of $`r_0`$ as the smoothing length) the biased density field is effectively smoothed at a larger scale than the unbiased. This is also apparent in Fig. 3.

The analysis is done by looking at the contours at a given filling factor and marking out paths (in a certain direction; ie., x or y) between contour crossings (Ryden, 1988; Ryden et al., 1989). Various different paths can be marked, in particular those measuring structures and voids. Structure paths are marked between a contour crossing where the density increases from below the threshold to above (an up-crossing) and a crossing where the density decreases from above the threshold to below (a down-crossing). Hence, a structure path is an “up-to-down” path. Conversely, void paths are marked between a down-crossing and an up-crossing; they are “down-to-up” paths. See Fig. 4 for examples of both types of paths.

Our analysis of these paths proceeds as follows. We start with individual path lengths in a given slice (at a given filling factor), $`l_{i,j}^T`$ and $`l_{i,j}^L`$, where superscripts denote path direction ($`T`$ = transverse, $`L`$ = line of sight), subscript $`i`$ denotes a path and subscript $`j`$ denotes a slice. We then calculate the mean of these path lengths in a single slice,

$$\overline{l}_j^T=\frac{1}{n_j^T}\underset{i}{}l_{i,j}^T,\overline{l}_j^L=\frac{1}{n_j^L}\underset{i}{}l_{i,j}^L,$$
(5)
where $`n_j^T`$ and $`n_j^L`$ are the number of paths in each direction for slice $`j`$ (ranges from about 70 at $`\varphi =0.01`$ to about 500 at $`\varphi =0.5`$). In addition to the mean, we have investigated the dispersion and rms of path lengths as possible probes of the effect. The mean was found to give the strongest signal.

Once the mean of path lengths in a slice is found, we calculate for each slice the ratio of mean path lengths in the line of sight to transverse directions,

$$r_j=\frac{\overline{l}_j^L}{\overline{l}_j^T}.$$
(6)
We then find the median value of this ratio for the ensemble of 576 slices previously described. (The median is used because the distribution of ratios is non-Gaussian.) For use with this method, we record only paths from every fourth column (in both directions) in order to reduce correlation between paths within a slice. For most slices this distance (four cells) is similar to or slightly larger than the correlation length, which typically is between three and four cells.

Once the median of the ratios of mean path lengths has been calculated, it is plotted against filling factor, $`\varphi `$. This type of plot is shown in Fig. 5 for real and redshift space, for both structure and void paths. In these plots, the error bars (denoted by thin lines) correspond to 68% confidence intervals of the median. This is the confidence in the median for a measurement including a single slice. For the volume (all slices taken together), the size of the errors are reduced by a factor of $`\sqrt{n}`$, where $`n`$ is the number of slices. For the void paths, empty columns in the simulation box present a problem, since they are essentially paths of indeterminate length (due to periodic boundary conditions). We have, therefore, chosen to assign empty columns a length of twice the box size. This gives us paths which are larger than the box (as they must be), but which are not unrealistically large. Inclusion of empty columns in this way has very little effect on the results.

A general characteristic of the redshift space plots in Fig. 5 (for both void and structure paths) is that at low filling factor the ratio is greater than one and at higher filling factors the ratio is less than one. (The $`\varphi `$ value where this transition occurs varies for the different models.) For the void paths, a ratio greater than one indicates that the mean of void size is greater in the line of sight direction (as was also found by Ryden & Melott (1996) and Schmidt et al. (2001)) as compared with the transverse direction. This is the result of the distance between widely separated structures becoming larger due to large scale coherent motion, which leads to an increase in the size of large voids (which are the majority at small filling factors). This, in turn, results in the shrinking or disappearance of small voids from their compression by larger voids (Melott, 1983). At higher filling factors, where the ratio dips below one, the majority of voids are small and hence the shrinking of small voids results in a shorter mean path length in the line of sight direction compared with that in the transverse direction.

For the structure paths, a ratio greater than one indicates that the mean of structure size is greater in the line of sight direction as compared with the transverse direction. This is caused by the Fingers of God effect due to small scale random motion, which lengthens the structures in the line of sight. At higher filling factors, where the ratio dips below one, the structures are larger and the compression effect dominates (since it is a large scale effect), compared to low filling factor where the structures are smaller and hence the small scale effect (Fingers of God) dominates.

We note that comparison of the redshift and real space plots in Fig. 5 clearly shows that this method detects the presence of the Bull’s Eye effect in redshift space. Our expectation of no detection of the effect in real space is confirmed since the ratios are at or near 1 (in contrast to redshift space ratios) for all models.

In this study, we wish to probe the large scale, coherent motion (which leads to compression of structures) since it is primarily this effect which, by physical reasoning, should be dependent upon the value of $`\mathrm{\Omega }_0`$, as described in Section 1. Conversely, small scale, random motion (which produces Fingers of God) is noise for our purposes. Therefore, further discussion of the analysis of path lengths will focus on the void paths since, as noted above, these paths are more sensitive to compression while structure paths are more strongly affected by Fingers of God.

Before we proceed we must comment on the peak which appears in the upper left panel of Fig. 5 (void paths in redshift space). The turn-over which occurs around filling factor 0.03 appears to be unphysical. We believe this to be so for two reasons: 1) Extremely low filling factor values correspond to a small percentage of the density above the threshold, yielding few structures to use as references for voids, leading to a small number of void-paths. 2) The scarcity of structures, as well as their small size, means that the majority of void paths present are large compared to the simulation box size (greater than half the box in length). For instance, at $`\varphi =0.01`$, just to the left of the peak, a typical slice has about 70 paths, 90% of which are greater than half the box. At the peak ($`\varphi =0.03`$), a typical slice has about 100 paths, 60% of which are large. At $`\varphi =0.05`$, just to the right of the peak, there are about 120 paths, 50% of which are large. Moving further to the right, at $`\varphi =0.15`$ a typical slice has about 280 paths, 10% of which are large.

It is apparent, then, that at small filling factors (less than about 0.03) any difference in path length between the line of sight and transverse directions is too small of a fraction of the total length to be significant. This conclusion is strengthened by the fact that there is no turn-over at low filling factors for the structure paths, where small $`\varphi `$ yields few paths, none of which are long (since structures are few and small). We conclude, therefore, that the turn over in the upper left panel of Fig. 5 is unrealistic and our analysis will only be applied to filling factors greater than 0.03.

In addition to detecting its presence, we wish to use the strength of the Bull’s Eye effect to estimate the value of the mass density. Therefore, we must determine whether the method at least distinguishes between widely separated values of $`\mathrm{\Omega }_0`$. Further, we look for the method to have a low sensitivity to bias. For a certain range of filling factor in the redshift space, void paths plot in Fig. 5 (upper left panel), it is apparent that high and low $`\mathrm{\Omega }_0`$ can be distinguished. In addition, in the range $`\varphi 0.1`$ to $`\varphi 0.15`$ biased and unbiased high $`\mathrm{\Omega }_0`$ are nearly identical. Therefore, in the filling factor range $`\varphi 0.05`$ to $`\varphi 0.2`$ the method both distinguishes between widely separated values of $`\mathrm{\Omega }_0`$ and is reasonably insensitive to bias. We have tested the method using a variety of power spectra and the quoted $`\varphi `$ range is consistently best.

In order to more precisely determine how well the method distinguishes between $`\mathrm{\Omega }_0`$ values, as well as how sensitive it is to bias, we have applied a Wilcoxon rank-sum test to the ensembles of ratios of mean path length. This test was used because it determines whether one set of data is systematically *not greater* in magnitude than another (for details, see Bluman (2001) and Lehmann (1998)). In this case, we expect the ratios of mean path length to be systematically larger for high $`\mathrm{\Omega }_0`$ (both biased and unbiased) compared to low $`\mathrm{\Omega }_0`$. As a test for sensitivity to bias, we compare between biased and unbiased high $`\mathrm{\Omega }_0`$.

We have applied this test to void paths in redshift space, at a given filling factor. We use the same ensembles of ratios used to produce the median plot in the upper left panel of Fig. 5. This gives us 576 ratios (for a given model), produced from every fourth slice of our simulation boxes. In choosing a value of the filling factor at which to apply the test, we must balance the desire for a strong signal with the need to be insensitive to bias. Therefore, we take $`\varphi 0.1`$ as a good choice since the signal is relatively large here and both the biased and unbiased high $`\mathrm{\Omega }_0`$ give similar results. This is true for both void and structure paths, indicating reliability of this choice of filling factor value.

In Table 1 we present results of comparing between values of $`\mathrm{\Omega }_0`$ for void paths in redshift space at filling factor 0.1 using the Wilcoxon rank-sum test. That is, we test to see how well the method distinguishes between models. For each comparison we give the probability that the ratios of mean path length for model A are *not* greater than those for model B. Therefore, small P indicates that the ratios for A *are* greater than those for B. We see, then, that the ratios for high $`\mathrm{\Omega }_0`$ (both biased and unbiased) are greater than for low $`\mathrm{\Omega }_0`$ with confidence approaching 100%. We also note that sensitivity to bias is not great, as there is only 88% confidence that the ratios for unbiased high $`\mathrm{\Omega }_0`$ are greater than those for biased high $`\mathrm{\Omega }_0`$.

As an additional check on these results, we have applied a Kolmogorov-Smirnov (K-S) test to the same data. The K-S test determines whether two sets of data are drawn from the same distribution (for details, see Lehmann (1998)). In this case, we find that low $`\mathrm{\Omega }_0`$ *cannot* be distinguished from both biased and unbiased high $`\mathrm{\Omega }_0`$ with less than $`10^5`$ % confidence. That is, low *can* be distinguished from high (both biased and unbiased) with confidence approaching 100%. Also, biased and unbiased high $`\mathrm{\Omega }_0`$ *cannot* be distinguished with about 33% confidence. That is, they *can* be distinguished with only about 67% confidence. These results are consistent with those described above for the Wilcoxon rank-sum test.

We have shown, therefore, that our path lengths method both detects the presence of the Bull’s Eye effect (as seen in Fig. 5) and distinguishes between values of $`\mathrm{\Omega }_0`$ with high confidence. In addition, it is reasonably insensitive to bias.

### 2.2 Percolation

A second method we have employed utilizes a percolation routine developed by Dominik & Shandarin (1992). Given a density threshold the routine denotes sites above the threshold as occupied and then checks to see whether each occupied site is part of a cluster. The criterion for a site becoming part of a cluster is that it be adjacent to another occupied site. Diagonal sites may also be connected if the mean density of four closest sites is above the threshold. The routine then finds the largest cluster and checks whether it spans the density slice in each (orthogonal) direction. The threshold is systematically varied to determine the highest density level at which the largest cluster spans in each direction. For our analysis, we then convert this threshold to a filling factor value. As discussed in Section 2.1, using filling factor reduces sensitivity to bias, which is true for this analysis as well. For a real space density field, the filling factors at which the largest cluster spans the slice in each orthogonal direction should be close together. This is due to the structures being isotropic. On the other hand, for a redshift space field the filling factors should be farther apart. This is because isotropy is destroyed by redshift distortions. Largest cluster spanning in redshift space occurs at a higher filling factor in the line of sight direction than in the transverse direction. The difference in filling factor depends on the strength of the redshift distortion and therefore on $`\mathrm{\Omega }_0`$.

Our analysis utilizes the difference in filling factor at which the largest cluster spans the density slice in the line of sight and transverse directions. Specifically, we take $`\mathrm{\Delta }\varphi =\varphi _L\varphi _T`$, where $`L`$ = line of sight and $`T`$ = transverse. From the ensemble of slices, we assemble lists of $`\mathrm{\Delta }\varphi `$ corresponding to different models in both real space and redshift space. We wish to make two comparisons between these lists of filling factor differences. First, to show a detection of the effect, we compare between $`\mathrm{\Delta }\varphi `$ in real space and in redshift space for a given model. Second, to determine how well the method distinguishes between values of $`\mathrm{\Omega }_0`$ (and how sensitive it is to bias), we compare $`\mathrm{\Delta }\varphi `$ in redshift space for different models.

In order to quantify these comparisons we have applied a Wilcoxon rank-sum test to the lists of filling factor differences. As in Section 2.1, this test was used because it determines whether one set of data is systematically *not greater* in magnitude than another. In this case, we expect the filling factor differences to be larger for redshift space compared to real space and for high $`\mathrm{\Omega }_0`$ compared to low $`\mathrm{\Omega }_0`$. This method is applied to the same ensemble of 576 slices as the previous method.

In Table 2 we present results comparing filling factor differences between real space and redshift space using the Wilcoxon rank-sum test. That is, we test to see how well the method detects the Bull’s Eye effect. We perform the test for each model and give probabilities that the filling factor differences for redshift space are *not* greater than those for real space. Therefore, small P indicates that the $`\mathrm{\Delta }\varphi `$ for redshift space *are* greater than those for real space. For each model we have confidence approaching 100% that the $`\mathrm{\Delta }\varphi `$ for redshift space are systematically greater than for real space. From these results, it is apparent that we have a strong detection of the effect by this method. The test indicates that redshift and real space can be strongly distinguished. Sensitivity to bias in this case appears to be low since the results for biased and unbiased high $`\mathrm{\Omega }_0`$ are the same.

In Table 3 we present results comparing models in redshift space using the Wilcoxon rank-sum test. That is, we test to see how well the method distinguishes *between* models. For each comparison we give the probability that the filling factor differences for model A are *not* greater than those for model B. Therefore, small P indicates that the $`\mathrm{\Delta }\varphi `$ for A *are* greater than those for B. The $`\mathrm{\Delta }\varphi `$ for high $`\mathrm{\Omega }_0`$ are greater than low $`\mathrm{\Omega }_0`$ with 99.99997% confidence. Conversely, the $`\mathrm{\Delta }\varphi `$ for biased high $`\mathrm{\Omega }_0`$ are *not* greater than low $`\mathrm{\Omega }_0`$ with 99.95% confidence. This result is contrary to our expectation that high $`\mathrm{\Omega }_0`$ will have $`\mathrm{\Delta }\varphi `$ greater than those for low $`\mathrm{\Omega }_0`$, regardless of bias. Similarly, the $`\mathrm{\Delta }\varphi `$ for unbiased high $`\mathrm{\Omega }_0`$ are greater than those for biased high $`\mathrm{\Omega }_0`$ with confidence approaching 100%. Hence, this method exhibits a high sensitivity to bias.

As an additional check on these results, we have applied a K-S test to the same data. In this case, the test shows that the three models *cannot* be distinguished from each other with less than $`10^5`$ % confidence. That is, all the models *can* be distinguished from each other with confidence approaching 100%. The same is true for the comparison between real space and redshift space (for a given model). The K-S test simply determines whether two data sets can be distinguished. Hence, the Wilcoxon rank-sum test is more informative for our purpose since we wish to determine the relative magnitude of the filling factor differences for various models.

We have shown that our percolation method both detects the presence of the Bull’s Eye effect and distinguishes between high and low $`\mathrm{\Omega }_0`$ with high confidence. This method is, however, strongly sensitive to bias, as shown by the comparison between biased high $`\mathrm{\Omega }_0`$ and unbiased low $`\mathrm{\Omega }_0`$. Sensitivity to bias is further seen in the strong distinction between biased and unbiased high $`\mathrm{\Omega }_0`$. It is therefore apparent that while this method may distinguish between values of $`\mathrm{\Omega }_0`$, it can be deceived by the presence of bias.

## 3 Discussion

We have found two distinct methods which quantitatively show the effect of redshift distortions known as the Bull’s Eye effect. These methods are also capable of distinguishing between the effect for high and low values of $`\mathrm{\Omega }_0`$ with high confidence, and the path lengths method is reasonably insensitive to bias.

Our results may be compared to measurements of redshift distortions using the galaxy-galaxy correlation function. While the correlation function does measure the compression of the galaxy distribution along the line of sight in redshift space, we have employed different methods which measure the same effect. The correlation function is sensitive to the combination $`\beta \mathrm{\Omega }_{0}^{}{}_{}{}^{0.6}/b`$ (Peacock et al., 2001; Hoyle, et al., 2001), while we expect our approach to be sensitive to $`\mathrm{\Omega }_0`$ *independent* of $`b`$ (as shown in Section 1). This expectation is borne out for our path lengths method, where we see little sensitivity to bias. The difference here could be that the correlation function is only sensitive to the fluctuation amplitude (not phases) and hence is affected by bias (which affects the amplitude), while our path length method apparently detects more than just amplitude.

The results of the path lengths method applied to void paths indicates that large voids are elongated in the line of sight direction in redshift space (see Section 2.1 and Fig. 5). This effect was also noted by Ryden & Melott (1996) and Schmidt et al. (2001). However, in both of these papers the effect was found to be relatively minor. In comparison, our results show a significant detection of the elongation of large voids along the line of sight direction.

It should be noted that for the percolation method, the filling factors at which largest cluster spanning occurs are typically around 0.45. We have seen with the path lengths method that biased and unbiased high $`\mathrm{\Omega }_0`$ diverge at these high filling factors. Hence, it may be reasonable to expect that the percolation method will be sensitive to bias as we have seen.

Overall, we have shown that it is possible to quantify the Bull’s Eye effect in N-body simulations. This indicates that we should be able to do the same in redshift surveys and thereby estimate the value of the mass density, $`\mathrm{\Omega }_0`$. In addition, since the path lengths method is relatively insensitive to bias, we can therefore give an estimate of $`\mathrm{\Omega }_0`$ without knowledge of the bias parameter, $`b`$, unlike methods which use redshift distortions to measure the distortion parameter, $`\beta `$.

Future work will concentrate on adapting our methods to apply to redshift survey data and then using these adapted methods to estimate $`\mathrm{\Omega }_0`$ for our universe.

Thanks to E. Praton for helpful comments on the manuscript. We are grateful for supercomputing support from the National Center for Supercomputing Applications. This work was supported by NSF grant AST-0070702 and the Madison and Lila Self Graduate Fellowship at the University of Kansas.


text:
# Knot Probability for Self-Avoiding Loops on a Cubic Lattice

## Abstract

We investigate the probability for appearance of knots in self–avoiding loops (SALs) on a cubic lattice. A set of $`N`$–step loops is generated by attempting to combine pairs of $`N/2`$–step self–avoiding walks constructed by a dimerization method. We demonstrate that our method produces unbiased samples of SALs, and study the knot formation probability as a function of loop size. Our results corroborate the conclusions of Yao et. al. with loops generated by a Monte Carlo methodyao .

Knots and links naturally appear in long polymers catenanes , and play a prominent role in biological systems and processesfrisch . Examples include chromosomes during cell division alberts , knots in bacterial DNA ringclosure , or knots in the native states of proteinstaylor . It can be shown rigorously that very long self–avoiding loops are always knottedsumners . However, the theoretical proofs do not provide the frequency of the knots, or the functional dependence of their frequency on loop size. Quantitative insight into this question was first provided by a numerical study of random walks on a lattice volog . Excluded volume effects, or self–avoiding (SA) interactions, are certainly crucial for correct description of polymersdegennesSC ; however, their incorporation into numerical studies is not simple. Earlier studies considered continuum models of self–avoiding loop (SAL) polymers with varying degrees of self–repulsion koniar ; deguchi , and demonstrated that with increasing number of monomers $`N`$, the fraction of unknotted loops decreases as $`\mathrm{e}^{N/N_o}`$. The characteristic size at which knots appear is surprisingly large: It increases from several hundred steps in the absence of self–avoidance, to hundreds of thousands for strongly SA polymers. Since the value of $`N`$ used in typical simulations does not exceed several thousands, for SAL one can assume that the probability of an unknotted configuration simply decays exponentially. The value of $`N_o`$ can then be extracted by noting that for $`NN_o`$ the probability of the knotted configurations is $`P_NN/N_o`$. A recent study by Yao et al.yao  investigated the knotting probability of SALs on a cubic lattice with $`N3000`$ and found $`N_o2.5\times 10^5`$. Our results corroborate the above study using a different approach to generating SALs.

Generating sufficiently large numbers of SALs has been the main obstacle to the study of statistics of knots in polymers. Static methods for creating SA polymers one at a time have the advantage of producing configurations that are independent of each other. If we do not require the two ends of a polymer to meet, several methods generate samples of “properly weighted” configurations: For lattice or for continuum models with “hard” potentials, this means that every configuration has the same weight. In the dimerization algorithmdimerization , an $`N`$-step SA walk (SAW) is created by generating two $`(N/2)`$-step SAWs and attempting to concatenate them. (If the concatenated walk is self–intersecting, it is discarded and the pair creation process is restarted.) The resulting SAWs are properly weighted, i.e. each has the same weight. While this method is very efficientmadrassokal , it is not well suited for generating SALs: To produce a properly weighted $`2N`$–step SAL, we could first generate two $`N`$–step SAWs, assume that they both start at the origin, and check whether they do not intersect and end at the same point creating a loop, discarding the pair if they do not form a proper SAL. However, the probability of two SAWs in $`d`$ dimensions accidentally terminating at the same point is proportional to $`R^dN^{d\nu }`$, where $`R^2`$ is the mean squared end-to-end distance of $`N`$-step SAWs, and in $`d=3`$ the swelling exponent is $`\nu =0.588`$ madrasslade . Full circles in the Fig. 1 depict this decay with $`N`$. In addition, the probability that two such walks do not intersect each other decays as $`N^{2(1\gamma )}`$, where $`\gamma =1.158`$ in $`d=3`$cara , as demonstrated by the diamonds in Fig. 1. Thus, successful formation of a three-dimensional SAL using this direct method is proportional to $`N^{2.08}`$, as depicted by open circles in Fig. 1.

To circumvent the problem of the rarity of cases where two SAWs starting from the origin terminate at the same point, Chen suggested chen  accumulating lists of SAWs. Whenever two SAWs in the list have the same end-to-end distance, they are rotated to make their endpoints coincide, thus producing a loop. Keeping the information on large numbers of SAWs is computationally memory–intensive. A more serious concern is that the ensemble of generated loops is biased: The solid line in Fig. 2 depicts the probability distribution of the end-to-end distance $`r`$ for 256-step SAWs; the dashed line corresponds to the subset of cases where a pair of SAWs has the same end–point. The latter is the correct ensemble for loops, while the former is the ensemble produced when accumulated SAWs are rotated and linked. This distinction in weight (bias) is well-known, and is properly corrected in Chen’s algorithm. However, as in all bias-corrected methods, while the expectation value of a desired quantity is correctly reproduced, the variance of this quantity increases with system size. Essentially, bias–corrected methods attempt to reconstruct one distribution from a tail of another distribution; the distributions moving further apart with size. Fortunately, as can be seen from Fig. 2, for moderate values of $`N`$ the distributions are not very different, and good results were obtained in Refs. koniar ; deguchi  by this method. Alternatively, Yao et al.yao  used a Monte Carlo pivoting methodmos  to generate many SAL configurations starting from an initial loop on a cubic lattice. While this method produces correctly weighted configurations, they may be statistically dependent. Cognizant of this difficulty, Yao et al. sampled the system at time intervals significantly exceeding the decay time of geometrical correlations, and verified that the times of appearance of knots behave like in a Poisson process.

In this work we employ a direct unbiased approach to generating $`(2N)`$–step SALs. We first use successive dimerizations to generate two $`N`$–step SAWs, and then check whether (if starting from the same origin) they form a SAL. As discussed before, the probability of successful pairing of two SAWs decays as $`N^{2.08}`$; this includes the rejection probability that two segments of a formed loop intersect. We find that the latter probability has only a weak dependence on the size of the loop. Indeed, the distribution of the distance between the origin and the 128th monomer in 256–step SALs, depicted by full circles in Fig. 2, practically coincides with the probability distribution (dashed line) of generating loops from two SAW segments, irrespective of their mutual intersections.

The small probability of two SAWs forming a loop can be enhanced 48-fold, by taking advantage of symmetries of the lattice to consider only SAWs whose end–point coordinates $`\stackrel{}{r}=(x,y,z)`$ satisfy the relations $`xyz0`$. This is achieved by generating a regular $`N`$-step SAW, and then performing the following transformations: If the coordinate $`x`$ of the end–point is negative, the walk is reflected with respect to $`yz`$ plane; and similarly for the other end–point coordinates. If the end–point has $`x<y`$ then $`x`$ and $`y`$ coordinates are interchanged along the whole walk; with corresponding interchanges for other pairs of axes. As a result of these transformations, we reduce the space in which the end–point can be located by a factor of 48, and thus increase the probability of forming a loop. In continuum, such reduction of space does not introduce any bias, since each SAW in the allowed subspace corresponds to 48 SAWs in the original space. This is, however, not true for SAWs ending on the boundaries of the allowed subspace; e.g. there are only 24 SAWs with $`x=y>z>0`$, and consequently loops created by SAWs ending at such points are under–represented by a factor of 2. In continuum, configurations with $`x=y`$ have zero measure, but on discrete lattices a finite fraction of loops have this property, thus, creating a slightly biased sample. Since this is a boundary effect, it decreases with the inverse linear size of the space considered, i.e. it will be proportional to $`1/R`$.

The direct method is quite efficient for generating small loops. The minimal (trefoil) on a cubic lattice consists of 24 stepsdiao , and the probability of its occurrence is extremely small. For $`N=64`$, we had to examine 510 million pairs of 32–step SAWs in order to find 1.4 million SALs, out of which 27 formed a trefoil knot of the type depicted in Fig. 3. At this point the probability to have a knot is $`0.000022\pm 0.000004`$. We used the Alexander polynomialalexander ; alextheor  (at the value of its argument equal to -1) to determine the presence and type of a knot. Since almost all observed knots were trefoils, and very few more complicated knots were encountered, this invariant sufficed for our purposes. The limiting factor in our simulations was computer time, and we could not go beyond $`N=1024`$, for which several months of CPU on a multiple processor computer were necessary to reach sufficient accuracy. At this value of $`N`$, the probability to form SALs dropped to $`10^5`$, even after the 48–fold enhancement of the sampling. While the correction to bias in sampling symmetric configurations was important for $`N=64`$, it became negligible for $`N=1024`$.

Figure 4 depicts our results for the probability $`p_N`$ of having a knot in a SAL of length $`N`$. The curvature in the results for this range of values of $`N`$ prevents a reliable extraction of $`N_o`$ by a linear fit. Our two highest data points compare well with the lowest data points of Yao et al.yao  that are included for comparison. While our results give a slightly larger probability than in Ref. yao , the difference is only one standard deviation, and is not significant at this level of accuracy. We also note that out of the 59 knots detected for $`N=512`$, only two were not trefoils, consistent with the results of Ref. yao .

In conclusion, we used the most direct method for generating an unbiased ensemble of SALs on a cubic lattice (up to symmetry factors). The method does not require large memory, but is very time consuming. While quite efficient for small and moderate sized SALs, it is limited to loops of around 1000 steps, at which the probability of forming a knot is around half a percent. Our results are consistent with those of Ref. yao , and provide an independent support of the conclusions from Monte Carlo sampling.

###### Acknowledgements.

This work was supported by US-Israel Binational Science Foundation (grant 1999-007), and by the National Science Foundation (grant DMR-01-18213).