R Reche

H Beust

J.-C Augereau

O Absil

On the observability of resonant structures in planetesimal disks due to planetary migration

Keywords: N-body simulations -Celestial mechanics -(Stars, ) planetary systems -(Stars, individual, ) Vega

Context. The observed clumpy structures in debris disks are commonly interpreted as particles trapped in mean-motion resonances with an unseen exo-planet. Populating the resonances requires a migrating process of either the particles (spiraling inward due to drag forces) or the planet (moving outward). Because the drag time-scale in resolved debris disks is generally long compared to the collisional time-scale, the planet migration scenario might be more likely, but this model has so far only been investigated for planets on circular orbits. Aims. We present a thorough study of the impact of a migrating planet on a planetesimal disk, by exploring a broad range of masses and eccentricities for the planet. We discuss the sensitivity of the structures generated in debris disks to the basic planet parameters. Methods. We perform many N-body numerical simulations, using the symplectic integrator SWIFT, taking into account the gravitational influence of the star and the planet on massless test particles. A constant migration rate is assumed for the planet.

Results. The effect of planetary migration on the trapping of particles in mean motion resonances is found to be very sensitive to the initial eccentricity of the planet and of the planetesimals. A planetary eccentricity as low as 0.05 is enough to smear out all the resonant structures, except for the most massive planets. The planetesimals also initially have to be on orbits with a mean eccentricity of less than than 0.1 in order to keep the resonant clumps visible. Conclusions. This numerical work extends previous analytical studies and provides a collection of disk images that may help in interpreting the observations of structures in debris disks. Overall, it shows that stringent conditions must be fulfilled to obtain observable resonant structures in debris disks. Theoretical models of the origin of planetary migration will therefore have to explain how planetary systems remain in a suitable configuration to reproduce the observed structures.

Introduction

