Hirsh Majid

Hassane Abouaïssa

Daniel Jolly

Gildas Morvan

Real-Time Dynamic Traffic Routing Using Variable Structure Control

Keywords: Macroscopic model, variable structure control, sliding mode control, dynamic traffic routing

In this paper, we present a new algorithm to deal with the real-time dynamic traffic routing (DTR) problem using a second order traffic flow model (METANET). We investigate variable structure control (VSC) as a high-speed switched feedback control resulting in sliding mode. The control objective is to minimize the difference of the travel times on alternate routes in a network setting. Accordingly, the congestion could be reduced/eliminated. The relevance of the proposed algorithm is demonstrated via a set of numerical simulations in the case of the same as well as different geometric conditions.

I. INTRODUCTION

Traffic control represents an efficient way to improve the freeway throughput and to ensure an efficient, safe and less polluting transportation of goods and persons [START_REF] Kostialos | Motorway network traffic control systems[END_REF]. It also contributes to a high reduction of direct and indirect costs.

Freeway traffic control can be achieved via a set of actions and measurements such as: dynamic speed limits, route guidance, dynamic assignment and dynamic traffic routing, ramp metering, etc.

Dynamic traffic assignment (DTA) and dynamic traffic routing (DTR) represent some of the most efficient solutions to steer congestion problems especially the no-recurrent ones. DTR/DTA or control of traffic diversion rests on the determination of time-dependent split variables at the diversion point in order to achieve a user-equilibrium traffic pattern [START_REF] Kachroo | Solution to the user equilibrium dynamic traffic routing problem using feedback linearization[END_REF] in a way that one route not be overloaded and the others underused. Accordingly, the congestion could be alleviated/eliminated. Several algorithms have been developed to solve this problem. Most of them rely on optimization approaches such as those proposed by [START_REF] Petta | Multiple users classes real time traffic assignement for on-line operations: a rolling horizon solution framework[END_REF]. As stated in [START_REF] Kachroo | System dynamics and feedback control problem formulations for real-time dynamic traffic routing[END_REF], the proposed optimization algorithm attempts to solve the DTR problem by optimizing the objective functions for the nominal model over the planning horizon and it is not adapted for on-line DTR control. Papageorgiou [START_REF] Papageorgiou | Dynamic modeling, assignment and route guidance traffic networks[END_REF] and Messmer & Papageorgiou [START_REF] Messmer | Route diversion control in motorway networks via nonlinear optimization[END_REF] have proposed algorithms for dynamic traffic assignment (DTA) based on the linear quadratic regulator, and nonlinear optimization techniques. Other techniques use expert systems to deal with the diversion problem [START_REF] Gupta | Development of prototype kbes in the management of non-recurrent congestion of the massachusetts turnpike[END_REF]. Liu & al. [8] have proposed a strategy based on a model reference adaptive control in order to guide the real-world traffic flow to evolve towards the desired states 1 Univ. Lille Nord France, F-59000 Lille, France. U-Artois, LGI2A, EA. 3926 Technoparc Futura, F-62400 Béthune, France.

(hassane.abouaissa, daniel.jolly, gildas.morvan)@univ-artois.fr. [START_REF] Kachroo | Real-time dynamic traffic routing based on fuzzy feedback control methodology[END_REF] have designed feedback and fuzzy control laws for an on-line diversion problem. Another algorithm for DTR and DTA problems based on the nonlinear H ∞ was proposed in [START_REF] Kachroo | Feedback control solutions to network level user-equilibrium real-time dynamic traffic assignement problems[END_REF]. See also [START_REF] Karimi | Integration of dynamic route guidance and freeway ramp metering using model predictive control[END_REF] and [START_REF] Pavlis | Simple decentralized feedback strategies for route guidance in traffic networks[END_REF] for more information about route guidance strategy.

In this paper, we exploit the variable structure control theory (VSC). The control algorithm consists of two parts; trajectory planning (open loop control) which is achieved using a specific class of complex systems called 'Differentially flat systems'. The trajectory tracking (closed loop control) is ensured by a high-speed switched feedback control resulting in sliding mode. In contrast of [START_REF] Abouaïssa | Flatness based control for real-time dynamic traffic routing problem[END_REF], the non-destination oriented METANET model [START_REF] Messemer | [END_REF] [START_REF] Messemer | Metanet: a macroscopic simulation program for motorway networks[END_REF] is used in the design of the control algorithm.

The paper is organized as follows: Section II presents the main principles of the DTR problem and its mathematical formulation. Section III recalls the main definitions of the control system. Section IV shows the control design methodology using the concept of flatness and sliding mode control. Section V provides some numerical simulations for a sample network. Section VI concludes the paper and outlines some tracks for further developments.

II. DYNAMIC TRAFFIC ROUTING FORMULATION

In this section, we first present a mathematical formulation which is used for the design of DTA/DTR flat controller. We expose a second order traffic flow model ( METANET).

For the sake of simplicity, we consider the case of two alternate routes divided into n 1 and n 2 sections, respectively, as depicted in Fig. 1.

The proposed control approach rests on the use of METANET model. The dynamic equation of the traffic density reads:

ρi,j (t) = 1 L i,j [q i,j-1 (t) -q i,j (t)] (1)
where (i, j)=(1, 1), (1, 2), . . . , (1, n 1 ), (2, 1), . . . , (2, n 2 ). L i,j is the section length and ρ i,j represents the traffic density at section j of the route i.

The relation between the traffic flow q i,j and the mean speed v i,j is:

q i,j (t) = ρ i,j (t)v i,j (t)λ i,j (2)
where λ i,j is the number of lanes. v i,j is the dynamic mean speed of section j on the route i.

The equation of dynamic mean speed is:

vi,j (t) = 1 τ (V ρi,j (t) -v i,j (t)) + 1 Li,j v i,j (t)(v i,j-1 (t) -v i,j (t)) -ν τ Li,j
ρi,j+1(t)-ρi,j (t) ρi,j (t)+κ

(3) V ρi,j is defined as a nonlinear expression called fundamental diagram mean speed (see [START_REF] May | Traffic Flow Fundamentals[END_REF]):

V ρi,j (t) = vf i,j exp - 1 a ρ i,j (t) ρ ci,j a (4)
ρ ci,j , vf i,j , a, ν, κ and τ are constant parameters which reflect particular characteristics of a given traffic system [START_REF] Messemer | [END_REF].

The control input α(t) ∈ [0, 1] is defined as an exogenous variable of the system [START_REF] Pavlis | Simple decentralized feedback strategies for route guidance in traffic networks[END_REF], it is a split rate that allows to reach a user equilibrium traffic pattern. For the case of two alternate routes, the equations of entry flow in each route are:

