Klaus Mølmer

Yann Le Coq

Signe Seidelin

Dispersive coupling between light and a rare-earth ion doped mechanical resonator

Keywords: numbers: 42.50.Wk, 42.50.Ct., 76.30.Kg

By spectrally hole burning an inhomogeneously broadened ensemble of ions while applying a controlled perturbation, one can obtain spectral holes that are functionalized for maximum sensitivity to different perturbations. We propose to use such hole burnt structures for the dispersive optical interaction with rare-earth ion dopants whose frequencies are sensitive to crystal strain due to the bending motion of a crystal cantilever. A quantitative analysis shows that good optical sensitivity to the bending motion is obtained if a magnetic field gradient is applied across the crystal during hole burning, and that the resulting opto-mechanical coupling strength is sufficient for observing quantum features such as zero point vibrations.

The ability to control and probe physical systems at the quantum mechanical level has been demonstrated in a variety of examples, ranging from single photons and atoms, over currents and voltages in electronic circuits to the motion of mechanical devices. Many of these systems are promising candidates for sensitive and high precision measurements and for transmission, processing and storage of quantum information, but one single system is typically not adequate for all of these functions. This has spurred the interest in so-called hybrid system [1] which combine physical components which are separately optimized for different tasks. Outside the technical challenge of handling physically very different systems in a single laboratory experiment, the mismatch of their physical properties (resonance excitation frequencies, spatial overlap, and coherence time scales) presents a main obstacle against the efficient transfer of quantum states between them. For instance, single atoms have microscopic dipole moments and interact in general only weakly with physical observables of mesoscopic quantum systems which occupy orders of magnitude larger spatial volumes. One successful remedy to this weak coupling is to use ensembles of many atomic particles with correspondingly increased coherent coupling strength. This is, indeed, the rationale behind the use of atomic ensembles for optical interfaces and memories and of spin ensembles in conjunction with superconducting circuits [1]. In this article, we address the application of rare-earth ion doped crystals for hybrid quantum technologies. Such crystals have found rich applications in quantum communication protocols where their strong inhomogeneous broadening is favorable for speed and bandwidth. One particular hybrid technology which holds promise for an efficient coupling between radically different degrees of freedom relies on opto-mechanical interactions [2,3]. Major achieve-ments such as the ability to prepare a mechanical oscillator in the quantum ground state [4,5] have spurred ambitious goals to prepare non-classical states of motion and use such systems in precision measurements and quantum information applications. While mechanical oscillators can be coupled via light beams to other systems, such as atomic ensembles [6,7] we propose a simpler setup in which the bending motion of an inorganic crystal is coupled to the optical transitions in an ensemble of rareearth ion dopants in the same crystal. The crystal strain generated by the harmonic motion of a bulk mechanical oscillation provides an intrinsically stable coupling to the dopant optical properties as, contrary to other coupling mechanisms to externally imposed fields, it circumvents instabilities arising from drifts in position of field source and oscillator. Weak strain coupling to single emitters has been previously observed [8][9][10], but using a rareearth ion ensemble, we benefit from their strong collective coupling as well as their narrow homogeneous linewidths and long coherence times, allowing unambiguous observations of quantum features.

We suggest to use spectral hole burning to prepare a transmission window with no optical absorption and obtain the opto-mechanical coupling through the offresonant, dispersive interaction with ions with transition frequencies outside the hole. As a concrete application, we will focus on Eu 3+ ions in a Y 2 SiO 5 host matrix as they exhibit record-narrow optical transitions among solid-state emitters. Despite its long coherence time [11], the 580 nm optical transition 7 F 0 → 5 D 0 in Eu 3+ ions is sensitive to crystal strain [12], making it an attractive candidate system for strain-coupled optomechanics. We investigate a crystal in the shape of a micro-meter scale cantilever whose bending produces a local strain-induced perturbation to the ions resonant frequencies. During bending, there will be a strain gradient across the cantilever, ranging from tensile to compressive strain (see fig. 1), and ions with identical frequencies in the unbent crystal (a) experience different frequency shifts in a bent cantilever (b). In the latter case, ions located in three different layers (subject to identical strain within each layer) experience different strain and the emission lines shift accordingly: if the oscillator is bent upwards, emitters on the top face experience a compression, whereas emitters on the bottom face are subject to a tensile strain. As the cantilever vibrates, the collective line shape periodically broadens and narrows, and the motion should be readily detected by an optical probe.

However, if static local perturbations of the emitters lead to additional inhomogeneous broadening with a very wide absorption line shape, the strain induced frequency shifts cannot be probed dispersively. To circumvent this, we propose to prepare the ion ensemble by burning of narrow spectral holes with the cantilever at rest, and subsequently use their optical response to the cantilever motion. Spectral holes are formed by resonant optical excitation and decay of the ions into other long-lived (dark) states, typically other hyperfine state within the electronic ground state manifold. Due to the narrow homogenous linewidths, one can scan the excitation laser over a finite frequency interval and burn a spectral hole which must be narrower than the ground state hyperfine splitting to avoid repopulating the hole by subsequent excitation of the dark ions. A weak optical field with a frequency at the center of the hole burnt interval will subsequently interact dispersively with the ions with frequencies outside the spectral hole, and we will in the following propose a method to make this dispersive coupling sensitive to the cantilever motion. As indicated by the example in fig. 1, bending of the crystal shifts the transition frequency by an amount proportional to the strain experienced by the individual ions. This implies that, to lowest order, the weak probe field is not sensitive to the bending, unless we can ensure an asymmetric response of the ion frequency distribution. Such an asymmetry can be prepared by a particular hole burning protocol as discussed in the following sections.

For our analysis we assume a spatially homogeneous ion distribution, and we also consider the frequency distribution due to the inhomogeneous broadening to be constant within a wide detuning range. We assume that the bending of the beam shifts the resonance frequency of all ions by a quantity proportional to the strain at their position in the crystal. The bending of the cantilever resonator with thickness e and length L is represented by a single collective coordinate X, measured as the vertical displacement of the cantilever tip shown in fig. 1. The ions are all located near the anchoring of the cantilever with spatial coordinates (x i , y i , z i ), where x i is in the same direction as X and runs from -e/2 to +e/2, and z i is along the beam axis, and takes values from 0 to L. From the Euler-Bernouilli theory of beams, we know that the local strain of the matrix around ion i is proportional to both X, x i and L-z i . We optically address only atoms near the base of the cantilever, and we therefore neglect the z i dependence, so the strain induced frequency shift reads kx i X, with k a constant depending on the beam geometry and physical properties of the doped material.

To enhance the position sensitivity of the cantilever, we prepare our system by applying conventional spectral hole burning while the unbent crystalline cantilever is subject to a magnetic gradient [15]. We assume a convenient choice of light polarization and magnetic field directions with reference to the crystalline axes so that all ions are shifted by the linear Zeeman effect such that the overall transition frequency of the i th ion at position (x i , y i , z i ) can be written as

ν i = ν 0 + δ 0 i + kx i X + gx i ∇B = ν 0 + δ i , (1)
where ν 0 is the chosen central frequency for our hole burnt structure, δ 0 i is the intrinsic inhomogeneous frequency shift of the i th ions (without the effect of cantilever motion or magnetic field), ∇B is the externally applied magnetic gradient and g the linear Zeeman effect sensitivity. We finally also assume that δ 0 i and x i are uncorrelated.

Our functionalized hole burning procedure consists of two steps:

1) Apply a static magnetic gradient ∇B 0 and scan the hole burning laser across the spectral range [ν 0 -∆, ν 0 + 3∆] within the inhomogeneous linewidth. The value of ∆ can be chosen arbitrarily, but must fulfill the condition ∆ ≥ gx max i ∇B 0 . This step is represented in fig. 2a. 2) Apply the reverse static magnetic gradient -∇B 0 and scan the hole burning laser across the interval [ν 0 -3∆, ν 0 + ∆]. This step is represented in fig. 2b.

As shown in fig. 2, after the magnetic gradient has been switched off, we obtain an inhomogeneous absorption profile of the ions with a transmission hole around the central frequency ν 0 , and with a width that depends on the location x i in the cantilever but is on average close to 6∆. We are then ready to probe and interact with the system by transmission of a weak laser beam with frequency ν 0 , at the center of the hole. In fig. 2c, we show how the frequency shift kx i X due to bending of the cantilever transforms the shape of the transmission window. In particular, we note that our hole burning procedure ensures a non-vanishing phase shift of the probe beam due to changes in X, because contributions from ions with positive x i are not canceled by the ions with negative x i as they are further detuned.

At a finite temperature, corrections to the hole shape arise because the thermal Brownian motion of the cantilever amplitude during hole burning exposes the whole inhomogeneous profile to strain distortions like the one shown in fig. 2d. These distortions lead to excitation and decay into dark states of those ions that are brought into resonance with the hole burning laser during their random Brownian motion. We estimate the consequences of this mechanism by assuming that ions with frequencies within kx i X of the spectral hole in fig. 2a are also transferred to dark states, where X is of same order of magnitude as the root-mean-square thermal excursions of X. Fig. 2d shows how the spectral hole is modified by the thermal Brownian motion under hole burning and effectively excludes ions at vertical position x i that are detuned with respect to ν 0 by an amount

δ ∈ [ L (x i ), R (x i )], where L (x i ) = -3∆ + gx i ∇B 0 - (+)kx i X and R (x i ) = 3∆ -gx i ∇B 0 + (-)kx i X for pos- itive (negative) values of x i .
To describe the interaction between the probe laser field and the hole burnt ensemble, we use a semi-classical approach based on the Maxwell-Bloch equations, which, on the one hand, accounts for the state of the individual atoms by solution of the quantum optical master equation and, on the other hand, describe the damping and dispersion of the Maxwell electric field amplitude due to the interaction with the mean atomic dipole density. Due to the spectral hole, the laser is far off-resonant from all atomic transitions, and we are in the dispersive and perturbative regime of linear excitation of the atoms.

The cantilever is a mechanical resonator of eigenfrequency ν M , for which we assume a small modulation amplitude. Under this assumption, and imposing the extra condition that ν 2 M |δ i | 2 for adiabatic following (see appendix), the phase shift of the probe field with wavelength λ and cross section A interacting with a two-level system with detuning δ i and radiative decay rate Γ is given by

∆ϕ i = - λ 2 Γ 8πA 1 δ i (2)
per atom.

The phase change ∆ϕ i of the classical (coherent state) field amplitude is equivalent to a phase change of the same value occurring on each single photon due to the presence of the ion. The combined state of the ion and the stream of photons thus accumulates a quantum phase θ that changes at a rate θ = R∆ϕ i , where R = IA ω0 (ω 0 = 2πν 0 ) is the rate of photons passing through A with intensity I. The ion is left in its ground state and the photons leave with unchanged frequency, and the time evolving phase due to the light-matter interaction is equivalent to the light induced AC Stark shift of the ion due to the laserfield,

V i = -θ = -R∆ϕ i = σ 0 I δ i , (3)
where σ 0 = λ 3 Γ 16π 2 c . The expression for V i has, indeed, the familiar form I/δ i , valid for large detunings. We note that we applied a simple two-level description and to evaluate the accurate value of the optical phase shift and atomic energy shift, one must specify the actual transition, the light polarization, and also the index of refraction of the host material.

Since the ion transition frequencies are modified by the crystal strain, their AC Stark shifts become functions of X as well as of the optical intensity, i.e., we obtain the desired optomechanical coupling. The total energy shift of all ions -due to the dispersive light matter interaction, is the sum of Eq.( 3) over all ions.

To evaluate the dependence of the total interaction energy on the collective cantilever coordinate X, we integrate the contributions over the thickness e of the cantilever and the corresponding distribution of δ from all ions outside the spectral hole:

V (X) = nσ 0 I 0 -e/2 dx L (x)+kxX -∞ dδ δ + ∞ R (x)+kxX dδ δ + e/2 0 dx L (x)+kxX -∞ dδ δ + ∞ R (x)+kxX dδ δ . (4)
Here, n is the density of emitters per unit of frequency and unit vertical distance, and the interaction limits which correspond to the edges of the spectral holes are expressed in terms of the functions defined previously in the text. By applying a first order Taylor expansion of 1/δ in kxX, the integrals can be readily evaluated and we obtain Figure 2: Schematics of the result of the hole burning procedure on the ground state ion distribution as a function of the ion position and "intrinsic" inhomogeneous frequency detuning δi. The lines denote the border between the zones with large |δi| values where the ions are left in the ground state, and the central zones where they transferred to the dark state during burning. In the first step of the hole burning procedure (a), the hole is burned in the [ν0 -∆, ν0 + 3∆] interval, indicated with dotted lines, while the magnetic gradient ∇B0 is on. The solid lines show the hole after the magnetic gradient is turned off, realizing the "right" edge of the final hole. The second step (b) realizes the "left" edge by hole burning the interval [ν0 -3∆, ν0 + ∆] while the magnetic gradient -∇B0 is on. The solid lines show the final hole after the 2 steps are completed and the magnetic gradient is turned off. Panel (c) shows the asymmetric effect of bending after burning both edges (in black, the unbent cantilever, and in grey dashed lines, the cantilever bent by a positive amount X). Panel (d) shows the effect on the hole due to thermal Brownian motion during the burning procedure: dotted lines correspond to burning at zero temperature and solid lines to burning with finite Brownian motion of the cantilever.

V (X) = 2kXnσ 0 I g∇B 0 e (g∇B 0 ) 2 -(k X) 2 - 3∆ (g∇B 0 + k X) 2 ln 1 + (g∇B 0 + k X)e 6∆ + 3∆ (-g∇B 0 + k X) 2 ln 1 + (-g∇B 0 + k X)e 6∆ . (5)
Note that for temperatures below 100 mK during hole burning, the Brownian motion is negligible, k X ∆, and the expression for V (X) simplifies, and a third order Taylor expansion in g∇B 0 e/(6∆) yields

V (X) = -kXnσ 0 I g∇B 0 e 3 54∆ 2 . ( 6
)
The interaction strength given by this expression is explicitly proportional to X and I as desired for the opto-mechanical coupling. Note that it is also proportional to the bias magnetic gradient ∇B 0 , applied during hole burning. The parameter ∇B 0 thus serves as an (adjustable) gain coefficient for the optomechanical coupling. To obtain a maximum coupling strength, it is desirable to match the bias magnetic gradient to the width of the spectral hole gx max i ∇B 0 = ∆. The coupling of the lightfield with the emitters causes an AC Stark shift that depends on the bending of the cantilever and hence results in an displacement X disp of the resonator equilibrium position given by

X disp = - dV (X) dX | X=0 /mω 2 M , (7)
where ω M = 2πν M is the angular frequency of oscillation of the resonator.

Conversely, the coupling between the light field and the emitters produces a phase shift on the laser. The same proportionality applies between the optical phasehift and interaction energy for the whole system as for a single ion, and we can thus directly write the phase shift of the optical beam as