Since the first direct imaging of a debris disk around β Pictoris by Smith & Terrile (1984), a dozen other optically thin dust disks have been spatially resolved around nearby main-sequence stars showing an infrared excess (Kalas et al. 2007;Schneider et al. 2006, and references there in). The images often reveal asymmetric structures and clumps, interpreted as the signature of gravitational perturbations. A planet immersed in a debris disk usually produces structures such as a gap along its orbit, by ejecting particles during close encounters, or density waves (e.g. a one-arm spiral), by modifying the precession rate of the dust particles (Wyatt 2005a). However, such structures cannot explain the observations of clumpy, non-axisymmetric disks [START_REF] Augereau | Extrasolar Planets: Today and Tomorrow[END_REF][START_REF] Meyer | Protostars and Planets V[END_REF], and resonant mechanisms with unseen planets have been proposed to account for the observed asymmetries (Ozernoy et al. 2000;Quillen & Thorndike 2002;Kuchner & Holman 2003;Wyatt 2003). A particle belongs to a mean motion resonance (MMR) when the particle to planet period ratio is a rational number, m : n with m and n integers. An MMR is located at a semi-major axis a given by a/a p = (m/n) 2/3 , where a p is the planet semi-major axis. In the Solar System, for example, about 15% of the known Kuiper Belt objects, including Pluto, are trapped in the 3:2 resonance with Neptune [START_REF] Chiang | Protostars and Planets V[END_REF]). The interesting property of MMRs for modeling Send offprint requests to: R. Reche e-mail: Remy.Reche@obs.ujf-grenoble.fr asymmetric disks is that, as explained for example in [START_REF] Murray | Solar System Dynamics (Solar System Dynamics[END_REF], resonant objects are not uniformly distributed in azimuth around a star: rather they gather at specific longitudes relative to the perturbing planet and subsequently form clumps. This arises from properties specific to MMRs as a given particle trapped in a MMR with a planet undergoes conjunctions with the planet at specific locations along its orbit. The particles tend to gather around the most stable orbital configurations that ensure that the conjunctions occur at the maximum relative distance. The clumps, which are the result of the collective effects of resonant particles, generally corotate with the planet (Kuchner & Holman 2003), while each of these resonant bodies has a different period from that of the planet (except for 1:1 resonant planetesimals): hence the motion of these density waves differs from the orbital motion of the resonant particles. However, MMRs are very thin radial structures that usually trap a small number of particles in a given disk. Therefore, any structure due to MMRs has a high chance of being totally hidden by the emission of the non-resonant particles, as illustrated in Fig. 1.

For clumps due to MMRs to be observed, the population of resonant particles must be significantly enhanced by an additional physical process. Two mechanisms can account for this: Poynting-Robertson (P-R) drag and planet migration. Dust particles that are too large to be ejected from the system by radiation pressure can spiral inward into the star due to P-R drag and to some other minor forces like stellar wind drag (e.g. Augereau & Beust 2006). In the course of its inward migration, a dust particle can be trapped into exterior MMRs with a planet, hence increasing the contrast of their asymmetric patterns (Kuchner & Holman 2003;Deller & Maddison 2005). Planet migration, on the other hand, involves particles of all sizes, except those ejected by radiation pressure. Many particles can be trapped in MMRs by a planet migrating outward in the disk. Each nonresonant particle crossing an MMR has a chance trapped and subsequently migrating, following the resonance (Wyatt 2003).

Several authors have studied either the effect of P-R drag, or of planet migration, on disk structures, using different methods (analytic, semi-analytic or numerical) and various planet parameters (mass and orbital eccentricity). A summary of the main previous studies is provided in Table 1. The P-R drag scenario has been extensively studied for a wide range of parameters, while the migrating planet scenario has been investigated only for a planet on a circular orbit by Wyatt (2003). It is thus important to better characterize the latter scenario in order to distinguish which of the two dominates the morphology of debris disks. Moreover, a number of studies (Lecavelier Des Etangs et al. 1996;Wyatt 2005b;Krivov et al. 2007) have shown that collisions may prevent MMRs from being populated by P-R drag since the collision timescale in massive debris disks might be much shorter than the P-R drag migration timescale.

Therefore, we propose in this paper to extensively study the planet migration scenario, using numerical modeling, by generating a synthetic catalog similar to what has been done for the P-R drag scenario using analytical (Kuchner & Holman 2003) or numerical studies (Deller & Maddison 2005). We extend the pioneering work done by Wyatt (2003) in studying the influence of the planet eccentricity on the visibility of the resonant patterns. In Section 3, we discuss the case of a planet migrating on a circular or low-eccentricity orbit. In Section 4, we extend this study to planets on orbits with eccentricities up to 0.7. This study is generalized to various migration rates and disk initial states in Section 5, and compared to previous studies in Section 6. The limitations of our approach are discussed in Section 7.

Numerical model

We consider a planetary system consisting of a star surrounded by a planet and a debris disk. We address the case of large particles, which are insensitive to pressure forces (radiation, stellar wind or gas pressure). The simulated disk is thus rather a planetesimal disk than a debris disk and we only consider gravitational forces. Importantly, we also do not take into account the gravitational interactions between planetesimals as they are negligible, nor mutual collisions. Dynamically speaking, the planetesimals are thus considered as test particles. A typical configuration for the simulations is a Vega-like central star (2.5 M ⊙ ) and a planet orbit with a 40 AU pericenter at the starting time. The initial planetesimal disk consists of 50 000 planetesimals distributed between 40 and 75 AU on circular orbits, with a surface density distribution proportional to r -1 . The disk midplane coincides with the orbital plane of the planet, and the inclinations of the planetesimals are randomly distributed within ±3 • .

In this model, the planet keeps a Keplerian orbit around the star, or migrates at a constant rate without modification of its eccentricity. This basic model is easy to implement and to analyze, and corresponds to the case described by Wyatt (2003). The goal of this paper is to extend this initial work to a wider range of planet eccentricities by a numerical study. We start by studying planets on low-eccentricity orbits (e < 0.1), and then extend our work to larger eccentricities.

To perform our simulations we have used the symplectic package SWIFT (Wisdom & Holman 1991;Levison & Duncan 1994), to which we have added planetary migration. To do this, we plugged in the Wyatt (2003) prescription. This method consists of adding an acceleration in the direction of the orbital motion of the planet, with an intensity equal to: vp = 0.5 ȧp GM * /a 3 p , where G is the gravitational constant, M * the stellar mass and ȧp the variation rate of the planet semi-major axis a p . This causes a change in the planet semi-major axis without modifying its eccentricity (for a planet on a low-eccentricity orbit; for a planet on a higher eccentricity orbit the change is not significant) or its inclination. We do not discuss here the origin of the migration. It can be due either to the migration of a large internal planet, or to the gravitational influence of the planetesimals themselves. The most important here is that we keep ȧp constant during each simulation, generally at 0.5 AU Myr -1 , to match the Wyatt (2003) model for the Vega disk. We have used the RMVS3 version of the SWIFT integrator, in order to have a better modeling of the close encounters between the planet and the planetesimals. In all the simulations, the system evolution is followed for 40 Myr.

Planets on low-eccentricity orbits

The scenario of a planet on a low-eccentricity orbit is the most studied case (see Table 1), for several reasons. First, planets were originally expected to be on almost circular orbits because, during the protoplanetary phase, circumstellar gas forces the planets to remain on very low eccentricity orbits. A planet on such an orbit therefore corresponds to the "standard scenario". Also, a low or zero planetary eccentricity simplifies an analytical analysis (Kuchner & Holman 2003). Nevertheless, it must be noted that many of the extrasolar planets detected so far have high eccen- tricities 1 , and we will therefore extend our study to high eccentricities in Section 4. Although this standard scenario has already been well studied, all its aspects have not yet been investigated. Wyatt (2003) studied the case of an outward migrating planet on a strictly circular orbit, while Kuchner & Holman (2003) and Deller & Maddison (2005) studied the case of planets on fixed loweccentricity orbits, considering only inward dust migration due to P-R drag. We propose in this section to numerically study a migrating planet on a circular or low-eccentricity orbit to search for possible differences with respect to previous studies. -External MMRs at large a/a p are weaker and less able to trap numerous planetesimals than closer MMRs (Kuchner & Holman 2003). First-order resonances (i.e. m = n + 1) with large m values are the closest and hence the strongest MMRs. However, MMRs near the planet are located close to each other and compete. A resonance overlap criterion (Wisdom 1980;Duncan et al. 1989) predicts that first-order resonances become completely chaotic when m > 0.45µ -2/7 + 1, where µ = m p /M * . This places a limit on the nearest resonance that can be populated. In the Solar System, the first completely overlapped MMR is the 17:18 for the Earth and the 4:5 in the case of Jupiter. - Wyatt (2003) showed that the planetary mass acts as a threshold for the probability of capture in first-order external MMRs: for a fixed migration rate, the probability drops quickly to 0 below a certain planetary mass, while above this mass the probability grows quickly to 1. For MMRs of higher order, the transition is less sharp. -An MMR that traps all the planetesimals at a given semimajor axis a r stops the growth of any other resonance that reaches this a r afterward. Resonances like 3:2 or 5:3 can thus trap a large number of planetesimals thanks to a large enough spatial separation, while resonances near the planet are too close to each other to trap large populations. This shows that 1 www.extrasolar.eu any modeling of structure generation with this process must be done globally, as the various MMRs compete with each other to be filled. -The long term evolution of planetesimals trapped in MMRs globally results in an increase of their eccentricity. If this eccentricity becomes high enough, planetesimals can become planet-crossing. Thanks to the resonance, they are nevertheless phase-protected against close encounters with the planet. But if the eccentricity grows too high, this phase-protection does not hold any longer, due to strong modulations of the angular velocity. Hence MMRs are limited in eccentricity (Sicardy et al. 1993), with e max = (2/5m p ) 1/2 . This concerns more specifically the MMRs that are close to the planet, as the planetesimals may easily become planet crossing. We thus expect the closest MMRs to lose planetesimals when they reach a given eccentricity limit. This is illustrated in Fig. 3 where the limits in eccentricity of the resonant populations are clearly visible.

General trends

In conclusion, only a few resonances, namely the most external of the first or second order resonances (4:3, 3:2, 5:3 and 2:1 resonance) capture most of the planetesimals. However, depending on the planet mass, the planetesimals that are not trapped in resonances can, or not, change the shape of the disk. In our simulations, planets with a mass above 1 Jupiter mass eject almost all the non-resonant planetesimals: the MMR structures therefore appear clearly. Below 1 Jupiter mass, planets cannot eject all the non-resonant planetesimals, which can then partly hide the MMR structures.

Circular orbits

In the circular orbit case (Fig. 2 and3), differences appear between the simulations, depending on the planet mass:

-For an Earth mass planet, non resonant planetesimals are still bound to the system and hide all the resonant structures. However, a small hole at the planet location can be observed. -For a Neptune mass planet, two clumps of equal density are generated by the 3:2 resonance and they are located in opposition with respect to the star. It also appears that the planet is not massive enough to capture many planetesimals in the 2:1 resonance and thus does not generate thin rings at large distances like more massive planets. Non-resonant planetesimals create a ring inside the planet orbit. -For a Saturn mass planet, the 2:1 resonant pattern appears in addition to the 3:2 one. The 2:1 MMR produces only one

Low-eccentricity orbits

Most of the structures discussed in Section 3.2 disappear as the planet eccentricity increases (Fig. 4). The libration amplitude of resonant planetesimals indeed increases, smoothing the density waves along the orbit. For eccentricities between 0.05 and 0.1, the disk looks like a ring with a hole at the location of the planet, or no longer shows structures for the lower mass planets. The rings are not only due to the 1:1 resonant planetesimals corotating with the planet, but are also populated by other major resonances (e.g., 2:1,3:2). Massive planets (last row of Fig. 4) are less sensitive to this effect because they can more efficiently eject planetesimals during close encounters, even if they are in the MMRs. Only resonant planetesimals with low libration amplitudes can survive and the disk remains structured as in the strictly circular case.

Planets on eccentric orbits

At eccentricities larger than 0.1, planets can significantly modify the disk geometry and produce a dynamically warm disk where all planetesimals may reach high eccentricities. Analytical de-velopments are more complex, but are still feasible as shown by Kuchner & Holman (2003) for resonant planetesimals and by Wyatt et al. (1999) for non-resonant planetesimals. So far, only the P-R drag migration scenario has been tested for such planetary orbits. We explore here the planet migration scenario in the case of eccentric orbits (Fig. 6 and7).

According to previous works, three phenomena should change in our simulations with respect to the low-eccentricity orbit case. First, close encounters between the planet and the planetesimals are more frequent, increasing the depletion rate of the disk. As the probability of ejection increases with the planetary mass, the more massive the planet, the more depleted the disk. Second, even the non-resonant planetesimals see their eccentricities rise significantly due to the gravitational perturbation of the planet. Finally, the trapping probability is also modified (Quillen 2006): increasing the planetary eccentricity decreases the trapping probability of the first order resonances but can increase it for higher order resonances. With a planet on an eccentric orbit, the resonant planetesimals are also not well protected against close encounters, which limits even more the number of planetesimals in MMRs. The migration thus cannot populate the MMRs enough to generate detectable patterns: for structures generated by a planet on a moderate or high eccentricity orbit, the MMRs, and thus the migration of the planet, is no longer an important factor whatever the planet mass as the non-resonant dynamics dominates the shape of the disk.

The depletion rate is so efficient that the disk is almost entirely depleted during our 40 Myr simulations, except in the case of an Earth mass or Neptune mass planet on moderate eccentric orbits (below 0.2 or 0.3). An Earth mass planet on a more eccentric orbit (above 0.5) can also produce transient collective non-resonant effects, spatially fixed, in the azimuthal distribution of the planetesimals (Fig. 9) before the disk is depleted, as explained below.

Non-resonant structures

The theoretical background for the dynamics of the planetesimals in our simulations is the restricted three-body problem, i.e., a problem where a mass-less test particle orbiting a star is perturbed by a planet orbiting the star on an unperturbed Keplerian orbit. We restrict ourselves to the planar case for simplicity. In

H 0 = - GM * 2a -Gm p        1 r -r p - r • r p r 3 p        , ( 1
)
where a is the osculating semi-major axis of the orbit of the planetesimal, M * is the mass of the star, m p that of the planet, G is the constant of gravitation, and r and r p are the heliocentric position vectors of the planetesimal and of the planet respectively. If the test particle is not locked into a mean-motion resonance with the planet, its secular motion is investigated by performing a double temporal averaging of H 0 over the orbital motions of the planet and of the planetesimal (Beust & Dutrey 2006). In this context, the semi-major axis a is a secular constant, as it is canonically coupled with the mean anomaly that has been removed from the Two examples for e p = 0.1 and e p = 0.5 are shown in Fig. 8. They both correspond to a = 1.3a p (a p is the semi-major axis of the planet). The whole dynamical problem obviously simply scales with a p , so that a only needs to be given in units of a p . The present case (a = 1.3a p ) corresponds to a planetesimal orbiting outside the planet's orbit, like those we are simulating. Note also that the topology of the Hamiltonian is independent of the mass m p of the planet, because the non-constant part of H 0 is The plots in Fig. 8 hold for a = 1.3a p , but for other values of a, we have similar plots. Conversely, the shape of the phase portrait depends critically on e p . We see that for a small e p , a planetesimal with a small initial e will keep e small for ever. For a large e p however, any planetesimal with a small initial e will be driven to high e values and ν ≃ 0. Starting from a population of planetesimals with negligible eccentricities, we end up after a certain delay with many highly eccentric planetesimals with their lines of apsides more or less aligned with that of the planet ν ≃ 0. This naturally generates a clump of planetesimals close to the apoastron of their orbits, as due to Kepler's second law, the planetesimals spend most of their time near apoastron. This is the origin of the clumps we obtain in our simulations with low-mass planets (Fig. 9).

So, why does this not hold for more massive planets? The secular dynamics described above is valid as long as the planetesimal does not undergo any close encounter with the planet. In the case of a close encounter, the orbit of the planetesimal is suddenly changed, and it is often ejected. Many regions in Fig. 8 correspond to a planet crossing orbit. The probability of having a close encounter with the planet within a given timespan is higher if the planet is more massive. It scales as m 2/3 p , because the Hill radius r H of the planet scales as m 1/3 p , and the encounter cross-section is expected to scale as r 2 H . The mass ratio between a 3 Jupiter mass planet and an Earth-sized planet is ∼ 1 000. We thus expect a planetesimal to undergo 100 times more en-Fig. 7. Same as Fig. 3, for similar planets, but on a moderate eccentricity orbit (e p = 0.3).See the electronic edition of the Journal for a color version of this figure .   counters with the first planet than with the second. Finally, with massive planets, most of the planetesimals are subject to a close encounter with the planet within the timespan of the simulation described in Fig. 6. This is why the disk appears so depleted at the end. Conversely, for low-mass planets, the close encounter probability is so low that many planetesimals keep following the secular dynamics until the end of the run. Therefore, they have enough time to generate a strong asymmetric clump.

We stress here that this clump is not due to any mean-motion resonance. There is thus no need for planet migration in this case and this may appear as an alternative scenario to mean-motion resonance for generating transient clumps. Nevertheless, even with an Earth mass planet, these clumps do not last as long as the resonant clumps. Planetesimals are not protected against close encounters with the planet, which finally deplete most of the disk after 35 Myr, in our simulations.

Generalization

As explained in the previous sections, we have run many simulations in order to sample correctly the parameter space of planetary eccentricity versus planetary mass. We can thus address the question of the visibility of asymmetric structures in a disk. As already mentioned, the shape of a disk is dominated either by resonant or by non-resonant planetesimals. The region where MMRs dominate the disk shape correspond to planets on low-eccentricity orbits. In this region, two situations can occur: planets can generate clear resonant patterns with several visible clumps in the disk (generally while on circular or very low ec- centric orbits) or produce smooth patterns, with only a hole at the planet location as the visible structure (generally while on an orbit with an eccentricity between 0.05 and 0.1). Concerning the non-resonant planetesimals, Earth mass planets on eccentric orbits can generate observable structures by secular perturbations. Outside these regions, the disks do not show any observable structures, when they are not totally depleted. The results for all these simulations are summarized in Table 2 and Fig. 10. Three main regions can be identified. In zones I and II, observable structures in the disk are generated by MMRs while in zone III, transient structures are generated by non-resonant mechanisms. In the remaining region, the disk does not show any structure. In zone I, MMRs create clumpy disks while in zone II they generate a smooth disk with a hole at the planet location. This figure also shows the fraction of planetesimals still bound to the system after 40 Myr (background color). However, this quantity is sensitive to several parameters (stellar mass, duration of the simulation, initial distribution of the planetesimals ...) while the limits of the three zones are quite independent.

However, we have assumed for these simulations a constant migration rate of 0.5 AU Myr -1 and a disk of planetesimals initially on circular orbits. With different assumptions the outcomes of the simulations could be changed: we have thus investigated these two parameters in order to discuss the robustness of our conclusions. For a given planetary mass and eccentricity, we expect the structures to change, as the trapping probability depends on the migration rate and on the planetesimal eccentricity. But, again, our main focus is to determinate if the resonant structures are visible or not. For instance, in the low-eccentricity orbit case, Neptune mass and Jupiter mass planets do not produce the same structures but they have the same sensitivity to the planet eccentricity. Here, we investigate if the migration rate or the initial planetesimal eccentricity change these conclusions. The results of these additional simulations are summarized in Tables 3 and4, in the same manner as in Table 2 for the nominal case. Table 3 corresponds to simulations with different migration rates but unexcited initial planetesimal disks. Table 4 describes simulations with initially excited disks.

Migration rate

In our work, the migration rate parameter was chosen to match the best fit obtained by Wyatt (2003) for the Vega debris disk to allow direct comparison. In his paper, Wyatt discussed the impact of the migration rate in the restrictive case of a planet on a circular orbit, and has shown that the trapping probability in an MMR increases with decreasing migration rates (or when the star is less massive).

We have thus performed simulations with a lower (0.05 AU Myr -1 ) and a higher (5 AU Myr -1 ) migration rate than previously (Table 3). Overall, our results are in good agreement with those of Wyatt (2003), even for eccentric orbits: the trapping probability increases when the migration rate decreases. For example, with a low migration rate, a Neptune mass planet traps planetesimals in the 2:1 MMR, while a Saturn mass planet populates the second libration center of the 2:1 resonance.

Although the disk resonant shape is modified because the populated MMRs change with the migration rate, the depen- dence on the mass and eccentricity of the planet remains unchanged: non-resonant planetesimals always dominate at moderate eccentricity and clear resonant patterns are only visible with a planet on a quasi circular orbit. We still note some subtle changes between the different migration rates. With a higher migration rate, low mass planets have less time to increase the eccentricities of the non-resonant planetesimals which perturb less the resonant structures. With a very low migration rate, planets have more time, during close encounters, to eject non-resonant planetesimals, and resonant planetesimals in the eccentric planet orbit case. Nevertheless, as these changes are small, it is possible to summarize our simulations in the (planet eccentricity; planet mass) parameter space and to discuss the visibility of resonant clumps using only these two parameters, as in Fig. 10.

Warm disks

We have assumed in our simulations that the planetesimal disk was initially dynamically cold (e disk = 0.0), but Wyatt (2003), while using the same hypothesis, pointed out that dynamically warm disks are an interesting alternative to be studied. We have thus extended our work in this direction, assuming disks with planetesimal eccentricities uniformly distributed between 0 and e limit (Table 4).

It appears that an initially warm disk does not differ much from a cold disk in the case of a planet on a moderate or high eccentricity orbit. The secular perturbations, or the close encounters, due to such a planet raise the planetesimal eccentricities on a timescale shorter than the migration time. A cold disk therefore becomes warm in the course of planetary migration. However, with a planet migrating on a circular orbit, the disk structures are different depending on the initial eccentricities of the planetesimals. A warm disk (e limit = 0.1) with a planet on a circular orbit is actually roughly equivalent to a cold disk with a planet on a low-eccentricity orbit (0.05 or 0.1): the MMRs trap many planetesimals but the clumps are smoothed by the large libration amplitude (Fig. 11). From our simulations, it appears that, for a Saturn mass or Jupiter mass planet on a circular orbit, an e limit of 0.1 is already too high to keep the resonant clumps visible. On the other hand, planets above 2 Jupiter masses have a large enough depletion rate to keep the resonant clumps visible with an e limit up to 0.2. These results show that the resonant structures are as sensitive to the planet eccentricity as to the planetesimals eccentricity.

Comparison with previous works

P-R drag scenario

For both circular and low-eccentricity orbits, Kuchner & Holman (2003) obtain much larger differences in the structures between low and high mass planets that those observed in our simulations: the P-R drag scenario seems to have a stronger dependence on the planetary mass than the migration scenario. The particles migrating inward due to P-R drag encounter first the most external resonances, but the capture probability for these resonances is low and the dust is more likely trapped in resonances that are closer to the planet. With massive planets, however, the probability of capture in distant resonances is high enough to trap a large number of dust particles in these outer disk regions and stop them before they can populate closer resonances. The MMRs populated by P-R drag, and thus the structures generated, are therefore not the same, depending on the planetary mass. Conversely, in the case of planetary migration, all the resonances can be populated at the same time, as planetesimals are initially present in the whole disk: as explained in Section 3, the MMRs actually populated are thus roughly the same, whatever the planet mass.

In practice, the differences between our results and those of Kuchner & Holman (2003) are more due to the initial planetesimal distribution in the disk than to the physical process. In their P-R drag scenario, the authors assume that the dust starts migrating inwards far away from the planet, outside the 2:1 resonance, while our planetary migration scenario assumes the planetesimals to be initially closer to the planet. If we started with planetesimals at a larger distance, we would have expected to see, as in the P-R drag scenario, more differences in the resonant structures for planets of different masses.

Planetary migration

Our results show that, with the forced planetary migration scenario, it is easy to distinguish a planet on a circular orbit from another on a low-eccentricity orbit, except for very low mass planets or very massive planets, because the resonant structures are drastically different. Constraining the planetary mass is more difficult than in the P-R drag scenario and only an order of magnitude can be expected. Wyatt (2003) used this scenario to reproduce Vega disk observations at submillimetric wavelengths (Holland et al. 1998). We must take into account that, with the large SCUBA PSF, only the two major clumps can be observed. However, this is enough to distinguish between our three planet mass examples, at least for a migration rate of 0.5 AU Myr -1 :

-It is possible to distinguish between a Saturn mass planet and a Neptune mass or Jupiter mass planet because it is the only one out of the three that produces asymmetric clumps in density. -We can also distinguish between a Neptune mass planet and a Jupiter one because, in the first case, the two clumps are in opposition with respect to the star, while, in the second case, the two clumps are separated in longitude by less than 180 • .

It is thus possible to obtain an estimate of the planetary mass. The situation is well summarized by Figure 11 of Wyatt (2003). It defines several regions in the (planetary mass, migration rate) parameter space that can be observationally distinguished from each other, but inside each region a wide range of planetary masses is possible.

Contradictions however appear between our simulation results and this previous study. The asymmetry in the emission of the two observed clumps was interpreted as the migration of a Neptune mass planet by Wyatt (2003). In his model, a Neptune can trap planetesimals in the 3:2 and 2:1 resonances and generate two asymmetric clumps, like a Saturn mass planet in our simulations. With our numerical model, we have found that a Neptune mass planet cannot trap planetesimals in the 2:1 resonance, but only in the 3:2 resonance: the two clumps are thus symmetric and cannot reproduce the Vega disk. A Neptune mass planet at a migration rate of about 0.5 AU Myr -1 lies at the sharp transition between a 0 and a 100% trapping probability (Wyatt 2003, Fig. 4a). A small change in the planetary mass or the migration rate in this configuration produces a large modification in the population of this resonance. As Wyatt (2003) uses a scaling law to predict the trapping probability, differences between our results may be explained by the approximation of this scaling law.

Nevertheless, the 2:1 resonance has an interesting behavior in the Neptune mass planet case: it perturbs all the planetesimals that cross it, but as soon as they reach an eccentricity of about 0.02 (in 0.25 Myr), most of them escape. While a Saturn mass planet (or a more massive one) cleans up the space between the initial and final position of the 2:1 resonance during the migration by trapping all the planetesimals, a Neptune mass planet only slightly rises the eccentricity of planetesimals entering the 2:1 resonance during the migration process (Fig. 12).

This phenomenon is better seen when the planetesimal trajectories are drawn in a semi-major axis, eccentricity diagram where the semi-major axis of the planetesimals are in units of that of the planet, in order to hide the migration, as in Fig. 13. All the planetesimals initially have roughly the same trajectory (they move from right to left because they do not migrate), but the small variations have a strong impact when the planetesimals cross the resonance. Some planetesimals remain trapped in the MMR, while others escape after being temporarily perturbed. However, within each of these two subgroups, nearly all planetesimals have similar behavior: the permanently trapped planetesimals have the same libration amplitude and the temporarily perturbed planetesimals escape the resonance roughly at the same eccentricity.

The width (in semi-major axis) of an MMR is proportional to the square root of the planetary mass (Beust & Morbidelli 1996). A Saturn mass planet therefore has wider resonances than a Neptune mass planet: its trapping probability is thus larger because it is less sensitive to the orbital parameters of the planetesimals which cross the MMR.

Discussion and limitations

Dust content and collisional activity

In this paper, we have only discussed the spatial distribution of the planetesimals. However, it is the dust produced by collisions between these planetesimals, and not the planetesimals themselves, that contributes to the emitted flux and therefore to the observed disk shape. One has to investigate the influence of the collisions on the dust spatial distribution in order to reproduce the observations. Collisions can have several consequences for the dust particles:

-As collisional cascades produce dust particles of all sizes, radiation pressure cannot be neglected any longer. This has a strong incidence on the populations of the MMRs. As explained by Wyatt (2006); Krivov et al. (2007), depending on their sizes, the dust particles can either stay in the resonance, leave the resonance but still stay bound to the system, or be blown out of the disk. Even in the case where the MMRs dominate the dynamics of the planetesimals, the disk can look rather smooth if the observations are sensitive to the smallest of the produced dust particles. -The collision rate is not the same in the whole disk: it is enhanced in the MMRs, as shown by Queck et al. (2007). These authors have even calculated that the average collision rate in an MMR is the highest for first-order resonances (by a factor about 2 times the non-resonant collision rate). This implies that at the wavelengths where the observed dust particles can still appear trapped in the MMRs (typically at 850µm), the resonant clumps are more visible than in our simulations. This also means that the 3:5 resonance, often populated in the simulations, is actually not so prominent because its clumps have a lower average collision rate than the first order resonances. -The destructive collisions produce fragments with different velocities from the parents bodies. Since these velocities are usually small compared to the orbital velocities, most of the fragments cannot leave the resonances. However, analytic results from Krivov et al. (2007) show that small particles produced by the collisional cascade are fast enough to escape the MMRs. Dust particles are mainly lost by this "velocity effect" before they become small enough to escape the MMR by radiation pressure.

In conclusion, predictions on the distribution of the smallest dust particles, observable in the near or mid-infrared wavebands, is difficult since it requires a proper description of the collisional activity and of the dynamics of particles influenced by radiation/wind forces. This is not the case for larger dust grains which are observable at submillimeter wavelengths: according to the previously quoted papers, these particles, created by mutual collisions between the planetesimals, stay in the MMRs. However, the most visible resonant clumps are not necessarily the most populated, because the collision rate depends on the order of the resonance.

Origin of planetary migration

The described simulations obtained with our model raise numerous questions about the interaction between a planetesimal disk and a planet. The structures generated are very sensitive to the eccentricity of the planet, and to a lesser extent to the migration speed. Moreover, some of the generated structures (especially those obtained for high eccentricity, low-mass planets) appear only after a fairly long time. A key issue to address in this question is the temporal evolution of the eccentricity (and the migration rate) and the combined evolution of the disk. A low mass, high-eccentricity planet can only generate the features described above (Fig. 9) if its eccentricity remains high for a sufficiently long time. This may appear unrealistic, as the planet eccentricity should decrease due to the interactions with the planetesimals. Another important question is the transition between circular and low-eccentricity orbits: we have shown in the previous section that resonant patterns change significantly when the eccentricity deviates from 0. What will be the resulting disk structures if the planet eccentricity undergoes periodic secular modulation between these two regions (like the giant planets of the Solar System) ? And how do the particles stay at an eccentricity lower than 0.1 all the time ?

These questions can be resolved by better modeling the origin and the evolution of planetary migration. In any realistic simulation, the orbit of the planet will be subject to secular evolution. There is even no need for other planets for this. The interaction with the disk particles themselves can be sufficient to significantly affect the planetary orbit. Thus taking into account the influence of the disk on the planets is necessary to derive realistic simulations. Several studies have already discussed the origin of planetary migration either by "planet-planetesimals" (Ida et al. 2000;Gomes et al. 2004) interactions or by "planetplanet" interactions [START_REF] Morbidelli | AAS/Division for Planetary Sciences Meeting Abstracts[END_REF]. Depending on the scenario and on the initial conditions, one can observe migration on low-eccentricity orbit or more chaotic migration after a short time on an eccentric orbit.

Conclusion

We have studied the problem of the presence of observable structures in planetesimal disks due to mean motion resonance with an unseen planet migrating outward in the disk. Using numerical simulations, we have explored a large range of parameters for the planet (mass and orbital eccentricity) and the disk (initial distribution of planetesimal eccentricities). In the case of a planet on a circular orbit migrating inside a dynamically cold disk, our results are in agreement with previous analytical studies.

In the cases not already addressed, namely planets on eccentric orbits or dynamically warm disks, we have found that the observability of resonant structures demands very specific orbital configurations. The clumps produced by MMRs with a planet on a circular orbit are smoothed in the case of a planet on an even moderately eccentric orbit. An eccentricity as low as 0.05 is enough to smooth all the resonant structures, except for the most massive planets. These results indicate that although trapping planetesimals in MMRs is an efficient mechanism to generate clumpy disks, stringent conditions must be fulfilled for this scenario to occur. Theoretical modeling of the origin of the planetary migration therefore will have to explain how planetary systems can remain under these conditions. Moreover, we only consider a planet migrating at a constant rate. A more realistic model with a variable, stochastic migration rate can reduce the population of resonances and thus their observability. A better model of planet migration thus should be developed in future studies.

Table 2. Summary of results for all simulations done in the present study with an initially unexcited disk (initial planetesimal eccentricities are equal to zero) and a standard migration rate of 0.5 AU Myr -1 . For each simulation, we list the migration rate, the planet mass and eccentricity, the fraction of surviving planetesimals at the end of the simulation (40 Myr) and the resulting disk shape, following the convention of Fig. 10. disk. We also thank Alessandro Morbidelli, Jens Rodmann and the anonymous referee for helpful comments on this paper. Most of the computations presented in this paper were performed at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI).