α(t)q e (t) = q 1,in (t) (1 -α)q e (t) = q 2,in (t) (5)
where q e (t) is the traffic demand. The control objective is to find the optimal split rate α in order to minimize the differences J(α) between the travel time T T of the two alternate routes. Kachroo and Özbay [START_REF] Kachroo | System dynamics and feedback control problem formulations for real-time dynamic traffic routing[END_REF] have formulated the DTR problem as follows: finding α 0 and the optimal α(t), which minimizes

J(α) = tf 0   n1 j=1 T T (ρ 1,j ) - n2 j=1 T T (ρ 2,j )   2 dt (6)
where T T (ρ 1,j ) and T T (ρ 2,j ) are the travel times function of section j on the route 1 and 2 respectively. t f is the final time.

In the general case, consider a traffic flow system with n alternate routes. The system is described by the same ordinary differential equation [START_REF] Kostialos | Motorway network traffic control systems[END_REF] where, (i, j) = ((1, 1), . . . , (1, n 1 ), (2, 1), . . . , (2, n 2 ), (n, 1), . . . , (n, n n )).

The problem consists then in finding a set of split variables, α n-1 , where, n-1 i=1 α i = 1, minimizes total travel time T T , (see [START_REF] Kachroo | System dynamics and feedback control problem formulations for real-time dynamic traffic routing[END_REF], for optimal formulation of the DTR problem for a general case)

         α 1 (t)q e (t) = q 1,in (t) α 2 (t)q e (t) = q 2,in (t) . . . (1 -α 1 -α 2 -, . . . , -α n-1
)q e (t) = q n,in (t) [START_REF] Gupta | Development of prototype kbes in the management of non-recurrent congestion of the massachusetts turnpike[END_REF] Note that for such a complex problem, an open loop control structure is not sufficient [START_REF] Kachroo | Feedback control solutions to network level user-equilibrium real-time dynamic traffic assignement problems[END_REF], it calls for a more robust feedback control one. In this context, we highlight the interest of an attractive approach based on the concept of differential flatness [START_REF] Fliess | Sur les systèmes non linéaires différentiellements plats[END_REF] [START_REF] Fliess | Flatness and defect of non-linear systems: Introductory theory and examples[END_REF]. The following section, recalls the main definitions and principles of differentially flat systems.

III. SLIDING MODE FLATNESS-BASED CONTROL SYSTEMS

In this paper, we exploited one of the most attractive methods that can be applied to a broad class of nonlinear systems resulting in controllers that are robust to modeling errors and unknown disturbances. We investigated variable structure control (VSC) as a high-speed switched feedback control resulting in sliding mode. The gains in each feedback path switch between two values according to a rule that depends on the value of the state at each instant. The purpose of the switching control law is to derive the non-linear system's state trajectory onto a prespecified surface in the state space and to maintain the system's state trajectory on this surface (switching surface) for subsequent time.

In standard sliding mode control, or first order sliding mode control (FOSMC) [START_REF]Sliding Mode Control in Engineering[END_REF] and [START_REF] Rashid | Robust State Estimation and Control of Highway Traffic Systems[END_REF] , the sliding surface is chosen so that it has a relative degree of one with respect to the control input. In such a case, the control input u is, for example, of the form; u(t) = u eq (t) + u d (t), where u eq is the continuous function and u d (t) = -K d sign(s(t)) -K p s(t) is the discontinuous function. K d specifies the speed of convergence of the closed-loop system in order to s(t) = 0. To combine a small switching gain with fast convergence, the discontinuous control term could be extended with a proportional feedback term K p s(t) [START_REF] Monsees | Discrete-Time Sliding Mode Control[END_REF].

u eq ensures that ṡ = 0. It is the equivalent control, it depends only on the switching surface s(t) and not on the control function u d [START_REF] Rashid | Robust State Estimation and Control of Highway Traffic Systems[END_REF]. Thus, it can be said that u eq introduces the trajectory planning (open loop control) of the output. In this context, the term u d should introduce the trajectory tracking (closed loop control) of the output.

In this paper, instead of u eq , the trajectory is made using a specific class of control system called 'differential flatness' .