∆ϕ = i ∆ϕ i = V (X)ω 0 IA , (8)
where the dependence on the intensity I cancels due to the explicit proportionality between V (X) and I in the dispersive coupling regime. This phase-shift can be detected interferometrically, which provides a readout of the position of the resonator. In particular, for the light intensity that leads to the displacement expressed in Eq. 7, the corresponding phase shift of the carrier optical frequency is V (X disp )ω 0 /IA.

The vibratory movements, on the other hand, can be detected on the transmitted laser field by observing spectral sidebands at the mechanical oscillator frequency ±ν M from the carrier. Thermal excitation of the cantilever motion yields a total integrated power of sideband phase fluctuations [16] equal to dV (X) dX | X=0 Xω 0 /IA, which converges to the value dV (X) dX | X=0 X zpf ω 0 /IA at zero temperature. To detect such sidebands requires their amplitude to exceed the shot noise limit for the given optical power and integration time (measurement bandwidth). The resolution of the phase shift is thus ultimately limited by photon absorption and transfer of ions to other long lived states, which gradually "overburns" the spectral hole and destroys it. The transition line in Eu:YSO, however, has a very small natural linewidth Γ and therefore, even for high optical power, overburning will not occur before a typical time τ = 4π/Γ, which sets the limit for the integration time. Increasing optical power within the integration time τ can lead to arbitrarily high resolution. For practical consideration, we choose here to limit the maximal optical power to a few mW, based on the fact that a typical standard commercial close-cycle cryocooler exhibit near 3 K a heat load capacity around 10mK/mW near the sample. Consequently, an optical power above a few mW (part of which may be absorbed during burning or probing phase) may therefore lead to excess heat load and/or thermal instability. Further improvement of the resolution can be obtained by repeating the experiment (erasing, re-imprinting and probing the spectral hole) and averaging the successive results.

As an example we consider a single-clamped cantilever with the dimensions 100 × 10 × 10µm 3 interacting with a laser beam traversing the cantilever near its fixed end for maximum strain as illustrated in fig. 1. We consider a cantilever which consists of Y 2 SiO 5 (Young Modulus of 135 GP) with an effective mass m eff = 1.1 × 10 -11 kg, and of which the first excited mechanical mode vibrates at ν M = 890 kHz. The cantilever contains a 0.1 % doping of Eu 3+ ions, with a 7 F 0 → 5 D 0 transition centered at 580 nm and with Γ = 2π × 122 Hz (and Γ * typically ≤ 2π × 1 kHz) at T=3 K. We choose a power of 1 mW and a hole width of 6∆=6 MHz to fulfill the condition for large detuned adiabatic following, see appendix. To calculate the proportionality constant k we use the values from reference [12]: -211.4 Hz/Pa for crystal site 1, which is the most sensitive of the two non-equivalent sites [17]. The sensitivity to magnetic field is g = 3.8 kHz/G [14]. In order to maximize the coupling while fulfilling the condition ∆ ≤ g∇B 0 e we choose ∇B 0 = 530 T/m, which seems possible to achieve in a cryogenic environment over short distances (mm or less) using small superconducting coils.

In this configuration, the static displacement of the tip of the resonator due to the light field amounts to X disp =0.4 pm and the corresponding phase shift of the laser (the carrier) equals 0.2 mrad. This shift is easily observable as, for the 1 mW laser power, the shot noise limited phase resolution is 0.45 µrad within the allowed detection time, before hole-overburning becomes non-negligible (approximately 16 ms for the 122 Hz linewidth). For comparison, a direct reflection of a 1 mW laser on the resonator would give rise to a much smaller static displacement (20 fm).

The amplitude of the Brownian motion at 3 K is X=0.2 pm, and the spectral sidebands contain an integrated phase of 0.11 mrad due to this thermal excitation. For an integration time equal to the inverse of the thermal linewidth (25 µs), the shot-noise limited phase-resolution is 14 µrad. The thermally excited sidebands are therefore readily observed, even within such short integration time.

Moreover, by increasing the integration time up to the maximum before over-burning, it is possible to observe and measure accurately the detailed shape and size of the sidebands. Zero-point motion of the resonator, averaged over the measurement, induces a small excess integrated phase of 0.4 µrad in the sidebands. As this value is close to the phase resolution achieved within the maximum integration time before hole-overburning, the shot noise limited resolution is therefore sufficient to observe the effect of the zero-point motion of the mechanical resonator. By using either dilution fridge or an active laser cooling mechanism (or a combination of these), the temperature can be lowered near a point where the thermal excitation does not dwarf the effect of zero point fluctuations. For example, at 30 mK, the zero point fluctuation induces approximately 10 -3 relative excess integrated phase over the effect of Brownian fluctuations alone. Such a deviation seems measurable, provided sufficient knowledge of the relevant parameters (temperature, Q factor,...). Several measurements at different temperatures may also be used to estimate the various relevant parameters with the necessary accuracy. Note that the resolution can be further increased by repeating the full hole imprinting and measurement sequence several times, or use optical repumpers to preserve the spectral hole.

A potentially perturbing effect arises due to fluctuations of the laser power. As the static displacement corresponding to 1mW of laser power is approximately 0.4 pm, this laser power must be stable to within ∼ 10 -4 to ensure a perturbation much smaller than the zero point fluctuations (∼ 1 fm), a power stability requirement well within reach of standard stabilization techniques. In addition the laser frequency must also be carefully controlled. The zero-point motion induces a frequency shift of the ions closest to the edge of the resonator of approximately 37 Hz. The probe laser must therefore have a frequency stability substantially better than that, which is well within reach of nowadays commercially available ultra-high finesse Fabry-Perot cavity stabilized lasers (which typically have sub-Hz frequency stability). Note that probing the Brownian motion alone at 3 K exhibits a much less stringent frequency stability requirement at the sub-kHz level.

In summary, we have proposed a hybrid scheme which provides an efficient optomechanical strain based coupling between a mechanical resonator and an ensemble of narrow linewidth ions. The coupling mechanism via a functionalized spectrally burnt hole does not suffer from inhomogeneous broadening of the emitters, nor from the inhomogeneity of the strain across the resonator during displacement. It is even possible to further increase the coupling by exploiting not just one but an ensemble of functionalized spectral holes within the inhomogeneous linewidth, which can be obtained by using a comb of optical frequencies. The narrow linewidth and dispersive nature of the coupling allows use of a relatively large optical power without destruction of the functionalized hole and observation of the cantilever motion at the quantum level. Moreover, the large back-action of the dispersive coupling on the resonator (more than one order of magnitude superior to that of simple light pressure and comparable to the amplitude of the Brownian motion at 3 K) holds promise for efficient resonator manipulation such as the implementation of active cooling protocols or quan-tum engineering.

Figure 1 :

 1 Figure 1: Schematics of a cantilever anchored in the foreground of the drawing, and with the far end (a) in its equilibrium position, and (b) bent upwards. The colors on the cantilever indicate different values of the strain: green corresponds to a non-strained material, blue and red to compressive and tensile strain, respectively. The white dashed lines indicate planes of equal strain, neglecting variations along zi.The bottom graphs illustrate how a single emission line is shifted due to the strain (arbitrary units).

We thank Stefan Kröll, Philippe Goldner and Sébastien Bize for discussions, and SS thanks Daniel Estève for help and discussions. KM acknowledges support from the Villum Foundation, YLC from the Ville de Paris Emergence Program and SS from la Région Rhône-Alpes (CMIRA Explora-pro). This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 712721.

APPENDIX

