Yann G Boucher
email: yann.boucher@enib.fr

Alberto Parini

Patrice Féron

Resonant modal conversion in a two-mode waveguide

Keywords: Mode Division Multiplexing, Optical Networks on Chip (ONoCs), Silicon photonics, Scattering parameters, Nonlinear resonators, Bistability

We characterize a system consisting of a two-mode waveguide coupled to a single-mode microring resonator possibly presenting a nonlinear response of Kerr type. By using the scattering parameter formalism extended to the multimode domain, we show that in the linear regime and for an ideally transparent medium, each resonance of the system can be exploited to perform complete even-to-odd (respectively, odd-to-even) modal conversion. Moreover, when the Kerr nonlinearity is effective, the microring enables a power-dependent modal switching mediated by phase bistability. Thanks to its mode-processing capabilities, this configuration is suitable to find application as a functional building-block in mode-division multiplexing (MDM) photonic integrated circuits.

Introduction

A deep research effort is currently carried out by academic and industrial laboratories toward the integration on a same technological substrate of both microelectronic and photonic devices. In fact, as energy issues affect progressively the scalability of high-performance multi-core microprocessor architectures, the implementation of optical links (Optical Network-on-Chip -ONoC) between computational cores is foreseen as a possible way to reduce the overall microprocessor power consumption while increasing at the same time the transmission bandwidth and data throughput (see [START_REF] Miller | Device Requirements for Optical Interconnects to Silicon Chips[END_REF], [START_REF] Bergman | Photonic Network-on-Chip Design[END_REF], [START_REF] Shacham | Photonic Networks-on-Chip for Future Generations of Chip Multiprocessors[END_REF] and references therein). The large majority of ONoCs solutions proposed in the literature rely for their operation on wavelength-division Multiplexing (WDM). WDM channels are easily routed, switched and multiplexed by means of functional devices, such as microring resonators [12], [START_REF] Sherwood-Droz | Optical 4 4 hitless silicon router for optical Networks-on-Chip (NoC)[END_REF], [START_REF] Parini | BER evaluation of a passive SOI WDM router[END_REF] or Mach-Zehnder interferometers [START_REF] Yang | Non-Blocking 4x4 Electro-Optic Silicon Switch for On-Chip Photonic Networks[END_REF], [START_REF] Yang | Reconfigurable nonblocking 4-port silicon thermo-optic optical router based on Mach-Zehnder optical switches[END_REF], which are relatively easy to integrate on silicon photonic chips. The on-chip exploitation of WDM-based communication strategies implies the availability of silicon-compatible laser sources, eventually multi-wavelength. Although some recent remarkable fabrications have confirmed the technological feasibility of silicon photonic transceivers including III-V bonded laser sources [START_REF] Zhang | 8 8 40 gbps fully integrated silicon photonic network on chip[END_REF], the integration of active cavities on silicon remains a challenging task, with the existing technological solutions providing contrasting performance in terms of available optical power, electric-to-optical power conversion efficiency and footprint. To reduce the number of sources or, in a complementary way, to optimize the use of the available ones, a possible solution relies on the exploitation of a further dimension inherent to optical communications, namely the modal order [START_REF] Peucheret | Signal processing for On-Chip Space Division Multiplexing[END_REF]. In mode-division multiplexing (MDM) techniques, communication channels are allocated on the same wavelength to the various propagation modes of a multi-mode waveguide. Therefore, used in conjunction with WDM, MDM can potentially increase the aggregate communication bandwidth of an optical link [START_REF] Luo | WDM-compatible mode-division multiplexing on a silicon chip[END_REF] and add a further degree of freedom for the design of on-chip integrated optical systems. The implementation of a complete MDM on-chip interconnection system requires the definition of a full range of new functional building-blocks able to selectively process optical signals with respect to the spatial (modal) dimension. Such functional blocks may include add-drop mode multiplexers [START_REF] Ding | On-chip two-mode division multiplexing using tapered directional coupler-based mode multiplexer and demultiplexer[END_REF] [START_REF] Yang | Silicon nitride three-mode division multiplexing and wavelength-division multiplexing using asymmetrical directional couplers and microring resonators[END_REF], modeselective routers [START_REF] Stern | On-chip mode-division multiplexing switch[END_REF] and mode-exchangers. These last, in particular, play a key role as they can enable transparent reconfigurability of the transmission paths within the mode-multiplexed network. Different possible schemes enabling mode-exchange have been recently experimentally demonstrated in literature [START_REF] Sun | Integrated switchable mode exchange for reconfigurable mode-multiplexing optical networks[END_REF], [START_REF] Ye | On-chip data exchange for mode division multiplexed signals[END_REF]. In this work we specifically focus on a configuration consisting of a two-mode waveguide coupled to a single-mode resonator eventually presenting a nonlinear response of the Kerr type. Here, the mode-conversion functionality relies on the presence of the side-coupled resonator that, by breaking the invariance of the system, allows a power exchange between the even and odd modes of the two-mode waveguide, which would otherwise be orthogonal. Moreover, the possible presence of a nonlinear cubic (Kerr) response in the microring entails a dependence of the mode-conversion effect on the signal power through a bistable dynamic [START_REF] Chen | Bistability and self-pulsation phenomena in silicon microring resonators based on nonlinear optical effects[END_REF], [START_REF] Rukhlenko | Analytical study of optical bistability in silicon ring resonators[END_REF], [START_REF] Dumeige | Dispersive tristability in microring resonators[END_REF], [START_REF] Sarid | Analysis of bistability in a ring-channel waveguide[END_REF]. A compact analytical expression for the transfer function of this configuration, providing also a full insight on its nonlinear behavior, is derived by observing that a multimode propagation can be interpreted as a coupling problem between single-mode fields. By exploiting this consideration, the scattering-parameter formalism, typically adopted for the dimensioning of single-mode devices and circuits, can be extended to also model multi-mode configurations. The article is organized as follows: In Section (2) the system configuration is presented, together with the definition of the relevant governing parameters. In Section (3), the formal equivalence between a two-mode waveguide and a system of two coupled single-mode waveguide is outlined. In Section (4) the interaction zone between the two-mode waveguide and the microring resonator is described in terms of a ternary single-mode directional coupler. In Section (5) the modal transfer function of the entire system is established in a compact matrix form. Finally, in Section (6), the investigation is extended to the nonlinear domain by introducing a Kerr response in the microring resonator. Conclusions and perspectives are drawn at the end.

2 System layout and governing parameters

The system under investigation consists of a two-mode waveguide (WG), coupled to a single-mode racetrack microring resonator (MRR) of total perimeter L, as depicted in Fig. 1. Evanescent coupling between the two-mode waveguide and the single-mode microring takes place along an interaction zone of total length d. For the sake of clarity, we limit the discussion to a single eigenstate of polarization, either TE or TM, so that each mode can be represented by an associated scalar value. The time dependence is taken as exp(+iωt) and z is the direction of propagation in the straight waveguide. The two-mode waveguide is dimensioned to guide a fundamental mode with even symmetry and a first-order mode with odd symmetry, both characterized by their respective propagation constants β even and β odd , with β even > β odd . We consider here a device designed to perform on-chip optical signal processing on a CMOS compatible platform, therefore we identify the even and odd modes as the T E 0 and the T E 1 modes of a multi-mode silicon on insulator (SOI) waveguide.

In order to derive for this system an analytical input-output relation under the form of a scattering (Jones-like) matrix, we exploit the equivalence that can be established between a two-mode waveguide and a set of two coupled single-mode waveguides [START_REF] Boucher | Analytical model for the coupling constant of a directional coupler in terms of slab waveguides[END_REF], as will be detailed in the next section.