Fig. 1 .

 1 Fig. 1. Example of a planetesimal disk without outward planet migration, nor inward P-R drag migration of the test particles, according to our numerical simulations as explained in Section 2. The star and planet locations, projected onto the orbital plane of the planet, are represented by large red points, and the planet orbit by a thin green line. The initial planetesimal disk consists of 50 000 planetesimals distributed between 40 and 75 AU, with the surface density distribution proportional to r -1 . Although some planetesimals are trapped in MMRs with the planet, they are not sufficiently numerous to generate spatial structures (besides the 1:1 MMR). See the electronic edition of the Journal for a color version of this figure.

Figures 2

 2 Figures 2 to 5 show examples of results obtained with our numerical model. It appears that, with a planet on a loweccentricity orbit, the planetesimals trapped in MMRs are numerous and dominate the shape of the disk. Four important factors must be taken into account to determine which resonances govern the aspect of the structures in the disk:

Fig. 2 .

 2 Fig. 2. Spatial distribution of planetesimals for a planet on a strictly circular orbit. The star and planet locations, projected onto the orbital plane of the planet, are represented by larger red points, and the planet orbit by a thin green line. The 4 rows correspond respectively to Earth mass, Neptune mass, Saturn mass and 3 Jupiter mass planets, from top to bottom. The 3 columns show the disk after 5, 15 and 40 Myr. The initial planetesimal disk consists of 50 000 planetesimals distributed between 40 and 75 AU, with the surface density distribution proportional to r -1 . See the electronic edition of the Journal for a color version of this figure.