The concept of flat systems was first introduced by Fliess and al [START_REF] Fliess | Sur les systèmes non linéaires différentiellements plats[END_REF] [18] more than a decade ago using the formalism of differential algebra (see [START_REF] Ph | Contribution l'étude des systèmes Différentiellement plat[END_REF] for a slightly different approach of differentially flat systems). This special class of non-linear control systems described by ordinary differential equations: differentially flat systems form a special class of nonlinear control systems for which systematic control methods are available once a flat-output is explicitly known.

The flatness-based control methods may be expected to play a very significant role in high technology applications in the next few years, similar to what happened for nonlinear control in the last decade [START_REF] Sira-Ramírez | Differentially flat systems[END_REF]. The main property of flat systems is that all the state and input variables can be expressed directly, without integration of any differential equation, in terms of the set of so-called "flat output" and a number of its time derivatives. More precisely, the entire system behavior is determined by the trajectory of a finite collection of quantities: flat outputs. This leads to a simple and elegant trajectory design. For a given system, the number of flat outputs is equal to the number of the system inputs.

For more details on flatness see also [START_REF] Fliess | Flatness and defect of non-linear systems: Introductory theory and examples[END_REF], [START_REF] Sira-Ramírez | Differentially flat systems[END_REF], [START_REF] Manish | Higher order sliding mode control of differentially flat systems[END_REF], [START_REF] Isidori | A sufficient condition for full linearization via dynamic state feedback[END_REF], [START_REF] Rouchon | Flatness based control of oscillators[END_REF] and [START_REF] Ph | Flat systems[END_REF].

IV. SLIDING MODE FLATNESS-BASED CONTROL FOR THE DTR PROBLEM

Here, we introduce a first order sliding mode flatnessbased control using the METANET model for dynamic traffic routing. In the freeway portion, there are two alternate routes with two sections in each route, see Fig. 2.

A. DTR Control for Two Alternate Routes With Two Sections

The equations of the system discretized in space are

ρ11 (t) = 1 L 11 [αq e -ρ 11 (t)v 11 (t)] (8)
vji (t) = 1 τ (V ρji (t) -v ji (t)) + 1 Lji v ji (t)(v ji-1 (t) -v ji (t)) -ν τ Lji ρji+1(t)-ρji(t) ρji(t)+κ (9)
ρ21 (t) = 1 L 21 [(1 -α)q e -ρ 21 (t)v 21 (t)] (10)
ρj2 (t) = 1 L j2 [ρ j1 (t)v j1 (t) -ρ j2 (t)v j2 (t)] (11)
V ρji (t) = vf ji exp - 1 a ρ ji (t) ρ cji a
The travel time function is obtained as follows:

y 1 (t) = L 11 v 11 (t) + L 12 v 12 (t) (12)
y 2 (t) = L 21 v 21 (t) + L 22 v 22 (t) (13)
The variable y is equal to the difference in the travel times the two sections.

y(t) = L 11 v 11 (t) + L 12 v 12 (t) - L 21 v 21 (t) + L 22 v 22 (t) (14)
The first time derivative leads to:

ẏ(t) = A + B + C + D + αE (
C1= -L 12 (v 12 ) 2     -σ 13 + v 12 - vf 12 σ 11 τ - vf 12 σ 15 (ρ 11 v 11 -ρ 12 v 12 ) L 12 ρc 12 σ 11 τ     C2= 1 v 12 v 11 - vf 11 σ 19 τ - v 12 - vf 12 σ 11 τ +σ13- ν(ρ 11 -ρ 12 ) L 11 τ (κ+ρ 11 ) C3= σ 17 (v 12 ) 2 v 12 - vf 12 σ 11 τ -σ13
D is composed of three terms (D 1 , D 2 and D 3 );

D1= L 22 (v 22 ) 2     -σ 14 + v 22 - vf 22 σ 12 τ - vf 22 σ 16 (ρ 21 v 21 -ρ 22 v 22 ) L 22 ρc 22 σ 12 τ     D2= -1 v 22 v 21 - vf 21 σ 20 τ - v 22 - vf 22 σ 12 τ +σ14- ν(ρ 21 -ρ 22 ) L 21 τ (κ+ρ 21 ) D3=- σ 18 (v 22 ) 2 v 22 - vf 22 σ 12 τ -σ14
E is composed of two terms. The first one concerns section 11 and the second belongs to section 21;

E11= L 11 qe (v 11 ) 2   vf 11 σ 7 L 11 ρc 11 σ 3 τ - ν L 11 L 11 τ (κ+ρ 11 ) + νσ 5 (L 11 ) 2 τσ 10   E21= L 21 qe (v 21 ) 2   vf 21 σ 8 L 21 ρc 21 σ 4 τ - ν L 21 L 21 τ (κ+ρ 21 ) + νσ 6 (L 21 ) 2 τσ 9   σ 3 = exp 1 a ( ρ11 ρc 11 ) a , σ 4 = exp 1 a ( ρ21 ρc 21
) a , 1) ,

σ 5 = ρ 11 -ρ 12 , σ 6 = ρ 21 -ρ 22 , σ 7 = ( ρ11 ρc 11 ) (a-
σ 8 = ( ρ21 ρc 21 ) (a-1) , σ 9 = (κ + ρ 21 ) 2 , σ 10 = (κ + ρ 11 ) 2 , σ 11 = exp 1 a ( ρ12 ρc 12 ) a , σ 12 = exp 1 a ( ρ22 ρc 22
) a ,

σ 13 = v12σ17 L12 , σ 14 = v22σ18 L22 , σ 15 = ( ρ12 ρc 12
) (a-1) ,

σ 16 = ( ρ22 ρc 22 ) (a-1) , σ 17 = v 11 -v 12 , σ 18 = v 21 -v 22 , σ 19 = exp( 1 a ( ρ11 ρc 11
) a ) and σ 20 = exp( 1 a ( ρ21 ρc 21

) a ).

The studied system is characterized by one input (control) variable, we have then one flat output F = y, which represents the difference in travel time function between the two routes. Thus the equation ( 15) can be rewritten as:

Ḟ (t) = A + B + C + D + αE (16)
from which

α(t) = Ḟ (t) -(A + B + C + D) E (17)
The expression of the state variable allows to choose a suitable trajectory of the travel time (the flat output). The equation of the control (input) variable allows to add additional constraints to this travel time trajectory. This means that all important properties of the system (see equation ( 15)) are contained in such a differential parametrization.

B. Trajectory Planning:

The equation ( 17), corresponds to an open loop control algorithm. In order to define the trajectory planning, a suitable desired trajectory F * has to be defined. According to the expression of the control variable in equation [START_REF] Fliess | Sur les systèmes non linéaires différentiellements plats[END_REF], this trajectory must have smooth derivatives up to order two. In order to reduce computational effort in real time situation, one can build this reference trajectory for the travel time (flat output) using a polynomial interpolation [START_REF] Sira-Ramírez | Differentially flat systems[END_REF] from the initial and final conditions of the travel time

(F (t i ) = F i , Ḟ (t i ) = 0) and (F (t f ) = F f , Ḟ (t f )) = 0).
This is accomplished by prescribing the following desired trajectory for the flat output F :

F * (t) =    F ti for t < t i F ti + (F tf -F ti )ϕ(t, t i , t f ) for t i ≤ t ≤ t f F tf for t > t f (18
) where ϕ(t, t i , t f ) is a polynomial function of time which exhibit a sufficient number of zero derivatives at times t i and t f . (For the polynomial calculation see e.g. [28] [23]):

ϕ(t) =      0 if t < t i 3 t-ti tf -ti 2 -2 t-ti tf -ti 3 if t i ≤ t ≤ t f 1 if t > t f (19
) Thus, by replacing the term Ḟ (t) in equation ( 17) by the term of desired trajectory Ḟ * (t), we obtain the nominal open loop control:

α * (t) = Ḟ * (t) -(A + B + C + D) E (20)

C. Trajectory Tracking:

Due to the parameter variations and disturbances in traffic flow and frequent changes in the traffic conditions, the open loop control is not sufficient. To ensure a steady state and reduce the influence of parameter variations, the traffic flow has to be operated in closed loop. Using the flatness-based open loop presented in equation ( 17), an additional feedback can be determined in order to achieve a desired dynamic behaviour and to compensate the external disturbances. Here, a first order sliding mode controller is used to feedback the system which asymptotically regulates the output towards the desired equilibrium position (see Fig. 3). Thus, the closed-loop control schema including the flatness-based loop can be obtained as follows:

ψ(t) = Ḟ * (t) -K d sign(s(t)) -K p s(t) (21)
where k d and k p are parameters that must be selected so as to satisfy the desired performances of the closed loop system and to ensure asymptotically stabilization of the input variable. The control law α now reads:

V. NUMERICAL SIMULATIONS For the numerical simulations, consider the freeway section depicted in Fig. 2. It shows two routes with the same length and geometric conditions. Each route is partitioned into 2 identical cells (each cell = 400m). At the end of each cell, one loop detector (sensor) is installed. The model parameters are depicted in table I. The control algorithm value is shown via variable message signs (VMS). For the sake of simplicity, we assume a full compliance of the drivers. The used data is collected between 6 AM -10 PM. The simulation time step is about 20 s. The simulations have been done in two different cases; when the two routes have the same geometric conditions as well as when the two routes have different geometric conditions (L11, L12 = 300 m). The traffic demand used in the simulations is shown in Fig. 4. Fig. 5, 6, 7 and 8 show the densities, mean speeds, control values and the travel times evolution in the case when the two routes have the same geometric conditions. From these

Fig. 2 .

 2 Fig. 2. Simulated freeway section

15 ) 6 L 5 L 11 L

 156511 The terms A, B, C, and D represent derivatives of sections 21, 11, 12 and 22 respectively. The term E consists of all the terms in which α does exist. It comes from sections 11 and 21. A is composed of three terms (A 1 , A 2 and A 3 ); 21 τ (κ+ρ 21 ) τ -ν(ρ 21 v 21 -ρ 22 v 22 ) L 22 L 21 τ (κ+ρ 21 ) of two terms (B 1 and B 2 ); 11 τ (κ+ρ 11 ) τ -ν(ρ 11 v 11 -ρ 12 v 12 ) L 12 L 11 τ (κ+ρ 11 ) 11 τ (κ+ρ 11 ) + νσ 5 (L 11 ) 2 τσ 10   C is composed of three terms (C 1 , C 2 and C 3 );