3 Equivalence between a two-mode waveguide and a symmetric single-mode coupler

In the framework of coupled-mode theory (CMT) [START_REF] Yariv | Optical Waves in Crystals[END_REF], [START_REF] Tamir | Guided-Wave Optoelectronics[END_REF], the periodic power exchange that takes place in a symmetric coupler between two single-mode waveguides separated by a distance g, such as the one schematized in Fig. 2(A), can be interpreted as a pattern resulting from the simultaneous propagation of the even and odd super-modes of the structure. The common propagation constant β of the two single-mode waveguides, and their mutual coupling strength χ, can be related to the propagation constant β even , and β odd of the two super-modes by the following exact relations:

β = β even + β odd 2 , (1)
χ = β even -β odd 2 . (2)
In the limiting case of a vanishing gap (g = 0), the even and odd super-modes of the coupler progressively converge toward the fundamental (even-symmetry) and first order (odd-symmetry) modes of a two-mode waveguide, as illustrated in Fig. 2(B). This allows us to say that a two-mode waveguide can also be thought of as a system of two tightly coupled single-mode waveguides, in which the gap distance has collapsed to zero [START_REF] Boucher | Analytical model for the coupling constant of a directional coupler in terms of slab waveguides[END_REF]. Moreover, from the knowledge of the propagation constants of the even and odd modes (the T E 0 and T E 1 ones in the SOI case), and by exploiting Eqs. ( 1) and (2), it is straightforward to establish an equivalent description for the two-mode waveguide involving only two single-mode structures, as will be detailed in the subsequent sections.

Scattering matrix model of the coupling section

By exploiting the formal equivalence outlined in the previous section, the coupling zone between ring and two-mode waveguide of the initial configuration, as depicted in Fig. 3(A), can be substituted with an asymmetric ternary coupler between three single-mode waveguides, as represented in Fig. 3(B). In this equivalent framework, W G 1 is the waveguide belonging to the microring (B) g=0

WG 1 + WG 2 (A) WG 1 WG 2 g>0
even supermode (β even ) odd supermode (β odd ) β β χ Fig. 2 The even (blue) and odd (red) supermodes of a symmetric coupler between two single-mode waveguides (A) characterized by the same propagation constant β and exchanging power with a coupling coefficient χ, converge towards the even and odd modes of a two-mode waveguide (B) as the gap g progressively vanishes. Conversely, any two mode waveguide can therefore be interpreted as a symmetric single-mode coupler with g = 0.

resonator, while W G 2 and W G 3 are the two new single-mode waveguides resulting from the decomposition of the original two-mode section. According to Eqs. ( 1) and ( 2), the common propagation constant of both W G 2 and W G 3 is β = (β even + β odd ) /2 and their mutual coupling coefficient is χ b = (β even -β odd ) /2, β even and β odd being the propagation constants of the even and odd modes of the initial two-mode waveguide, respectively. Moreover, for the sake of simplicity, we assume perfect phase matching between the microring and the decomposed single-mode waveguides, so W G 1 shares the same propagation constant β. The interaction between the microring resonator (W G 1 ) and W G 2 is accounted for by a second coupling coefficient χ a .

We further assume that coupling takes place only between neighboring waveguides (i.e. no direct coupling between W G 1 and W G 3 ) 1 . We define F 1 , F 2 and F 3 as the waves propagating in the three waveguides, with:

F 1 (0) = a 1 (3a) F 2 (0) = a 2 (3b) F 3 (0) = a 3 , (3c)
and

F 1 (d) = b 4 (4a) F 2 (d) = b 5 (4b) F 3 (d) = b 6 , (4c)
so that a n (respectively b n ) represents the complex amplitude of the wave entering into (respectively coming out of) the n-th port of the ternary coupler of length d. The propagation of F n waves within the ternary coupler (assuming the W G 1 ring waveguide to be open) obeys the following matrix equation:

i ∂ ∂z   F 1 F 2 F 3   =   β χ a 0 χ a β χ b 0 χ b β     F 1 F 2 F 3   . (5)
It is more convenient to work in the framework of the slowly-varying envelopes

A n , with F n = A n exp(-iβz) and to define two new parameters χ = χ 2 a + χ 2 b
and α, such as:

χ a = χ sin (α) (6a)
χ b = χ cos (α) . (6b)
Since the coupling strength between the microring and the two-mode waveguide is weaker than the coupling strength between the two propagating modes inside the two-mode waveguide (in the even-odd supermodes framework), we can assume χ a ≤ χ b with α ∈ [0, π/4]. Eq. ( 5) can be rewritten under the following form:

i ∂ ∂z   A 1 A 2 A 3   = χ   0 sin α 0 sin α 0 cos α 0 cos α 0     A 1 A 2 A 3   (7)
Taking Eq. ( 7) as a starting point, we follow a standard diagonalization procedure of the operator to derive its eigenvalues and eigenvectors. After some

1 By means of simple analytical evaluations made on slab waveguides, one can quickly verify the validity of this assumption, as long as W G 1 remains reasonably far from the system made by W G 2 + W G 3 . In fact, given the transverse power distribution of the modes in this configuration, the overlap integral between W G 1 and W G 3 (which is roughly proportional to the coupling constant) appears to be vanishingly smaller than the one between W G 1 and W G 2 . calculations, we obtain eventually the following scattering matrix representation, which relates the complex field amplitudes a n and b n of the ternary coupler on the initial (z = 0) and final (z = d) sections, respectively:

  b 4 b 5 b 6   = e -iβd   ρ -i κ X -i κ τ -i κ ′ X -i κ ′ ρ ′     a 1 a 2 a 3   , (8)
where we define the following parameters:

θ = χd, (9a)
ρ = sin 2 α cos θ + cos 2 α, (9b)
κ = sin α sin θ, (9c) X = sin α cos α (cos θ -1), (9d)
τ = cos θ, (9e)
ρ ′ = cos 2 α cos θ + sin 2 α, (9f)
κ ′ = cos α sin θ, (9g)
Beside this, it can also be checked that 1 -ρ = ρ ′ -τ = sin 2 α (1 -cos θ), which will be useful in the calculations of Section 5. It must be observed that an equivalent matrix expression for the coupler could have been derived starting from a different basis set, namely the one explicitly dealing with even and odd mode amplitudes in the two-mode waveguide. However, our reference system made of single mode waveguides makes the ternary coupler simpler to model, as detailed in Appendix (A) .

Modal transfer function of the system

By operating a substitution of the two-mode waveguide with its equivalent single-mode representation, we have derived the scattering matrix (8) which describes the mode-amplitude evolution within the ternary coupler, assuming the W G 1 waveguide to be open. A matrix expression accounting also for the presence of the microring (W G 1 closed on itself) can be derived by observing that:

a 1 = b 4 γe -iβ(L-d) , ( 10
)
where γ is the propagation loss factor over the distance (L -d).

By introducing Eq. ( 10) in Eq. ( 8), we obtain the following (2 × 2) matrix expression:

b 5 b 6 = A -i B -i B C a 2 a 3 , (11)
which links directly the (a 1 , a 2 ) and (b 5 , b 6 ) mode amplitudes. The elements in 11 are:

A = e -iβd τ -κ 2 γe -iβL 1 -ρ γe -iβL , (12a)
B = e -iβd κ ′ + κX γe -iβL 1 -ρ γe -iβL , (12b)
C = e -iβd ρ ′ + X 2 γe -iβL 1 -ρ γe -iβL . ( 12c
)
For the sake of clarity, we stress again that Eq. ( 11) relates the field amplitude on the z = 0 and z = d sections of the ternary coupler, when W G 1 is closed on itself (closed loop) in the framework where the two-mode waveguide is interpreted as a system of two tightly coupled single-mode waveguides. In order to switch back to the initial even-odd mode amplitude representation, the following rotation matrices can be straightforwardly applied to Eq. ( 11):

b even b odd = √ 2 2 1 1 -1 1 b 5 b 6 , (13a)
a 2 a 3 = √ 2 2 1 -1 1 1 a even a odd , (13b)
leading to the matrix expression:

b even b odd = [J] a even a odd = 1 2 A + C -2i B -A + C -A + C A + C + 2i B a even a odd , (14)
which relates the even/odd mode amplitudes between the z = 0 and z = d sections of the device. The explicit form for the elements of [J] reads as:

J 11 = e -iβd 2 ρ ′ + τ -2i κ ′ + γe -iβL D (X -i κ) 2 , (15a)
J 12 = J 21 = e -iβd 2 ρ ′ -τ + γe -iβL D (X 2 + κ 2 ) , (15b)
J 22 = e -iβd 2 ρ ′ + τ + 2i κ ′ + γe -iβL D (X + i κ) 2 , (15c)
where D = 1 -ρ γe -iβL is a common denominator factor characterizing the microring. As the elements of the matrix J depend on the spectral parameter β, the degree of mode mixing produced by the system depends intrinsically on the working wavelength. In the specific case of a transparent microring (γ = 1) and for wavelengths corresponding to one of the resonances of the microring [exp(-iβL) = 1], since D = 1 -ρ the elements of J reduce to:

J 11 = J 22 = 0, (16a)
J 12 = J 21 = e -iβd . (16b)
The matrix J has therefore an anti-diagonal structure, thus enabling a complete even-odd (odd-even) cross-coupling between the modes as they transit through the system. More explicitly, an even (odd) mode incident on the input port (section z = 0) is transmitted to the output port (section z = d) as an odd (even) one. Conversely, on the anti-resonance wavelengths, for which [exp(-iβL) = -1], the vanishing terms of the matrix J are the anti-diagonal ones (J 12 = J 21 = 0) and any incident modal state is transmitted unchanged as a result of the diagonal shape of J. In Fig. ( 4) are presented the evolutions of the b even (thick blue line) and b odd (thin red line) output mode amplitudes calculated directly through Eq. ( 14) as a function of the parameter βL in the range [-3π, 3π] in the case of a pure even mode excitation on the input port (a even = 1, a odd = 0). As expected, on each of the microring resonances [where (βL = 2qπ), with q integer], the even mode power on the input port is completely converted to the odd mode on the output port (b odd = 1) while on the anti-resonances [where (βL = (2q + 1)π)] no conversion takes place. Moreover, between resonance and anti-resonance conditions, an intermediate degree of mode mixing is achievable. We also observe that, following the assumption of a transparent media, the total power carried by even and odd modes is conserved through the passage within the system

(|a even | 2 + |a odd | 2 = |b even | 2 + |b odd | 2 = const)
, the only change concerns the modal repartition.

In view of the forthcoming nonlinear modelling, it is useful to derive a closed expression for the field amplitude b 4 inside the resonator as a function of the mode amplitudes a 2 and a 3 :

b 4 = e -iβd -i κ a 2 + X a 3 1 -ρ γe -iβL , (17)
or, alternatively, as a function of the even/odd input mode amplitudes:

b 4 = a even √ 2 2 e -iβd X -i κ 1 -ρ γe -iβL + a odd √ 2 2 e -iβd X + i κ 1 -ρ γe -iβL , ( 18
)
The power carried by the system depends on the possible coherence between the even and odd modes. If these are mutually incoherent, then the internal power is reduced to the mere sum of the even and odd contributions. Moreover, as a result of the phase quadrature between factors X and iκ, both even and odd modes are equally coupled into the ring.

6 Nonlinear regime

Phase bistability and mode conversion

In the previous sections of this work, we have derived a matrix expression relating the amplitudes of the even/odd modes at the input-output ports of the device and we have shown that, under particular conditions, this system can act as an ideal mode converter. We now extend our discussion by considering a microring resonator presenting a nonlinear response of the Kerr type, as schematically shown in Fig. 5. We assume (i) the nonlinear zone perfectly matched at both ends, so that the whole system remains unidirectional and (ii) that the coupling zone of length d is still linear. For the sake of clarity, we consider a strictly even input mode on the left port of the device (a odd = 0), and a purely loss-less microring (γ = 1). The Kerr medium introduces an additional nonlinear phase shift which is proportional to the field power |b 4 | 2 circulating in the microring. It is useful to introduce the reduced (dimensionless) quantity:

Y int = |b 4 | 2 P N ORM , (19)
normalized with respect to a term P N ORM chosen in such a way that the total internal phase shift ϕ in the microring (comprehensive of both linear and nonlinear contributions) could be expressed as:

ϕ = ϕ lin + 2π Y int . (20)
By defining Y e = |a even | 2 /P N ORM as the reduced input power for the (even) excitation, the power within the ring Y int reads as:

Y int = Y e K 0 A (ϕ) , (21)
where K 0 is a coupling factor that depends only on the coupler characteristics:

K 0 = |X -iκ| 2 2 (1 -ρ) 2 = 1 2 1 + ρ 1 -ρ , (22)
and A (ϕ) is a normalized comb-like Airy function

A (ϕ) = 1 1 + m sin 2 (ϕ/2) , (23)
whose spectral selectivity depends on the parameter m m = 4ρ

(1 -ρ) 2 . ( 24
)
We would like to emphasize that, for a given reduced input power Y e , neither the microring internal power Y int nor the total phase shift ϕ are known quantities, this being the nonlinear nature of the problem. To determine Y int and ϕ in a self-consistent way, we need to suitably rewrite the expressions ( 20) and ( 21) in order to put in evidence their simultaneous dependence on the unknown parameter ϕ:

Y int Y e = K 0 A (ϕ) , (25a)
Y int Y e = ϕ -ϕ lin 2π Y e . ( 25b
)
The right-hand term of Eq. (25a) has the form of a comb-like Airy function, while the right-hand term of Eq. (25b) describes a family of straight lines whose slope is inversely proportional to Y e . By equating (25a) and (25b) we get the following expression in the unknown total phase shift ϕ:

K 0 1 + m sin 2 (ϕ/2) = ϕ -ϕ lin 2πY e , (26)
whose solution(s) as a function of Y e can be found either graphically [START_REF] Boyd | Nonlinear Optics[END_REF] or numerically. Moreover, if the microring resonances are sharp enough to be assimilated to Lorentzian functions, these solutions can be also derived analytically as the roots of a third order algebraic equation (see Appendix B). The system behavior as a function of the reduced input power Y e can be assessed by referring to Fig. 6, where the left-hand side (LHS) and the right-hand side (RHS) of Eq. ( 26) are plotted separately against the unknown total phase shift ϕ. Intersection points between RHS and LHS represent the graphical solutions of Eq. ( 26) and their projections of on the x-axis give directly the total phase shift values associated to the different input powers Y e . For lower input power (Y e = Y 1 ), the RHS straight line intersects the LHS curve only in (a) and its projection on the x-axis identifies the total phase shift of the microring in the quasi linear regime (ϕ = ϕ lin ). When the input power Y e is increased (Y e = Y 2 ), the slope of the RHS line decreases progressively until intersecting the LHS curve simultaneously in (b) and (c). In this conditions Eq. ( 26) admits 26). When the reduced input power Ye follows the path highlighted by the red dashed arrows, the system evolves from the working point (a) (microring off-resonance, no mode conversion) to the working point (c) (microring on resonance, complete mode conversion).