Fig. 3 .

 3 Fig. 3. Same simulations as in Fig. 2, with the planetesimals represented in a (semi-major axis, eccentricity) plane. The plotting conventions are the same as in Fig. 2.See the electronic edition of the Journal for a color version of this figure.

Fig. 4 .

 4 Fig. 4. Same as Fig. 2, for similar planets, but on a low-eccentricity orbit (e p = 0.05). See the electronic edition of the Journal for a color version of this figure.

Fig. 5 .

 5 Fig. 5. . Same as Fig. 3, for similar planets, but on a low-eccentricity orbit (e p = 0.05). See the electronic edition of the Journal for a color version of this figure.

  Hamiltonian by the averaging process. The secular Hamiltonian of the planar problem turns out to have only one degree of freedom. It depends for instance only on the eccentricity of the planetesimal e and of its longitude of periastron ν with respect to that of the planet. For any given fixed values for a and for the planet eccentricity e p , we can draw a phase portrait (i.e., level curves of Hamiltonian) of the dynamics in an (ν, e) plane.

Fig. 6 .

 6 Fig. 6. Same as Fig. 2, for similar planets, but on a moderate eccentricity orbit (e p = 0.3). See the electronic edition of the Journal for a color version of this figure.

Fig. 8 .

 8 Fig. 8. Phase portraits (Hamiltonian level curves in an (ν, e) plane) of the secular non-resonant dynamics in the planar restricted three body problem, for a fixed semi-major axis a = 1.3a p , and planet eccentricities e p = 0.5 (left plot) and e p = 0.1 (right plot). The red line separates regions where the planetesimal orbit does not cross that of the planet (low eccentricities for e p = 0.1 and a small island around ν = 0 for e p = 0.5) from regions where both orbits cross each other.See the electronic edition of the Journal for a color version of this figure.