Fig. 3 .

 3 Fig. 3. Closed loop control structure

α

  (t) = Ḟ * (t) -Kdsign(s(t)) -Kps(t) E -A + B + C + D E

figures, it is clear that the densities and mean speeds of the system in each reciprocal cell are equal, it means that the control algorithm distribute the traffic demand between two routes somehow one route not be overloaded and the other underused. Fig. 9, 10, 11 and 12 show the densities, mean speeds, control values and the travel times evolution in the case when the two routes have some different in theirs geometric conditions (L11, L12 = 300 m, L21, L22 = 400 m). We can realize the control algorithm distribute always the traffic demand between the two routes equally. Contrary, because of different length of the two routes, the mean speed of the two routes are different. In this paper, a first order sliding mode flatness-based control algorithm was designed in order to control the incoming flow into two alternate routes which have the same origin and destination. The objective of the control algorithm is to minimize the difference between the travel times of two alternate routes, accordingly, to optimize the flow on the freeway portion. The control algorithm was applied using the METANET model which is a nonlinear traffic flow model. The results show the relevance of the control algorithm in the cases of the same and different geometric conditions of the two alternate routes. These results inspire us to extend this work to deal with more complex network and to design an integrated control in the future.


Paul Amouroux
email: pamourou@yahoo.fr

Patricia Estrada

Christian Gonzalez

Patrich Cerpa

Francisco Urra

Jorge Pérez-Schultheiss

Arthropods of Bahía Exploradores (Chile: Aysén) as bioindicators of human activities: Identifying diversity and abundance

published or not. The documents may come

Introduction

Arthropods are commonly used as bioindicators for evaluating the quality of habitats due to their diversity and/or abundance rapidly evolve in response to environment changes. Bahía Exploradores (Aysén Region, northwest Patagonia, latitude 46 °S) is a remote and nearly pristine valley of glaciated southern Chile. The arthropods of this area can illustrate the human colonization within the valley, but they are little known. The objective of this study was to characterize the diversity and abundance of arthropods in two areas with contrasted levels of human impact in order to identify the most adequate taxon.

The first site is the human-environment observatory (OHM) Patagonia where the human presence dates back to more than 50 years. The second site is the locality 'Teresa' where the human settlement is contemporary.

Materials and Methods

When One week during January 2019, 2020 and 2021

Where Two contrasted human activities sites: OHM Patagonia and Teresa (Fig. 1).

Results and discussion

Conclusion

This survey provided information on the diversity of this Chilean southern region and specimens were deposited in public national collections (MNHN, IEUMCE). Coleoptera could be one of the most useful terrestrial taxon for future investigation on human impacts on the environment. Nevertheless, taxa with aquatic larvae could also be useful to investigate river/water quality. However, sampling effort and methods need to be standardized to future comparison. Moreover, surveys must to be done in other seasons (months) to increase the species detection. Who CRUSTACEA: Isopoda (Fig. 4) and Amphipoda ARACHNIDA: Opiliones (Fig. 5) INSECTA: Orthoptera, Hemiptera, Coleoptera, Lepidoptera, Hymenoptera and Diptera.

A total of 981 individuals were collected in the OHM Patagonia and 429 in Teresa.

The most abundant insect orders were Coleoptera (36%), Diptera (15%) and Hymenoptera (12%) (Fig. 6). In general, Teresa appeared less abundant and diverse than OHM Patagonia (Table 1). The order Coleoptera was the most diverse and abundant taxon. However, the results is highly related to the sampling effort, the climatic conditions, and the techniques used. @Pérez-Schultheiss @Pérez-Schultheiss @Amouroux @Amouroux

Fig. 1 :

 1 Fig. 1: Localization of the two sampled sites. Modified from Google Earth.

Fig. 6 :

 6 Fig. 6: Treemap of main order abundances.

  Fig. 2: Malaise trap Fig. 3: Pitfall trap

Acknowledgments

This work was cofunded by the Labex DRIIHM, French programme "Investissements d'Avenir" (ANR-11-LABX-0010), and by the ANID -FONDECYT postdoctorado #3180344


S Galland

G Udry

H Montagnier

I Beust

M Desort

A.-M Lagrange

F Galland

S Udry

G Montagnier

H Beust

I Boisse

X Bonfils

F Bouchy

X Delfosse

A Eggenberger

D Ehrenreich

T Forveille

G Hébrard

B Loeillet

C Lovis

M Mayor

N Meunier

C Moutou

F Pepe

C Perrier

F Pont

D Queloz

N C Santos

D Ségransan

A Vidal-Madjar

Extrasolar planets and brown dwarfs around A-F type stars , VII. θ Cygni radial velocity variations: planets or stellar phenomenon?

Keywords: . VII. theta Cygni radial velocity variations, planets or stellar phenomenon? radial velocities -stars, early-type -stars, planetary systems -stars, individual, θ Cygni

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

Introduction

