Vincent Laquerbe
email: vincent.laquerbe@isae.fr

Romain Pascaud

Thierry Callegari

Laurent Liard

Olivier Pascal

Modeling and Experimental Measurements of a Tunable Microstrip Resonator using Plasma Discharges

Keywords: Cold plasma, tunable resonator, microstrip, DC voltage, 3D modeling

In this paper, we suggest the use of a cold plasma as tunable material inside a microstrip resonant cavity. Plasma dielectric constant can indeed be moved to values below 1 to tune its resonant frequency. DC plasma analysis were conducted and integrated into classic electromagnetic solvers to investigate tuning abilities. Numerical simulations are consistent with experimental results and make this original tuning techniques viable for high power applications.

I. INTRODUCTION

Due to a growing demand for multiband system, frequency reconfigurability is a fundamental need for today's microwave circuits. Frequency tuning can be achieved using integrated lumped elements (PIN diodes, varactors, ...) or tunable materials (liquid crystals, ferroelectrics, ...), [START_REF] Petosa | An overview of tuning techniques for frequency-agile antennas[END_REF]. However, these solutions usually concede power handling issue.

In this context, integrated plasma discharges in microwave circuit offer a promising alternative. Plasmas can indeed withstand massive microwave power. Besides, their microwave response can be controlled to emulate different values of dielectric constant (smaller than 1 and negative). Such a versatility can offer many applications, for example power limitation [START_REF] Pascaud | Low insertion loss microplasma-based limiter integrated into a microstrip bandpass filter[END_REF], radiation [START_REF] Rayner | Physical characteristics of plasma antennas[END_REF], stealth [START_REF] Sakai | Plasma metamaterials as cloaking and nonlinear media[END_REF], or frequency tuning, which is considered here. Frequency tuning using plasma discharges has been proposed by considering plasma as a lumped microwave component [START_REF] Semnani | Plasma-enabled tuning of a resonant rf circuit[END_REF] in planar circuits. In this paper, we suggest the use of plasma as a tunable material integrated into a microstrip resonator. Provided adequate control, the plasma dielectric constant can be tuned from unity to zero and enable frequency agility of a microstrip resonator.

We computed numerical simulations to characterize the nonuniform plasma dielectric permittivity. A prototype has been developed to validate the concept of a plasma based tunable microstrip resonator. Finally, S parameters of the whole resonator were simulated and compared with experimental measurements.

II. THEORETICAL APPROACH

A. Theoretical Background

A plasma is a globally neutral ionized gas. In this study, we specifically consider non-magnetized low pressure nonequilibrium plasma (cold plasma), which means only electrons have enough energy to ionize the gas and the degree of A cold plasma can be represented by its relative permittivity ε p which is derived from the Drude model [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF]:

ε p (ω) = 1 - ω 2 p ω(ω + iν) , ω 2 p = e 2 n e m e ε 0 (1)
where ν, ω p , e, n e , m e , and ε 0 are the electron collision frequency which represents the electromagnetic losses, the plasma angular frequency, the Coulomb charge, the electron density, the mass of electron and the vacuum permittivity, respectively. From (1), it is clear that when ω p is slightly larger than ω the real part of ε p takes values between 0 and 1.

B. Design of a plasma based tunable microstrip resonator

As mentioned previously, the main principle of our tunable microstrip resonator consists in modifying the dielectric constant between its conductors. We therefore consider a microwave system composed of a 1 cm radius circular microstrip resonator whose gas-filled cavity height is 813 μm. A plasma is then ignited between the patch and the ground plane to change the dielectric permittivity inside the cavity (typically 0 < ε p < 1) and achieve frequency tuning. Note that a superstrate was added to uphold the patch above the cavity.

To minimize losses, ν must be much smaller than ω in (1). Considering for example argon, the gas pressure must be reduced to P = 2 Torr ( 2.67 mbar) in order to obtain an electron collision frequency ν ≈ 3 GHz (using BOLSIG+, [START_REF] Hagelaar | Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[END_REF]) which is ten times smaller than ω. At this pressure, DC plasmas reach medium electron densities (10 9 cm -3 < n e < 10 12 cm -3 ). Plasma ignition however requires around 800 V at 2 Torr in a 813 μm cavity, when the two microstrip conductors acts as electrodes. Hence, we suggest the solution depicted in Fig. 1 to lower the breakdown voltage.

An extra 1 cm thick step-shaped cavity and a hollow cathode were added below the microwave circuit to increase the inter electrode distance to 1.1 cm. In this case, the plasma supply current is delivered through the microstrip line and the cathode (see Fig. 1) which reduces the breakdown voltage to values between 200 V and 300 V. We also notice from Fig. 1 that a 0.5 cm radius hole is drilled through the cathode plane. The additional cavity central hole is gradually drilled from 0.5 to 1 cm to fit the resonator radius and enhance the plasma diffusion. Likewise, a slotted grid is designed to enable plasma diffusion above the microstrip ground plane. It is composed of a central hole of diameter 2 mm and four radial layers of holes. The distance between two holes is 1 mm and the angular opening of the outer holes was limited to 20 • not to excite extra microwave resonances on the ground plane. Two outer slots were also added to the patch to lower its natural frequency to 4 GHz. Given the gas pressure and the geometry from Fig. 1, plasma dielectric permittivity is indeed more flexible at these frequencies.

III. PLASMA AND MICROWAVE NUMERICAL MODELING

In the light of the difference between electron and ion mobility, a homogeneous plasma region cannot be stable: thus a sheath is necessarily created between any plasma and its surrounding medium to balance this difference. In our case, the plasma contact area is complex (mainly due to the slotted ground plane) which suggests a strong non uniform electron density distribution inside the cavity.

To understand the impact of a non uniform plasma inside the resonator cavity, simulations based on fluid dynamics, electrostatics, and Maxwellian particle distribution were conducted using our in-house software GDSim (Glow Discharge SIMulation) developed by Boeuf et al. [START_REF] Boeuf | Two-dimensional model of a capacitively coupled rf discharge and comparisons with experiments in the gaseous electronics conference reference reactor[END_REF]. The derived profiles of n e were then included in a full-wave electromagnetic simulation with commercial software Ansys HFSS.

A. Plasma ignition and diffusion numerical modeling

The simplified simulated geometry comprised the microstrip resonator electrode, the microwave cavity, the slotted ground plane, the step-shaped cavity, and the cathode. The superstrate was not included as it do not have direct contact with the plasma. Axial symmetry was assumed and the mesh surrounding the slotted grid was denser (30 μm compared to 400 μm elsewhere) to best represent the plasma diffusion and avoid overestimating the sheath thickness. The series resistance R was set to 680 Ω to allow large current variations, and thereby higher n e . Fig. 2 exhibits several electron density distributions inside both the plasma and RF cavities for different values of U DC . As expected, the electron density is not uniform between the ground plane and the microstrip line. Higher excitation voltage leads to higher electron density and thus larger expected tuning ability. We can also notice that high density plasma hardly reaches the outer holes of the ground plane; this phenomenon was also observed experimentally.

B. Full-wave simulation including plasma model

From the 2D-axisymmetric models from Fig. 2, several isocontours of electron density have been extracted. Each contour has been extrapolated to fit the non-axisymmetric ground plane (basically a local sweep around each hole) and imported into Ansys HFSS.

As the number of simulated plasma layers increases, the required mesh becomes denser and leads to very large computational resources. The plasma model was hence restricted to a five layers geometry which is depicted in Fig. 3.

Between two isocontours the plasma electron density is assumed uniform. Plasma electrical properties have been computed using [START_REF] Petosa | An overview of tuning techniques for frequency-agile antennas[END_REF] to define new materials with dispersive dielectric constant and bulk conductivity to represent plasma losses. Full-wave analysis from 3.5 GHz to 5 GHz have been conducted.

IV. NUMERICAL AND EXPERIMENTAL RESULTS

A prototype has been designed to confirm the tunability of a plasma based microstrip resonator. The whole system is presented in Figure 4. A bias tee has been mounted outside the glass jar to bring DC current to the discharge cavity without overloading the Vector Network Analyser (VNA). U DC was generated using floating potential and the ground plane was connected to earth through the VNA.

The simulated and measured S parameters are presented in Fig. 5 for different values of I DC . We first see that the measured resonant frequency without plasma, when I DC = 0 mA, is different from the simulated one which is found to be f mes = 4.26 GHz. This deviation is likely due to experimental issues such as substrates misalignments or small variations of their dielectric constant. Though, highly excited plasma, when I DC = 147 mA, leads to a 231.1 MHz measured frequency shift, which corresponds to a tuning range of 5.43 %. Note that, in this case, the measured relative insertion loss was strongly reduced by -3.53 dB at resonance.

To better characterize the frequency tuning range, different plasma configurations were considered for U DC varying from 0 V to 550 V. Fig. 6 summarizes the simulated and measured frequency shift of the plasma based microstrip resonator and the relative insertion loss at the resonant frequency for different values of input current. In both cases, a roughly linear current-dependent behaviour is observed and confirms that a simple tuning control can be achieved. However, the simulated tuning range seems limited, compared to measurements, and hardly reaches few percent. This mitigation is mainly attributed to the plasma model. First, we think that numerical results can be improved by considering more plasma layers. Besides, the 3D extrapolation likely underestimates the electron diffusion in regions right above the ground plane (between holes) which narrows the resulting plasma dielectric permittivity range inside the cavity and thus limits the resonator tuning ability. To the best of our knowledge, it is still the first time that realistic 3D plasma models were included in full-wave simulations of integrated microwave circuit.

V. CONCLUSION

In this paper, we have implemented a plasma discharge inside a microstrip resonant cavity to tune its frequency. DC plasma simulations were held to understand and characterize the electron density diffusion through a slotted plane. The resulting non uniform plasma models were used to investigate tuning abilities of such a resonator with commercial electromagnetic software. Numerical predictions are consistent with experimental results. A 5.43 % tuning range was experimentally achieved. We believe this value can be improved by designing more sophisticated diffusion grid or considering more efficient gases, such as neon or xenon.