Fig. 9 .

 9 Fig. 9. An Earth mass planet on a very eccentric orbit (e p = 0.6). The plotting conventions are the same that Fig. 1. The disk structure does not rotate with the planet and is spatially fixed. See the electronic edition of the Journal for a color version of this figure.

Fig. 10 .

 10 Fig. 10. Overview of the simulation outputs for planetesimal disks, in a (planet mass, planet eccentricity) plane. All simulations shown in this figure concern a planet with a constant migration rate of 0.5 AU Myr -1 . Color scale indicates the fraction of planetesimals still bound after 40 Myr. The fraction of surviving planetesimals is linearly interpolated between the simulations of the Table 2. In zones I and II, observable structures in the disk are generated by MMRs while in zone III, transient structures are generated by non-resonant mechanisms. In the remaining region, the disk does not show any structure. In zone I, MMRs create clumpy disks while in zone II they generate a smooth disk with a hole at the planet location. See the electronic edition of the Journal for a color version of this figure.

Fig. 11 .

 11 Fig. 11. A Saturn mass planet on a circular orbit migrating outward a dynamically warm disk. Initial eccentricities of planetesimals are uniformly distributed between 0 and 0.2. The plotting conventions are the same as in Fig. 2 and 3. See the electronic edition of the Journal for a color version of this figure.

  -1 . b In Jovian mass. c Fraction of surviving planetesimals at the end of the simulation, i.e. 40 Myr. d As in Fig. 10.

  -1 . b In Jovian mass. c Fraction of surviving planetesimals at the end of the simulation, i.e. 4 Myr for 5 AU Myr -1 migration rate and 200 Myr for 0.05 AU Myr -1 . d As in Fig. 10.