Radial-velocity (RV) surveys have led to the detection of more than 300 planets during the past decade 1 . These surveys mainly focus on solar and later-type main-sequence (hereafter MS) stars ( > ∼ F7) that exhibit numerous lines with low rotational broadening, making them ideal targets for classical velocimetry. However, it is crucial to understand how planetary systems form over a wide variety of parent stars and to know, in particular, if there is a correlation between (i) the planet masses and the parent star masses as predicted, for instance, by Kennedy & Kenyon (2008), in order to constrain formation models such as those from Ida & Lin (2005) (but also [START_REF] Boss | [END_REF]and Laughlin et al. 2004 for M dwarfs); and/or between (ii) the planet existence and the parent star masses. For massive stars, surveys of subgiant/giant stars have started to provide first information on planets at orbital distances typically greater than 0.7 AU (e.g., Johnson et al. 2006Johnson et al. , 2007;;Hatzes et al. 2005;Niedzielski et al. 2007;Lovis & Mayor 2007;Sato et al. 2008). Closer separations have to be investigated by observing massive main-sequence stars. In this framework, we have developed a tool dedicated to searching for planets around early (A-F) type stars. The method allowing measurement of the RV of rapid rotators is described by Galland et al. (2005a, hereafter Paper I).

In 2005 we started two surveys dedicated to searching for extrasolar planets and brown dwarfs around a volume-limited sample of A-F main-sequence stars i) with the ELODIE fiberfed echelle spectrograph (Baranne et al. 1996) mounted on the 1.93-m telescope at the Observatoire de Haute-Provence (OHP, France) in the northern hemisphere; and ii) with the HARPS spectrograph (Pepe et al. 2002) installed on the 3.6-m ESO telescope at La Silla Observatory (Chile) in the southern hemisphere. In 2006, the ELODIE spectrograph was replaced by SOPHIE [START_REF] Bouchy | SOPHIE: the Successor of the Spectrograph ELODIE for Extrasolar Planet Search and Characterization[END_REF]. We detected a planet around an F6V star with ELODIE (Galland et al. 2005b, Paper II) and a brown dwarf around an A9V star (Galland et al. 2006b, Paper IV), and with HARPS a two-planet system around an F6IV-V star (Desort et al. 2008, Paper V). We also derived the first statistics of planet existence around A-F stars thanks to our HARPS survey (Lagrange et al. 2009, Paper VI).

We present and analyse in this paper the RV variations of θ Cygni. Section 2 provides the stellar properties and the various data obtained on this object. In Sect. 3, we discuss the origin of the observed RV variations.

Stellar characteristics and measurements

The star θ Cygni (HD 185395, HIP 96441, HR 7469) is a M 1 = 1.38 ± 0.05 M star, with an age estimated at 1.5 +0.6 -0.7 Gyr (Nordström et al. 2004), and located at 18.33 ± 0.05 pc from the Sun (ESA 1997;van Leeuwen 2007). We took its rotational velocity v sin i, effective temperature T eff , and surface gravity log g from Erspamer & North (2003) and Gray et al. (2003) (values in Table 1). We assumed a spectral type F4V, commonly attributed to this star as, e.g., in the Bright Star Catalogue [START_REF] Hoffleit | Bright Star Catalogue[END_REF] or in the HIPPARCOS catalogue (ESA 1997). The wavelength range is 3850-6800 Å for ELODIE and 3872-6943 Å for SOPHIE. Typical exposure times were 15 and 3 min for ELODIE and SOPHIE, respectively, leading to a S /N of ∼200. The exposures were performed with simultaneousthorium spectra to follow and correct for the possible drift of the instrument due to local temperature/pressure variations (whose impact shows a standard deviation of 2.5 m s -1 ). With SOPHIE we used the high-resolution (R ≈ 75 000) mode. [mas yr -1 ] -262.99 † The proper motion is affected from the orbital motion that we discuss in Sect. 2.2. Note: photometric and astrometric data from the HIPPARCOS catalogue (ESA 1997;van Leeuwen 2007); spectroscopic data from Nordström et al. (2004) and Erspamer & North (2003).

Radial velocity variations

The radial velocities (Fig. 1) were measured using a dedicated tool (SAFIR) described in Paper I and based on the Fourier interspectrum method developed in Chelli (2000). The uncertainty associated with ELODIE data is 9 m s -1 on average, consistent with the value obtained from our simulations (see Paper I). In the case of SOPHIE data, the uncertainty is 5 m s -1 on average (taking the photon noise and instrument stability into account).

As θ Cygni has a relatively low projected rotational-velocity (v sin i = 7 km s -1 ), we could also measure the RV using a Gaussian adjustement to the cross-correlation function (CCF). The results obtained by the two different methods are found to be consistent.

ELODIE and SOPHIE RV data show (Fig. 1) a quasi-periodic signal with peak-to-peak amplitude of about 220 m s -1 , much larger than the uncertainties. A drift in the RV curve moreover seems to be present over the whole data set, which could be attributed to a stellar companion (Sect. 2.2). We finally note that the amplitude of the RV variations could also be slightly variable.

We used the CLEAN algorithm [START_REF] Högbom | [END_REF], applied to Lomb-Scargle periodograms to derive the periodograms of the radial velocities measured with both instruments (Fig. 2). This algorithm removes the aliases associated with temporal sampling of the data, and it deconvolves the window function iteratively from the initial "dirty" periodogram to produce the resulting cleaned periodogram.

In the case of ELODIE data, the periodogram shows one peak at a period of 128 ± 5 days. The uncertainty is evaluated with the full width at half maximum of the highest peak. In the case of SOPHIE data, the peak corresponds to a period of 158 ± 10 days.

Line profile variations

SOPHIE Lomb-Scargle periodograms of the BVSs and curvatures (defined as in Hatzes 1996) are presented in Fig. 3, together with false-alarm probabilities (FAP, Kürster et al. 1997). In the case of the SOPHIE data, a peak is seen at approximately 140 days, i.e., not very different from the one measured in the Figure 4 presents the bisector velocity spans (BVSs) as a function of RVs for ELODIE and SOPHIE spectra. The amplitude of the BVS variations is quite small: 50 m s -1 , much smaller than the amplitude of the RV variations. No clear correlation is seem between the BVS and RV variations. It shows that the spectra are mainly shifted in radial velocity without significant changes in the line shape.

Stellar jitter

Finally, to quantify its short-term variations, we monitored θ Cygni for 1 h consecutively (as we did for HD 60532, Paper V). High-frequency variations are due to stellar phenomena and produce a noise (jitter) that has to be taken into account in the analysis of longer period variations. In October 2008, 46 consecutive spectra were then taken with SOPHIE under average observation conditions (airmass below 1.2, S /N 160). The resulting RVs and associated bisector velocity spans are presented in Fig. 5. It appears that the short-term variations can account for an RV amplitude of 30 m s -1 (σ rv = 6.4 m s -1 ) and that the total bisector velocity span amplitude over the whole SOPHIE data (Fig. 4, bottom) set can only be explained with those short-term variations. In the following attempt to analyse the high-amplitude RV variations, we adopt increased RV uncertainties (at least ±6.4 m s -1 , fixed to that value) to take this stellar jitter into account.