The AC Stark shift of the atomic ground state level and the phase shift of the transmitted field are both governed by the response of the two-level atomic system to the laser field, and in particular, the mean optical coherence induced in the ions. We denote the lower and upper levels |g and |e , and we assume that we only excite the ions very weakly, i.e., the ground state population ρ gg = 1 to first order in the Rabi frequency Ω = DE/ , where D is the dipole matrix element and E denotes the driving field amplitude. This leads to the equation of motion for the optical coherence

where the decoherence rate Γ may include both radiative (∝ D 2 ) and non-radiative contributions, and where we have defined δ r = 2πδ as the detuning, measured in radians per second.

The steady state solution reads, ρ st eg = (iΩ/2)/(Γ/2iδ r )

-Ω/(2δ r ), and the resulting macroscopic polarization of the doped cantilever (with an ion density n) is then given by P = nDρ eg -n(Γ/δ r )E, leading to the optical phase shift and energy V ∝ P • E analyzed in the article.

The purpose of this appendix is to analyze the case where X performs oscillatory motion, X(t) = X 0 cos(ω M t), and the strain induced coupling leads to a temporally modulated detuning, δ r (t

). Due to the long lifetime of the ion coherence and excited state, we cannot merely assume that the coherence ρ eg attains the stationary solution ρ st eg , evaluated at the time dependent value of the detuning. We can assume, however, that the ions reach a periodic steady state, and for weak modulation, we make the Ansatz,

Applying this Ansatz in (9), and isolating terms ∝ e ±iω M t leads to the relations, a

ω M +δ a 0 where δ = δ r + iΓ/2. Inserting these results in the equation for the non-oscilating terms in the periodic solutions to (9), we obtain, a 0 -Ω/(2δ r ), which holds to first order in the small modulation amplitude of the detuning ε.

Collecting terms, the periodically varying coherence writes,

and the coherence, indeed, follows the stationary solution adiabatically in this coupling regime.

When ω M is comparable to or larger than |δ r |, we see a more complicated dependence of the real (dispersive) and imaginary (absorptive) part of the coherence, and also a weaker dependence on the position in the high frequency limit, where the coherence is unable to follow the time dependent steady state.


Fabien Mathieu
email: fabien.mathieu@nokia-bell-labs.com

Kleinberg's Grid Reloaded

Keywords: 1998 ACM Subject Classification C.2.2 Routing Protocols, C.2.4 Distributed Systems Keywords and phrases Small-World Routing, Kleinberg's Grid, Simulation, Rejection Sampling Digital Object Identifier 10.4230/LIPIcs.OPODIS.2016.21

diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In a prescient 1929 novel called Láncszemek (in English, Chains), Karinthy imagines that any two people can be connected by a small chain of personal links, using no more than five intermediaries [10].

Years later, Milgram validates the concept by conducting real-life experiments. He asks volunteers to transmit a letter to an acquaintance with the objective to reach a target destination across the United States [START_REF] Milgram | The small world problem[END_REF][START_REF] Travers | An experimental study of the small world problem[END_REF]. While not all messages arrive, successful attempts reach destination after six hops in average, popularizing the notion of six degrees of separation.

Yet, for a long time, no theoretical model could explain why and how this kind of smallworld routing works. One of the first and most famous attempts to provide such a model is due to Jon Kleinberg [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF]. He proposes to abstract the social network by a grid augmented with shortcuts. If the shortcuts follow a heavy tail distribution with a specific exponent, then a simple greedy routing can reach any destination in a short time (O(log 2 (n)) hops). On the other hand, if the exponent is wrong, then the time to reach destination becomes Ω(n α ) for some α. This seminal work has led to multiple studies from both the theoretical and empirical social systems communities.

making any concession about the assumptions made in the original model, which is kept untouched.

Fueled by the capacity to obtain quick and accurate results even for very large grids, we give a fresh look on Kleinberg's grid, through three independent small studies. First, we show that the model is in practice more robust than expected: for grids of given size there is quite a large range of exponents that grant short routing paths. Then we observe that the lower bounds proposed by Kleinberg in [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF] are not tight and suggest new bounds. Finally, we compare Kleinberg's grid to Milgram's experiment, and observe that when the grid parameters are correctly tuned, the performance of greedy routing is consistent with the six degrees of separation phenomenon.

Roadmap

Section 2 presents the original augmented grid model introduced by Kleinberg and the greedy routing algorithm. A brief overview of the main existing theoretical and experimental studies is provided, with a strong emphasis on the techniques that can be used for the numerical evaluation of Kleinberg's model.

In Section 3, we give our algorithm for estimating the performance of greedy routing. We explain the principle of dynamic rejection sampling and detail why it allows to perfectly emulate Kleinberg's grid with the same speed that can be achieved by toroidal approximations. We also give a performance evaluation of the simulator based on our solution. For readers interested in looking under the hood, a fully working code (written in Julia) is given in Appendix A.

To show the algorithm benefits, we propose in Section 4 three small studies that investigate Kleinberg's model from three distinct perspectives: robustness of greedy routing with respect to the shortcut distribution (Section 4.1); tightness of the existing theoretical bounds (Section 4.2); emulation of Milgram's experiment within Kleinberg's model (Section 4.3).

Model and Related Work

We present here the model and notation introduced by Kleinberg in [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF][START_REF] Kleinberg | Navigation in a small world[END_REF], some key results, and a brief overview of the subsequent work on the matter.

Kleinberg's Grid

In [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF], Kleinberg considers a model of directed random graph G(n, r, p, q), where n, p, q are positive integers and r is a non-negative real number. A graph instance is built from a square lattice of n × n nodes with Manhattan distance d:

if u = (i, j) and v = (k, l), then d(u, v) = |i -j| + |k -l|.
d represents some natural proximity (geographic, social, . . . ) between nodes. Each node has some local neighbors and q long range neighbors. The local neighbors of a node u are the nodes v such that d(u, v) ≤ p. The q long range neighbors of u, also called shortcuts, are drawn independently and identically as follows: the probability that a given long edge starting from u arrives in v is proportional to (d(u, v)) -r . The problem of decentralized routing in a G(n, r, p, q) instance consists in delivering a message from node u to node v in a hop-by-hop basis. At each step, the message bearer needs to choose the next hop among its neighbors. The decision can only use the lattice coordinates of the neighbors and destination. The main example of decentralized algorithm is the greedy routing, where at each step, the current node chooses the neighbor that is closest to destination based on d (in case of ties, an arbitrary breaking rule is used).

The main metric to analyze the performance of a decentralized algorithm is the expected delivery time, which is the expected number of hops to transmit a message between two nodes chosen uniformly at random in the graph.

This paper focuses on studying the performance of the greedy algorithm. Unless stated otherwise, we assume p = q = 1 (each node has up to four local neighbors and one shortcut). Let e r (n) be the expected delivery time of the greedy algorithm in G(n, r, 1, 1).

Theoretical Results

The main theoretical results for the genuine model are provided in the original papers [START_REF] Kleinberg | Navigation in a small world[END_REF][START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF], where Kleinberg proves the following:

e 2 (n) = O(log 2 (n));
for 0 ≤ r < 2, the expected delivery time of any decentralized algorithm is Ω(n (2-r)/3 ); for r > 2, the expected delivery time of any decentralized algorithm is Ω(n (r-2)/(r-1) ). Kleinberg's results are often interpreted as follows: short paths are easy to find only in the case r = 2. The fact that only one value of r asymptotically works is sometimes seen as the sign that Kleinberg's model is not robust enough to explain the small-world routing proposed by Karinthy and experimented by Milgram. However, as briefly discussed by Kleinberg in [START_REF] Easley | The small-world phenomenon[END_REF], there is in fact some margin if one considers a grid of given n. This tolerance will be investigated in more details in Section 4.1.