The system parameters are the same as in Fig. 6 two possible solutions which entail a bistable behavior. A further input power increase (Y e = Y 3 ) makes the system evolve toward further solution points (d, e) thus tracing a complete hysteresis loop. It is now important to observe that the intersection point (c) lies exactly on the peak of the RHS function, and that the total phase shift associated to this solution is ϕ = 2π. In this specific working point the microring is in a resonance condition and the matrix J has, consequently, a diagonal structure which enables a complete mode conversion. A possible functionalization of such a system is clearly a nonlinear mode switch (or modal commutator). In fact, if the linear working point is initially set to an anti-resonance of the microring (ϕ = ϕ lin = π), as the input power Y e of the mode is increased, the system is led from an initial condition where no mode conversion takes place (J in diagonal form), to a condition of complete mode conversion (J in anti-diagonal form). More explicitly, an even (odd) mode tuned to an anti-resonance in the linear limit, would switch to an odd (even) mode as the system is progressively driven through its hysteresis cycle. The bistable nature of the system is clearly revealed in Fig. 7, where the values of the total phase shift ϕ (blue solid line) are traced for increasing values of the reduced (even-mode) input power Y e . On the same curve, also the solution points (a),(b),(c),(d) and (e) previously discussed in Fig. 6 are reported.

As one can note, by following with the input power Y e the path highlighted by the red dashed line, the system can be led from the point (a), corresponding to a condition of mode decoupling (no mode conversion), to the point (c) where the mode conversion is complete, thus confirming the potential applicability of this configuration as a power-controlled mode switcher/converter.

Relation between Q-factor and mode-conversion power level

For a practical implementation of the device, it is useful to assess the dependence between the microring quality factor Q, which finally quantifies the maximum possible bandwidth of the signal to be processed, and the input power level Y e of the point (c) where mode conversion takes place. By observing Fig. [START_REF] Dumeige | Dispersive tristability in microring resonators[END_REF], one can note that the (c) point has Cartesian coordinates (x, y) = (2qπ, K 0 ) with q integer, and that on each of the microring resonances (for which ϕ = 2qπ), the value of the corresponding reduced input power can be directly derived by simple inversion of Eq. ( 26):

Y e = 2qπ -ϕ lin K 0 = (4qπ -2ϕ lin ) 1 -ρ 1 + ρ , (27)
in which the explicit dependence from the ρ parameter is obtained through Eq. ( 22).

As long as m is high enough to have well defined Lorentzian lineshapes, also the quality factor Q of the resonator is linked to the coupling coefficient ρ through the following relation:

Q = ϕ 0 2 √ ρ 1 -ρ = qπ √ ρ 1 -ρ , (28)
where ϕ 0 = 2qπ is the actual phase of the resonance peak chosen as working point.

It should now be remembered that ρ represents the fraction of field which is reinjected in the ring at each loop (Eq. 8), therefore, an increase (decrease) of ρ produces an increase (decrease) of the K 0 and Q parameters via Eqs. ( 22),(28) with a corresponding decrease (increase) through Eq. ( 27) of the Y e reduced point corresponding to the (c) point, as graphically evidenced in Fig. [START_REF] Matsko | Optical hyperparametric oscillations in a whispering-gallery-mode resonator: Threshold and phase diffusion[END_REF]. As the ρ coefficient approaches progressively to unity, the quality factor Q grows towards infinity with a consequent enhancement of the nonlinear effect at the expenses of a reduced available bandwidth for signal transmission. This dynamic is therefore typical of resonant systems encompassing third order (Kerr) effects, as already evidenced in literature at both analytical and experimental levels [START_REF] Rasoloniaina | Analysis of third-order nonlinearity effects in very high-q wgm resonator cavity ringdown spec-Yann G[END_REF], [START_REF] Treussart | Evidence for intrinsic kerr bistability of high-q microsphere resonators in superfluid helium[END_REF], [START_REF] Matsko | Optical hyperparametric oscillations in a whispering-gallery-mode resonator: Threshold and phase diffusion[END_REF].

Conclusions

In 8 Evolution of microring quality factor Q (red solid line) and reduced power level Ye (blue solid line) enabling mode conversion, as a function of the coupling parameter ρ. The parameters of the system are the same as in Fig. 6.

and nonlinear Kerr regimes. The presence of the microring breaks the intrinsic even/odd symmetry of the two-mode waveguide, thus enabling mode mixing. In the loss-less case, energy can be totally transferred from even to odd modes and vice versa, provided to operate on a resonance wavelength of the microring. The system can therefore act as a mode converter. We have also derived the analytical expression for the internal field of the microring and, when Kerr nonlinearity is effective, we have demonstrated mode-switching associated with ring phase bistability. One key hypothesis of our model relies on the existence of a phase matching condition between the microring propagation constant β and the average value (β even + β odd )/2 of the propagation constants of the even and odd modes in the two-mode waveguide. This condition gives rise to a particularly simple and compact form for the ternary coupler scattering matrix. When relaxing this constraint on β, a scattering matrix representation for the ternary coupler can be still derived analytically, but the detailed calculation proves to be too cumbersome to be presented in the frame of this work. Nevertheless, the main symmetries of the matrix are preserved, so that the main conclusions remain essentially the same, except for a slight shift of the resonant frequencies of the system.

A Considerations on the choice of the modal basis

In Section 4, the matrix expression for the coupler has been obtained by treating the two-mode waveguide (TMWG) as a system of mutually coupled single-mode waveguides (SMWG), thus obtaining a description of the field evolutions in the fully single-mode basis {e 1 , e 2 , e 3 }, where en denotes the unit (normalized) vector associated to the mode propagating in the n-th single-mode channel. A possible alternative basis is the one using {e 1 , eeven, e odd }, in which eeven (resp. e odd ) denotes the even (resp. odd) amplitudes of the modes in the two-mode waveguide (see Fig. 3). Within this alternative hybrid single/multimode framework, the evolution equation takes the following form:

i ∂ ∂z   F 1 Feven F odd   =   β χ 1e χ 1o χ 1e βeven 0 χ 1o 0 β odd     F 1 Feven F odd   , (29)
where β denotes, as previously, the propagation constant of the upper SMWG, and βeven (resp. β odd ) stands for the propagation constant of the even (resp. odd) mode in the lower TMWG. As far as coupling constants are concerned, we need an expression for both χ 1e (between mode 1 and mode even) and χ 1o (between mode 1 and mode odd). This means two degrees of freedom for a complete description of the coupler. Besides, thanks to their intrinsic orthogonality, no direct coupling between the even and odd modes needs to be accounted for.

On the other hand, within the framework exploited in Section (4), the TMWG is treated as a system of two mutually coupled SMWG interacting with the upper SMWG, the coupling being always limited to the nearest neighbour(s). For the sake of algebraic simplicity, the three SMWG are assumed identical, sharing therefore the same propagation constant β. The upper (resp. lower) coupling constant is denoted χa (resp. χ b ). In this basis, according to Eq. ( 5):

i ∂ ∂z   F 1 F 2 F 3   =   β χa 0 χa β χ b 0 χ b β     F 1 F 2 F 3   . (30)
The complex amplitudes of the modes can be transformed from one basis into the other by means of the following rotation matrices:

  F 1 Feven F odd   =   1 0 0 0 √ 2/2 √ 2/2 0 - √ 2/2 √ 2/2     F 1 F 2 F 3   ,   F 1 F 2 F 3   =   1 0 0 0 √ 2/2 - √ 2/2 0 √ 2/2 √ 2/2     F 1 Feven F odd   .