Fig. 12 .

 12 Fig. 12. Capture in the 2:1 resonance for a Neptune mass planet (top), en Earth mass planet (middle) and a Saturn mass planet (bottom), all on a circular orbit, after 5 Myrs. Dashed lines show the position of the resonance at the beginning and at the end of the simulations. In those simulations, an Earth mass planet does not capture at all planetesimals while a Saturn mass planet traps all of them. A Neptune mass planet only traps a fraction of the planetesimals but gives a small kick in eccentricity for the others. The plotting conventions are the same as in Fig. 3. See the electronic edition of the Journal for a color version of this figure.

Fig. 13 .

 13 Fig. 13. Projection on the (semi-major axis, eccentricity) plane of the trajectory of 40 planetesimals near the 2:1 resonance of a Neptune mass planet, migrating on a circular orbit. The semimajor axis is in units of the planet semi-major axis. Green lines are for permanently trapped planetesimals while the red lines are for the temporarily perturbed planetesimals. See the electronic edition of the Journal for a color version of this figure.

Table 1 .

 1 Summary of recent papers on particles trapping in MMRs with a planet. Planetary mass in Jovian mass divided by stellar mass in solar mass.

	Authors	Method	Planet parameters	Migration Notes
			Mass ratio a eccentricity
	Wyatt (2003)	semi-analytic 0.0003 to 3	0	planet	Forced migration
	Wyatt (2006)	semi-analytic			none	This work extends the previous one to smaller
						particles sensitive to radiation pressure.
	Kuchner & Holman (2003)	analytic	0.005 to 15	0 to 0.6	particles	Only resonant particles trapped during migra-
						tion due to P-R drag.
	Deller & Maddison (2005)	numerical	0.01 to 3	0 to 0.7	particles	Particles migrate due to P-R drag and solar
						wind.
	Moro-Martín & Malhotra (2002)	numerical	0.05	0	particles	Study of the Kuiper Belt.
	This paper	numerical	0.001 to 3	0 to 0.7	planet	Only planetesimals disk, forced migration

a

  Mig. rate a Mass b Ecc. Surv. planetesimals c Disk shape d Mig. rate a Mass b Ecc. Surv. planetesimals c Disk shape d

	0.5	0.0035	0.0	100%	None	0.5	0.05	0.0	100%	I
			0.01	100%	None			0.01	100%	I
			0.05	100%	None			0.05	100%	II
			0.1	100%	None			0.1	100%	None
			0.2	100%	None			0.2	50%	None
			0.3	100%	None			0.3	10%	None
			0.4	35%	III			0.4	5%	None
			0.5	10%	III			0.5	10%	None
			0.6	5%	III			0.6	5%	None
			0.7	5%	III			0.7	5%	None
	0.5	0.33	0.0	75%	I	0.5	1	0.0	70%	I
			0.01	70%	I			0.01	70%	I
			0.05	85%	II			0.05	55%	II
			0.1	75%	None			0.1	25%	II
			0.2	25%	None			0.2	10%	None
			0.3	10%	None			0.3	10%	None
			0.4	10%	None			0.4	10%	None
			0.5	10%	None			0.5	5%	None
			0.6	5%	None			0.6	5%	None
			0.7	5%	None			0.7	0%	None
	0.5	2	0.0	65%	I	0.5	3	0.0	55%	I
			0.01	60%	I			0.01	50%	I
			0.05	30%	I			0.05	20%	I
			0.1	15%	II			0.1	10%	I
			0.2	5%	None			0.2	0%	None
			0.3	5%	None			0.3	0%	None
			0.4	5%	None			0.4	0%	None
			0.5	5%	None			0.5	0%	None
			0.6	0%	None			0.6	0%	None
			0.7	0%	None			0.7	0%	None
	a In AU Myr

Table 3 .

 3 Same as Table 2, but for different migration rates. Mig. rate a Mass b Ecc. Surv. planetesimals c Disk shape d Mig. rate a Mass b Ecc. Surv. planetesimals c

	Disk shape d

Table 4 .

 4 Same as Table2but for simulations with initially excited disks. The maximum initial eccentricity of the planetesimals is mentioned for all simulations. The migration rate is 0.5 AU Myr -1 for all simulations. Max. ecc. a Mass b Ecc. Surv. planetesimals c Disk shape d Max. ecc. a Mass b Ecc. Surv. planetesimals c Disk shape d Fraction of surviving planetesimals at the end of the simulation, i.e. 40 Myrs.

	0.1	0.035	0.0	100%	None	0.1	0.05	0.0	100%	II
			0.01	100%	None			0.01	100%	None
			0.05	100%	None			0.05	100%	None
			0.1	100%	None			0.1	100%	None
			0.2	100%	None			0.2	55%	None
	0.1	0.33	0.0	70%	II	0.1	1	0.0	60%	I
			0.01	70%	II			0.01	55%	II
			0.05	85%	II			0.05	50%	II
			0.1	70%	None			0.1	25%	II
			0.2	20%	None			0.2	5%	None
	0.1	2	0.0	50%	I	0.1	3	0.0	40%	I
			0.01	45%	I			0.01	40%	I
			0.05	25%	I			0.05	15%	I
			0.1	15%	II			0.1	10%	I
			0.2	5%	None			0.2	0%	None
	0.2	2	0.0	35%	I	0.2	3	0.0	30%	I
			0.01	35%	I			0.01	25%	I
			0.05	20%	I			0.05	10%	I
			0.1	10%	I			0.1	5%	I
			0.2	5%	None			0.2	0%	None
	a For planetesimals.

b In Jovian mass. c d As in Fig. 10.

Acknowledgements. We are grateful to Philippe Thebault, Alexander V. Krivov and Martina Queck for enlightening discussions about dust collisions in debris


F Petit

Y Phulpin

M Saguan

Ph Dessante
email: philippe.dessante@supelec.fr

A CONTRIBUTION OF EXPERIMENTAL ECONOMICS TOWARD CHARACTERIZATION OF THE USE OF MARKET POWER IN OLIGOPOLISTIC ELECTRICITY MARKETS

Despite the numerous researches about imperfect competition, the market power remains difficult to quantify using traditional economics methods. In this paper, we propose an experimental economics design and outline some ways of analysis of its results toward characterization of the use of market power. A simple system with two regions and a limited interconnection transfer capacity allocated by an implicit auction is studied. Depending on the experiments two or three subjects share equitably the production capacity in one region, while the production capacity is equitably shared among 5 subjects leading to a more competitive situation in the second one. In both regions, we observe a market price that is different from the theoretical results allowing a quantification of the use of market power. Results are also analyzed based on a characterization of the subjects' behaviour. Further the impact of subjects' behaviour on the market price evolution is described.

Introduction

The ongoing liberalisation process of electricity markets has caused an important evolution of the industry structures over the world. A major goal of the reforms was to achieve economic efficiency and low electricity prices by introducing competition in former regional electricity markets. As way of example, a trend toward development of transmission interconnection capacities in order to reach a higher degree of competition in regional electricity markets [START_REF] Boucher | Towards a common European electricity market[END_REF]] has been observed in Europe. In this context, interconnection capacities play not only the physical role of connector between regional transmission systems but also the role of threatening competitors and therefore promoting competition.

The electricity market integration requires well designed rules in order to allow efficiently scarce transmission capacities among economic agents. Two main methods have thus been applied in Europe to cope with limited interconnection capacities between countries: i) implicit auctions, where energy and transmission are priced together and ii) explicit auctions, where energy and transmission are priced separately [Ehrenman-Smeers (2005)]. These methods are efficient if the "perfect competition" assumption is respected. However, recent experience has shown that this prerequisite is usually not true.

In fact the recent merging of generation utilities has reinforced the oligopolistic context remaining from former regional monopolies. Most of European regional electricity markets have thus reached a high level of market concentration with for instance monopoly in France, Belgium and Greece, duopoly in Spain, Finland, and Sweden or triopoly in United Kingdom [START_REF] Glachant | A Competitive Fringe in the Shadow of a State Monopoly: the Case of France[END_REF], [START_REF] Newbery | Electricity liberalisation in Britain: the quest for a satisfactory wholesale market design[END_REF][START_REF] Von Der Fehr | The Nordic market: signs of Stress[END_REF], [START_REF] Crampes | The Spanish Electricity Industry: Analysis and Prospects[END_REF]]. There is a large consensus among economists to state that these types of oligopoly induce an important market power for its participants. This notion of market power could be defined as the ability to alter profitably prices away from competitive levels [START_REF] Mas-Collel | Microeconomic Theory[END_REF]]. The high market concentration is thus a characteristic of imperfect competition and a serious threat for the economic efficiency of electricity markets. There has been an increasing research activity to study the interaction between market power in national electricity market, limited interconnection capacities and allocation methods.