Imaging data

We observed θ Cygni at high angular resolution and high contrast with the adaptive-optics (AO) instrument PUEO (Rigaut et al. 1998) mounted on the 3.6-m Canada-France-Hawaii Telescope (CFHT, USA). The near-infrared camera used, KIR (Doyon et al. 1998), has a field of view of ∼35 × 35 , with a scale of ∼0.035 per pixel. We performed non saturated exposures, as well as 30-s saturated ones, to investigate the vicinity of the star at a deeper level. Special care was taken to ensure that the non saturated exposures could be used as references for accurately measuring the possible companions positions relative to the central star and as references for measuring the photometric contrast between the star and the possible companions.

The star θ Cygni was observed in June 2004, September 2005, and November 2007. The log of observations is given in Table 2. A classical reduction was performed using the software ECLIPSE [START_REF] Devillard | The Messenger[END_REF]. A candidate companion (CC) is seen in non saturated images (see Fig. 6) with the narrow bandwidth Fe ii filter (λ 0 = 1.644 μm, Δλ = 0.015 μm).

A deconvolution algorithm using the method described in Véran (1997) was applied to derive the contrast and angular separation (ρ) between the star and the companion.

Figure 7 shows the relative positions of the CC between 2004 and 2007. Clearly, the CC is not a background star, but it is bound to θ Cygni. Given the HIPPARCOS distance, 18.6 pc, we derive a projected separation of 46.5 AU between the two objects, thus a minimum period of roughly 230 years, assuming a circular orbit. Moreover, we see in Fig. 7 that the orbit of the companion is -still very partially -resolved over a three-year period of observation.

The measured contrast between θ Cygni and its companion θ Cygni B is 4.6 ± 0.1 mag in H band; in the K band, the measured contrast is 4.5 ± 0.1 mag. Given the star's apparent magnitudes, provided by Skrutskie et al. (2006), and distance (see above), we deduce H and K absolute magnitudes of 7.0 ± 0.1 and 6.7 ± 0.1 mag, respectively. Using the BCA98 evolutionary models (Baraffe et al. 1998), and assuming any age above ∼100 Myr, we deduce a mass m 2 0.35 M for the companion. The evolutionary effects are negligible. Using the empirical relation given by Delfosse et al. (2000) for M H 7, we find a comparable mass m 2 0.33 M for the companion.

The source θ Cygni was classified as a double star [START_REF] Dommanget | [END_REF]. The Washington double star (WDS) data (Hartkopf & Mason 2001) indicate a visual companion detected several times since 1889, with a magnitude of ∼12, i.e., comparable to the visual magnitude expected from a 0.35-M star at θ Cygni distance. The relative position of this object varied between 1892 and 1968 (see Table 2), and indicates that this companion is bound to θ Cygni. The companion was not detected by HIPPARCOS because the contrast with θ Cygni was too high. It is reasonable to say that the WDS companion and the one found with PUEO is the same. Then we can track its motion over more than a century (Fig. 8), but its orbit is still very incomplete.

Origins of the observed RV variations

We investigate hereafter different possible origins for the RV variations: stellar phenomenon (spots, pulsations) and planets. Beforehand, we estimate the possible impact of the stellar companion on the spectroscopic data.

Impact of the stellar companion on the RV data

The star θ Cygni B could a priori play a role on the measured radial velocities in two ways. Firstly, given their separation ( 2 ) and given usual and the entrance width of the optical fiber, the spectra of the two stars are superimposed in the ELODIE or SOPHIE data (fiber entrance of 2 and 3 , respectively). This could introduce a bias in the measurements. However, the contrast of 4.6 mag in H band translates into a contrast of 7.9 mag (a flux ratio of ∼1500) in the V band, hence the negligible signal of the secondary in our spectra. We expect a radial-velocity effect below 1 m s -1 . If very active, it could still produce a weak Hα signature superimposed on the θ Cygni Hα line. This does not affect our results, as this line is not considered in the RV measurements.

With the classical cross-correlation technique, the potential pollution of the spectrum by a stellar companion can be tested using various masks. We therefore checked that the RV amplitudes remain identical when using various masks. Also, the RVs are identical if we use either the red part or the blue part of the spectra to measure them. This confirms that the spectrum of the companion has no impact on the measured radial velocities.

Secondly, the stellar companion of course induces radial velocity variations in the primary. Given the companion properties and assuming the system seen edge-on, we estimate the maximum drift possibly induced on the primary star by plotting the acceleration projected on the line of sight z, with respect to the true separation r (Fig. 9)

d 2 z dt 2 = Gm 2 r 2 cos arcsin ρ r ,
where G is the gravitational constant, m 2 = 0.35 M the mass of the secondary, and ρ = 46 AU the closest projected separation measured. We find a maximum drift of 12 m s -1 yr -1 , which would lead to a maximum drift of 60 m s -1 over 5 years. In fact, a similar drift has to be included in our fit of the radial-velocity curve (see Sect. 3.3). The star θ Cygni B could then explain such a drift, but of course the observed periodic radial-velocity variations with an amplitude larger than 150 m s -1 are not explained by the presence of this stellar companion.

Stellar phenomenon

Stellar spots

We saw in the previous section that there is no correlation between the star RV variations and the BVS variations. Given the star's projected rotational-velocity, the instrument resolution, and according to the study presented in Desort et al. (2007), we can definitely conclude that the observed RV variations are not due to stellar spots. Indeed, were this the case, a correlation between the BVS and the RV would be observed. Typically, given the star's properties, one or more spots on an inclined star would be needed to reproduce a periodic signal, and a linear correlation coefficient -0.5 between the bisector velocity spans and the RVs. Given the observed RV amplitude, the amplitude of the BVS would therefore be much higher than what is actually observed. Also, in such a case, one would expect significant photometric variations. A single spot producing such an RV variation would induce a photometric amplitude between 5 and 30 mmag, depending mainly on the star inclination and the spot location (Desort et al. 2007). The photometry given by HIPPARCOS (ESA 1997) is constant with a scatter of only 0.004 mag. We recognize, however, that the HIPPARCOS data were not recorded simultaneously with the spectroscopic ones, so this photometric argument is certainly weaker than the absence of correlation between RV and BVS variations. No clear emission in the core of the Ca ii lines is observed (see Fig. 10 for the Ca ii K line) which excludes a high level of activity.