(31) As a result, we get eventually:

i ∂ ∂z   F 1 Feven F odd   =   β χa/ √ 2 -χa/ √ 2 χa/ √ 2 β + χ b 0 -χa/ √ 2 0 β -χ b     F 1 Feven F odd   . (32)
Eq. (32) should be now compared with Eq. (30). Note that the three modes {e 1 , eeven, e odd } do not propagate at the same velocity, and that the coupling between e 1 and eeven has the same magnitude (but not the same sign) as that between e 1 and e odd . Not only do we recover the (rather obvious) fact that βeven = β + χ b and βeven = β -χ b : by identification, we also determine a compact expression for the coupling constants χ 1e = -χ 1o = χa/ √ 2. In other words, the approach dealing with single mode waveguides is simpler, as one less degree of freedom is required. Moreover, a multimode propagation problem looks more easily tractable and understandable in terms of coupled single-mode channels, thus giving an alternative picture with an immediate physical insight.

B Intersection of a Lorentzian function with a straight line

If the microring resonances are sharp enough to be assimilated to Lorentzian functions, then, with a suitable variable substitution, Eq. ( 26) can be recast into the following equivalent Yann G. Boucher et al.

form: 1 x 2 + 1 = x + a K ⇔ x 2 + 1 (x + a) = K, (33)
modelling the intersection of a Lorentzian curve with a straight line of horizontal intercept (anchor point) a and its inverse slope K, as depicted in Fig. 9. Eq. ( 33) is polynomial equation of the third-order, which is algebraically solvable by means of the Cardano's formula, and whose roots are parametrized by a and K. By operating the substitutions: b = (a/3), z = x + b, p ′ = -b 2 + (1/3) and q ′ = b 3 + b -(K/2), Eq. ( 33) can be recast into its reduced cubic form:

z 3 + 3p ′ z + 2q ′ = 0, (34)
which allows to investigate the nature of the roots more easily. More specifically, Eq. (33) admits three real (physical) solutions only if the factor:

R = (p ′ ) 3 + (q ′ ) 2 = K ′ 2 -2K ′ b b 2 + 1 + 18b 2 + 1 27 , K ′ = K/2 (35)
is strictly negative (R < 0). The condition on the negative sign for R is fulfilled only if p ′ < 0, that is if a > √ 3. Therefore, in order to have R < 0 and thus three real solutions, the anchor point a of the straight line must be located well outside the Lorentzian line shape. This first condition on a being assured, the factor R is negative for values of the reduced slope K ′ lying within the two roots of the polynomial (35). The explicit values for these roots can be evaluated through the reduced discriminant:

∆ ′ = b 2 b 2 + 1 - 18b 2 + 1 27 , (36)
which is a quantity unconditionally positive since b > (1/3). The two roots K ′ 1,2 of (35) are finally:

K ′ 1,2 = b b 2 + 1 ± √ ∆ ′ .
To summarize, for a given detuning a > (3), the bistability range of the system described by Eq. ( 33) corresponds to the zone where the reduced slope K ′ lies between the two roots K 1,2 of (35), where R < 0.

Fig. 1

 1 Fig. 1 Schematic representation of the system under investigation, consisting of a twomode waveguide (W G) and a single-mode racetrack microring resonator (M RR) of total perimeter L, which are evanescent coupled over a distance d.

Fig. 3

 3 Fig.3The initial configuration of the system (A) is described by the even and odd mode amplitudes at the input (aeven, a odd ) and output (beven, b odd ) ports. After substitution of the two-mode waveguide with an equivalent coupler between two single-mode ones, the interaction zone highlighted by the green-dotted rectangle is modelled as a ternary directional coupler between three single-mode waveguides (W G 1 ,W G 2 and W G 3 ) which share the same propagation constant β. Coupling is supposed to take place only between neighboring waveguides with strength χa between W G 1 and W G 2 and strength χ b between W G 2 and W G 3 .

2 Fig. 4 5

 24 Fig.4Even (blue thick solid line) and odd (red thin solid line) output mode evolution versus βL, for system parameters χa = χ b = 0.6, d = 1, L = 20, aeven = 1 and a odd = 0 (i.e. pure even excitation on the input port). The dashed black line represents the total (conserved) modal power.

Fig. 6

 6 Fig. 6 RHS (red solid line) and LHS (green solid line) of Eq. (26) traced against the nonlinear phase shift ϕ. Here the parameters of the system are: χa=χ b =0.6, d = 1, L = 20, aeven = 1 and a odd = 0.