While we focus here on the original model, let us give a brief, non-exhaustive, overview of the subsequent extensions that have been proposed since. Most proposals refine the model by considering other graph models or other decentralized routing algorithms. New graph models are for example variants of the original model (studying grid dimension or the number of shortcuts per node [START_REF] Barrière | Efficient routing in networks with long range contacts[END_REF][START_REF] Easley | The small-world phenomenon[END_REF][START_REF] Fraigniaud | Greedy routing in small-world networks with power-law degrees[END_REF]), graphs inspired by peer-to-peer overlay networks [START_REF] Singh Manku | Know thy neighbor's neighbor: The power of lookahead in randomized P2P networks[END_REF], or arbitrary graphs augmented with shortcuts [START_REF] Fraigniaud | Universal Augmentation Schemes for Network Navigability: Overcoming the (n)-Barrier[END_REF][START_REF] Fraigniaud | On the searchability of small-world networks with arbitrary underlying structure[END_REF]. Other proposals of routing algorithms usually try to enhance the performance of the greedy one by granting the current node additional knowledge of the topology [START_REF] Fraigniaud | Eclecticism shrinks even small worlds[END_REF][START_REF] Singh Manku | Know thy neighbor's neighbor: The power of lookahead in randomized P2P networks[END_REF][START_REF] Martel | Analyzing kleinberg's (and other) small-world models[END_REF].

A large part of the work above aims at improving the O(log 2 (n)) bound of the greedy routing. For example, in the small-world percolation model, a variant of Kleinberg's grid with O(log(n)) shortcuts per node, greedy routing performs in O(log(n)) [START_REF] Singh Manku | Know thy neighbor's neighbor: The power of lookahead in randomized P2P networks[END_REF].

Experimental Studies

Many empirical studies have been made to study how routing works in real-life social networks and the possible relation with Kleinberg's model (see for example [START_REF] Liben-Nowell | Geographic routing in social networks[END_REF][START_REF] Easley | The small-world phenomenon[END_REF] and the references within). On the other hand, numerical evaluations of the theoretical models are more limited to the best of our knowledge. Such evaluations are usually performed by averaging R runs of the routing algorithm considered.

In [START_REF] Kleinberg | Navigation in a small world[END_REF], Kleinberg computes e r (n) for n = 20, 000 and r ∈ [0, 2.5], using 1,000 runs per estimate. However, he uses a torus instead of a regular grid (this will be discussed later in the paper). In [START_REF] Singh Manku | Know thy neighbor's neighbor: The power of lookahead in randomized P2P networks[END_REF], networks of size up to 2 24 (corresponding to n = 2 12 in the grid) are investigated using 150 runs per estimate.

Closer to our work, Athanassopoulos et al. propose a study centered on numerical evaluation that looks on Kleinberg's model and some variants [START_REF] Athanassopoulos | An experimental study of greedy routing algorithms[END_REF]. For the former, they differ from the original model by having fixed source and destination nodes. They compute e r (n) for values of n up to 3,000 and r ∈ {0, 1, 2, 3}, using 900 runs per estimate.

To compare with, in the present paper, we consider values of n up to 2 24 ≈ 16, 000, 000 and r ∈ [0, 3], with at least 10,000 runs per estimate. To explain such a gain, we first need O P O D I S 2 0 1 6

21:4

Kleinberg's Grid Reloaded to introduce the issue of shortcuts computation.

Drawing shortcuts

As stated in [START_REF] Athanassopoulos | An experimental study of greedy routing algorithms[END_REF], the main computational bottleneck for simulating Kleinberg's model comes from the shortcuts.

There are n2 shortcuts in the grid (assuming q = 1); When one wants to a shortcut, any of the n 2 -1 other nodes can be chosen with non-null probability. This can be made by inverse transform sampling, with a cost Ω(n 2 );

The shortcut distribution depends on the node u considered, even if one uses relative coordinates. For example, a corner node will have i + 1 neighbors at distance i for 1 ≤ i < n, against 4i neighbors for inner nodes (as long as the ball of radius i stays inside the grid). This means that, up to symmetry, each node has a unique shortcut distribution 1 . This prevents from mutualising shortcuts drawings between nodes.

In the end, building shortcuts as described above for each of the R runs has a time complexity Ω(Rn 4 ), which is unacceptable if one wants to evaluate e r (n) on large grids.

The first issue is easy to address: as observed in [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF][START_REF] Singh Manku | Know thy neighbor's neighbor: The power of lookahead in randomized P2P networks[END_REF][START_REF] Athanassopoulos | An experimental study of greedy routing algorithms[END_REF], we can use the Principle of deferred decision [START_REF] Mitzenmacher | Probability and Computing: Randomized Algorithms and Probabilistic Analysis[END_REF] and compute the shortcuts on-the-fly as the path is built, because they are drawn independently and a node is never used twice in a given path. This reduces the complexity to Ω(Rn 2 e r (n)).

Torus approximation

To lower the complexity even more, one can approximate the grid by the torus. This is the approach adopted in [START_REF] Kleinberg | Navigation in a small world[END_REF][START_REF] Singh Manku | Know thy neighbor's neighbor: The power of lookahead in randomized P2P networks[END_REF]. The toroidal topology brings two major features compared to a flat grid:

The distribution of the relative position of the shortcut does not depend on the originating node. This enables to draw shortcuts in advance (in bulk); There is a strong radial symmetry, allowing to draw a "radius" and an "angle" separately.

To illustrate the gain of using a torus instead of a grid, consider the drawing of k shortcuts from k distinct nodes. In a grid, if one uses inverse transform sampling for each corresponding distribution, the cost is Ω(n 2 k). In the torus, one can compute the probabilities to be at distance i for i between 1 and n (the maximal distance in the torus), draw k radii, then choose for each drawn radius i a node uniformly chosen among those at distance i. Assuming drawing a float uniformly distributed over [0, 1) can be made in O(1), the main bottleneck is the drawing of radii. Using bulk inverse transform sampling, it can be performed in O(n + k log(k)), by sorting k random floats, matching then against the cumulative distribution of radii and reverse sorting the result.

Fast Estimation of Expected Delivery Time

We now describe our approach for computing e r (n) in the flat grid with the same complexity than for the torus approximation.

u #1 #2 #3
Figure 1 Main idea of the dynamic rejection sampling approach (n = 4).

Dynamic rejection sampling for drawing shortcuts

In order to keep low computational complexity without making any approximation of the model, we propose to draw a shortcut of a node u as follows:

1. We embed the actual grid G (we use here G to refer to the lattice nodes of G(n, r, p, q)) in a virtual lattice B u made of points inside a ball of radius 2(n -1). Note that the radius chosen ensures that G is included in B u no matter the location of u. 2. We draw a node inside B u such that the probability to pick up a node v is proportional to (d(u, v)) -r . This can be done in two steps (radius and angle):