We looked for possible long-term stellar variations using classical Hα, β, γ indicators. Very faint Hα variations are detected, but they are not correlated with the RV variations. As the Hα indexes are moreover quite sensitive to pollution by the thorium lamp, we cannot attribute them to a stellar origin. Finally, the stellar rotational period is less than 7 days according to its v sin i and assuming a stellar radius typical for this type of star, a spot or a set of spots cannot explain RV variations with periods of one hundred days or more. We can therefore safely conclude that spots are most probably not responsible for the observed RV variations.

Granulation

Another source of variation could be the attenuation or suppression of the convective blueshift due to the presence of plages (e.g., Deming & Plymate 1994;Marquez et al. 1996). The expected convective blueshift expected for θ Cygni would be in the range 400-800 m s -1 depending on the line (see Gray 2009 for a F5IV-V, which is the closest to θ Cygni in his sample). The observed variations could then be due to a cyclic variation with a plage filling factor varying between 0 and 30%, typically in order to produce the 220 m s -1 variation. However, such a variation of the convection properties should lead to a strong variation in the bisectors and in the Ca ii index, which is not the case.

Stellar pulsations

Generally, pulsations induce line-profile variations that strongly affect the bisector velocity span (Paper VI). Besides, the timescale of the observed radial-velocity variations (≥100 days) is far larger than the ones of pulsations known for this type of main-sequence stars. It is in fact more characteristic of giant stars variability. Variability periods of a few days are observed in the case of the pulsating γ Doradus stars (Mathias et al. 2004). Moreover, if we integrate the RVs measured between a minimum and a maximum of the amplitude for half a period, we end up with a total displacement close to the stellar radius, which, if even possible, would lead to detectable photometric variations.

It is therefore unlikely that classical pulsations are responsible for the observed RV variations. However, the presence of a peak at about 150 days in the BVS periodogram and in the RV periodogram indicates that the period of the RV variations is linked in some way to the low amplitude line shape variations. Such a situation has never been reported to our knowledge and is indeed quite puzzling. We cannot at this stage exclude that we could be facing a new type of stellar variability, undetected so far because of a lack of long-term, very-precise RV monitoring of main-sequence stars.

Planet(s) around θ Cygni

The θ Cygni RV curve shows quasi-periodic variations with a period of ∼130-150 days, together with a variable positive trend. Moreover, the RV curve seems to be modulated in amplitude.

We now try to see whether these variations can be attributed to a planet or a planetary system. We first tried to fit the whole data set (with uncertainties set to ±6.4 m s -1 ) assuming a single planet and allowing a drift. We failed to correctly reproduce the observed RVs. Figure 11 shows an example of a fit with a planet on a circular orbit (plus a drift). The planet that produces such an RV curve has a mass of 2.3 M Jup and is on a 0.6 AU circular orbit (Table 3). The algorithm ends up with a period of 155 days hung on the SOPHIE data set, but fails to properly fit the ELODIE data set, as if there was period/phase change with time (clearly visible as structure in the residuals near the range [3600-3800] days). The additional drift needed is 16 m s -1 yr -1 , slightly higher than the maximum value that the binary companion would probably produce. We note that a system that fails to fit satisfactorilly the whole set of data (taken over more than 5 years) would allow the data to be fitted if they were limited to one or two consecutive periods.

We then tried to simultaneously fit the whole set of data, assuming the presence of several planets and using a genetic algorithm search. No stable solution was found. The only satisfactory fits are achieved by unstable systems with orbits that cross each other. We show an example in Fig. 12 of a fit obtained with a two-companion plus drift model. The residuals are still very high, the fit is not improved compared to the one-companion plus drift model (residual rms 35 m s -1 versus 39 m s -1 , and the drift that we get is approximately the same: 17 m s -1 yr -1 versus 16 m s -1 yr -1 ). Moreover, these kinds of configurations with massive planets on such close orbits are not dynamically stable.

The system is then obviously more complex than just consisting of one, two, or even three planets plus a drift. We then explored more exotic configurations:

instead of harbouring one single planet, the system could consist of a binary planet system orbiting the star, much like the Pluto-Charo system, but with higher masses. However, the radial-velocity signal generated by this configuration would be very close to the one generated by a single planet [m s -1 ] 38.7 reduced χ 2 6.1 having the total mass. The only departure from the pure Keplerian system would be related to the secular perturbations of the orbit due to the binary nature of the system. The associated secular period would be at least several hundred primary orbital periods, i.e., much longer than our observation time span; hence, we should mainly detect the primary orbital signal with good accuracy; one could also think of 2 co-orbiting planets locked in 1:1 mean-motion resonance. Such a configuration has been observed in the system of Saturn satellites. The two satellites orbit Saturn on the same orbit. In the rotating frame, they have synchronized horseshoe-like libration motions that prevent them from colliding. The less massive one has the largest amplitude motion (Yoder et al. 1983). The associated libration period is here again a few hundred orbital periods of the primaries. Over a shorter time span, the radial-velocity signal of the whole system will mainly consist of the sum of the individual signals of the two planets. If we assume that the two planets have zero eccentricities, the radial velocity signal will have the form A cos(nt) + B cos(nt + φ), where n is the common mean motion of the two planets, A and B are amplitudes related to the masses of the planets, and φ is a phase shift that depends on their current mutual configuration. This can be rewritten as C cos(nt + ψ) where C and ψ are new amplitudes and phases that depend on A, B, and φ. This is equivalent to the signal generated by one single planet. This remains true if the planets have small eccentricities, and if they have larger eccentricities, the mutual system is not stable. Here again, the only departure from this signal will be the temporal variation of C and ψ that is related to the mutual libration motion of the two bodies. Hence we should not expect to see changes before several decades. The signal detected over our observation timespan should not be different from that of a single planet.

a 1 sin i [10 -3 AU] 0.99 f (m) [ 1 0 -9 M ] 5 . 5 M 1 [M ] 1 . 3 8 m 2 sin i [M Jup ] 2 . 2 9 a [AU] 0.63 drift [m s -1 yr -1 ] 1 6±

Concluding remarks

The radial velocities obtained with ELODIE and SOPHIE on θ Cygni are quasi-periodically variable, with a 150 day period.

We investigated different possible origins to these RV variations, either star or planet related.

Given knowledge today of stellar activity, using several criteria as usual in this type of study, we are unable to attribute these variations to the star itself (spots, pulsations); however, the peak in the periodogram of the BVSs at 140 days, i.e., close to the period of the RV variations prevents us from totally excluding a stellar origin.

We tried to fit the data with one planet orbiting at less than 1 AU and with a mass of a few Jupiter masses, also taking the impact of the observed companion star into account. It appears that such an hypothesis allows us to fit only part of the data, recorded over a limited number of periods, but fails to satisfactorily fit the whole set of data (taken over more than 5 years). More complex systems were investigated, but no convincing result was obtained. More observations and detailed studies of gravitational interaction between the 2 planets are needed to understand this sytem.