Fig. 7 (

 7 Fig. 7 (Blue solid line) Evolution of the microring total phase shift ϕ against the (purely even) reduced input power Ye. The points of the curve are obtained by numerically solving Eq. (26). When the reduced input power Ye follows the path highlighted by the red dashed arrows, the system evolves from the working point (a) (microring off-resonance, no mode conversion) to the working point (c) (microring on resonance, complete mode conversion). The system parameters are the same as in Fig.6

2 Fig. 9

 29 Fig. 9 Intersection of a straight line of variable slope K and anchor point a with a Lorentzian curve. The vertical blue dashed line correspond to the linear regime, whereas the two green-dashed lines delimit the range of bistability inside which two or three distinct real solutions may exist.

Acknowledgements This work is supported by the Labex CominLabs (ANR-10-LABX-07-01) through the 3D-Optical-ManyCores project: (http://www.3d-opt-many-cores.cominlabs.ueb.eu/). The authors would like to thank Prof. Christophe Peucheret for helpful discussions.


Vania Conan

Jérémie Leguay

Timur Friedman

The heterogeneity of inter-contact time distributions: its importance for routing in delay tolerant networks

Prior work on routing in delay tolerant networks (DTNs) has commonly made the assumption that each pair of nodes shares the same inter-contact time distribution as every other pair. The main argument in this paper is that researchers should also be looking at heterogeneous inter-contact time distributions. We demonstrate the presence of such heterogeneity in the often-used Dartmouth Wi-Fi data set. We also show that DTN routing can benefit from knowing these distributions. We first introduce a new stochastic model focusing on the intercontact time distributions between all pairs of nodes, which we validate on real connectivity patterns. We then analytically derive the mean delivery time for a bundle of information traversing the network for simple single copy routing schemes. The purpose is to examine the theoretic impact of heterogeneous inter-contact time distributions. Finally, we show that we can exploit this user diversity to improve routing performance.

Introduction

In the kind of delay tolerant networks (DTNs) [START_REF] Fall | A delay-tolerant network architecture for challenged internets[END_REF] that we consider in this paper, nodes are mobile and have wireless networking capabilities. They are able to communicate with each other only when they are within transmission range. The network suffers from frequent connectivity disruptions, making the topology only intermittently and partially connected. This means that there is a very low probability that an end-to-end path exists between a given pair of nodes at a given time. Such DTNs can be consider in ad hoc networking when connectivity is very low (e.g. in tactical military communications), in transportation systems as in the DieselNet project [START_REF]UMassDieselNet: A Bus-based Disruption Tolerant Network[END_REF] or in Pocket Switched Networks (PSN) [START_REF] Chaintreau | Impact of human mobility on the design of opportunistic forwarding algorithms[END_REF] which are formed by devices that people carry everyday (cell phones, PDAs, music players). In all these contexts, end-to-end paths can exist temporarily, or may sometimes never exist, with only partial paths emerging. This paper addresses the extreme case, where only temporal paths exist. We call such networks temporal DTNs, or t-DTNs. When a node in a t-DTN receives a "bundle" of information from a neighboring node, it keeps it until it meets another node which provides an opportunity to relay the bundle. The bundle is transfered from one node to another instantly and this transfer is atomic.

Prior work on routing in t-DTNs has commonly made the assumption that each pair of nodes shares the same inter-contact time distribution as every other pair. The main argument in this paper is that researchers should also be looking at cases in which inter-contact time distributions are heterogeneous. Chaintreau et al. [START_REF] Chaintreau | Impact of human mobility on the design of opportunistic forwarding algorithms[END_REF] posit that there might be heterogeneity, but we show it and characterize it. We also show how exponential distributions can be composed to yield the heavytailed distributions that Chaintreau et al. observed. As we shall see, the heterogeneity that we highlight allows us to usefully extend the work of Spyropoulos et al. [START_REF] Spyropoulos | Single-copy routing in intermittently connected mobile networks[END_REF][START_REF] Spyropoulos | Multi-copy routing in intermittently connected mobile networks[END_REF], which analyzes numerous routing schemes for t-DTNs, but that uses mobility models that yield homogeneous distributions.

We show, on the well known Dartmouth Wi-Fi data set [START_REF] Henderson | The changing usage of a mature campus-wide wireless network[END_REF], that despite the existence of a heavy-tailed distribution when inter-contact times are considered in the aggregate, a large portion of the node pairs present inter-contact time distributions that can be well fitted by an exponential distribution. We found these distributions to be heterogeneous, with a wide variation in exponents.

We also provide the first formal analysis of the impact of heterogeneous exponential inter-contact time distributions on simple single-copy routing schemes. We show that routing strategies can benefit, in terms of delay, from this heterogeneity, and in particular from knowing these distributions. A node can choose among possible relay nodes based upon their expectations for meeting other relays or the destination.

Inter-contact time model

This section presents the model we use to analytically derive the delay expectations for the routing protocols we study later in this paper.

Exponential t-DTNs

We consider a network composed of n nodes. Let's first look at the inter-contact time between two individual nodes (i, j): t 1 i j < t 2 i j < t 3 i j < ... are the successive instants at which a contact between i and j occurs.

∆t k i j = t k+1 i j -t k i j (1)
is the inter-contact time between the k th and (k + 1) th contact instants.

We assume that the ∆t k i j are samples from independent and identically distributed random variables that follow an exponential law with parameter λ i j , which we note τ i j = exponential(λ i j ). The mean inter-contact time between i and j is thus given by E[τ i j ] = 1/λ i j .

In the overall network, all n nodes are supposed to behave independently, so that the n(n -1)/2 pairwise inter-contact times τ i j are independent exponential processes with differ- ent parameters. The τ i j family of processes is symmetric and ∀i, τ ii = 0. The exponential t-DTN is thus entirely and uniquely characterized by the n(n -1)/2 strictly positive real parameters λ i j .

Assumptions

The model focuses on the temporal dynamics of node connectivity in a DTN. In this way it provides a common framework to analyse different DTNs. In particular it applies very well to social networks for which the position of nodes at a given time is not of primary importance. We believe this abstraction helps focus on the inherent characteristics of intermittent connectivity in DTNs.

Characterizing inter-contact time behavior helps abstract away from the spatial information that is essential in the analysis of mobile ad hoc networks. There is no reference to geographic, localisation or any other such spatial information. There is also no reference to air interface parameters, quality of or contention on the links, etc. Node mobility is not explicitly modelled: only its aggregated impact on the inter-contact time is taken into account.

The model makes a stationarity hypothesis with respect to node inter-contact time distributions. In other words, nodes behaviors are assumed to change on a slower scale than bundle exchanges. We also suppose that nodes have infinite capacity in bandwidth and storage.

Another key hypothesis is that contacts (and thus bundle transfers) are assumed to be instantaneous. In the model, pairwise contacts do not overlap: in contrast to the mobile ad hoc network cases, no partial routes (involving more than two nodes) exist at any given time.

In this respect, the results with the proposed model are upper bounds, but, as we will see, still provide valuable information and insight on how to route bundles in DTNs. We leave refinements of the model for future work.

Fitting the model

In the t-DTN model just elaborated, we assume that the inter-contact time distribution for each pair of nodes is exponential. The main reason is that it will allow us to go beyond asymptotic results and provide explicit formulas for the bundle delivery time, and other parameters, of different routing protocols. In this section we look at real data to evaluate how reasonable this hypothesis might be.

Experimental data set

To validate the hypothesis, we use real data from the Wi-Fi access network of Dartmouth College [START_REF] Henderson | The changing usage of a mature campus-wide wireless network[END_REF]. These data track users' sessions in the wireless network by showing the time at which a node associates or dissociates from an access point. We use the subset of data pre-processed by Song et al. for their prior work [START_REF] Song | Evaluating location predictors with extensive Wi-Fi mobility data[END_REF] on mobility prediction.

As we describe in prior work [START_REF] Leguay | Evaluating mobility pattern space routing for DTNs[END_REF], we must select from the data, and make some assumptions, in order to constitute a useful DTN data set. We take the subset of users who are present in the network every day between January 26 th 2004 and March 11 th 2004, a class period during which we expect nodes' activity to be fairly stationary. This data set contains 834 users, or nodes. Then, we assume that two nodes are in contact if they are present at the same time at the same access point (AP). Finally, we filter these data to remove the well known ping-pong effect. Indeed, wireless nodes, even non-mobile, can oscillate at a high frequency between two APs. To counter this, we filter all the inter-contact times below 1,800 seconds. Note that defining better filtering methods, albeit challenging, would be of interest for the community. As this is not the purpose of this work, We choose here the threshold that Yoon et al. [START_REF] Yoon | Building realistic mobility models from coarse-grained traces[END_REF] used for the same purpose. We use this new data set for the remainder of this paper.

The Wi-Fi scenario may be not a perfect fit for interactions between nodes in t-DTNs. Indeed, in opposition to always-on devices carried by humans, Wi-Fi nodes are typically turned off, transported, and then turned on again, thus missing potential contacts en route. However, the size, quality, and public availability of the data set make it nonetheless one of the best resources for this kind of study. Jones et al. [START_REF] Jones | Practical routing in delaytolerant networks[END_REF] and Chaintreau et al. [START_REF] Chaintreau | Impact of human mobility on the design of opportunistic forwarding algorithms[END_REF] recently used these traces in a similar way.

Exponential inter-contacts

Fig. 1 shows the distribution of E[τ i j ], the expected intercontact time for the pair of nodes (i, j). This plot has been computed over all the 28, 490 source-destination pairs that experienced an average inter-contact time lower than one week within the two months period that we considered in Dartmouth data. We can see that the distributions are heterogeneous with expectations varying over three orders of magnitude. The average E[τ] is 11.6 hours with a variance of 7.1 hours. Then, we test for whether the inter-contact process between any two nodes can be modelled by an exponential process with a parameter λ = 1/E[τ]. We use the Cramer-Smirnov-Von-Mises [START_REF] Eadie | Statistical Methods in Experimental Physics[END_REF] hypothesis test. For each pair (i, j), we compare the cumulative distribution I i j N for the N inter-contacts observed and the hypothesis function whose cumulative distribution is F i j (x) = 1exp(-λ i j x). We also compare I i j N with that of a power law distribution. Note that we only perform the computation for pairs that show a sufficient level of connectivity by having a mean inter-contact time lower than one week and that have more than 20 contacts. We identify 8,402 pairs to be exponentially distributed and 28 with a power law which makes respectively, 62.3% and 0.2% of the 13, 482 pairs that we retain for the test.

From these observations, it seems clearly more reasonable, in this data set, to model pairwise inter-contact time distributions as exponential rather than power law since a large number of pairs have shown inter-contact times exponentially distributed. Despite that few are power law distributed, we conjecture that the rest of pairs might follow distributions that are a mix of exponential and power law distributions. As we have examined only one data set, albeit an often-used one, we cannot draw many conclusions about what will be revealed elsewhere. It is reasonable to expect that other mobility traces in campus environments will show similar characteristics. However, it is surprising that a memoryless process seems to be at work in such a high proportion of node pairs in an environment in which one would expect some temporal correlations. We hope this will be a spur to study these distributions in other data sets. Traces from the Haggle project [START_REF] Chaintreau | Impact of human mobility on the design of opportunistic forwarding algorithms[END_REF] or the ones of the Reality Mining project [START_REF] Eagle | Social serendipity: Mobilizing social software[END_REF] might be considered. We let this study for future work.

Power laws

Chaintreau et al. [START_REF] Chaintreau | Impact of human mobility on the design of opportunistic forwarding algorithms[END_REF] observed that aggregated inter-contact times follow power laws in a number of DTN traces (also including one based on the Dartmouth data). Fig. 2 shows that, for our data set, the cumulative distribution of aggregated inter-contact times also follows a power law of the form f (x) = cx δ , with exponent δ = -0.16 and scale parameter c = 3.45. Let's now consider what happens for pure exponential t-DTNs where all pairwise inter-contact time distributions are exponentially distributed. Under which conditions do the aggregated inter-contact time distributions follow a power law, or is the pairwise exponential assumption too strong to yield a power law in the aggregate?

Let Θ be the aggregated inter-contact time for all pairs of nodes, and let p(λ ) be the probability distribution of the λ parameters:

P(Θ > t) = ∞ λ =0 e -λ t p(λ )dλ (2)
What eqn. 2 says is that, for exponential t-DTNs, the aggregated inter-contact time distribution is fully characterized by the distributions of the λ parameters, and thus of the E[τ i j ] matrix. More precisely, the tail cumulative distribution of the aggregated inter-contact times is given by the Laplace transform of the distribution p of the λ parameters.

A Pareto law of the form ( b t+b ) α , t ≥ 0, with shape parameter α > 0 and scale parameter b > 0, is observed if and only if the λ follow a gamma distribution p(λ

) = λ α-1 b α e -bλ Γ (α)
, λ ≥ 0.

To verify this on the data set we proceed in the following way: we estimate parameters α and b from the cumulative distribution of the λ parameters for pairs that were shown to follow an exponential behavior (the ones that pass the Cramer hypothesis test). We find b = 113, 766.9 and α = 2.26. Fig. 3(a) shows the estimated cumulative gamma distribution g(x) with the experimental lambda cumulative distribution for all pairs that have shown to be exponential. Then, we plot in Fig. 3(b) the corresponding power-law h(t) with cumulative distribution of aggregated inter-contact times. As one can see, the two experimental curves fit the theoretical curves. What this result shows is that when one considers an exponential t-DTN, we can regain the power law behavior for the aggregated inter-contacts when the distribution of the parameters is a gamma, which is the case in the data we used when considering the subset of pairs that have inter-contact times exponentially distributed.

Single copy routing strategies

Having defined a stochastic model that is realistic for the data set under study, we now examine different simple single copy routing strategies. We derive analytical formulas that we will use to study the impact of heterogeneous λ i j parameters on routing.

In all routing strategies, we consider that nodes know all pairwise mean inter-contact times for all nodes in the network, i.e., each node knows the λ i j matrix. This knowledge could be diffused through an epidemic type of routing, or learned by each node from past contacts.

Wait scheme

Under the Wait routing strategy, the source node s waits until it meets d, the destination, to deliver the bundle in one hop.

If the bundle is injected at time t, its delivery time is equal to R t sd , the remaining inter-contact time before the next contact between nodes s and d. The memoryless nature of exponentials implies that R t sd also follows an exponential distribution with the same parameter. The mean expected delivery is thus given by:

E[D w sd ] = 1/λ sd (3)
This straightforward result gives an upper bound on the delivery time that a routing strategy should meet, since the Wait strategy is the most rudimentary one hop single copy scheme.

MED

The Minimum Expected Delay (MED) routing strategy was first introduced by Jain et al. [START_REF] Jain | Routing in a delay tolerant network[END_REF]. This strategy, similar to source routing, defines which path the bundle will follow from s to d, that is, the ordered list of intermediate relay nodes it will have to go through. The list is chosen to provide minimum expected end-to-end delay.

If a path is given by the following ordered list of nodes r 0 = s < r 1 < r 2 < r 3 < ... < r n-1 < r n = d, and relaying occurs at time instants t 1 < t 2 < ... < t n , the total delivery time along path (s, r 1 , r 2 , ..., r n-1 , d) is given by the remaining intercontact time after each relaying instant, that is:

D med s,r 1 ,r 2 ,...,r n-1 ,d = R t 1 sr 1 + R t 2 r 1 r 2 + ... + R t n r n-1 d (4)
Using the fact that E[R r i r j ] = 1/λ i j , the expected delivery time along the path is thus given by: E[D med s,r 1 ,r 2 ,...,r n-1 ,d ] = 1/λ sr 1 + 1/λ r 1 r 2 + ...1/λ r n-1 d (5) Finding the optimal path thus amounts to finding a lowestweight path between nodes s and d in a graph in which the weight on each link (i, j) is defined as 1/λ i j . Dijkstra's algorithm can be used.

Spray and Wait routing

The Spray and Wait strategy was first introduced by Grossglauser and Tse [START_REF] Grossglauser | Mobility increases the capacity of ad-hoc wireless networks[END_REF], and is designed to take advantage of opportunistic contacts. It consists of two steps. First the source node uses the first nodes encountered as relays to the destination. This is the "spraying" step. A relay node then uses the "wait" strategy to relay the bundle, i.e. it waits until it meets the destination to deliver the bundle. Here, we study the case where only one relay is used, which we designate 1-SW.

Let us first consider the spraying step. The bundle is injected at source s at time instant t. The first node r it encounters may be any of the n -1 other nodes d, r 1 , r 2 , ..., r n-2 and the time X it takes to meet this first node is the infinum of the inter-contact times with all other nodes:

X = in f (R t sd , R t sr 1 , ..., R t sr n-2 ) (6)
Since all R t sr i are independent exponentials with parameters λ sr i , we have (see [15, p.328]):

-The random index r of the first node encountered is independent of the first encounter time X -X is exponentially distributed, with parameter:

Λ s = λ sd + ∑ n-2 i=1 (λ sr i ) -Pr(First node encountered is r) = λ sr Λ s
This means that we can represent the spraying step as independently identifying the encountered node (with probability λ sr Λ s ) and adding an exponential waiting time with parameter Λ s .

Two cases may arise: either the first node encountered r equals d, and s delivers the bundle, or r = d and node r waits to meet node d to deliver the bundle.

The delivery time Z d , when node d is encountered first is thus given by:

E[Z d ] = 1 Λ s (7)
The delivery time Z r along path r, i.e., conditioned on using node r as a relay, is thus the sum of the first encounter time X and the remaining delivery time between nodes r and d, and thus:

E[Z r ] = 1 Λ s + 1 λ rd (8)
The total delivery time Z is computed by conditioning on all possible first encountered nodes d, r 1 , r 2 , ..., r n-2 , events whose probabilities are given by λ sr Λ s :

E[Z] = λ sd Λ s E[Z d ] + n-2 ∑ i=1 ( λ sr i Λ s E[Z r i ]) (9)
After simplification, we can state that: in a network composed of n nodes, 1-SW delivers a bundle from source s to destination d with mean delivery time given by:

E[D 1-sw sd ] = (1 + ∑ r =s,r =d λ sr λ rd ) ∑ r =s λ sr (10)

Comparing routing protocols

This section looks at routing performance of the protocols we considered in the presence of heterogeneity in inter-contact time distributions.

In this context, we present 1-SW * , a variation of 1-SW. Instead of spraying its bundle to the first node that it encounters, the source node s sprays only to nodes in a subset R. We call this a 1-SW R scheme. We define 1-SW * to be a 1-SW R scheme which uses a subset R that minimizes E[D 1-sw R sd ]. Following the same line of reasoning as in Sec. 4.3, and defining 1/λ dd = 0, one finds that the expected delivery time is given by:

E[D 1-sw R sd ] = (1 + ∑ r∈R λ sr λ rd ) ∑ r∈R λ sr ( 11
)
We performed simulations using Dartmouth traces (see Sec. 3) to study how the algorithms behave in the case of heterogeneous connectivity. We simulate the following protocols: Wait and 1-SW which are naive schemes, and, 1-SW * and MED that are designed to take advantage of heterogeneity. We slightly modified 1-SW, to better compare it with 1-SW * : a node i is a potential relay only if λ id > 0, i.e., if it has a chance of meeting the destination. In MED, we authorized intermediary relays to directly transfer bundles to the destination whenever met.

We choose at random 100 different source destination pairs (s, d) and replay the contacts between the 835 nodes present in the data to see how for each pair a bundle, generated at the beginning of the two months period, is delivered.

λ values used for route selection in 1-SW * and MED, and to determine theoretical delays of 1-SW, 1-SW * and MED have been computed over the data filtered due to the ping-pong effect (see Sec.3). However, the contacts replayed in simulations were that in the original traces as it does not impact the results and because filtering was only of interest for modelling. 1.0 ±0.0 1-SW 87.3 ±3.0 23.0 ±0.9 22.7 ±2.7 15.3 ±0.9

2.0 ±0.1 1-SW * 86.9 ±2.0 18.8 ±0. [START_REF] Spyropoulos | Single-copy routing in intermittently connected mobile networks[END_REF] 15.4 ±0.9 13.0 ±1. [START_REF] Fall | A delay-tolerant network architecture for challenged internets[END_REF] 2.0 ±0.1 MED 87.9 ±2.2 20.9 ±0.7 18.0 ±1.1 1.3 ±0. [START_REF]UMassDieselNet: A Bus-based Disruption Tolerant Network[END_REF] 7.2 ±0.2

Table 1 Simulation results with Dartmouth data.

Table 1 presents the simulation results averaged over 5 runs with the 90% confidence levels that are obtained using the Student t distribution. It presents, for each of the protocols, the average delivery ratio, the average delay ("A delay") and the median delay ("M delay") computed over the delivered bundles, the average theoretical delay over all the bundles generated (infinite delay is assumed to be the length of the simulated period, i.e. 45 days), and the average hop count, also obtained on delivered bundles.

The major result in Table 1 is that schemes that make use of heterogeneity of inter-contact times (1-SW * and MED) perform better, either in delivery ratio or delay, than the ones that do not exploit it (Wait and 1-SW). Wait only delivers 11.2% of bundles because most of the source, destination pairs selected at random satisfy λ sd = 0 (e.g. they never met). 1- SW, 1-SW * and MED achieve almost the same delivery ratios with respectively 87.3%, 86.9% and 87.9%. About 13% of the bundles were thus not delivered. In terms of delay, among these last three protocols, 1-SW plots the highest with a mean of 23.0 and a median of 22.7 in days, 1-SW * the lowest with a mean of 18.8 and a median of 15.4. MED appears to be in the mid-range with a mean of 20.9 and a median of 18.0 in days.

The difference between the modified 1-SW and 1-SW * gives a further insight on the type of heterogeneity that should be considered. The modified 1-SW is a one hop strategy that uses only true relays to the destination: relay nodes in 1-SW must meet both the source and the destination. The scheme is not completely ignorant of heterogeneity, as it exploits binary connectivity information, the fact that not all nodes meet one another. 1-SW * goes beyond that and differentiates between neighboring nodes based on the quantitative expected intercontact time. The fact that 1-SW * ourperforms the modified 1-SW thus indicates that routing actually benefits from the quantitative inter-contact time heterogeneity, and not just from node connectivity.

Table 1 shows a descrepancy between the theoretical and the experimental delays. This can be first explained by the presence of node pairs that do not have an exponential behavior. This is particularly true for 1-SW * . In this case the computation of expected delays on mean inter-contact times misses possible inter-dependencies of node contacts.

Simulation artifacts also come into play. The routing simulation is carried out on a limited time scale. The λ values are computed over the entire data set in a prior pass, so a relay node may meet the destination for the last time before having met the source for the first time. This pre-computation being not realistic, we could have used on-line predictive or learning methods. However, as they are challenging to define, we let this study for future work and intend here to provide early validation results to motivate research in the domain. The fact that 1-SW * delivers slightly less bundles than 1-SW is clearly due to this artefact. Indeed, because in 1-SW, the source transfers more rapidly the bundle to a relay, we have a lower probability to contacts between that relay and the destination. Also, MED suffers from the same simulation artefact, it would have delivered 100% of bundles otherwise.

Through these simulations, we validate the natural feeling that we should take into account the heterogeneity of intercontact time distributions in the design of routing solutions for t-DTNs. Furthermore, because 1-SW * , which is only a twohop protocol, achieves better performance than MED by delivering bundles with lower delays and a lower impact on network resources (MED delivers bundles in 7.2 hops in average while 1-SW * uses only 2), we expect promising future work inspired from 1-SW * to be done. The opportunistic nature of 1-SW * is the main reason of this superiority over MED, in which bundles follow a strict sequence of relays, in a network which is not a perfect exponential t-DTN.

Conclusion

We have first shown that, in a widely-used t-DTN data set, distributions of inter-contact times are heterogeneous. As a consequence, one has to take it into account while modeling. Second, we have validated the insight that considering heterogeneity in routing improves performance. We presented a simple routing strategy, 1-SW * , adapted from the Spray and Wait scheme, which is capable of using only a subset of relays to improve routing performance, measured in term of average delay.

Clearly, our work, based as it is upon one data set, will benefit from validation against others as mentioned in Sec. 3. Work also needs to be done to examine why a memoryless model fits so many node pairs in an environment in which one would expect to find more temporal correlations. Finally, formal studies and validations should be conducted with more elaborate schemes, in terms of number of copies distributed or in terms of the number of hops traversed.

Fig. 1

 1 Fig. 1 Distribution of E[τ].

Fig. 2

 2 Fig. 2 Distribution of inter-contacts in logarithmic scale.

Fig. 3

 3 Fig. 3 Distributions with exponential pairs.

  Wait 11.2 ±0.9 19.8 ±3.7 16.3 ±10.3 41.3 ±0.5

	delivery A delay M delay th. delay hop count
	ratio (%) (days) (days)	(days)	(hops)

Acknowledgments

We gratefully acknowledge David Kotz for enabling our use of wireless trace data from the CRAWDAD archive at Dartmouth College. We also thank Augustin Chaintreau for his valuable comments.