For the radius, we notice that the probability to draw a node at distance i is proportional to i 1-r , so we pick an integer between 1 and 2(n -1) such that the probability to draw i is to i 1-r / 2(n-2) k=1 k 1-r . For the angle, pick an integer uniformly chosen between 1 and 4i 3. This determines a unique point v among the 4i points at distance i from u in the virtual lattice, chosen with a probability proportional to (d(u, v)) -r . If v belongs to the actual grid, it becomes the shortcut, otherwise we try again (back to step #2).

This technique, illustrated in Figure 1, is inspired by the rejection sampling method [START_REF] John Von Neumann | Various techniques used in connection with random digits[END_REF]. By construction, it gives the correct distribution: the node v that it eventually returns is in the actual grid and has been drawn with a probability proportional to (d(u, v)) -r .

We call this dynamic rejection sampling because the sampled distribution changes with the current node u. Considering u as a relative center, the actual grid G moves with u and acts like an acceptance mask. On the other hand, the distribution over the virtual lattice B u remains constant. This enables to draw batches of relative shortcuts that can be used over multiple runs, exactly like for the torus approximation.

The only possible drawback of this approach is the number of attempts required to draw a correct shortcut. Luckily, this number is contained. ◮ Lemma 1. The probability that a node drawn in B u belongs to G is at least 1 8 .

Proof. We will prove that

v∈G\{u} (d(u, v)) -r v∈Bu\{u} (d(u, v)) -r > 1 8 .
We use the fact that the probability decreases with the distance combined with some geometric arguments. Let G c be a n × n lattice that has u as one of its corner. Let H u the (2n -1) × (2n -1) lattice centered in u.

In terms of probability of drawing a node in G \ {u}, the worst case is when u is at some corner: there is a bijection

f from G to G c such that for all v ∈ G \ {u}, d(u, f (v)) ≥ d(u, v).
Such a bijection can be obtained by splitting G into G ∩ G c and three other sub-lattices that are flipped over G ∩ G c (see Figure 2). This gives

v∈G\{u} (d(u, v)) -r ≥ v∈G\{u} (d(u, f (v))) -r = v∈Gc\{u} (d(u, v)) -r .
Then we observe that the four possible lattices G c obtained depending on the corner occupied by u fully cover H u . In fact, axis nodes are covered redundantly. This gives

v∈Hu\{u} (d(u, v)) -r < 4 v∈Gc\{u} (d(u, v)) -r .
Lastly, if one folds B u \ H u back into H u like the corners of a sheet of paper, we get a strict injection g from

B u \ H u to H u \ {u} (the diagonal nodes of H u are not covered). Moreover, for all v ∈ B u \ H u , d(u, v) ≥ d(u, h(v)). This gives v∈Bu\Hu (d(u, v)) -r ≤ v∈Bu\Hu (d(u, h(v))) -r < v∈Hu\{u} (d(u, v)) -r .
This concludes the proof, as we get

v∈G\{u} (d(u, v)) -r v∈Bu\{u} (d(u, v)) -r ≥ v∈Gc\{u} (d(u, v)) -r v∈Hu\{u} (d(u, v)) -r + v∈Bu\Hu (d(u, v)) -r > 1 8
.

◭

Fabien

Remarks

When r = 0 (uniform shortcut distribution), the bound 1 8 is asymptotically tight: the success probability is exactly the ratio between the number of nodes in G \ {u} and B u \ {u}, which is

n 2 -1 4(n-1)(2n-1) -→ n→+∞ 1 8
. On the other hand, as r grows, the probability mass gets more and more concentrated around u (hence in G), so we should expect better performance (cf Section 3.2). The dynamic rejection sampling approach can be used in other variants of Kleinberg's model, like for other dimensions or when the number of shortcuts per node is a random variable (like in [START_REF] Fraigniaud | Greedy routing in small-world networks with power-law degrees[END_REF]). The only requirement is the existence of some root distribution (like the distribution over B u here) that can be carved to match any of the possible distributions with a simple acceptance test. Only the nodes from H u may belong to G, so nodes from B u \ G are always sampled for nothing. For example, in Figure 1, this represents 36 nodes over 84. By drawing shortcuts in H u \ {u} instead of B u \ {u}, we could increase the success rate lower bound to 1/4. However, this would make the algorithm implementation more complex (the number of nodes at distance i is not always 4i), which is in practice not worth the factor 2 improvement of the lower bound. This may not be true for a higher dimension β. Adapting the proof from Lemma 1, we observe that using a ball of radius β(n -1) will lead to a bound β!(2β) -β , while a grid of side 2n -1 will lead to 2 -β . The two bounds are asymptotically tight for r = 0. In that case the grid approach is β β β! more efficient than the ball approach (this grows exponentially with β).

Performance Evaluation

We implemented our algorithm in Julia 0.4.5. A working code example is provided in Appendix A. Simulations were executed on a low end device (a Dell tablet with 4 Gb RAM and Intel M-5Y10 processor), which was largely sufficient for getting fast and accurate results on very large grids, thanks to the dynamic rejection sampling.

Unless said otherwise, e r (n) is obtained by averaging R = 10, 000 runs. consider powers of two ranging from 2 7 (about 16,000 nodes) to 2 24 (about 280 trillions nodes). We focus on r ∈ [0, 3].

Regarding the choice of the bulk size k (we assume here the use of bulk inverse transform sampling), the average number of shortcuts to draw over the R runs is Re r (n). This leads to an average total cost for drawing shortcuts in O(⌈ Rer (n) k ⌉(n+k log(n))). k can be optimized if e r (n) is known, but this requires bootstrapping the estimation of e r (n). Yet, we can remark that for n fixed and R large enough, we should choose k of the same order of magnitude than n, which ensures an average cost per shortcut in O(log(n)). This is the choice we made in our code, which gives a complexity in O(max(Re r (n), n) log(n)). This is not efficient for n ≫ Re r (n) (the bulk is then over-sized, so lot of unused shortcuts are drawn), but this seldom happens with our settings. Figure 3 presents the time needed to compute e r (n) as a function of n for r ∈ {1, 2, 5 2 }. We observe running times ranging from seconds to a few minutes. To compare with, in [START_REF] Athanassopoulos | An experimental study of greedy routing algorithms[END_REF], which is to the best of our knowledge the only work disclosing computation times, one single run takes about 4 seconds for n = 400. With a similar time budget, our implementation averages 10,000 runs for n = 2 17 ≈ 130, 000 (r = 1 or r = 2.5), or up to n = 2 24 ≈ 16, 000, 000 for the optimal exponent r = 2. In other words, we are several orders of magnitude faster.

We also looked at the cost of the dynamic rejection sampling approach in terms of shortcuts drawn outside of G. Figure 4 shows the success frequency of the sampling as a function of r for n = 2 14 (the actual value of n has no significant impact as long as it is large enough). We verify Lemma 1: the success rate is always at least 1/8 = 0.125, which is a tight bound for r = 0. For r = 1, the rate is about 0.29, and it climbs to 0.86 for r = 2. Failed shortcuts become negligible (rate greater than 0.99) for r ≥ 2.5. Overall, Figure 4 shows that the cost of drawing some shortcuts outside G is a small compared to the benefits offered by drawing shortcuts from a unique distribution.

Remark

In terms of success rate, Figure 4 shows that r = 2.5 is much more efficient than r = 1, but running times tend to be slightly longer for r = 2.5 (cf Figure 3). The reason is that e 1 (n) is lower than e 2.5 (n), which overcompensates the success rate difference.

Applications

Given the tremendous amount of strong theoretical results on small-world routing, one can question the interest of proposing a simulator (even a fast one!).

In this Section, we prove the interest of numerical evaluation through three (almost) independent small studies.

Efficient enough exponents

In [START_REF] Kleinberg | Navigation in a small world[END_REF], Kleinberg gives an estimation of e r (20, 000) using a torus approximation (cf Section 2.3.2). A few years later, he discussed in more details the results, observing that [START_REF] Easley | The small-world phenomenon[END_REF]: e r (n) stays quite similar for r ∈ [1. [START_REF] Easley | The small-world phenomenon[END_REF][START_REF] Barrière | Efficient routing in networks with long range contacts[END_REF]; the best value of r is actually slightly lower than 2.

We believe that these observations are very important as they show that routing in Kleinberg's grid is more robust than predicted by theory: it is efficient as long as r is close enough to 2.

Using our simulator, we can perform the same experiment. As shown by Figure 5, the results are quite similar to the ones observed in [START_REF] Kleinberg | Navigation in a small world[END_REF].

Yet, there is a small but essential difference between the two experiments: Kleinberg approximated his grid by a torus while we stay true to the original model. Why do we use the word essential? Both shapes have the same asymptotic behavior (proofs in [START_REF] Kleinberg | Navigation in a small world[END_REF] are straightforward to adapt), so why should we care? It seems to us that practical robustness is an essential feature if one wants to accept Kleinberg's grid as a reasonable model for routing in social networks. To the best of our knowledge, no theoretical work gives quantitative

ropt(n) rmin(e2(n)) rmin(2e2(n))
Figure 6 Reasonable values of r results on this robustness, so we need to rely on numerical evaluation. But when Kleinberg uses a torus approximation, we can not rule out that the observed robustness is a by-side effect of the toroidal topology. The observation of a similar phenomenon on a flat grid discards this hypothesis. In fact, it suggests (without proving) that the robustness with respect to the exponent for grids of finite size may be a general phenomenon.

We propose now to investigate this robustness in deeper details. We have evaluated the following values, which outline, for a given n, the values of r that can be considered reasonable for performing greedy routing:

The value of r that minimizes e r (n), denoted r opt ;

The smallest value of r such that e r (n) ≤ e 2 (n), denoted r min (e 2 (n));

The smallest and largest values of r such that e r (n) ≤ 2e 2 (n), denoted r min (2e 2 (n)) and r max (2e 2 (n)) respectively.

The results are displayed in Figure 6. All values but r opt are computed by bisection. For r opt , we use a Golden section search [START_REF] William | Minimization or maximization of functions[END_REF]. Finding a minimum requires more accuracy, so the search of r opt is set to use R = 1, 000, 000 runs per estimation. Luckily, as the computation operates by design through near-optimal values, we can increase the accuracy with reasonable running times.

Besides confirming that r = 2 is asymptotically the optimal value, Figure 6 shows that the range of reasonable values for finite grids is quite comfortable. For example, considering the range of values where e r (n) is less than twice e 2 (n), we observe that:

For n ≤ 2 11 (less than four million nodes), any r between 0 and 2.35 works;

For n ≤ 2 14 (less than 270 million nodes), the range is between 0.85 and 2.26. Even for n = 2 24 (about 280 trillions nodes), all values of r between 1.58 and 2.16 can still be considered efficient enough.

Asymptotic Behavior

Our simulator can be used to verify the theoretical bounds proposed in [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF]. For example, Figure 7 shows e r (n) for r equal to 1, 2, and 5 2 .

21:11 As predicted, e r (n) seems to behave like log 2 (n) for r = 2 and like n α for the two other cases. Yet, the exponents found differ from the ones proposed in [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF]. For both r = 1 and r = 2.5, we observe α = 1 2 , while the lower bound is 1 3 . Intrigued by the difference, we want to compute α as a function of r.

However, we see in Figure 7 that a log 2 (n) curve appears to have a positive slope in a logarithmic scale, even for large values of n. This may distort our estimations. To control the possible impact of this distortion, we estimate the exponent at two distinct scales:

n ∈ [2 15 , 2 20 ], using the estimation α 1 := log 2 (er (2 20 ))-log 2 (er (2 15 )) 5

;

n ∈ [2 20 , 2 24 ], using the estimation α 2 := log 2 (er (2 24 ))-log 2 (er (2 20 ))

4

. The Results are displayed in Figure 8. The range of r was extended to [0, 4] due to the observation of a new transition for r = 3

As feared, our estimations do not indicate 0 for r = 2, but the fact that α 2 is closer to 0 than α 1 confirms that this is likely caused by the log 2 (n) factor. Moreover, we observe that α 1 and α 2 only differ for r ≈ 2 and r ≈ 3, suggesting that both estimates should be accurate except for these critical values.

Based on Figures 7 and8, it seems that the bounds proposed by Kleinberg in [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF] are only tight for r = 0 and r = 2. Formally proving more accurate bounds is beyond the scope O P O D I S 2 0 1 6 and spirit of this paper, but we can conjecture new bounds hinted by our simulations:

For 0 ≤ r < 2, we have e r (n) = Θ(n 2-r 3-r ); For 2 < r < 3, we have e r (n) = Θ(n r-2 ); For r > 3, we have e r (n) = Θ(n). This conjectured bounds are consistent with the one proposed in [START_REF] Barrière | Efficient routing in networks with long range contacts[END_REF], where it is proved that for the 1-dimensional ring, we have:

For 0 ≤ r < 1, e r (n) = Ω(n

1-r 2-r ); For r = 1, e r (n) = Θ(log 2 (n)) For 1 < r < 2, e r (n) = O(n r-1 ); For r = 2, e r (n) = O(n log(log(n)) log(n) ) For r > 2, e r (n) = O(n).
It is likely that the proofs in [START_REF] Barrière | Efficient routing in networks with long range contacts[END_REF] can be adapted to the 2-dimensional grid, but some additional work may be necessary to demonstrate the bounds' tightness. Moreover, our estimates may have missed some slower-than-polynomial variations. For example, the logarithms in the bounds for r = 2 in [START_REF] Barrière | Efficient routing in networks with long range contacts[END_REF] may have a counterpart in the grid for r = 3. This would explain why the estimates are not sharp around that critical value (like for the case r = 2 and the term in log 2 (n)).

Remark

For r ≥ 3.5, we have observed that e r (n) ≈ 2 3 n, which is the expected delivery path in absence of shortcuts: for these high values of r, almost all shortcuts link to an immediate neighbor of the current node, which make them useless, and the rare exceptions are so short that they make no noticeable difference.

Six degrees of separation

What makes Milgram's experiments [START_REF] Milgram | The small world problem[END_REF][START_REF] Travers | An experimental study of the small world problem[END_REF] so famous is the surprising observation that only a few hops are required to transmit messages in social networks, a phenomenon called six degrees of separation in popular culture.

Yet, the values of e r (n) observed until now are quite far for the magic number six. For example, in Figure 5, the lowest value is 140. In addition to the asymptotic lack of robustness with respect to the exponent (discussed in Section 4.1), this may partially motivate the amount of work accomplished to increase the realism of the model and the performance of the routing algorithm.

In our opinion, a good model should be as simple as possible while being able to accurately predict the phenomenon that needs to be explained. We propose here to answer a simple question: is the genuine model of greedy routing in G(n, r, p, q) a good model? Sure, it is quite simple. To discuss its accuracy, we need to tune the four parameters n, r, p and q to fit the conditions of Milgram's experiments as honestly as we can. Size The experiments of Milgram were conducted in the United States, with a population of about 200,000,000 at the late sixties. All inhabitants were not susceptible to participate to the experiments: under-aged, undergraduate or disadvantaged people may be considered as de facto excluded from the experiments. Taking that into consideration, the correct n is probably somewhere in the range [5000,14000]. We propose to set n = 8, 500, which corresponds to about 72,000,000 potential subjects. Exponent In [START_REF] Easley | The small-world phenomenon[END_REF], Kleinberg investigates how to relate the r-harmonic distribution with reallife observations. He surveys multiple social experiments and discusses the correspondence with the exponent of his model, which gives estimates of r between 1.75 and 2.2. Neighborhood The default value p = q = 1 means that there are no more than five "acquaintances" per node. This is quite small compared to what is observed in real-life social networks. For example, the famous Dunbar's number, which estimates the number of active relationships, is 150 [START_REF] Dunbar | Neocortex size as a constraint on group size in primates[END_REF]. More recent studies seem to indicate that the average number of acquaintances is larger, ranging from 250 to 1500 (see [START_REF]Ithiel de Sola Pool and Manfred Kochen[END_REF][START_REF] Tyler | How many people do you know?: Efficiently estimating personal network size[END_REF][START_REF] Wellman | Is Dunbar's number up?[END_REF] and references within). We propose to set p and q so that the neighborhood size 2p(p + 1) + q is about 600, the value reported in [START_REF] Tyler | How many people do you know?: Efficiently estimating personal network size[END_REF]. Regarding the partition between local links (p) and shortcuts (q), we consider three typical scenarios: p = 1, q = 600 (shortcut scenario: the neighborhood is almost exclusively made of shortcuts, and local links are only here to ensure the termination of greedy routing). p = 10, q = 380 (balanced scenario). p = 15, q = 120 (local scenario, with a value of q not too far from Dunbar's number).

Having set all parameters, we can evaluate the performance of greedy routing. The results are displayed in Figure 9. We observe that the expected delivery time roughly stands between five and six for a wide range of exponents.

r ∈ [1.4, 2.3] for the shortcut scenario. r ∈ [1.3, 2.3] for the balanced scenario. r ∈ [1. [START_REF]Ithiel de Sola Pool and Manfred Kochen[END_REF][START_REF] Barrière | Efficient routing in networks with long range contacts[END_REF] for the local scenario. Except for the local scenario, which leads to slightly higher routing times for r > 2, the six degrees of separation are achieved for all values of r that are consistent with the observations surveyed in [START_REF] Easley | The small-world phenomenon[END_REF]. This allows to answer our question: the augmented grid proposed by Kleinberg is indeed a good model to explain the six degrees of separation phenomenon.

O P O D I S

Conclusion

We proposed an algorithm to evaluate the performance of greedy routing in Kleinberg's grid. Fueled by a dynamic rejection sampling approach, the simulator based on our solution performs several orders of magnitude faster than previous attempts. It allowed us to investigate greedy routing under multiple perspective.

We noted that the performance of greedy routing is less sensitive to the choice of the exponent r than predicted by the asymptotic behavior, even for very large grids. We observed that the bounds proposed in [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF] are not tight except for r = 0 and r = 2. We conjectured that the tight bounds are the 2-dimensional equivalent of bounds proposed in [START_REF] Barrière | Efficient routing in networks with long range contacts[END_REF] for the 1-dimensional ring.

We claimed that the model proposed by Kleinberg in [START_REF] Kleinberg | The small-world phenomenon: An algorithmic perspective[END_REF][START_REF] Kleinberg | Navigation in a small world[END_REF] is a good model for the six degrees of separation, in the sense that it is very simple and accurate.

Our simulator is intended as a tool to suggest and evaluate theoretical results, and possibly to build another bridge between theoretical and empirical study of social systems. We hope that it will be useful for researchers from both communities. In a future work, we plan to make our simulator more generic so it can handle other types of graph and augmenting schemes.

A

Code for Expected Delivery Time (tested on Julia Version 0.4.5) using S t a t s B a s e # For u s i n g J u l i a b u i l t -i n sample f u n c t i o n

Figure 2

 2 Figure 2Graphical representation of G, Gc and Hu (n = 4). The sub-lattices used to built a bijection between G and Gc (cf proof of Lemma 1) are represented with distinct shapes and gray levels.

2 Figure 3

 23 Figure 3 Computing er(n) using R = 10, 000 runs

Figure 4

 4 Figure 4 Frequency of shortcuts drawn in Bu that belongs to G observed during a computation of er(n) (n = 2 14 )

Figure 5

 5 Figure 5 Expected delivery time for n = 20, 000

Figure 7 Figure 8

 78 Figure 7 Expected delivery time for different values of r

rFigure 9

 9 Figure 9 Performance of greedy routing for parameters inspired by Milgram's experiments.

function s h o

  r t c u t s _ b u l k ( n , probas , b u l k _ s i z e ) r a d i i = sample ( 1 : ( 2 * n -2) , probas , b u l k _ s i z e ) s h o r t c u t s = Tuple { Int , I n t } [ ] f o r i = 1 : b u l k _ s i z e r a d i u s = r a d i i [ i ] a n g l e = f l o o r (4 * r a d i u s * rand ()) -2 * r a d i u s push ! ( s h o r t c u t s , ( ( r a d i u s -abs ( a n g l e ) ) , ( s i g n ( a n g l e ) * ( r a d i u s -abs ( r a d i u s -abs ( a n g l e ) ) ) ) ) ) end return s h o r t c u t s end # E s t i m a t e s the Expected D e l i v e r y Time f o r G( n , r , p , q ) o v e r R runs function edt ( n , r , p , q , R) b u l k _ s i z e = n probas = w e i g h t s ( 1 . / ( 1 : ( 2 * n -2 ) ) . ^( r -1)) s h o r t c u t s = s h o r t c u t s _ b u l k ( n , probas , b u l k _ s i z e ) s t e p s = 0 f o r i = 1 :R # s : s t a r t / c u r r e n t node ; a : t a r g e t node ; d : d i s t a n c e t o t a r g e t s_x , s_y , a_x , a_y = t u p l e ( rand ( 0 : ( n -1) , 4 ) . . . ) d = abs ( s_x -a_x ) + abs ( s_y -a_y ) while d>0 sh_x , sh_y = -1, -1 # sh w i l l be b e s t s h o r t c u t node d_s = 2 * n # d_s w i l l be d i s t a n c e from sh t o a f o r j = 1 : q # Draw q s h o r t c u t s ch_x , ch_y = -1, -1 # ch w i l l be c u r r e n t s h o r t c u t c_s = 2 * n # c_s w i l l be d i s t a n c e from ch t o a # Dynamic r e j e c t i o n s ampling while ( ch_x < 0 | | ch_x >= n | | ch_y < 0 | | ch_y >= n ) r_x , r_y = pop ! ( s h o r t c u t s ) ch_x , ch_y = s_x + r_x , s_y + r_y i f is empty ( s h o r t c u t s ) s h o r t c u t s = s h o r t c u t s _ b u l k ( n , probas , b u l k _ s i z e ) end end c_s = abs ( a_x -ch_x ) + abs ( a_y -ch_y ) i f c_s < d_s # maintain b e s t s h o r t c u t found d_s = c_s sh_x , sh_y = ch_x , ch_y end end i f d_s < d-p # Follow s h o r t c u t i f e f f i c i e n t s_x , s_y = sh_x , sh_y d = d_s e l s e # Follow l o c a l l i n k s d = dp delta_x = min ( p , abs ( a_x -s_x ) ) delta_y = p -delta_x s_x += delta_x * s i g n ( a_x -s_x ) s_y += delta_y * s i g n ( a_y -s_y ) end s t e p s += 1 end end s t e p s /= R return s t e p s end

To take advantage of symmetry, one can consider the isometric group of a square grid, which can be built with the quarter-turn and flip operations. However, its size is 8, so even using symmetry, there are at least n

8 distinct (non-isomorphic) distributions.