Numerous economic studies have already covered the theme of quantification of the use of market power in presence of transmission constraints. A traditional theoretical approach mainly based on Game Theory has been used to study imperfect competitive set up in electricity markets with transmission constraints [START_REF] Cardell | Market Power and strategic Interaction in electricity networks[END_REF], [START_REF] Borenstein | The competitive effects of transmission capacity in a deregulated electricity industry[END_REF], [START_REF] Hobbs | Linear Complementarity Models of Nash-Cournot Competition in Bilateral and POOLCO Power Markets[END_REF], [START_REF] Ralph | EPECs as models for electricity markets[END_REF]]. According to most of game theory approaches, it is supposed that all market participants have perfect information and perfect rationality for computing the market (Nash) equilibrium and that they agree in playing it. However, there is still no satisfactory method for quantifying the exercise of market power in real systems, nor to evaluate the market power for any market designs. The main drawbacks of the theoretical approach are the lack of modelling bounded rationality and the lack of consideration of idiosyncrasies of market designs. A (fairly new) research method that solves these drawbacks is the experimental approach.

Experimental economics shows some properties that are decisive for a better understanding of the use of market power in electricity market. The experimental technique, which allows taking into account human decision making in a complex situation, is described in [START_REF] Smith | Economics in the Laboratory[END_REF]]. As long as a set of conditions is respected, for instance if the experiment is simple enough for players to understand their impact and for economists to analyze results properly, experimental economics can provide useful information in seven cases listed by Smith; two reasons out of the seven justify the use of this method in this study, namely to compare environments and to establish empirical regularities as a basis for a new theory.

In this paper an experimental approach is applied for measuring the exercise of market power in a simple system with two regions and a limited interconnection transfer capacity allocated by an implicit auction. Results from experiences are compared with theoretical results allowing a first quantification of the use of market power. Then, we propose a definition of some types of market behaviours, namely marginal-bidder, follower and leader, and we analyze the market price and its evolution depending on the behaviour of players. An accurate characterization of players' behaviour is indeed necessary for developing adequate regulation of new electricity markets, which could also be evaluated through experimental economics.

This paper is organized as follows: Section 2 describes the experimental design and the theoretical results that should be obtained according to classic game theory. Section 3 outlines the experimental process. A complete overview of the results is then presented in section 4. Finally conclusions and propositions for further works are presented in section 5.

Experimental design and assumptions

In this section, we set the characteristics of the particular case that was studied and experimented. We start by defining the market structure, then we draw up the market design and theoretical results are finally presented.

Market structure

The market structure is presented in Figure 1. Producers and consumers are located within two zones, which are most often linked by an interconnection 3 . Subjects control only producers, while consumers' actions are determined by a computer according to a predefined demand curve. The available production is 150 MWh in each zone. This amount is equitably divided among 5 producers in the 1 st zone and 2 (experimental design A) or 3 (experimental design B) producers in the 2 nd zone. Production costs depend on the zone. For each producer, half of his production has a low marginal cost (10 and 15 €/MWh in zone 1 and 2, respectively) and the other half has a high marginal cost (20 and 25 €/MWh in zone 1 and 2, respectively). All numerical details are provided in Table I The base case is experimental design A. It features a more competitive zone in which 5 producers own a production capacity of 30 MWh at a relatively low cost, and an oligopolistic zone in which only 2 producers can produce up to 75 MWh at a higher cost. Demand profiles are chosen in order to control the producers' ability to use market power unilaterally or not; additionally, they are drawn so that the line be congested under competitive assumption in phase 2, but not in phase 1. The ability to exercise market power is affected in the experimental design B by splitting of the 2 producers in zone 2 into 3.

Market design

The interconnection capacity is allocated following the method known as implicit auction (or market coupling mechanism [ETSO-EuroPEX (2004)]). In this context, each zone has a local market, where only producers that are physically located in the same zone can bid through a uniform price sealed auction.

Once offers have been submitted, aggregated bid and ask curves are computed on each market to calculate the two equilibrium prices PA and PB as represented in Figure 3 When the price PB is higher than PA, a quantity ∆Q is exchanged from zone A toward zone B. This export is equivalent to an additional ask ∆Q at any price in market A while the corresponding import is equivalent to an additional bid ∆Q at any price in market B. This quantity ∆Q is progressively raised from zero until one of the following cases is reached:

-1st case: no congestion:

As we can see in figure 4, the quantity exchanged ∆Q = QA* -QA = QB* -QB is smaller than the Available Transfer Capacity of the interconnection A-B: there is thus only one price P* = PA* = PB* for both markets. In both cases, once equilibrium prices have been computed, bids submitted in market X are:

-fully accepted when the offered price is below PX*;

-fully rejected when the offered price is above PX*;

-partially accepted when the offered price equals PX*.

Theoretical results

Competitive and (static) Cournot equilibrium are computed for each experimental design and each phase following [START_REF] Hobbs | Linear Complementarity Models of Nash-Cournot Competition in Bilateral and POOLCO Power Markets[END_REF]. Tables 3 and4 show results corresponding to experimental designs A and B respectively. Given that the only difference between them is market structure ("concentration") in zone 2, results corresponding to the perfect competition case are equivalent.

In these cases the difference between phase 1 and 2 is the occurrence of congestion. In phase 1, both zones have the same price as there is no congestion. In phase 2, price in zone 2 is higher than price in zone 1; the interconnection is congested.

In Cournot case, market structure is important and this can be seen comparing results from experimental design A and B4 . Cournot prices of experimental design A corresponding to more concentrated market structure in zone 2 are higher than prices of experimental design B with a less concentrated market structure. Moreover, the price obtained when only one player exerts his market power is determined; "NP" means that no producer has market power on this zone, namely that no player can unilaterally and profitably raise the market price above the competitive one. It can be noticed that in phase 2, any producer in zone 1 is pivotal (as long as the market price is below 60 €/MWh): under these circumstances, unilateral exercise of market power in zone 1 leads to a leveling of prices in both zones.

Experimental process

9 experiences were run with different settings and subjects from various origins, as shown in Table 5.

The experimental economics laboratory

All experiments were led in the laboratory for experimental economics in Supelec, provided with fifteen isolated computers so that no subject can communicate with each other. Subjects pass their offers and receive results through a web interface represented in Figure 6. All computations and data storage are performed on a server, allowing any numerical analysis of the results ex-post.

Instructions

Every experience is foregone by thirty minutes of instructions; the market structure is described and the market design is detailed. Subjects are also taught how to use the interface. Two demonstration periods close the instruction phase in order to ensure that subjects are familiar with the interface.

Subjects are given the phase parameters (demand curves and Available Transfer Capacities) and they are told when each phase begins. They know the total amount of production capacity available per zone (150 MWh) but they are not told others' costs and capacities.

Experimental progress

The experiment itself lasts from one hour to one hour and a half depending on the number of periods to run. During first periods, subjects have 60 seconds to submit their bids; this term is progressively decreased while subjects are getting trained; however time is never binding during the experience.

Subjects can submit and modify as many bids as many times as they want; nonetheless offered quantities and prices must be integers not only for computational purposes, but also to accelerate price convergence. In addition to information on the market structure, subjects have access to their last three periods' profits, their aggregated profits, results for each previous period bids (fully accepted, partially accepted or rejected) and the price evolution on each market (zone 1 or 2) from the beginning of the experience. There is no transaction cost and production costs are only due for sold units.

First results

Comparison of results with theory