Fig. 1 .

 1 Fig. 1. Exploded view of the considered system with DC electrical connections between electrodes. Superstrate edges were round-cut to solder SMA connector at the end of the microstrip

3 Fig. 2 .Fig. 3 .

 323 Fig. 2. Simulated electron density for different values of supply voltage, (2a) U DC = 400 V, (2b) U DC = 450 V, (2c) U DC = 500 V, (2b) U DC = 550 V. Views were zoomed for clarity and do not represent the entire simulated region.

Fig. 4 .

 4 Fig. 4. Picture of the resonator inside the glass jar with plasma (top view).

Fig. 5 .

 5 Fig. 5. Simulated (solid line) and measured (dotted line) S parameters for different values of supply voltage U DC .

Fig. 6 .

 6 Fig. 6. Simulated and measured (6a) frequency shift and (6b) relative insertion loss obtained for different values of current.

ACKNOWLEDGMENT

The authors would like to thank the French Defense Agency (DGA) for their financial support.


Ido Levy

Candice Forrester

Haiming Deng

Manuel Roldan-Gutierrez

Martha R Mccartney

David J Smith

Christophe Testelin

Lia Krusin-Elbaum

Maria C Tamargo

Manuel Goldan

Compositional Control and Optimization of Molecular Beam Epitaxial Growth of (Sb 2 Te 3 ) 1-x (MnSb 2 Te 4 ) x Magnetic Topological Insulators

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction:

Incorporation of magnetic elements into 3D topological insulators (TIs) such as Bi2Te3 and Sb2Te3 has garnered a great deal of interest in recent years due to the potential to achieve quantum anomalous Hall effect (QAHE) [START_REF] Yu | Quantized Anomalous Hall Effect in Magnetic Topological Insulators[END_REF] . QAHE was first demonstrated by incorporating Cr as a dopant in thin (Bi,Sb)2Te3 films [START_REF] Chang | Experimental Observation of the Quantum Anomalous Hall Effect in a Magnetic Topological Insulator[END_REF] . However, defects caused by the incorporation of the Cr atoms into the TI crystal lattice resulted in the effect being only evident at sub-Kelvin temperatures [START_REF] Chang | Thin Films of Magnetically Doped Topological Insulator with Carrier-Independent Long-Range Ferromagnetic Order[END_REF][START_REF] Checkelsky | Trajectory of the Anomalous Hall Effect Towards the Quantized State in a Ferromagnetic Topological Insulator[END_REF][START_REF] Mogi | Magnetic Modulation Doping in Topological Insulators Toward Higher-Temperature Quantum Anomalous Hall Effect[END_REF][START_REF] Ou | Enhancing the Quantum Anomalous Hall Effect by Magnetic Codoping in a Topological Insulator[END_REF] . Recently it was discovered that crystal growth of TIs, specifically Bi2Te3 and Sb2Te3, with Mn leads to the formation of a new crystal phase, MnBi2Te4 and MnSb2Te4, respectively, having a septuple layer (SL) structure rather than the typical quintuple layer (QL) structure of the non-magnetic TI [START_REF] Lee | Crystal Structure, Properties and Nanostructuring of a new Layered Chalcogenide Semiconductor[END_REF] . The resulting materials are intrinsically magnetic rather than magnetically doped and may lead to fewer defects. A single SL is ferromagnetic, with the magnetic spins aligned out of plane. However, when stacked to form a bulk structure, the SLs couple antiferromagnetically and are not conducive to the observation of QAHE, except in the ultrathin layers of odd number (5-7) of SLs [START_REF] Deng | Quantum Anomalous Hall Effect in Intrinsic Magnetic Topological Insulator MnBi2Te4[END_REF] . One approach used to resolve this problem was to separate the magnetic SLs layers with a few non-magnetic QLs, enabling ferromagnetic coupling between the non-adjacent SLs. The QAHE was recently demonstrated in a 1:4 MnBi2Te4:Bi2Te4 system at several Kelvin [START_REF] Deng | High-Temperature Quantum Anomalous Hall Regime in a MnBi2Te4/Bi2Te3 Superlattice[END_REF] . Another approach recently shown to produce ferromagnetic coupling of the SLs in MnSb2Te4 to introduce excess Mn into the pure SL structure. It has been proposed that MnSb antisite defects induce ferromagnetic coupling between SLs [START_REF] Riberolles | Evolution of Magnetic Interactions in Sb-Substituted MnBi2Te4[END_REF][START_REF] Liu | Site Mixing for Engineering Magnetic Topological Insulators[END_REF] .

The growth of MnSb2Te4, or of mixed SL:QL layered structures, has been reported for bulk growth techniques as well as epitaxial growth techniques [START_REF] Wimmer | Mn-Rich MnSb2Te4: A Topological Insulator with Magnetic Gap Closing at High Curie Temperatures of 45-50[END_REF][START_REF] Li | Glassy Magnetic Ground State in Layered Compound MnSb2Te4[END_REF][START_REF] Murakami | Realization of Interlayer Ferromagnetic Interaction in MnSb2Te4 Toward the Magnetic Weyl Semimetal State[END_REF] . In most cases, the formation of SLs occurs by self-assembly depending on the amount of Mn provided during growth. Typically, to grow these materials, all the elements are mixed simultaneously in order to form the crystal, and the resulting mix of SLs and QLs depends on the amount of Mn added. Although this is a convenient method to grow these materials due to its simplicity, it has some important drawbacks, particularly the difficulty to control the formation of a specific structure or composition reproducibly and consistently throughout the crystal lattice. Molecular beam epitaxy (MBE) with its layer-by-layer growth mechanism promises to provide certain benefits for control of the desired structures. However, little work has yet been done to systematically investigate and understand the MBE growth parameters that affect the structural properties of these materials. Such an understanding is essential for the intentional growth of precise layered structures with accurate compositional control and the desired magnetic properties.

In this work we report the growth of (Sb2Te3)x(MnSb2Te4)y structures by MBE. Variation of the Mn flux during growth enables control of the SL to QL ratio in the layer, although competing processes, such as Mn dopant incorporation in the QL and excess Mn in the SLs result in irreproducible and less than accurate control. Moreover, the distribution of the SLs throughout the crystal is random so that the local properties of the material vary. Modified growth techniques, involving the implementation of an annealing step in the early stages of growth, and increasing the growth temperature are explored and the results show improved crystal quality and a preferential incorporation of Mn into SLs resulting from these changes. A mechanism is proposed to explain and understand the observed improvements in the control of growth and the materials properties.

Experimental:

All samples were grown in a Riber 2300P system with base pressure of 3-5x10 -10 Torr, equipped with reflection high-energy electron diffraction (RHEED) for in-situ growth monitoring, on epi-ready c-plane (0001) sapphire substrates. The substrates were heated in vacuum to 600°C for 1 hour prior to growth to remove impurities from the surface. High purity 6N Antimony (Sb), Tellurium (Te) and 5N8 Manganese (Mn) fluxes were achieved by a Riber double zone cell for Sb and single zone Knudsen cells for Mn and Te. The fluxes were measured by the beam equivalent pressure (BEP) read by a Riber ion gauge placed in the position of the substrate prior to growth. The Mn BEP fraction, BEP(Mn)/[BEP(Mn) + BEP(Sb)], proportional to the Mn flux fraction, was used to control the Mn content during growth, and was varied between 0.00 and 0.11. Growth was performed under excess Te flux. The samples were grown by MBE via a 2-step growth method consisting of an initial deposition of a thin (usually 3-5 nm) low temperature buffer (LTB) Sb2Te3 layer, at 200°C for ~3-5 minutes. The LTB layer growth was then stopped, while keeping constant Te flux at the sample surface, and the substrate temperature was raised to the desired growth temperature of 245 -255°C. A LTB growth has been shown to improve nucleation in MBE growth of TIs [START_REF] Bansal | Epitaxial Growth of Topological Insulator Bi2Se3 Film on Si(111) with Atomically Sharp Interface[END_REF][START_REF] Bansal | Thickness-Independent Transport Channels in Topological Insulator Bi2Se3 Thin Films[END_REF] . Next, the Mn containing layer was grown for 1-2 hours, during which the Mn, Sb and Te sources were simultaneously impingent on the sample surface. The Te:Sb BEP ratios were kept between 20-30, ensuring excess Te during growth. All samples used similar Sb flux, and the composition was varied by adjusting the Mn cell temperature to change the Mn flux. The growth rates for the samples varied between 0.4 and 1.2 nm/min depending on the amount of Mn incorporated.

All sample were characterized by a variety of post growth techniques. High resolution X-Ray diffraction (HR-XRD) measurements were performed using a Bruker D8 Discover diffractometer with a da-Vinci configuration and a Cu Kα1(1.5418 Å) source. Atomic force microscopy (AFM) images were obtained using Bruker Dimension FastScan AFM with a FastScan-A silicon probe. 2b-e) were performed (EAG Laboratories) using a Hitachi HD-2700 Spherical Aberration Corrected Scanning-TEM system coupled with energy dispersive X-ray spectroscopy (EDX). Further high-resolution STEM imaging (Figures 2a and6) was carried out using a probe-corrected JEOL ARM-200F operated at 200 keV. Hall transport measurements were performed in a 14 T Quantum Design physical property measurement system (PPMS) in 1 mTorr (at low temperature) of He gas.