With the data available, we are then not able to conclude on the origin of these puzzling RV variations. If a planet origin is confirmed, then θ Cygni, with a spectral type of F4V, would be the earliest main-sequence star hosting planets found so far. Moreover, it would be one of the few low-metallicity stars hosting planets. Its planetary system would not be simple, and it would in particular include strongly interacting planets.

If a stellar origin were confirmed, then this would show that, unexpectedly, some main-sequence stars, not classified as active by the usual criteria (e.g., Ca ii H&K indexes) or pulsating, may undergo intrinsic variations that produce quasi periodical, largeamplitude, and long-period (more than 100 days) RV variations, with at the same time, low levels of line shape deformations (hence small-amplitude BVS variations). Such situations have not been considered so far in the analysis of RV variations, and would need to be considered in future searches for such longperiod planets, both in terms of observational strategy and data analysis.

We can at least conclude that θ Cygni is an individual complex system that deserves many more observations to be understood, and it may also serve as an example for other searches.

  Description of the data Between 2003 and 2006, we recorded 91 high S /N spectra of θ Cygni with ELODIE and, between November 2006 and December 2008, we recorded 162 spectra with SOPHIE.

Fig. 1 .

 1 Fig. 1. Radial velocities of θ Cygni obtained with ELODIE (left) and SOPHIE (right).

Fig. 2 .

 2 Fig. 2. ELODIE (top) and SOPHIE (bottom) CLEAN periodograms of the radial velocities. On each data set, one peak is observed, either at 128 ± 5 days or at 158 ± 10 days, respectively.

Fig. 4 .

 4 Fig. 4. Bisector velocity span versus RVs for ELODIE (top) and SOPHIE (bottom) data, showing that there is no correlation between line profiles and radial velocities. The dashed lines on SOPHIE data show the minimum effect of short-term variations, see Sect. 2.1.4.

Fig. 5 .

 5 Fig. 5. Top: RVs for 1-h SOPHIE data. Bottom: bisector velocity spans versus RVs for 1-h SOPHIE data.

Fig. 6 .

 6 Fig. 6. High contrast and angular resolution image acquired with the adaptive optics system PUEO installed at CFHT: a companion is clearly visible.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Temporal evolution of the separation of the stellar companion of θ Cygni with respect to the central star. Filled symbols represent the position of the stellar companion relative to the main star for the three AO observations (filled circle: June 2004; filled square: September 2005; filled triangle: November 2007), whereas empty symbols show the positions at the same dates in case of a background star. The curve shows the path that would have followed the candidate companion if it was a background star, taking the star proper-and parallactic-motion into account.

Fig. 9 .

 9 Fig. 9. Possible drift induced by the stellar companion given the unprojected true separation r to the primary star.

Fig. 10 .

 10 Fig. 10. No emission is observed in the Ca ii K line for the θ Cygni spectra. This is the average spectrum of all the SOPHIE spectra used, after recentring by the RV variations measured.

Fig. 11 .

 11 Fig. 11. ELODIE and SOPHIE radial velocities and orbital solution for θ Cygni, with one planet (2.3 M Jup at ∼0.6 AU, circular orbit) and a drift. The residuals to the fitted orbital solution are displayed below.

Fig. 12 .

 12 Fig.12. ELODIE and SOPHIE radial velocities and orbital solution for θ Cygni, with two planets (2.1 M Jup at ∼0.6 AU and 0.8 M Jup at ∼0.5 AU on circular orbits) and a drift. The residuals to the fitted orbital solution are displayed below.

Table 1 .

 1 θ Cygni stellar properties.

	Parameter		θ Cygni
	Spectral type		F4V
	v sin i	[km s -1 ]	7
	V		4.49
	B -V		0.395 ± 0.015
	π	[mas]	54.54 ± 0.15
	Distance	[pc]	18.33 ± 0.05
	M V		3.14
	[Fe/H]		-0.08
	T eff	[K]	6745
	log g		4.2
	M 1 Age ppm (α) †	[M ] [Gyr] [mas yr -1 ]	1 .38 ± 0.05 1.5 +0.6 -0.7 -8.15
	ppm (δ) †

Table 2 .

 2 Log of WDS observations between 1889 and 1968 followed by the AO observations using PUEO at CFHT between June 2004 and November 2007.

	Date	ρ	θ	Contrast Band
		[ ]	[ • ]	[ m a g ]
	1889.37	3.62	43.9
	1892.38	3.79	47.0
	1898.46	3.37	49.2
	1898.63	3.71	46.9
	1958.58	3.42	51.7
	1968.72	2.92	59.9
	2004-06-28 2.510 ± 0.021 67.37 ± 0.48 4.6 ± 0.1	Feii
	2007-11-16 2.369 ± 0.005 69.02 ± 0.11 4.6 ± 0.1	Feii
	In the case of WDS data the uncertainties are unknown.

Table 3 .

 3 ELODIE/SOPHIE best orbital solution for θ Cygni, considering one planet and a drift.

	Parameter		θ Cygni b
	P	[days]	154.5 ± 0.4
	T 0	[JD-2 450 000]	4016 ± 1
	e		0 (fixed)
	ω	[deg]	0
	K	[m s -1 ]	7 0 ± 4
	N meas		253
	σ O-C

A comprehensive list of known exoplanets is available at http://exoplanet.eu.

Acknowledgements. We acknowledge support from the French CNRS and the support from the Agence Nationale de la Recherche (ANR grant NT05-4_44463). We are grateful to the Observatoire de Haute-Provence (OHP) and the CFHT for their help during the observations, and to the Programme National de Planétologie (PNP, INSU). These results have made use of the SIMBAD database, operated at the CDS, Strasbourg, France. They also make use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. We also thank Gérard Zins and Sylvain Cètre for their help in implementing the SAFIR interface. X.B. acknowledges support from the Fundação para a Ciência e a Tecnologia (Portugal) in the form of a fellowship (reference SFRH/BPD/21710/2005) and a programme (reference PTDC/CTE-AST/72685/2006), as well as the Gulbenkian Foundation for funding through the "Programa de Estímulo à Investigação". N.C.S. would like to acknowledge the support from Fundação para a Ciência e a Tecnologia (Portugal) in the form of a grant (references POCI/CTE-AST/56453/2004 and PPCDT/CTE-AST/56453/2004), and through programme Ciência 2007 (C2007-CAUP-FCT/136/2006).

Based on observations made with the ELODIE and SOPHIE spectrographs at the Observatoire de Haute-Provence (CNRS, France) and with the PUEO adaptive optics system at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France,