The market price obtained in each zone at the end of each phase is reported in the tables below; when no convergence was reached, the price is the average of the last 20 periods and is indicated by a * . As it is explained in the following part, a market price of 21 instead of 20 €/MWh is still close to the competitive price. Thus, in the first phase, the market is often close to perfect competition; in the second phase, results are quite similar to the unilateral use of market power. In accordance with the results predicted in part 2.3, market prices in both zones stay rather close, particularly in phase 2.

However, the market prices are far from (static) Cournot's model expectations, and it is necessary to study individual behaviours in order to understand this divergence.

A characterization of subjects' behaviour

In order to characterize subjects' behaviour, we propose a classification based on offer functions. Experimental economics offer the possibility to observe the exercise of market power by comparing bids with production costs. As way of example, figure 8 represents the percentage of quantity offered over all periods and experiences for a given markup (equals to offered price minus production cost); all withdrawn units are represented as having a markup of 100 €/MWh. This figure tends to show a massive exercise of market power. However, a more accurate study of the way subjects bid indicates that, since the price can never be below 20 €/MWh, the bids associated with the cheaper part of their production is not significant as far as subjects' behaviour is concerned: the corresponding price can range from their production cost (10 or 15 €/MWh) to the competitive price (20 or 25 €/MWh) without any particular reason nor any consequence on the market price. Thus, only the second half of their available units has been taken into account in figure 9.

Since the offered price must be an integer, a markup equal to 1 is considered as a competitive behaviour linked to experimental design: the subject's profits are similar if his bid is accepted at a price of 20 €/MWh or if it is rejected, it would be zero in both cases. As a consequence he may bid just above his marginal cost. Even under this condition, hardly half of the more expensive units are offered at a competitive price on zone 1, and less than a third on zone 2.

Markup for the most expensive units

Figure 9: Price markup by zone

To build their offer functions, experiences show that subjects do not only refer to their marginal cost, but also to the previous period market price. Thus figure 10 represents the distribution function of the difference between the offered price and the previous market price, without taking into account competitive bids.

Price difference for the more expensive units and non-competitive bids This figure displays two different behaviours for bids which are not offered at their marginal cost. First, the price can be chosen above the previous market price, or even withdrawn as this is the case for almost half of the non-competitive bids on zone 2; this behaviour, and particularly capacity withdrawal, is an attempt to raise prices and is a feature of a price-marking, or leading, strategy.

Second, half of the non-competitive bids on zone 1 is offered at the previous market price minus 1 or 2 €/MWh. While the competitive behaviour is a price-taking strategy ensuring the player he will sell whenever it is profitable for him and the leader behaviour is a price-making strategy allowing the player to raise market price at a risk of not selling his power, this third one can be qualified as a follower behaviour: by bidding just below the previous price, the subject aims at preventing a price drop whilst limiting the risk for his offer to be rejected. This strategy is a clear incentive for a leader to raise the price punctually even if it is non-profitable in short term: the price decrease will be slow enough for him to earn more during the next periods than what he loses by exercising market power. Subjects can then reach an implicit collusion similar to the one observed in a repeated Stackelberg game, in which the follower plays functions of the leader's previous period decision. An even stronger implicit collusion was observed during the experiment A6 in which the two producers of the second zone implicitly agreed to bid at the monopoly price.

Impact of players' behaviour on the market price and its evolution

An experiment involving a leader and several followers results in prices significantly higher than the competitive price; moreover, price does not converge, but a dynamic equilibrium is found: price variation is a succession of peaks followed by slow decreases (over five to ten periods).

In particular, when the interconnection is such that all the excess production in zone 1 can be exported into zone 2, congestion leads to a raise of the market price in zone 1 up to its level in zone 2. Simultaneously, as soon as the interconnection is congested, a lack of competition on zone 2 results in a price increase, from which a vicious circle may ensue.

Conclusion

Although the experimental design that is presented in this paper is extremely simple, it contributes to a better understanding of the use of market power. This analysis is supported by a classification of the subjects' bids. Three main types of bids could indeed be clearly identified. This possible quantification of each type of bid is an interesting vector toward understanding the evolution of the market price. Future works are expected to depict more precisely the evolution of market price depending on the quantity of bids in each of the three classes.

As its results are significantly different than those obtained using more classical theories, experimental economics could also be useful in order to evaluate the real efficiency of regulation designs.

Figure

  Figure 1: Market Structure

  Figure 2: Demand profiles

Figure 3 :

 3 Figure 3: Market clearing mechanism

Figure 4 :

 4 Figure 4: Market clearing by non congested interconnection -2nd case: congestion: As presented in figure 5, the quantity exchanged ∆Q = QA* -QA = QB* -QB is equal to the Available Transfer Capacity of the A-B interconnection: there is one price for each market and (PB* -PA*) x ATC is the congestion rent.

Figure 5 :

 5 Figure 5: Market clearing by congested interconnection

Figure 6 :

 6 Figure 6: Web interface for experiment subjects

  Figure 8: Price markup for all the experiences

Table 1 : Numerical data beside producers

 1 .Each experiment counts two or four phases; during which, the demand bids and the exchange capacity are maintained constant. These parameters are exposed in Table2.

	Zone	Experimental design	Unit type	Marginal production cost (€/MWh)	Available quantity (MWh)
	1	A and B	Low cost High cost	10 20	15 15
	2	A B		Low cost High cost Low cost High cost	15 25 15 25	38 37 25 25
	Experimental design reference	Phase	Duration (number of periods)	Zone 1	Demand curve	Zone 2
	A and B		1 2	30 30		Base load Peak load

Table 2 : Numerical data beside demands and interconnection

 2

	ZONE 1	Maximum exchange: 30 MWh	ZONE 2 2/3 producers
	5 producers		150 MWh
	150 MWh

3

The way this interconnection is run is described in part 2.2: Market design.

Table 3 : Price equilibrium for experimental design A

 3

		Phase 1 (no congestion)	Phase 2 (congestion)
		Price zone 1	Price zone 2	Price zone 1	Price zone 2
	Perfect competition	20,0	20,0	20,0	25,0
	Cournot	43,6	44,0	47,3	52,0
	Unilateral market power	NP	25,0	25,0	NP
		Phase 2 (no congestion)	Phase 2 (congestion)
		Price zone 1	Price zone 2	Price zone 1	Price zone 2
	Perfect Competition	20,0	20,0	20,0	25,0
	Cournot	40,8	40,8	47,3	50,0
	Unilateral market power	NP	25,0	25,0	NP

Table 4 : Price equilibrium for experimental design B

 4

Table 5 : Distribution of subject by experimental design and by origins

 5

	Designation	Experimental design	Subjects origin	Number of experiences led	Number of period per experience
	A1, A2, A3		Students in engineering	3	60
			PhD students & faculty
	A4		members (non-	1	60
		A	economists)
	A5		PhD students (economists)	1	60
	A6, A7		Students in engineering	2	40 5
	B1		Students in engineering	1	60
	B2	B	PhD students & faculty members (non-	1	60
			economists)

Table 6 : Results for experimental design A

 6

		Phase 1 (no congestion)	Phase 2 (congestion)
		Price zone 1	Price zone 2	Price zone 1	Price zone 2
	Perfect competition	20,0	20,0	20,0	25,0
	Cournot	40,8	40,8	47,3	50,0
	Unilateral market power	20,0	25,0	25,0	25,0
	B1	21	21	31,9 *	36,2 *
	B2	20,8 *	20,9 *	24,3 *	25,6 *

Table 7 : Results for experimental design B

 7

Difference between offered price and previous market price

			50,0%
			45,0%
			40,0%
			35,0%
			30,0%
	Frequency		25,0%						Zone 1 Zone 2
			20,0%
			15,0%
			10,0%
			5,0%
			0,0%
	-60	-40	-20	0	20	40	60	80	100
			Offer -previous market price (€/MWh)
	Figure 10:

F. Petit and M. Saguan are with GRJM research group at University of Paris XI. An up-dated version of this article will be available at the GRJM web site (www.grjm.net).

The authors want to thank the "Conseil Régional d'Île de France" for their support and for funding the experimental facilities.

It is important to note that Cournot equilibrium computations have been done under the assumption that players do not try to "game" transmission pricing[START_REF] Ralph | EPECs as models for electricity markets[END_REF]].