Scanning transmission electron microscope (STEM) images (Figures

Electrical contacts in the vdP configuration were made with indium bonded on the edge of the thin film.

Results and Discussion: Growth was monitored in-situ by reflection high energy electron diffraction (RHEED). Typical patterns at the end of the growth, and atomic force microscopy (AFM) surface images are shown in Figure 1. The initial RHEED pattern after LTB and during initiation of growth often exhibited some polycrystalline rings or spotty features which typically disappeared or became less prominent as the growth proceeded. Figures 1a-c show the RHEED images for three samples grown with different Mn BEP ratios, taken at the end of the growth. Most samples, regardless of composition, showed a streaky pattern, such as the one shown in Figure 1a; however, a few samples exhibited additional features suggesting rough surfaces and some disorder, as illustrated in Figures 1b-c. A corresponding variation in surface roughness is also visible in the AFM images presented in Figures 1d-f. Scanning transmission electron microscopy (STEM) observations were performed on some of the samples, primarily using high-angle annular-dark-field (HAADF) imaging because of its high-Z sensitivity. Figure 2 compares HAADF STEM cross sections of three samples having different Mn levels. The low magnification STEM images presented in Figures 2a-c enabled the analysis of large regions of each sample and estimates for the relative amount of septuple layers (%SL). The results illustrate the formation of SLs and QLs in different proportions, ranging from 90% SLs in Figure 2a to 13% SLs in Figure 2c. All of the STEM images show well ordered, highly crystalline structures and no obvious defects. The sample of Figure 2a, which was grown using a Mn BEP fraction of 0.11, has a crystal structure that consists almost completely of SLs, 90% calculated from the STEM image, and has a 27.6% Mn atomic percent, as obtained from energy dispersive X-ray (EDX). This result implies that, in addition to the Mn incorporated to form the SL structure, a significant amount of excess Mn is also incorporated as an impurity in the SLs and the QLs in that sample. Figure 2b shows the STEM image for another sample that contains 71% SLs. Interestingly, that sample was grown with a similar Mn BEP ratio as the sample in Figure 2a (0.10) and has similar total Mn content of 27.4% as obtained from EDX. The lower number of SLs in this sample suggests that an even larger fraction of the Mn was incorporated as impurity atoms (Mn in Sb and Te sites), and specifically, that the QLs may also contain large quantities of substitutional Mn. Figure 2c shows the STEM image for a third sample, which contains only 13% SLs and has a Mn content of 8.1% (from EDX). As for the previous two samples, based on the Mn content from EDX, a larger amount of SLs (over 50%) would be expected if all the Mn atoms were incorporated to form stoichiometric SLs. The presence of excess Mn throughout some regions of the sample can be seen in Figures 2d and2e, which show the Mn EDX cross sectional scan overlaid on the STEM image in two different regions of the sample of Figure 2c. In the region depicted in Figure 2d the Mn is mostly confined to the SLs, whereas the Mn atoms in Figure 2e are more randomly distributed throughout the crystal lattice and occupy other sites, while also forming SLs. 2a and2c.

All of the samples were characterized by high resolution X-ray diffraction (HR-XRD). A selection of the HR-XRD scans for samples ranging in Mn BEP fraction from 0.02 to 0.11 are shown in Figure 3a, where the sample presented in Figure 2a is the one corresponding to the top XRD scan. The expected position of the peaks for pure Sb2Te3 (bottom) and pure MnSb2Te4 (top), are also indicated in the figure. The XRD scans are ordered by the Mn content of the samples, from the lowest Mn content at the bottom to the highest at the top. The Mn BEP ratio used during growth is given on the right of each scan. It is evident that as the content of the Mn increases, there is a gradual change in the HR-XRD from that of a QL-based Sb2Te3 crystal structure to a MnSb2Te4 (SL-based) structure. The Sb2Te3 peaks corresponding to the (003), (006) and (0015) planes show a clear shift as the Mn fraction increases, evolving into the (003), (009) and (0021) peaks of MnSb2Te4, respectively, as can be seen by the red arrows superimposed on the plots. At intermediate Mn flux fractions (0.03-0.06), two weak peaks appear at around 2θ-ω ~10 -15°, as previously reported by others [START_REF] Lapano | Adsorption-Controlled Growth of MnTe(Bi2Te3)n by Molecular Beam Epitaxy Exhibiting Stoichiometry-Controlled Magnetism[END_REF] . Those authors proposed that this feature indicates a change of structure from Sb2Te3 to close to MnBi4Te7. At the high Mn fraction levels (0.09-0.11), a single peak that corresponds closely to MnSb2Te4 (006) and (0018) planes appears. To further illustrate the gradual shift in the XRD peaks, Figure 3b shows an expanded view of the region around the (0015) plane of Sb2Te3 for a few selected samples. A clear trend is visible in the shift of the peak position as the Mn content increases, ranging from the Sb2Te3 (0015) peak position to the MnSb2Te4 (0021) peak position.

Figure 3c shows the relationship between the shift of the HR XRD (0015) peak and the %SL calculated from the STEM images. The calculated %SLs in the structure measured from the STEM images plotted as a function of the peak position shows a very good linear correlation.

Thus, the shift in the XRD peak position is a direct measure of the %SL in the structure, and the linear relationship obtained, as given by Equation ( 1), can be used to estimate the %SL in all of our samples based on the HR-XRD (0015) peak position. While this shift correlates with the increased presence of SLs in the samples, it is not an indication of the total Mn content. This is exemplified by the results for the sample shown in Figure 2a, which although having nearly the same peak position as that reported for pure (stoichiometric) MnSb2Te4, its Mn concentration according to EDX is almost twice as high as that of stoichiometric MnSb2Te4 (14.3 % Mn). Thus, we conclude that our samples, grown under the conditions described above, are Mn-rich and likely contain a fair amount of disorder in the form of Mn impurity concentration. In an attempt to reduce the amount of apparent disorder in the crystals, a modified growth procedure was developed. Two changes were implemented. First, an annealing step was added after an initial thin portion of the layer was grown (from here on referred to as pre-anneal step).

Annealing has been used often in MBE and other growth techniques to improve crystalline quality [START_REF] Lee | Defect Reduction by Thermal Annealing of GaAs Layers Grown by Molecular Beam Epitaxy on Si Substrates[END_REF][START_REF] Guo | Single Domain Bi2Se3 Films Grown on InP(111)A by Molecular Beam Epitaxy[END_REF] . Using an anneal step after a short growth to promote SLs formation has been recently reported by others [START_REF] Gong | Experimental Realization of an Intrinsic Magnetic Topological Insulator Chin[END_REF] . In our experiments, the sample growth began the same way as before, with the thin Sb2Te3 LTB followed by an interruption of the growth and the initiation of growth of the Mn-containing layer at the desired growth temperature (Tg). After 15 minutes, the growth is again stopped, and the sample temperature is raised to 300°C for 10 minutes under a Te flux (pre-anneal step). The temperature is then lowered back to Tg and growth of the Mn containing layer is resumed for the desired total growth time.

As a first assessment of the effect of the anneal step, we examined the evolution of the RHEED pattern. The RHEED image of the sample after 15 minutes of growth, but before the preannealing step (Figure 4a), often shows noticeable rings, suggesting some polycrystalline growth. After completing the 10 min. annealing step, the RHEED pattern typically improves (Figure 4b), showing less evidence of polycrystalline growth. Finally, at the end of the full layer growth, the RHEED improves significantly (Figure 4c). More importantly, implementation of the pre-annealing step had other effects on the sample that is grown over the pre-annealed thin layer.

Figure 4d shows a comparison of the (0015) peak of the XRD scans of four samples, two grown with the pre-anneal step and two grown without. All samples in Figure 4d were grown using a similar Mn flux fraction, Tg and Te overpressure. The figure shows a clear shift toward larger angle in the (0015) XRD peaks for the two samples grown with the pre-anneal step, suggesting that larger %SLs form in the samples grown over the pre-annealed thin layer. The calculated %SLs for each sample of 4d, summarized in Table 1 as samples 1 through 4, increase from 40-50 %SLs to 70-90 %SLs with the pre-annealed layer (samples 3 and 4), indicating that more of the Mn is incorporated in the form of SLs when the layer growth occurs over a thin pre-annealed layer, consequently reducing the incorporation of substitutional Mn into QL and SL. In some samples, in addition to the pre-anneal step, the Tg was also increased to 265-270°C.

The (0015) peak of the HR-XRD scans of two samples grown with a pre-anneal step and the same value of Mn BEP, but with two different Tg values of 250°C and 265°C, are presented in Figure 4e. A further increase in the %SLs was observed for the sample grown at the higher growth temperature of 265°C indicated by the shift of the (0015) peak. The %SLs for these two samples, calculated from the XRD peak position, are given Table 1, as samples 5 and 6. Sample 6, grown at 265°C, shows a further increase of the %SLs bringing its composition up to 100% SLs.

In order to understand the effect of the pre-anneal step, an experiment was performed in which a comparison was made between a thin layer grown for 15 minutes and annealed for 10 minutes at 300°C, and another thin layer that was also grown for 15 minutes but not annealed. The HR-XRD of the two samples are shown in Figure 4f. The stronger and sharper XRD peaks for the annealed sample suggest improved crystallinity. However, the positions of the (0015) peak for the two samples, identified by a dashed line in the figure, are the same, indicating the same %SL for both samples. The AFM images of the surfaces of the two samples (Figures 4g-h) show that the surface of the annealed sample (Figure 4h) is smoother (rms = 3.0nm) than that of the sample (Figure 4g) that was not annealed (rms = 4.2nm). The data of Figures 4f-h suggests that the preanneal step improves the crystalline quality of the annealed thin layer and makes its surface smoother, but it does not change the %SL of the thin annealed layer. Thus, we conclude that the increased %SLs obtained in the layer grown over the pre-annealed thin layer must be explained by a different mechanism not evident from our experiment. The relationship between the Mn BEP fraction and the %SL of all the samples grown is illustrated in Figure 5a, which shows the samples grown using the original growth conditions as black open circles, the samples grown using the pre-anneal step as red open circles and those grown using the pre-anneal step as well as a higher Tg as full red circles. Within the range of MBE growth temperatures used here (250-270°C), it can be assumed that all the Mn incorporates in the layer (i.e., Mn will not re-evaporate from the growth surface). For the samples grown under the original growth conditions (black circles), there is a reasonable correlation between the Mn BEP ratio and the XRD peak shift (or %SLs), but a large amount of scatter is present in the data. The scatter is consistent with our assessment that Mn is incorporated into the crystal in different ways, and that a significant fraction of Mn is being incorporated as substitutional dopants without modifying the crystal structure (i.e., without forming SLs).

We also observe that the samples grown with the modified growth conditions shown by the points in red, all have higher %SLs for the same Mn BEP ratio. This is consistent with our conclusion that the pre-anneal step promotes the formation of SLs in the subsequent layer, rather than its incorporation in the lattice as a dopant. Within the red data points, the samples grown at a higher Tg, illustrated by the filled red circles, typically result in even higher %SLs for the same Mn BEP during growth.

Our observations are summarized schematically in Figure 5b. The Mn BEP ratio used during growth is illustrated by the green cone on the right, increasing from values of 0.00 to 0.11 in our growths. At very low Mn BEP ratios, only Sb2Te3 doped with Mn (Sb2Te3:Mn) can form. As more Mn is used during growth, (as the Mn BEP ratio increases), SLs can be formed resulting in mixed QL and SL structures. The implementation of a pre-annealing step results in SLs formed at lower Mn BEP ratio. Based on these observations we propose the following tentative mechanism. We assume that a critical quantity of excess Mn atoms must be present at the surface of the sample in order to form a SL rather than incorporate as impurities in the QLs. During growth, the Mn incorporates into the QL as impurity atoms and some excess Mn accumulates on the growing surface until there is a sufficient quantity required to form a SL. We propose that the pre-annealing step facilitates the accumulation of Mn on the surface, thus favoring the formation of SLs at lower Mn BEP ratios as compared to when no pre-annealing step is performed. Once this excess Mn is produced, it is more likely for an excess of Mn to remain on the surface as the layers grow, promoting the formation of SLs throughout the growth.

The effect of increasing the growth temperature is also illustrated in Figure 5b. An increase in the %SL formation for the same Mn BEP fraction is also observed by the use of a higher Tg. The reason for this enhancement might also be related to the ease of Mn accumulation on the surface 20 at a higher temperature, since at a higher temperature Sb may desorb from the surface leading to higher Mn concentration. However, the growth temperature effect may also be due to the higher formation energy of MnSb2Te4 21 , which may be favored at higher Tg. Magnetotransport measurements on our samples revealed that all samples having at least a few %SLs exhibited ferromagnetic behavior, as demonstrated by the presence of a hysteresis loop near ~0T magnetic field (B). This result is counter to the expected anti-ferromagnetic behavior predicted by DFT calculations for pure MnSb2Te4 [START_REF] Yan | Evolution of structural, magnetic, and transport properties in MnBi2-xSbxTe4[END_REF] . A typical Hall resistance as a function of magnetic field for a sample that has 72% SLs is shown in Figure 6a. The STEM image, in Figure 6b, confirms the %SLs, and shows no significant defects. The Mn content in this sample is ~21% (as measured by EDX), which is about 50% higher than the expected value for ideal (stoichiometric) MnSb2Te4 crystal with 100% SLs. This indicates that, as found previously in the samples of Figure 2, there is significant excess Mn in the crystal lattice. Ferromagnetism in Mnrich MnSb2Te4 has been recently reported and is attributed to Mn-Sb site exchange (illustrated in Figure 6c) [START_REF] Riberolles | Evolution of Magnetic Interactions in Sb-Substituted MnBi2Te4[END_REF] .

In order to determine the location of the excess Mn in the crystal in our samples, intensity profile measurements of the HAADF images, which are sensitive to atomic number, were performed. The result, depicted in the inset of Figure 6b, shows a combination of SLs evident by a sequence of seven intensity peaks, one for each atomic plane, and QLs, evident by the presence of five peaks corresponding to the five atomic planes. The profile is aligned with the section of the STEM image where the profile data was collected for ease of interpretation. The center peak of each SL or QL unit is marked by orange or black arrows, respectively. The SLs have a drop in intensity in the peak corresponding to the middle layer due to the small mass of the Mn atoms, suggesting predominantly Mn in that layer, as expected. However, the peaks corresponding to Sb and Te in the SLs show some variations in intensity, which possibly can be attributed to Mn site exchange. By contrast, the QLs show uniformly high intensity for all the peaks, as expected from the closely similar atomic mass of Te and Sb. Further observation of intensity variations in the STEM image also suggests Mn substitution in some Sb sites and Sb substitution in some Mn sites of the SL, examples of which are illustrated by the large orange and green circles in the expanded region of the STEM image shown in Figure 6d. Thus, we conclude that the ferromagnetic behavior of our MnSb2Te4 could also be due to Mn-Sb site exchange. A detailed study of the magnetic properties of our samples is currently underway and is beyond the scope of this paper.

Conclusions:

In this work, we have investigated the growth by MBE of (Sb2Te3)x(MnSb2Te4)y structures on sapphire substrates. Increasing the flux fraction of Mn (or Mn BEP fraction) during growth results in a gradual change of the composition of the structure from all Sb2Te3 QLs (y = 0) to all SLs (x = 0). The presence of significantly more Mn in the samples than what is accounted for by the percent (%) of SLs, suggests that Mn also incorporates as a dopant forming Sb2Te3:Mn and MnSb2Te4:Mn. We conclude that under our initial growth conditions, controlling the Mn flux ratio is not sufficient to accurately determine the resulting (Sb2Te3)x(MnSb2Te4)y structure, and that a large degree of excess Mn is incorporated in the layers. We also analyzed the evolution of the XRD scans of the layers as the Mn content is increased. The results show a gradual shift of several peaks from the Sb2Te3 position to the MnSb2Te4 position. We show that this shift can be used to calculate the %SLs in the resulting crystal.

By modifying the growth conditions with the incorporation of a pre-annealed thin layer at the beginning of growth, a significant increase of %SL formation for the same Mn flux fraction used during growth was observed in the overall layer. Thus, with the pre-annealed layer, the incorporation of Mn as a structural element to form SLs, is favored. We propose a mechanism consistent with our observations based on the assumption that a local excess of Mn is needed at the growth front for the formation of SLs during growth. We suggest that the formation of a Mn-

Figure 1 :

 1 Figure 1: RHEED images (a-c) of three samples and their corresponding 2μm x 2μm AFM (d-f)

Figure 2 :

 2 Figure 2: HAADF STEM images of three samples grown with different Mn BEP ratios, 0.11,

Figure 3 :

 3 Figure 3: (a) HR-XRD measurements of samples grown with increased Mn BEP Flux ratios. Peak positions for Sb2Te3 and MnSb2Te4 13 are noted as black lines below and above the

Figure 4 :

 4 Figure 4: (a-c) RHEED images of a sample surface before (a), and after (b), pre-annealing and at

Figure 5 :

 5 Figure 5: (a) Peak position of all samples (and %SLs), as a function of the Mn BEP ratio.

Figure 6 :

 6 Figure 6: Analysis of a sample with 72% SLs; (a) Hysteresis loop of the Hall resistance Rxy as a

Table 1 :

 1 Calculated %SLs for samples grown under different growth conditions.

	Sample Pre-annealed	Mn	BEP	Tg	% SLs
		Layer	ratio		(C)
	Effect of pre-annealed layer (Figure 4a)
	1	No	0.060		255	51
	2	No	0.060		260	41
	3	Yes	0.060		255	72
	4	Yes	0.065		255	92
	Effect of increased growth temperature (Figure 4b)
	5	Yes	0.090		250	79
	6	Yes	0.090		265	100

Acknowledgements

This work was supported by NSF Grant Nos. DMR-1420634 (MRSEC PAS 3 ), DMR-2011738 (PAQM), HRD-1547830 (CREST IDEALS) and HRD-2112550 (Phase II CREST IDEALS). The authors would like to acknowledge the Nanofabrication Facility of the CUNY Advanced Science Research Center (ASRC) for instrument use, scientific and technical assistance.

*mtamargo@ccny.cuny.edu

rich growth surface is facilitated by the pre-annealing step. A further increase in the %SLs for the same Mn flux fraction was observed by also increasing the growth temperature.

Hall resistance as a function of magnetic field measurements of our samples show that the samples are ferromagnetic, which may be in part related to the excess Mn incorporation in our layers. We present an HAADF intensity profile measurement to show evidence of substitution of Mn in Sb sites, as well as Sb in Mn sites. Mn/Sb anti-sites have been reported to lead to ferromagnetic behavior. The demonstrated control of the composition of Mn containing Sb2Te3 structures and the evidence for reduced excess Mn in samples grown under the modified growth conditions explored provide essential insights for the controlled growth by MBE of these highly promising and technologically relevant materials.


Dieter Fiems
email: dieter.fiems@ugent.be

Eitan Altman

Markov-Modulated Stochastic Recursive Equations with Applications to Delay-Tolerant Networks

Keywords: Stochastic recursive equations, Multitype branching processes, Random environment, Immigration, Delay tolerrant Ad-hoc Networks

We investigate Markov-modulated stochastic recursive equations. This class includes multi-type branching processes with immigration as well as linear stochastic equations. Conditions are established for the existence of a stationary solution and expressions for the first two moments of this solution are found. Furthermore, the transient characteristics of the stochastic recursion are investigated: we obtain the first two moments of the transient solution as well. Finally, to illustrate our approach, the results are applied to the performance evaluation of packet forwarding in delay-tolerant mobile ad-hoc networks.

Equations Stochastiques Recursives à Environement Aléatoire Markovien avec Applications aux Résaux Ad Hoc Tolérant les Délais

Résumé : Cet article est consacré à l'étude d'équations récursives stochastiques modulées par une chaîne de Markov. Cette classe d'équations contient entre autre le processus de branchement multitype avec immigration et les équations stochastiques linéaires. Des conditions sont établies pour l'existence d'une solution stationnaire et pour cette solution, on déduit des expressions pour les deux premiers moments. Ensuite, les caractéristiques transitoires de la récursion stochastique sont étudiées et les deux premiers moments de la solution transitoire sont établis. En conclusion, pour illustrer notre méthode, les résultats sont appliqués à l'évaluation de la performance de l'expédition de paquets dans les réseaux ad hoc mobiles tolérants aux retardements.

Mots-clés :

Equations stochastiques recursives, Processus de branchement Multitypes, Environement aléatoire, Immigration, Réseau ad-hoc tolérant les délais 1 Introduction In this contribution, we investigate a stochastic recursive equation of a form related to both branching processes and linear stochastic equations which operates in a random environment. Both steady-state and transient behaviour are investigated and, in either case, expressions for the first two moments are established. Our results are then applied in the context of packet forwarding in delay-tolerant mobile ad-hoc networks.

The object of this paper is equations of the type

X n+1 = A n (X n ) + B n ,
where X n ∈ R M are non-negative column vectors. B n ∈ R M , which we call the "immigration" term, is assumed to be stationary ergodic and A n are a sequence of independent and identically distributed (i.i.d.) stochastic processes. Moreover, the processes A n are assumed to have a divisibility property: if for some k, x = x 0 + x 1 + ... + x k then A n (x) can be represented as

A n (x) = k i=0 A (i) n (x i ) where { A (i)
n } i=0,1,2,...,k are identically distributed with the same distribution as A n (•), but they need not be independent. An example of an integer-valued process A(x) that is divisible is the number of arrivals at an M/G/1 queue during the service time of x customers. A non-descrete example is the time to serve all the customers that arrive at an M/G/1 queue during some time x.

The above stochastic recursive equations include various types of stochastic processes: stochastic difference equations (also called stochastic auto-regressive processes), branching processes in a discrete state space and branching processes in a continuous state space. These processes have in common that the expectation of the state, E[X n ], has the same linear dynamics: there is some matrix A

such that E[X n+1 ] = A E[X n ] + E[B n ].
The main objective of this work is the introduction of a random environment Y n in the above class of recursive equations (which we assume for simplicity to be a function of a state of some finite state Markov chain). More precisely, we introduce equations of the type

X n+1 = A n (X n , Y n ) + B n (Y n ).
where B n is the immigration process, and Y n is a Markov chain that can be used to produce a correlation between B n and A n as well as to remove the independence assumption on A n .

The theoretical foundations are then applied to various models of delaytolerant ad-hoc networks whose dynamics can be described by the stochastic recursive equations under investigation. In particular, the impact of the random environment on the performance of these networks is investigated.

Before proceeding to our main results, we survey some related literature. A special case of our framework are branching processes with a random environment. These have been well studied, both with and without immigration; see the survey [START_REF] Athreya | Handbook of Statistics[END_REF] and the references therein. For example, conditions are presented for the extinction when the random environment is stationary ergodic. Further, the stability, strong law of large numbers and central limit theorems for multitype branching processes with immigration in a random environment have been studied in [START_REF] Key | Limiting distributions and regeneration times for multitype branching processes with immigration in a random environment[END_REF][START_REF] Roitershtein | A note on multitype branching processes with immigration in a random environment[END_REF].These processes find applications in very diverse fields, including biological systems and queueing theory. For example, McNamara et. Al [START_REF] Mcnamara | Optimality models in behavioral biology[END_REF] consider an asexual species with non-overlapping generations. Individuals born some a year, reach maturity and reproduce one year later and then die. The number of individuals of the different genotypes in the consecutive years constitute a multi-type branching process. A numerical example is presented for which the expected population size has the same growth rate with and without a random environment. However, the probability that the population gets extinct is different in the two cases, being equal to one in the case of the random environment.

Prime examples in queueing theory where branching processes with immigration play a major role, include infinite server queues [START_REF] Eliazar | On the discrete-time G/GI/∞ queue[END_REF], processor sharing queues [START_REF] Grishechkin | On a relation between processor sharing queues and Crump-Mode-Jagers branching processes[END_REF][START_REF] Queija | Processor-Sharing Models for Integrated-Services Networks[END_REF], as well as various polling systems. Resing [START_REF] Resing | Polling systems and multi-type branching processes[END_REF] already demonstrated that the numbers of customers in the different queues of a polling system at polling instants are described by a multi-type branching process with immigration. Similarly, station times -the time the server remains with a particular queue -at polling instants are described by such a branching process as well, albeit with a continuous state space [START_REF] Altman | Expected waiting time in polling systems with correlated vacations[END_REF].

Many of the above queueing models have natural extensions that involve random environments. A polling example that can be modelled in our framework is studied in [START_REF] Macphee | A Markov chain model of a polling system with parameter regeneration[END_REF]. Every time one of the queues empties, some parameters can change at random. Another application is a polling system in which the polling order of the queues is determined according to a Markov chain: if the nth station being polled is station i, then the probability that the next one to be polled is j is given by the transition probability p ij . The infinite server queue with random environment has been studied recently in [START_REF] Auria | M/M/∞ queues in semi-markovian random environment[END_REF][START_REF] Falin | The M/M/∞ queue in a random environment[END_REF]. These authors assume a framework of i.i.d. exponentially distributed interarrivals and i.i.d. exponentially distributed service times. Our approach allows, in contrast, to obtain explicit expressions for the first and second moments in the more general setting of general stationary ergodic arrivals and general independent bounded service time, with Markov random environment.

The theory that we develop allows us to obtain the two first moments of the stochastic recursive equations under very general assumptions on the immigration process: it is only assumed to be stationary ergodic.

The remainder of this paper is organised as follows. In the next section, the stochastic recursive equation is introduced and some notation is established. Section 3 is concerned with conditions for the existence of a stationary solution and expressions for the first two moments of this stationary solution. Transient characteristics are found in Section 4. With the results established, we investigate various applications of our framework in the context of delay-tolerant mobile ad-hoc networks in Section 5. Finally conclusions are drawn in Section 6.

Stochastic model

We consider the sequence of random column vectors X n ∈ R M , adhering to the stochastic recursive equation,

X n+1 = A n (X n , Y n ) + B n (Y n ) , n ∈ Z . (1)

INRIA

Here Y = {Y n } denotes a Markov chain, taking values on a finite state-space Θ = {1, 2, . . . , N } whereas A n and B n denote random vector-valued functions with domain R M ×Θ and Θ, respectively. The functions A n (•, i) are independent random variables for all i ∈ Θ, n ∈ Z and further adhere to the following assumptions.

For each y ∈ Θ, A n (•, y) has a divisibility property. Let x = x 1 + x 2 + . . . + x k ∈ R M , then A n (x, y) has the following representation,

A n (x, y) = k i=1 Â(i) n (x i , y) , (2)
whereby Â(i) n (•, y), i = 1, . . . , k, are identically distributed, but not necessarily independent, with the same distribution as A n (•, y).

For each y ∈ Θ, A n (•, y) is linear in the mean, E[A n (x, y)] = A (n) y x , x ∈ R M , y ∈ Θ , n ∈ Z . (3)
Here {A

y , y ∈ Θ, n ∈ Z} is a set of fixed M × M matrices. Further, for each y ∈ Θ, the correlation matrix of A n (x, y) is linear in xx ′ and in x. For all x = [x 1 , . . . , x M ] ∈ R M , we have the following representation,

E[A n (x, y)A ′ n (x, y)] = F (n) y (xx ′ ) + M j=1 x j Γ (n) y,j , (4)
y ∈ Θ, n ∈ Z . For each y ∈ Θ and n ∈ Z, F

(n) y is a linear operator that maps M × M non-negative definite matrices on M × M non-negative definite matrices and satisfies

F (n) y (0) = 0. Further, {Γ (n) y,j , y, j ∈ Θ, n ∈ Z} is a set of fixed M × M matrices.
We mention some processes A n (•, y) for which the assumptions (2) to (4) hold. A first example constitutes the class of linear stochastic recursive equations,

X n+1 = A n+1 (Y n )X n + B n (Y n ).
Here A n+1 (y) constitutes a sequence of independent and identically distributed random matrices for all y ∈ Θ. Multi-type branching processes are a second example, either with a continuous or a discrete state-space. Examples of a branching processes with discrete state spaces are developed in Section 5. For an example of branching process with a continuous state space, the reader is referred to [START_REF] Altman | Expected waiting time in polling systems with correlated vacations[END_REF]. Finally, we note that any combination of the previous processes also complies with assumptions (2) to [START_REF] Eliazar | On the discrete-time G/GI/∞ queue[END_REF].

Before proceeding to our main results, we introduce some additional notation. Let p

(n) ij = Pr[Y n+1 = j|Y n = i] denote the transition probability of the Markov chain Y k at time n (i, j ∈ Θ) and let P (n) = [p (n)
ij ] denote the corresponding transition matrix. The probability that the Markov chain is in state k at slot n is denoted by π

(n) k = Pr[Y n = k].
For the immigration process B n , the following notation is introduced for the first and second order moments, b

(n) i = E[B n (i)] , B (m,n) ij = E[B m (i)B n (j)]
Finally, the following block matrices and block vector are defined to simplify further notation,

Â(n) =       A (n) 1 p (n) 11 A (n) 2 p (n) 21 . . . A (n) N p (n) N 1 A (n) 1 p (n) 12 A (n) 2 p (n) 22 . . . A (n) N p (n) N 2 . . . . . . . . . A (n) 1 p (n) 1N A (n) 2 p (n) 2N . . . A (n) N p (n) N N       , (5) b
(n) = k∈Θ π (n) k       p (n) k1 b (n) k p (n) k2 b (n) k . . . p (n) kN b (n) k       , (6)
and,

B(m,n) = j∈Θ π (m) j       B (m,n) j1 p (m) j1 B (m,n) j2 p (m) j1 . . . B (m,n) jN p (m) j1 B (m,n) j1 p (m) j2 B (m,n) j2 p (m) j2 . . . B (m,n) jN p (m) j2 . . . . . . B (m,n) j1 p (m) jN B (m,n) j2 p (m) jN . . . B (m,n) jN p (m) jN       . (7)

Stationary analysis

In this section, we investigate the existence of a stationary solution of the recursion (1) in a stationary ergodic framework, and obtain expressions for the first two moments of this solution. For this, we make the following additional assumptions.

The process {B n , n ∈ Z} is stationary ergodic.

The Markov chain Y n is ergodic.

The processes A n are independent and identically distributed.

In view of these assumptions, we may simplify notation as follows:

A (n) i = A i , b (n) i = b i , b(n) = b, B (m,n) ij = B (n-m) ij , p (n) ij = p ij , P (n) = P , π (n) k = π k , Â(n) = Â and B(m,n) = B(n-m) .
For ease of notation, the following operator will prove useful: for any x ∈ R M , let k i=n A i (x, Y i ) = x for k < n whereas, for k ≥ n, this operator is defined by the following recursion,

k i=n A i (x, Y i ) = A k k-1 i=n A i (x, Y i ), Y k .
The operator above can be applied analogously on Â(l) i for each l. We now state the stability theorem.

INRIA Theorem 1 Assume that (i) b i < ∞ component-wise for all i ∈ Θ; and (ii) that all the eigenvalues of the matrix  are within the open unit disk. Then, there exist a unique stationary solution X * n , distributed like,

X * n = d ∞ j=0 n-1 i=n-j Â(n-j) i (B n-j-1 (Y n-j-1 ), Y i ) , (8)
for n ∈ Z. The sum on the right side of the former expression converges absolutely almost surely. Furthermore, one can construct a probability space such that lim n→∞ X n -X * n = 0, almost surely, for any initial value X 0 . Proof. The proof follows a standard Loynes scheme. We define on the same probability space, the sequence of processes

X [l] = {X [l] n ; n ∈ Z, n ≥ -l}, with initial state zero X [l]
-l = 0 and governed by the recursive equation ( 1). One easily verifies that for fixed n, the sequence

X [l]
n is monotone increasing (componentwise) in l. Hence, the limit lim l→∞ X

[l] n X *
n is well defined; the right-hand side of (8) follows by consecutively applying the recursion (1), taking into account the divisibility (2).

Since, X * n is measurable on the tail σ-algebra of the stationary ergodic sequence

{(A n (•, Y n ), B n (Y n )); n ∈ Z}, it
is either finite, almost surely, or infinite, almost surely. In view of the recursion (1), we find,

E[X [l] n+1 1{Y n+1 = j}] = i∈Θ A i E[X [l] n 1{Y n = i}]p ij + b i π i p ij ,
for all j ∈ Θ. Here 1{} is the indicator function which equals 1 if its argument is true and 0 if this is not the case. The former system of equations can be expressed in matrix notation as,

µ [l] n+1 = µ [l] n + b . µ [l]
n denotes a (block) column vector with elements µ

[l] n (j) E[X [l] n 1{Y n = j}] (j ∈ Θ) .

In view of the former expression, one easily verifies by induction that µ

[l] n ≤ (I -Â) -1b , component-wise. This then implies that E[X [l]
n ] is bounded above by a finite constant for all l ≥ 0, n ≥ -l which in turn implies that E[X * n ] is bounded above. Since X * n is either finite, almost surely or infinite almost surely, we conclude that X * n is finite, almost surely. Consider now, the process {X n ; n ∈ N}, governed by the recursion (1), with initial state X 0 . Consecutively applying the recursion yields,

X n = n-1 j=0 n-1 i=n-j Â(n-j) i (B n-j-1 (Y n-j-1 ), Y i ) + n-1 i=0 Â(0) i (X 0 , Y i ) ,
for n ∈ N + . Hence, we have,

X n -X * n = n-1 i=0 Â(0) i (X 0 , Y i ) - n-1 i=0 Â(0) i (X * 0 , Y i ) , (9)
for n ∈ N + . We now show that both terms on the right-hand side of ( 9) converge to 0 almost surely. Let μn (x) denote the block column vector with elements,

E n-1 i=0 Â(0) i (x, Y i ) 1{Y n = l} , l ∈ Θ .
By conditioning on the state of the Markov chain Y and in view of the linearity in the mean (3), we find,

μn (x) = Â μn-1 (x) = Ân μ0 (x) , n ∈ N + ,
where μ0 (x) denotes a block column vector whose jth element equals π j x. Since  has no eigenvalues outside the open unit disk, we have lim n→∞ μn (x) = 0 for all x ∈ R M , which implies,

lim n→∞ E n-1 i=0 Â(0) i (x, Y i ) = 0 .
Since, the sequence of random variables

n-1 i=0 Â(0) i (x, Y i
) is non-negative, application of Fatou's Lemma yields that this sequence converges to 0 almost surely. In view of [START_REF] Key | Limiting distributions and regeneration times for multitype branching processes with immigration in a random environment[END_REF], this shows that X n -X * n converges to 0 component-wise, which concludes the proof. ⋄ With the stability conditions established, we now focus on expressions for the first and second moments of X * 0 , conditioned on the state of the Markov chain Y 0 . Let µ, the conditional first moment vector, be the block column vector with elements µ i E[X * 0 1{Y 0 = i}], i ∈ Θ. Analogously, let Ω, the conditional second moment matrix, be the block column vector with elements

Ω i E[X * 0 (X * 0 ) ′ 1{Y 0 = i}], i ∈ Θ.
The following theorem provides expressions for these vectors.

Theorem 2 Assume that the stability conditions of Theorem 1 are satisfied. The conditional first moment vector is then given by,

µ = (I -Â) -1 b . ( 10
)
Under the additional assumption that the second moments of B 0 (i) are finite, i ∈ Θ, the elements Ω i of the conditional second moment matrix of X * 0 are the unique solution of the system of equations,

Ω l = k∈Θ F k (Ω k ) + M j=1 µ (j) k Γ (j) k + B (0) kk π k + A k Λ k + Λ ′ k A ′ k p kl , l ∈ Θ, where Λ k denotes the kth diagonal (block) element of ∞ j=0
Âj B(j+1) and with µ

(j)
k the jth element of µ k .

INRIA

Proof. Taking the expectations in (1) immediately yields the following system of equations,

µ l = E[X * 1 1{Y 1 = l}] = k∈Θ E[(A 0 (X * 0 , k) + B 0 (k))1{Y 0 = k, Y 1 = l}] = k∈Θ A k µ k p kl + k∈Θ b k p kl π k , l ∈ Θ .
In block matrix notation, this system of equations is equivalent to, µ = µ + b , which implies [START_REF] Macphee | A Markov chain model of a polling system with parameter regeneration[END_REF]. Notice that I- is non-singular since  has only eigenvalues within the unit disk.

For the second moment, we first focus on the set of matrices Φ kl (n),

Φ kl (n) = E[X * 0 B n (l) ′ 1{Y 0 = k}] , k, l ∈ Θ, n ∈ N . ( 11
)
In view of the recursion (1) this matrix satisfies,

Φ kl (n) = E[X * 1 B n+1 (l) ′ 1{Y 1 = k}] = j∈Θ E[(A 0 (X * 0 , j) + B 0 (j))B n+1 (l) ′ 1{Y 0 = j}]p jk = j∈Θ A j Φ jl (n + 1)p jk + B (n+1) jl π j p jk ,
for k, l ∈ Θ, i ∈ N. In matrix notation, this system of equations reads,

Φ(n) = ÂΦ(n + 1) + B(n+1) = ∞ j=0 Âj B(n+j+1) , (12)
where Φ(i) denotes the block matrix with elements Φ kl (i), k, l ∈ Θ. Notice that the sum in [START_REF] Queija | Processor-Sharing Models for Integrated-Services Networks[END_REF] converges to a finite-valued matrix since (i) the finiteness of the second moments of B 0 (i) (i ∈ Θ) implies that the elements of B(n) are uniformly bounded and since (ii) the eigenvalues of  are within the unit disk.

With the expression of Φ(i) at hand, and in view of the recursion (1), we immediately find,

Ω l = E[X * 1 (X * 1 ) ′ 1{Y 1 = l}] = k∈Θ E[(A 0 (X * 0 , k) + B 0 (k))(A 0 (X * 0 , k) + B 0 (k)) ′ 1{Y 0 = k, Y 1 = l}] = k∈Θ   F k (Ω k ) + M j=1 µ (j) k Γ (j) k   p kl + k∈Θ B (0) kk π k p kl + k∈Θ A k Φ kk (0)p kl + k∈Θ Φ kk (0) ′ A ′ k p kl .
Here µ

(j) i
denotes the jth element of the vector µ i . This completes the proof. ⋄

RR n 6872

Transient analysis

The previous section focused on the existence and the moments of the steady state process. We now consider its transient analysis. Here, no additional assumptions are required but the assumption that the initial vector X 0 and the initial Markov state Y 0 are independent from the processes A n and B n . Let µ(n) denote the conditional first moment vector at slot n. That is, µ(n) is the column vector whose kth (block) element equals

µ k (n) = E[X n 1{Y n = l}].
Analogously, let Ω(n) denote the conditional second order moment (block) column vector whose kth element equals

Ω k (n) = E[X n X ′ n 1{Y n = l}].
Expressions for these block vectors are given in the Theorem below.

Theorem 3 The first moment vector at slot n + 1 can be obtained recursively by,

µ(n + 1) = Â(n) µ(n) + b(n) ( 13
)
The second moment vector at slot n+1 can be obtained by the following recursion

Ω l (n + 1) = k∈Θ   F (n) k (Ω k (n)) + M j=1 µ (j) k (n)Γ (n) k,j   p (n) kl + k∈Θ B (n,n) kk π (n) k p (n) kl + k∈Θ A (n) k Φ (n,n) k,k p (n) kl + k∈Θ (Φ (n,n) k,k ) ′ (A (n) k ) ′ p (n) kl , (14)
with

Φ (n,n) k,k
the kth diagonal (block) element of

Φ (n,n) = Â(n-1) Â(n-2) • • • Â(0) Φ (0,n) + n-1 k=0 Â(n-1) Â(n-2) • • • Â(k+1) B(k) , (15)
and where Φ (0,n) is a block matrix with elements

Φ (0,n) k,l = µ k (0)(b (n) l ) ′ , k, l ∈ Θ. Proof.
Taking the expectations in (1) immediately yields the following system of equations,

µ l (n + 1) = k∈Θ A (n) k µ k (n)p (n) kl + k∈Θ b (n) k p (n) kl π (n) k , l ∈ Θ .
which yields [START_REF] Resing | Polling systems and multi-type branching processes[END_REF]. Proceeding analogously, we find [START_REF] Roitershtein | A note on multitype branching processes with immigration in a random environment[END_REF] for Ω i (n),

Ω l (n + 1) = E[X n+1 (X n+1 ) ′ 1{Y n+1 = l}] = k∈Θ   F (n) k (Ω k (n)) + M j=1 µ (j) k (n)Γ (n) k,j   p (n) kl + k∈Θ B (n,n) kk π (n) k p (n) kl + k∈Θ A (n) k Φ (n,n) k,k p (n) kl + k∈Θ (Φ (n,n) k,k ) ′ (A (n) k ) ′ p (n) kl , INRIA with, Φ (m,n) k,l = E[X m B n (l) ′ 1{Y m = k}] .
For the latter expectation, we further have the following recursion,

Φ (m,n) k,l = j∈Θ E[(A m-1 (X m-1 , j) + B m-1 (j)) B n (l) ′ 1{Y m-1 = j, Y m = k}] = j∈Θ A (m-1) j Φ (m-1,n) j,l p (m-1) jk + B (m-1,n) j,l π (m-1) j p (m-1) jk , (16)
whereas the independence of the initial state yields,

Φ (0,n) k,l = E[X 0 B n (l) ′ 1{Y 0 = k}] = µ k (0)(b (n) l ) ′ .
The system of equations ( 16), can be rewritten in block matrix notation,

Φ (m,n) = Â(m-1) Φ (m-1,n) + B(m-1)
Solving the recursion then yields [START_REF] Small | Resource and performance tradeoffs in delay-tolerant wireless networks[END_REF], which completes the proof. ⋄

Applications in Delay Tolerant Mobile Ad-hoc Networks

With the theory established, we now focus on some applications in the context of packet-forwarding in delay tolerant mobile ad-hoc networks. The first two examples concern models with a fixed number of nodes. The last two examples are concerned with models where the number of nodes vary during time.

Fixed number of mobile nodes

As a first example, we consider a network that consists of N mobile nodes. Some fixed node wishes to send a packet to a destination node. As connectivity is assumed to be low, the source makes use of the mobility of other mobiles that serve as relays. Whenever the source is within the transmission range of another node, it transmits a packet to that node. Whenever a node with a copy of a packet is within the transmission range of the destination then it transmits the packet to it.

Remark 1

In this example as well as in all the following ones, we restrict to the two-hop routing scheme [START_REF] Garetto | On the effectiveness of the 2-hop routing strategy in mobile ad hoc networks[END_REF]; we do not consider epidemic routing, i.e. a relay node that receives a packet from the source does not relay it further to other intermediate nodes.

Time is discrete and it is assumed that at each time n, each node has a probability p (n) i ≥ p > 0 to meet each other node; this probability also depends on the state i ∈ Θ of a modulating Markov chain. The transition matrix of this chain at time n is denoted by

P (n) = [p (n) ij ] i,j∈Θ .

RR n 6872

The validity of (a continuous time version of) this model without the random environment has been discussed in [START_REF] Groenevelt | Message delay in MANETs[END_REF], and its accuracy has been shown for a number of mobility models (Random Walker, Random Direction, Random Waypoint). The random environment enables us to further capture the fluctuations in time of the connectivity probability. For example, reduced connectivity probabilities may occur during rain storms which result in worse radio channel conditions.

Let ξ n (ℓ, j) denote the indicator that equals 1 if nodes ℓ and j meet at time n and is otherwise 0. Let ξ (i) n n,i be i.i.d. copies of ξ n Let X n be the number of nodes with a copy of the packet at time n, then the following recursion can be established,

X n+1 = X n + N -Xn i=1 ξ (i) n .
The framework established in the preceding sections cannot be used directly for this recursion. A simple change of variables, Xn = N -X N , however yields a recursion which is captured by our framework,

Xn+1 = Xn - Xn i=1 ξ (i) n = A n ( Xn ) + B n .
where

A n ( Xn ) = Xn - Xn i=1 ξ (i) n , B n = 0 .
Notice that Xn denotes the number of nodes that do not have a copy of the packet at time n. In view of the recursion for Xn , we obviously find B(n) = 0 and b(n) = 0 for all n. Further, we have,

A (n) i = 1 -p (n) i , F (n) i (x 2 ) = (1 -p (n) i ) 2 x 2 , Γ (n) i,1 = (1 -p (n) i )p (n) i , for i ∈ Θ.
The steady state behaviour is trivial, Xn converges to 0 and X n converges to N . Hence, we here focus on the transient behaviour. For this, we assume that the Markov chain is a time-homogeneous 2-state chain; state 1 (state 2) corresponding to low (high) interference levels. Let α and β denote the transition probabilities from state 1 to state 2 and from state 2 to state 1, respectively and let p i denote the probability that a node receives the packet in state i (i = 1, 2). For ease of notation, let σ = (1 -β)/(2 -α -β) denote the (long-term) fraction of slots that the interference level is low (the chain is in state 1) and let the interference time τ = 1/(1 -α) + 1/(1 -β) denote the mean length of a high and a low interference period.

In Figure 1, the time-evolution of the mean and the variance of the number of nodes that have the packet is depicted for different parameter sets; for all sets, we have N = 200 nodes, p 1 = 0.5%, p 2 = 10% and at time 0 the interference level is low. Further, following transition probabilities are chosen: α = 17/18 and β = 1/2 for set 1, α = 89/90 and β = 9/10 for set 2 and α = 449/450 and β = 49/50 for set 3. The fraction of time that the interference level is INRIA high is equal for all three sets. However, in comparison with set 1, interference periods are on average 5 times (25 times) longer for set 2 (set 3). It is readily observed that the lengths of the interference periods have a huge impact on the performance. Longer periods yield a slower, more variable spreading of the packet among the nodes.

Packet discarding

We now move to some models where the steady state behaviour is non-trivial.

In order to avoid packets to remain forever at nodes, it has been suggested to use expiration timers [START_REF] Small | Resource and performance tradeoffs in delay-tolerant wireless networks[END_REF] for packets. Assume that each node uses a geometrically distributed initial value for the timer, with parameter q. The timer is initiated when a packet is received for relaying and the packet is discarded when the timer expires.

Retaining the assumptions and notation of the previous example, the number of nodes with the packet (excluding the source node) at consecutive slots are related as follows,

X n+1 = Xn i=1 ζ (i) n + N -Xn i=1 ξ (i) n .
Here, ζ (i) n are i.i.d. indicators that equal 0 if the ith copy of a packet is discarded at time n. Note that for this model, the source node is not included in X to ensure that the source does not discard the packet. As in the previous example, this model is not directly covered by our framework. However, we may rewrite the equation as follows,

X n+1 = Xn i=1 ζ (i) n + Yn i=1 ξ (i) n , Y n+1 = N - Xn i=1 ζ (i) n - Yn i=1 ξ (i) n .
In vector notation this set of equations can be written as follows,

X n+1 Y n+1 = A n X n Y n + B n , with, A n x y = x i=1 ζ (i) n + y i=1 ξ (i) n 1 -1 , B n = 0 N .
This shows that our framework is applicable. We get the following matrices,

A (n) i = q p (n) i -q -p (n) i , Γ (i) 1,n = q(1 -q) 1 -1 -1 1 , Γ (i) 2,n = p (n) i (1 -p (n) i ) 1 -1 -1 1 ,
and,

F (i) n x 2 xy xy y 2 = x 2 q 2 + 2xyqp (n) i + y 2 (p (n) i ) 2 1 -1 -1 1 .
Figure 2 depicts the mean and variance of the number of nodes that have the packet (in steady state) vs. the mean packet discarding time T = 1 -1/q for the same parameter sets as in Figure 1. Obviously, as nodes keep the packets longer, the mean number of nodes that have the packet increases. In contrast, the variance first increases with the discarding time and then decreases again. As in the previous example, correlation in the random environment negatively affects the mean number of nodes with the packet while the variance increases.

Variations in the total number of mobiles

As before, some fixed source node wishes to send a packet to a destination node, thereby making use of the mobility of other mobiles that serve as relays. Whenever the source is within the transmission range of another node, it transmits the packet to that node. Whenever a node with a copy of the packet is within the transmission range of the destination, it transmits the packet to it. Further, we here adopt the two-state Markov environment of the previous examples, modelling the variations of the interference level.

In contrast to the previous examples, the number of nodes does not remain constant. New nodes arrive and nodes depart from the system. Let W n denote the number of nodes that have the packet at time slot n and let Z n denote the Fiems & Altman number of nodes that do not have the file. We have the following recursion,

W n+1 = Wn j=1 ζ (j) n,1 + Zn j=1 ζ (j) n,2 ν (j) n Z n+1 = Zn j=1 ζ (j) n,2 (1 -ν (j) n ) + B n (17)
Here

ζ (j)
n,1 is the indicator that the jth node that has the packet leaves the system at slot n, ζ

(j)
n,2 is the indicator that the jth node that does not have the packet leaves the system at slot n and ν (j) n is the indicator that the jth node that does not have the packet, receives the packet at slot n. Finally, B n denotes the number of new nodes that arrive during slot n.

We make the following assumptions. The indicators ζ

(j)
n,1 and ζ

(j)
n,1 constitute doubly indexed sequences of independent Bernoulli distributed random variables; let q denote the probability that a node leaves the system. A node then remains in the system for T = 1/(1 -q) slots on average. Further, the indicators ν (i) n are Bernoulli distributed random variables whose distributions depend on the interference level (the state of the Markov chain) during slot n. Let p 1 (p 2 ) denote the probability that a node receives the packet if the Markov chain is in state 1 (state 2). Finally, the sequence of new nodes B n is stationary ergodic. Under these assumptions, the recursion (17) clearly fits the framework; let X n [W k , Z k ] ′ . Adhering to the notation of the framework of Sections 2 to 4, the following matrices and vectors characterise the recursion,

P = α 1 -α 1 -β β , A i = q qp i 0 q(1 -p i ) , Γ (i) 1 =
q(1 -q) 0 0 0 ,

Γ (i) 2 = qp i (1 -qp i ) -q 2 (1 -p i )p i -q 2 (1 -p i )p i q(1 -p i )(1 -q(1 -p i )) ,
and,

F (i) (xx ′ ) = x 2 1 q 2 1 0 0 0 + x 1 x 2 q 2 2p i 1 -p i 1 -p i 0 + x 2 2 q 2 p 2 i (1 -p i )p i (1 -p i )p i (1 -p i ) 2 , with x = [x 1 , x 2 ] ′ and i = 1, 2.
We focus on the steady-state behaviour. Figure 3 depicts the mean number of nodes E[W ] that have the packet and the mean number of nodes E[Z] that do not have the packet in steady-state. In Figure 3(a), these means are plotted vs. the mean residence time T of the nodes for different values of the interference time τ . The mean number of nodes in the system is equal to 50; the mean number of arrivals E[B] in a slot is reduced for increasing mean residence times. If there is no interference (state 1), a node receives the packet with probability INRIA p 1 = 0.1 whereas no transmission is possible during periods of high interference (p 2 = 0). For all curves, the interference level is high during σ = 90% of the slots. It is readily observed that the mean residence time of a node has a considerable impact on E[W ]. Obviously, if nodes remain longer, they carry the packet for a longer time which explains the increase in the mean number of nodes that carry the packet. Further, increasing interference times yield decreasing E[W ]. This is confirmed by Figure 3 In accordance with equation ( 10), the mean value does not depend on the second order moments of the number of arrivals during the consecutive slots. The variance however does depend on these moments.

E[W], E[Z] τ E[W] E[Z] E[B]=1 E[B]=2 E[B]=3 (
In Figures 4 and5, the variance of the number of nodes with and without the packet are depicted versus the mean residence time (a) and versus the mean interference time (b) for the same parameters as in Figure 3. Additionally, the number of new arrivals are assumed to constitute a series of Poisson distributed random variables; this means that the variance of the number of new arrivals equals the mean number of new arrivals. These variables constitute an independent sequence in Figure 4 whereas their autocorrelation function ρ(n) has a geometric decay -ρ(n) = 1/2 n -in Figure 5. Amongst others, discrete autoregressive processes of order 1 have such an autocorrelation function. From these figures, it is readily observed that the variance of the number of nodes with (without) the packet is heavily affected by the mean residence time, the mean interference time and by the correlation of the number of new arrivals. For increasing values of the residence time, the variance of the number of nodes with (without) the packet first increases and then again decreases. A number of different effects are at play here. (i) Longer residence times imply that there are fewer arrivals which remain longer (recall that the mean number of nodes in the system is constant). As such, longer residence times may or may not increase the variance of the number of nodes in the system: on the one hand, the variance of the residence times increases, on the other hand the variance of the (iii) Increasing the interference time means that there are longer periods with packet transmissions as well as longer periods without transmissions. Longer interference times thus imply higher variances.

Mobility of the source and the nodes

As a more advanced application of our framework, we now consider the following scenario. Again, some fixed source node wishes to send a packet to a destination node, thereby making use of the mobility of other mobiles that serve as relays. Whenever the source is within the transmission range of another node, it transmits the packet to that node. Whenever a node with a copy of a packet is within the transmission range of the destination, it transmits the packet to it. The source node moves according to a random walk through a spatial grid, for example the spatial grid depicted in Figure 6. In each of the regions of the grid, new nodes arrive according to a stationary ergodic process which then travel through the grid until they leave. If a node is in the same region as the source, the node receives the packet with a fixed (possibly region-dependent) probability. In the remainder, let N denote the number of regions.

Let X n (k) denote the number of nodes in region k at time n that have the packet and let Z n (k) denote the number of nodes that do not have the packet. Further let X n and Z n denote the column vectors with elements X n (k) and Z n (k), respectively. Let Y n denote the region where the source node resides at time n and let B n (k) denote the number of new nodes that arrive in region k at time n; B n is the column vector with elements B n (k). We then have the following recursion,

X n+1 = N i=1 Xn(i) j=1 ζ (i,j) n,1 + N i=1 Zn(i) j=1 ζ (i,j) n,2 ν (i,j) n Z n+1 = N i=1 Zn(i) j=1 ζ (i,j) n,2 (1 -ν (i,j) n ) + B n (18)
Here

ζ (i,j)
n,1 is a column vector of indicators; its kth element is the indicator that the jth node in region i that has the packet at time n moves to region k. The indicator vector ζ (i,j) n,2 is defined likewise. Its kth element is the indicator that the jth node in region i that does not have the packet at time n moves to region k. Further, ν (i,j) n denotes the indicator that the jth node in region i that does not have the packet at time n, receives the packet. Notice that some of the packets may leave the grid as not all packets necessarily move to any of the regions.

We make the following assumptions on these indicators. The indicators ν

X n+1 Z n+1 = A n X n Z n , Y n + 0 B n , (19)
where the sequence A n is independent and identically distributed, adhering assumptions (2) to (4) such that our framework is applicable.

To limit the number of parameters involved, we make the following assumptions. All nodes move in the grid depicted in Figure 6. Each node remains in the same region during T = 1/(1 -α) time intervals on average and then moves to any of the neighbouring regions with probability 1/6. T is referred to as the mean residence time of the nodes. Analogously, the source node remains in the same region for T slots on average and then moves to any of the neighbouring regions with equal probability. The source node never leaves the grid. A node that does not have the packet which is in the same region as the source node, receives the packet with probability p. Under these assumptions, the system parameters are determined in the appendix.

In Figure 7 the mean number of nodes with (a) and without (b) the packet in the different regions is depicted vs. the mean residence time T in a region for different values of the transmission probability p. The number of new arrivals in the different regions scales with the residence times of the nodes: E B (i) = 50/T for i = 1, 2, . . . , 7 such that the total number of nodes in the different regions remains constant. First, notice that by symmetry, the characteristics of regions 2 to 7 are the same. Further, it is clear that longer residence times imply that more nodes receive the packet. Clearly, nodes do not only remain longer in a region but also longer in the grid. Hence, the probability that they receive the packet increases. Since the number of new arrivals is scaled with the residence times, the mean number of nodes without the packet decreases for increasing mean residence times; see Figure 7(b).

Conclusions

This paper provides expressions for the first two moments -both transient moments and steady state moments -for stochastic recursive equations which encompasses both linear stochastic recursive equations and multi-type branching processes with immigration in a random environment. The immigration term in these recursions is taken to be stationary ergodic, whereas the random environment is taken to be Markovian. With the theory established, various examples in the context of delay-tolerant ad-hoc networks are developed: we find transient moments for a fixed number of mobile nodes in a Markovian environment and steady state moments for a fixed number of nodes with packet discarding after the exploration of an expiration timer, for a variable number of mobile nodes in a Markovian environment and for a variable number of mobile nodes moving in a grid.

A Moments for the mobility application

In this appendix, we calculate the moments of the branching process of the mobility application of Section 5.4. Under the assumptions of Section 5.4, the transition matrix of the Markov chain Y is given by the matrix P below.

P =                         
α 3 0 0 0 α 3 α                         
Further, we introduce the auxiliary sub-stochastic matrix Q, its entries being the probabilities that a node moves from one region to the other,

Q =                         
                        

RR n 6872

Given the matrix Q above, it is easily shown that the 14 × 14 matrices Âk which correspond with the process A n () of equation ( 19), are given by,

A k = Q QS k 0 Q(I -S k ) ( 20
)
where S k denotes a 7 × 7 matrix of zeroes except for its k diagonal element which equals p.

For the second moment, it is necessary to determine the linear operators F (i) and the matrices Γ i,j , i = 1, . . . , 7, j = 1, . . . , 14. These immediately follow from the following expressions for the second order (cross) moments,

Figure 1 :

 1 Figure 1: Mean and variance of the number of nodes that have the packet as a function of time for various parameter settings.

Figure 2 :

 2 Figure 2: Mean and variance of the number of nodes that have the packet as a function of the mean discarding time for various parameter settings.

  (b) where E[Z] and E[W ] are depicted vs. the interference time τ for various values of the mean number of arrivals E[B]. The same parameter setting is assumed as in Figure 3(a).

Figure 3 :

 3 Figure 3: Mean number of nodes that have the packet and of the number of nodes that do not have the packet vs. the mean residence time for various values of the mean interference time (a) and vs. the mean interference time for various values of the mean number of arrivals in a slot (b).

Figure 4 :

 4 Figure 4: Variance of the number of nodes that have the packet and of the number of nodes that do not have the packet vs. the mean residence time for various values of the mean interference time (a) and vs. the mean interference time for various values of the mean number of arrivals in a slot (b). There is no arrival correlation.

Figure 5 :

 5 Figure 5: Variance of the number of nodes that have the packet and of the number of nodes that do not have the packet vs. the mean residence time for various values of the mean interference time (a) and vs. the mean interference time for various values of the mean number of arrivals in a slot (b). The autocorrelation function of the new arrivals decays geometrically.

Figure 6 :

 6 Figure 6: Spatial grid of the nodes

are

  Bernoulli distributed random variables whose distributions only depends on the position Y n of the source node. Further, a node moves from one region to another with a fixed probability. Hence, the vectors ζ (i,j) n,1 and ζ (i,j) n,1 constitute triply indexed sequences of independent and identically distributed random vectors. Finally, new nodes arrive according to a stationary ergodic process. We

Figure 7 :

 7 Figure 7: Mean number of nodes with the packet (a) and without the packet (b) vs. the mean residence time in a region

RR n 6872

Fiems & Altman