UPE Vincent Aguiléra

Lvmt

A DEPARTURE TIME CHOICE MODEL FOR DYNAMIC ASSIGNMENT ON INTERURBAN NETWORKS

INTRODUCTION

The increasing interest for time-varying tolling policies in transportation networks has highlighted the need for dynamic traffic assignment models in which users are allowed to choose both their route from origin to destination and their departure time.

In the case of a single route with a bounded capacity, [START_REF] Vickrey | Congestion theory and transport investment[END_REF] derived the departure rate as the outcome of individual cost minimization. In Vickrey's model, for a given departure time, individual costs are modelled as continuous, uniminimal and convex functions of the schedule delay, i.e. as a penalty for early or late arrivals at destination. Since then a number of extensions have been proposed, notably to deal with complex networks, multiple origins and multiple destinations. Most of those extensions (see for instance [START_REF] Heydecker | Analysis of Dynamic Traffic Equilibrium with Departure Time Choice[END_REF] or de Palma and [START_REF] De Palma | Modelling and Evaluation of Road Pricing in Paris[END_REF] are tailored for studying equilibrium in metropolitan networks, in which commuting travel patterns are assumed: the preferred arrival time to a destination is often taken from a normal distribution around a peak hour, and individual schedule costs essentially follow Vickrey's model.

Unfortunately, those assumptions do not hold in the case of large interurban networks. In this case, the economic behaviour of a trip-maker with respect to its preferred arrival time varies widely from one to the other. Some travellers may want to be at destination early in the morning, some others late in the evening, while some others require to be at destination absolutely before a given hour, etc. Moreover, a significant percentage of the travellers may wish to avoid peak days, by rescheduling their departure day to the day before or to the day after, if incited to do so. This article presents extensions made to a dynamic traffic assignment model to encompass those needs. It has the following outline. Section 1 contains a brief review of the literature. It helps in highlighting the distinctive features of our model, which is introduced in section 2. The two main algorithms used to compute a user equilibrium with departure choice are provided in section 3. Some results from numerical studies are presented in section 4, before concluding.

LITERATURE REVIEW

The seminal paper on trip scheduling [START_REF] Vickrey | Congestion theory and transport investment[END_REF] essentially deals with morning commute. A fixed number of commuters is considered, travelling from a single origin to a single destination, through a single route. Congestion occurs at a bottleneck of bounded capacity. Users are modelled as microeconomic agents minimizing a cost function that depends on travel time and schedule delay. In the simplest version of the model, Vickrey considers homogeneous users that all have the same preferred arrival time and the same cost function, linear in travel time and piece wise linear convex in schedule delay. The main results from this model are recalled in subsection 1.1. Since then, several extensions have been proposed. Some of them are the topic of subsection 1.2. However, as shown in subsection 1.3, some empirical evidences from field data indicate that those various extensions may not fulfil the requirements of departure choice modelling in interurban networks.

Main results from the bottleneck model

A set of commuters wish to reach a CBD accessible by one route with bounded capacity. Each user is characterized by a preferred arrival instant * h , and assesses the decision of departing at an instant h using a cost function of the form:

( ) ( ) + - + + = ) ( . ) ( . ) ( . ) ( h d h d h t h G γ β α
(1) where:

• ) (h t is the travel time when departing at h,

• * ) ( ) ( h h t h h d - + =
is the delay incurred at arrival when departing at h,

• α is the value of time of the commuter,

• β [resp. γ ] it the marginal cost of arriving earlier [resp. later] than preferred,

• ( ) -and ( ) + stand for the positive and negative part of the delay.

There are two standard ways of describing the set of users. Either one considers a finite number of categories of commuters, differing by their * h (and also possibly their value of time, value of arriving late or early); or one considers that users have preferred arrival instants distributed among a set of possible values. In the later case, the "S-shape" assumption is made: there is a single interval during which the density of commuters exceeds capacity. This assumption makes the model analytically tractable, and induces a travel time pattern similar to the one with a unique * h shared by all commuters.

A typical equilibrium situation is depicted in Fig. 1 for a bottleneck of capacity K. The set of users is modelled by a cumulative distribution * X over the set of their preferred arrival instants. Their choices of departure instants are given by the cumulative distribution + X . The bottleneck model allows one to compute the cumulative distribution - X of users at the exit of the bottleneck. Fig. 1 can be interpreted as such: the horizontal difference between + X and - X (i.e ) (h t ) gives the amount of time needed to traverse the bottleneck, when entering the bottleneck at instant h. The horizontal difference between - X and * X (i.e. ) (h l ) gives the delay at arrival. The traversal time ) (h t is a piecewise function with only two admissible slopes and a single maximum. It increases at the beginning of the congestion period, when users are arriving earlier than preferred. When ) (h t is decreasing, users are arriving later than preferred. Note that the simple form of the delay cost function (i.e. the two last terms of Eq.1) implies the piece wise linear shape of the travel time function. Under more general assumptions, it would be smoother.

Extensions to the bottleneck model

Assumptions made in Vickrey's model are indeed well tailored for morning commutes: the simple congestion model makes sense for trips toward a small CBD area; a unique preferred arrival time is adapted for work starting time; finally a convex schedule delay cost function does fit with job schedule constraints. Since then, the transportation community has investigated the field in two main directions. Some works, mainly from transport economists, have focused on users heterogeneity. Others have proposed extensions to whole networks.

Departure time choice with users heterogeneity

Heterogeneity in preferred arrival times can be addressed either in discretely fashion, by allowing only a finite set of preferred arrival times, or continuously using a distribution. The finite case has been studied extensively by [START_REF] Lindsey | Existence, uniqueness and trip cost function properties of user equilibrium in the bottleneck model with multiple user classes[END_REF] while the continuous case was first treated by [START_REF] Hendrickson | Schedule delay and departure time choice in a determinist model[END_REF]. Heterogeneity pertaining to the costs of travel time and of schedule delay has been studied, among others, by [START_REF] De Palma | A structural model of peak period congestion: a traffic bottleneck with elastic demand[END_REF][START_REF] Van Der Zijpp | Multiclass continuous-time equilibrium model for departure time choice on single-bottleneck network[END_REF][START_REF] Van Der Zijpp | Multiclass continuous-time equilibrium model for departure time choice on single-bottleneck network[END_REF]. Other extensions include the modelling of stochastic demand and capacity, multiple routes or elastic demand (see [START_REF] De Palma | Recent developments in the bottleneck model[END_REF] for a review). When users are at equilibrium, the bottleneck model predicts a congestion pattern with a single peak in travel time. In the numerous extensions, the resulting congestion patterns are very similar to the homogeneous case. When considering a finite number of preferred arrival instants, there is a limited number of peaks in travel time (at most one per preferred arrival instant) and a spontaneous segregation among users is observed. Commuters with different preferred arrival instants depart at different instants [START_REF] Lindsey | Existence, uniqueness and trip cost function properties of user equilibrium in the bottleneck model with multiple user classes[END_REF]. The case where users' preferences are distributed over an interval has received less attention. Papers in this line mainly considered "S-shape" distribution (See Smith,1984 and[START_REF] Daganzo | The uniqueness of a time-dependent equilibrium distribution of arrivals at a single bottleneck[END_REF]. As exposed in section 1.1 this case is practically equivalent to the one with a single preferred arrival time and produces exactly the same travel time pattern.

Departure time choice on networks

The bottleneck model has been extended to networks, in an attempt to produce operational planning models. The computation of the user equilibrium in such a context is known as the dynamical traffic assignment problem with departure time choice. [START_REF] Friesz | A Variational Inequality Formulation of the Dynamic Network User Equilibrium Problem Operations[END_REF] first proposed a formulation of the user equilibrium with both route and departure time choice. His model considers users spread between origin-destination (o-d) pairs, with a unique preferred arrival instant by o-d. Since then most of the models proposed in the literature rely on Friesz's original paradigm (see for instance Wie et al, 2001). Rather than focusing on users heterogeneity, this part of the literature has made considerable efforts to improve congestion representation by integrating sophisticated traffic models.

Empirical evidences and their practical implications for departure time choice modelling.

Two assumptions underlie the bottleneck model and its various extensions: (i) preferred arrival times are taken from a discrete set of values and (ii) delay cost functions are convex. As discussed in this subsection, those two assumptions appear not to be appropriate means of modelling economic preferences of inter-urban trip makers.

Let us first have a look at the travel time patterns and flow rates observed in interurban trips. Fig. 1(a) shows the variations of travel time between Lyon and Montpellier, two cities of southern France, on a holiday departure day. The pattern is quite far from the single peak predicted by the bottleneck model: at least two peak periods can be observed, along with significant variations elsewhere. The flow rate on the same o-d pair is plotted in Fig. 1(b). One can observe significant flow rates during the whole day. This is clearly not consistent with a single preferred arrival instant. Another interesting point is the diversity of inter urban travellers. As opposed to the morning peak where the road traffic is mainly composed of commuters, inter urban trips have a wide variety of purposes, inducing significant differences in value of time and delay cost functions. In the same order of idea, a significant part of the traffic is composed of heavy vehicles, which has important consequences on congestion modelling.

Finally the results of a survey organized in 2008 by three French motorway operators show that, during summer holydays, trip makers can be divided into two categories: some can absolutely not afford arriving later than scheduled (e.g. because they need catching the key for their rental). Some others are far less constrained at arrival and are even ready to reschedule their departure day, if they can benefit from lower congestion and toll fares. This last point is especially interesting as it shows that in inter urban context the convexity of delay cost functions can no longer be assumed. Indeed in this case the cost of the delay does not necessarily decrease as the arrival time gets closer to the preferred schedule. A traveller considering to leave one day in advance to avoid traffic jams will not necessarily consider arriving at 2 a.m. a better option.

To sum it up, empirical observations show that, for inter urban trips, a departure time choice model should differ from the classical "bottleneck-like" approach, with respect to the three following requirements:

1. A high level of heterogeneity regarding both preferred schedules (several preferred arrival times per o-d pair) and economic characteristics (value of times and schedule delay cost functions) is required. 2. Multi class congestion modelling is to be considered. 3. Users should be able to choice their day of departure as well as the time of the day.

THE DYNAMIC ASSIGNMENT PROBLEM WITH DEPARTURE TIME CHOICE.

In order to meet the modelling challenges exposed just below, we present an extension based on the LADTA 1 model proposed by [START_REF] Leurent | On network assignment and supply demand equilibrium: an analysis framework and a simple dynamic model[END_REF]. The dynamic assignment with departure time choice is modelled in a supply demand framework where the supply is a network of bottlenecks and the demand a set of microeconomic agents. Those agents are characterized by economic preferences (value of time and delay cost function), temporal preferences (preferred arrival time) and physical characteristics (vehicle type). The network is subject to congestion and hence to a given demand the supply model associates a set of arc traversal time functions.

Similarly the demand reacts to the supply by adjusting the time varying flows at the entrance of each route of the network according to the level of congestion. The existence of solutions of such assignment problems has been established recently by [START_REF] Meunier | Existence results for dynamic congestion games[END_REF] in a continuous game theory framework.

In its original version, LADTA deals with multi-class dynamic traffic assignment considering a fixed demand: the demand is expressed as a dynamic o-d matrix that solely depends on departure instants from the origins. In the sequel, this problem is called the user optimum assignment with fixed departure time. It is stated as a fixed point problem in subsection 2.1, along with some helpful notations. Next, subsection 2.2 states the user optimum assignment with departure time choice problem in a similar way.

User optimum assignment with fixed departure time

Let us model a road network as a directed graph

) , ( A N G =
, where N is a finite set of nodes and ( )

N N A × ⊆
is the set of arcs in G . A route r in G is a finite, non empty, sequence of connected arcs. R denotes the set of routes in G . Users of the transportation network are modelled as a continuum of microeconomic agents, where each user belongs to a user class u . U denotes the set of user classes. To each route r in R is associated a route traversal time function

( ) h t u r ,
, where h is a departure instant from the head node of the route, taken in a continuous set of instants H . Also, a route traversal cost function ( )

h c u r , is associated to each route r in R . If ( ) d o d o i ≠ = , ,
is a distinguished pair of nodes in N N × , then R R i ⊆ denotes the set of routes starting at o and ending at d .

( )

h x u i,
is the density of users wishing to depart from o at instant h , and

( ) ( ) ∫ < ′ ′ ′ = h h u i u i h d h x h X , ,
is the cumulated distribution of such users.

Route choice

Given a fixed set

I of o-d pairs, a demand { } u i U I X X , , =
, and a set of route traversal cost functions

{ } u r U R c c , , = , a route choice RC is the definition of a set { } u r i U R I X X , , , ,
= of distributions of users on routes per o-d pair such that, for all ( )

u h i ,
, triples, the two following equations hold:

( ) ( ) h X h X u i R r u r i i , , , = ∑ ∈ (2) ( ) ( ) ( ) { } h c h c h x u r R r u r u r i i , , , , min 0 ∈ = ⇒ ≠ (3)
Eq. ( 2) expresses that, at every departure instant, the demand of users of class u on the o-d pair i is distributed among the routes of this o-d pair. Eq. ( 3) expresses that only routes of minimal cost are chosen. Using a compact notation:

( ) U R U I U R I c X X , , , , , RC = Traffic flowing Knowing a distribution U R I X , ,
of users on routes, a traffic flowing TF is a function that returns a set of routes traversal cost and time functions. Using a compact notation:

( ) ( ) U R I U R U R X t c , , , , TF , =
Details concerning the computation of a traffic flowing are omitted here. The reader interested by the actual implementation of traffic flowing in LADTA can refer to (Aguiléra and [START_REF] Bibliography Aguiléra | On large size problems of dynamic network assignment and traffic equilibrium : computational principles and application to Paris road network[END_REF].

User optimum assignment with fixed departure time Given a demand U I X , , a user optimum assignment with fixed departure time

( ) U I X , UOA is a route choice U R I X , , ~
such that: x , wishes to minimize its total traversal cost, he has to choose a departure instant h that minimizes

( ) ( ) ( ) { } U R I U I U R I U I X X X X , ,
( ) ( ) ( ) ( ) * , , , , * h h t h D h c h r G u r u u r h u - + + =

Departure choice

( ) ( ) U R U R U I U R I t c X X , , * , , , , , DC = (5)
User optimum assignment with departure time choice Given a demand * ,U I X , a user optimum assignment with departure time choice ( )

* , * UOA U I X is a route choice U R I X , , ~
such that:

( ) ( ) { } ) TF( , DC ŨOA , , * , , , * , * U R I U I U R I U I X X X X ∈ = (6)

ALGORITHMS.

For a given o-d pair, an exact computation of the departure choice as defined by Eq.5 involves managing the possibly infinite set of routes on this o-d (by definition, each route is finite, but looping routes are not excluded). However, it appears reasonable in practice to concentrate on a finite subset of efficient routes. The departure choice problem is then casted down to the departure time choice on a single route, for which an algorithm is detailed in subsection 3.1. An algorithm for the combined route choice and departure time choice is then proposed in subsection 3.2.

Departure time choice on a single route

Let r be a route and assume we are given the traversal time r t and cost r c for every departure instant h on that route. Now suppose that a demand on that route is expressed as a distribution * r X of users, as a function of preferred arrival instants * h . All users in * r X share a same delay cost function. This section presents an algorithm, called DTCR, which computes the distribution ( )

h X r
at the entrance of the route for every departure instant h and consistent with the pair ( r t ; r c ). The basic idea of the algorithm is as follows: for each * h of a sample of preferred arrival instants, find the set of departure instants that lead to a total traversal cost that is "close" to the minimum. Then, spread uniformly over this set the volume of users with preferred arrival time close to * h . • G δ : a small positive number.

Algorithm

• n : a positive integer. Output • ) (h X r
: a cumulated distribution of users for every departure instant from the origin of r Begin

Let * k h , n k K 0 = be a discrete sampling of [ ] * * ; b a h h . 0 ← r X For n k K 1 = 1. Let * ,k r X ∆ be the number of users in the interval [ ] * * 1 ; k k h h - ∫ - ← ∆ * * 1 ) ( * * , k k h h r k r dh h x X 2. Let * h be the midpoint of [ ] * * 1 ; k k h h - 2 * 1 * * - + = k k h h h 3. Compute the total cost function G such that : ( ) ( ) ( ) ( ) * h h t h D h c h G r r - + + = 4. Let g be the minimum value of G ( ) { } h G g H h∈ ← min

Compute the set *

H where G is close to its minimum

( ) { } G g h G H h H δ < - ∈ ← , * 6. Compute the function ) ( , h X k r ∆ such that the users * ,k r X ∆ are uniformly distributed on * H . 7. Add k r X , ∆ to r X k r r r X X X , ∆ + ←

End For End

Combined route choice and departure time choice

Algorithm DTCR solves a departure choice problem with no route choice and a fixed supply. We are now going to take advantage of it in order to solve the combined route and departure time choice problem. At each iteration k of the algorithm, we are given a fixed departure time o-d matrix

] 1 [ , - k U I X
, and we first find the user optimal assignment with fixed departure time, thus obtaining for every instants and on every o-d pair an optimal route together with its traversal time and cost ( )

] [ , ] [ , , k u i k u i t c
. Using DTCR it is then possible to compute the optimal departure time decision of the users with respect to ( )

] [ , ] [ , , k u i k u i t c
and to derive a new fixed departure time o-d matrix

] [ , k U I X . Algorithm UOA* Inputs • * ,U I X : a set of cumulated distributions over [ ] * * ; b a h h Output • U R I X , ,
: a user optimum assignment with departure time choice Begin ( )

] 0 [ , ] 0 [ , ] 0 [ , , DTCR U I U I U I t c X ← 1 ← k Do 1.
Compute a user optimum assignment with fixed departure time from the demand

] 1 [ , - k U I X ( ) ( ) ] 1 [ , ] [ , ] [ , ] [ , , UOA , , - ← k U I k U I k U I k U R I X t c X

For each o-d pair i and user class u , compute the auxiliary demand

] [ , k u i Y using ( ) ] [ , ] [ , , k u i k u i t c
as inputs of the DTCR algorithm

( ) ] [ , ] [ , ] [ , , DTCR k U I k U I k U I t c Y ← 3. Increase k and compute k w 1 + ← k k k w k 1 = 4. Compute ] [ , k U I X ( ) ] [ , ] 1 [ , ] [ , 1 k U I k k U I k k U I Y w X w X + - ← - While ] [ , k U I X differs significantly from ] 1 [ , - k U I X ( ) ] [ , , , UOA k U I U R I X X ← End 4. NUMERICAL EXAMPLES.
The two algorithms exposed in section 3 have been implemented within the LTK 2 , an implementation of LADTA that allow for computing dynamic traffic assignments on large networks (Aguiléra and [START_REF] Bibliography Aguiléra | On large size problems of dynamic network assignment and traffic equilibrium : computational principles and application to Paris road network[END_REF]. In this section, the results coming from two test cases are presented. The first one, called SR91, is a small network with one origin, one destination, two routes and two user classes. Albeit simple, it illustrates combined route and departure time choice. It is also a well known case study for transportation economists and as such most of the data regarding both scheduling cost functions and value of time are available in the literature. The second example, called VDR, is an application to a large network.

It has to be clear to the reader that what follows in subsections 4.1 and 4.2 has to be understood as numerical experiments. Their sole purpose is illustrative, and the conclusions, figures and charts presented therein have no particular meaning outside the scope of this paper.

The SR91 example

The State Route 91 is located in the orange county (California, USA) and was faced to an important congestion problem at the beginning of the 90's. The road is connecting a residential zone to a labour pool. Before 1995, it had a capacity of 8000 pcu/h, and four lanes in each direction. In 1995 two lanes were added in each direction. The two additional lanes are equipped with time varying tolls, in both directions. The tolls are used as a congestion management tool to alleviate traffic on peak hours by spreading demand. Yet this management is complex: a variation in the toll fare can induce both trip rescheduling and rerouting. Commuters are highly heterogeneous regarding their arrival time preference [START_REF] Sautter | Accessibilité et Modulation tarifaire[END_REF], so the travel time pattern at equilibrium is likely to have a non trivial form.

User classes

Two user classes, subscripted r and p , are considered. Both of them have a delay cost function D under the classic V-shape form:

( ) ( ) + - + = d d d D . . ) ( γ β
They differ only by their value of time α , their unit cost of arriving late, proportional to β and their unit cost of arriving early, proportional to γ . For user class p ,

Demand and supply

The distributions of desired arrival instants for each user class are also taken from [START_REF] Sautter | Accessibilité et Modulation tarifaire[END_REF] and are based on historical data from the road operator Cofiroute, in charge of the SR91 since 1995. The capacities of two routes are set to 8000 pcu/hour for the free route, and 2500 pcu/hour for the two additional lanes.

Scenarios

Two scenarios have been simulated. The first is called untolled. The amount of the toll fare on the two additional lanes was set to 0, and the user equilibrium with departure time choice has been computed using the LTK. The second scenario is called tolled. The time-varying toll fare was set equal to the curve plotted in Fig. 2

(d).

A comparative study allows us to discuss the net effect of tolling.

Results

The results of each route and each scenario are plotted in Fig. 2(a-c) . In both cases the results are consistent with the bottleneck model. In the untolled scenario, the travel time pattern shows a double peak with slopes corresponding to the one induced by the ratios α γ and α β . The travel time maxima are obtained when delays are close to 0. Users are using indifferently the two routes and the travel time is rigorously the same. In the tolled scenario, the two users types are segregated by the toll ; only users of the category r uses the tolled route. The travel time pattern is affected by the time-varying toll. Let us note that the tolling scheme globally reduces travel times. This is achieved by two mechanisms. On the tolled route the toll tends to spread the traffic, thus reducing congestion. On the untolled route, the reduction is less pronounced and is due to a small reduction in traffic resulting of the departure of the users of category r.

Table 1 exemplifies the results of an hypothetical socio-economic analysis based on the results of the simulation. The benefit of the toll is driven by the time savings made by the class r, which is natural as they have the entire benefit of the faster tolled lanes. Yet the increase in traversal costs of class p is compensated by a decrease in delay costs. This is not straightforward but a closer inspection reveals that as class r uses only very marginally the free route in the tolled scenario, the class p has relatively more capacity than in the untolled scenario. However as there is no toll on their route they tend to schedule their trips regarding to congestion costs, so the benefit of this additional capacity results mainly in a decrease in delay cost.

Untolled

The VDR example

During summer holydays, a significant part of the trans-european road traffic is concentrated in the Vallée du Rhône (VDR) area. Tourists coming from northern Europe (including Belgium, the Netherlands, Germany and Great Britain), travel across France to reach (or return from) southern countries (e.g Italy and Spain), meeting on their way people from the Paris area. The situation is depicted in Fig. 2 The map on the left hand side (Fig. 2 During summer Saturdays, traffic conditions on motorways are usually very bad, especially on the A7, because of high levels of congestion. To better operate the network, motorway operators may be interested in studying time varying tolling strategies. Results presented in the sequel illustrate the ability of our model to handle such kind of studies on large networks.

Input data

Most of the data was provided to us by courtesy of companies of the Vinci Group (ASF, APRR and Cofiroute) operating the motorway network in the area of interest. The network comprises 2404 arcs and 939 nodes. We knew for each arc its capacity, free flow travel time for passenger cars, and travel price for passenger cars. The demand was expressed for 628 o-d pairs. The simulated day was July the 14 th , 2007. Two user classes, named vl and vls, were considered. Those two classes correspond to passenger cars1 . They share most of their characteristics (same free flow travel time, same toll prices, etc). They are distinguished only by their delay cost functions, plotted in Fig. 3(a). For user class vl, the penalty of arriving later than scheduled grows linearly at a very high rate, while the penalty of arriving sooner grows at a lower rate. The delay cost function of user class vls is a little bit more complex. Around 0, it has a classical V-shape, except that the cost of a early arrival grows faster than the cost of a late arrival. Between 6 and 18 the cost of the delay is infinite. Around 24, the shape of the delay cost function is similar than around 0, except that it is shifted up by an amount that correspond to the cost of rescheduling the departure to the day after. The values of the delay cost evolve similarly around 24. This expresses the cost of rescheduling the departure to the day before.

We were provided with a departure choice demand matrix * The monetary price u a p , includes the amount of the toll fare. The existing toll fares were made time dependant by multiplying them by a time dependant factor, which is plotted in Fig. 3(b). This factor is greater than one (i.e. the fare is higher than usual) between 5 and 17. Its is lower than 1 between 20 and 34, meaning that the amount of the fare is lower than usual between 8 p.m of the simulated day and 10 a.m of the day after.

Calibration

Most of the motorways of the network under study are equipped with closed toll systems. As a consequence, we had at our disposal an accurate time-dependant o-d matrix for the simulated day, built from the toll stations records. This allowed for a fine grain calibration of the model, by adjusting its parameters until simulated traffic flows computed by a (fixed demand) traffic assignment match well traffic counts data, for a significant percentage of motorway sections.

Results

Once the model calibrated, we could run a user optimum assignment with departure time choice and time-varying tolls. Results were compared to the fixed-demand, constant-toll assignment. The difference between the fixed-demand and the computed demand is synthesized in Fig. 3(d). The two plots are the sum of the cumulative flows at departure, for all o-d pairs. The variations of the demand computed with departure time choice clearly indicates that a significant part of users have rescheduled their departure day to the day after. This is confirmed when observing the evolution of the pattern of congestion during the simulated day, Fig. 3(c). The plotted value is an estimate of the total time spent in congestion: at a given instant during the day, the input flow rate of each arc is multiplied by the time spent in queue on this arc; the sum for all arcs is plotted. In the fixed-demand scenario, congestion starts soon after 0 a.m. and vanishes around 9 p.m, with a significant peak between 9 a.m and 12. With departure time choice, a fraction of the demand has rescheduled its departure to the day after. This leads to a congestion period that last half longer than in the fixed demand scenario.

CONCLUSION

We have presented a model for dynamic traffic assignment with departure choice in interurban networks. While most of the models concerned with departure choice modelling are designed with commuting patterns in mind, it has been shown that the requirements for an appropriate modelling of interurban trips departure choice are substantially different. Using a continuous, rather than discrete, distribution of preferred arrival instants seems to better reflect the heterogeneity of users' preferences at arrival. The convexity of delay cost functions has to be relaxed, since it does not allow for some desired behaviours (e.g. day to day departure rescheduling). Semi empirical algorithms have been proposed, implemented, and illustrated by two numerical examples. The results are encouraging and show some of the benefits practitioners can take from using our approach when studying congestion management schemes. A substantial amount of work is still needed to properly state the appropriate numerical solution techniques.

Figure 1 :

 1 Figure 1: The bottleneck model.

  Figure 1. (a)Travel time pattern between Narbonne and Orange on the 14th of July 2007. (b) o-d flows between Narbonne and Orange (courtesy of ASF).

.

  of users of class u on the oUsing a compact notation:

  . r t ) : the route traversal cost (resp. time) function of r . instants.• D : a delay cost function.

  et Small (2001) study on SR91. p γ and p β are the ratios obtained by Small (1982) in its study for San Francisco. In dollars, the values are 87 class r , along the line of the technical reportfrom[START_REF] Sautter | Accessibilité et Modulation tarifaire[END_REF], it is assumed that the ratios

Figure 2 :

 2 Equilibrium values of the SR91 under two scenarios.

  (a)) shows the location of the VDR area and the structure of traffic flows from foreign countries. The network of main roads in the French network is mapped in Fig.2(b), along with the set of o-d pairs this example is concerned with. The main axis in the VDR area in the A7 motorway, located between Lyon (LY in Fig.2(b)) and Orange (OR in Fig.2(b)). The distance between those two cities is around 200km. The VDR area The VDR area motorways and national roads o-d pair (a) the VDR area (b) od-pairs under study.

Figure 2 :

 2 The major road network in France.

α

  traversal cost. The number of user classes, the user cost functions and the departure choice demand matrix were provide as such. They have been inferred for the simulated day from a large survey conducted by the motorway operators in year 2008.The traversal cost of an arc a for a user of class u entering the arc at an instant h is a linear combination of the monetary price u a p , and of the arc traversal time is the value of time for the user class u . The user classes vl and vls are such that vl α is by far lower than vls α .

Figure 3 :

 3 Inputs and outputs for the VDR example.

. User optimum assignment with departure time choice

												,		=	, ŨOA , RC ∈	,	,	TF	(4)
	2.2Delay cost function
	A delay cost function D is a positive and continuous mapping between delays and cost units, such that ( ) +∞ → d D when +∞ → d . In practice, D is often such
	that	D	) 0 ( =	0	, but this is not a requirement.
					h	+	t	r	,	u	h	-	h	*	=	0	. More generally, if u D is a delay cost function
	associated to the user class u , and if u i

Optimal routes and departure instants

When travelling from o to d , if a user u i x , wants to reach the destination d at a given arrival instant * h , following a given route r in i R , he has to choose a departure instant h such that ( )

Table 1 : Results from a hypothetical socio-economic analysis (values are in M$).

 1

	Tolled	Difference

The traffic of heavy vehicles is ignored since traffic regulation rules forbid truck traffic during some of the most congested days in summer.

ACKNOWLEDGMENTS

Special thanks are due to Jean Delons, Head of the Economy and Traffic Department at Cofiroute for all the support he and his team provided.


Antoine Bourgeade
email: antoine.bourgeade@cea.fr

Candice Mézel
email: mezel@celia.u-bordeaux1.fr

Olivier Saut
email: olivier.saut@math.u-bordeaux1.fr

Modeling the early ionization of dielectrics by ultrashort laser pulses

Keywords: Nonlinear optics, Modeling multiphoton ionization, Ultrashort laser pulses, Plasma Physics

In this paper, we present a model for propagation of intense and ultrashort laser pulses ionizing dielectrics. We consider early ionization so that this process is sufficiently weak to avoid requiring a complete description of the ionization process (e.g. the use of kinetic equations which are very expensive from a computational point of view). As the intensity of the field is small, one photon ionization is neglected. Ionization may only occur through multi-photonic and collisional ionization. The conduction band is discretized and multiple-rate equations are written for electron densities. The wavefield evolves through Maxwell equations.

Introduction

In this paper, we present a mathematical model adequate to describe the interaction of dielectrics with ultrashort laser pulses. This interaction involves several physical mechanisms. Provided that the laser intensity is high enough (about tens of TW/cm 2 ), electrons are generated by multiphoton ionization (MPI). They reach the conduction band where they get enough energy (more than the band gap value) to perform electronic collisional ionization. One can see [START_REF] Bergé | Ultrashort filaments of light in weakly ionized, optically transparent media[END_REF][START_REF] Couairon | Femtosecond filamentation in transparent media[END_REF] for the coupling of MPI with a propagating laser field (using the envelope approximation).

There are two other ways to perform ionization: the tunnel ionization, often competing with MPI in which the electrons cross the coulombian barrier instead of jumping over, and the above threshold ionization (ATI) in which the electrons get more than the number of photons required for MPI. These two processes are not described in this paper but could easily be included in the proposed method.

In this paper, we couple the propagation of the pulse described by Maxwell equations and a multigroup approach inspired from [START_REF] Rethfeld | Ultrafast dynamics of nonequilibrium electrons in metals under femtosecond laser irradiation[END_REF]. Thus we believe that the accounting of the ionization process will be more precise than for instance in [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF] where a temperature model of ionization is coupled with Maxwell equations. As we do not make the envelope hypothesis, the interaction will be more precisely described than in [START_REF] Bh Christensen | Modeling ultrashort-pulse laser ablation of dielectric materials[END_REF] where a multigroup model is coupled with an envelope equation for the propagation.

Femtosecond laser pulses can now be used for precise machining at nanometer scales [START_REF] White | Singlepulse ultrafast-laser machining of high aspect nano-holes at the surface of SiO 2[END_REF]. Indeed these pulses produce limited collateral damages in terms of stress waves, thermal conduction or melting. Ultrashort beams also permit to observe very short phenomena (chemical reaction in gases. . . ). Their enormous power allows them to be used for material ablation or nanomachining [START_REF] Joglekar | A study of the deterministic character of optical damage by femtosecond laser pulses and applications to nanomachining[END_REF][START_REF] Qiu | Femtosecond laser-induced microfeatures in glasses and their applications[END_REF] but does not prevent their application to biological materials [START_REF] Kennedy | A first-order model for computation of laser-induced breakdown thresholds in ocular and aqueous media[END_REF].

There are many approaches to describe this interaction. One can find models using kinetic approaches based on Fokker-Plank equations [START_REF] Apostolova | Modeling of laser-induced breakdown in dielectrics with subpicosecond pulses[END_REF][START_REF] Vatsya | Modeling of laser-induced avalanche in dielectrics[END_REF] or even physical models with an absorption term. However, these approaches are very computationally expensive. We are looking at modeling early ionization where only a small part of the material is ionized. This model should be numerically efficient and adequate to describe multi-photonic ionization. We are not really interested in rendering strong ionization.

On this topic of weak ionization, one can already find several models. For instance in [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF], free electrons are considered as a unique variable. The wave-field is driven by Maxwell equations. Two ionization models are introduced : multi-photonic and collisional ionizations. The collisional model is based on a temperature to compute collisional ionization rate. This assumes that the electrons have reached a thermal equilibrium, which may not be the case.

In [START_REF] Rethfeld | Free-electron generation in laser-irradiated dielectrics[END_REF], a multi-group approach is used through multiple rate equations. However, the evolution of the electric field is not accounted for. The ionization equations are solved in the time domain at a given spatial location for a given profile of the laser pulse. Furthermore, some probabilities describing one photon absorption by an electron are difficult to obtain or estimate experimentally.

In this paper, we have combined both approaches. For this matter, several assumptions are made. The displacement of free electrons is neglected as we are interested in weak ionization and short pulses. Finally, the electrons ionizing by the collision process turn back to the lowest energy level of the conduction band. The ionization model is coupled with a Maxwell equations solver. This allows us to obtain a model adapted to the early ionization of dielectrics by ultrashort laser pulses.

The outline of this paper is as follows. Section 1 describes the mathematical model. Section 2 describes the numerical methods and algorithms that were used to discretize this model. Finally, Section 3 presents some numerical experiments highlighting the accuracy of our model.

Theoretical model

In our model, we consider several phenomena affecting electrons and neutrals. It includes multiphoton and collisional ionization, ohmic heating of free electrons by the laser pulses, radiation (excitons) and collisions between electrons and electrons/neutrals.

Valence band

Conduction band E n e r g y

Wave-field ionization

Collisional ionization

Figure 1: Different phenomena involved in generation of free electron gas.

In this paper, we consider a pulse (of wavelength 800 nm) propagating in dielectric (silica). We assume that there is only one electron on the highest level of the valence band and that this level is not realimented. This is justified as strong ionization will not be considered in this work. The extension to more complicated materials is straightforward. Accounting for the whole valence band is the subject of a forthcoming work.

The laser pulse will evolve through Maxwell equations. For the free electrons, we use multiple rate equations (as in [START_REF] Rethfeld | Free-electron generation in laser-irradiated dielectrics[END_REF]) and a renormalization procedure to ensure energy conservation and render ohmic heating. The polarization of the material is given by a Bloch equation [START_REF] Bourgeade | Numerical methods for the bidimensional Maxwell-Bloch equations in nonlinear crystals[END_REF] or is directly computed from the number of neutrals [START_REF] Bourgeade | Dynamic Load Balancing Computation of Pulses Propagating in a Nonlinear Medium[END_REF].

Light description

To model the propagation of the laser in two and three dimensions, we use the fully time-dependent Maxwell equations. The laser pulse is described by the pair (D, H) evolving through Maxwell equations

       µ 0 ∂ t H = -∇ × E, ∂ t D = ∇ × H -J, ∇ • D = 0, ∇ • H = 0, (1)
where we have D = ǫE + P. The two relations ∇ • D = 0 and ∇ • H = 0 do not play any role in the sequel and are satisfied by the numerical schemes.

The interaction of light with matter occurs through the polarization P and current J. The current J will be obtained from the densities and velocities of free electrons. We will see later on how to calculate these terms. We choose to present two different approaches for the computation of these terms in order to have a better control on the results.

Free electrons

Free electrons are described by their densities according to their energy. As in [START_REF] Rethfeld | Free-electron generation in laser-irradiated dielectrics[END_REF], we discretize the conduction band in N p levels whose mean energies (starting from the fundamental level i = 0) are

E i = (i + 1 2 ) ω L , 0 ≤ i ≤ k, (2)
where ω L is the photon energy of laser light and k = [ Ecrit ω L + 1] = N p -1, E crit being the critical energy of ionization.

Remark 1 In fact, as the ohmic heating of the electrons is accounted for separately, there is no need to consider energy levels separated by the energy of a photon. The conduction band can be discretized differently as shown in Sec. 4.3 In the sequel, in our applications, we take N p = 6 (as we are considering fused silica but the extension to more energy levels is straightforward). The density of free electrons in the energy level i at time t and point x will be denoted by n i (t, x). In particular, the density of electrons above the critical energy for impact ionization is denoted by n k .

For the time being, we consider three processes acting on the free electrons (see Fig. 1). The first one, multi-photon ionization depends on the laser pulse and the number of neutrals. The second one, collisional ionization, depends on the number of free electrons having reached the ionization energy and of the number of neutrals. Finally, the electrons also lose energy by radiation.

Collecting our assumptions, we can write them down as evolution equations. As in [START_REF] Rethfeld | Free-electron generation in laser-irradiated dielectrics[END_REF], we have to describe the evolution of each group density. The densities evolve through

   ∂ t n 0 = α ion N -α des n 0 + 2α col N N0 n k , ∂ t n i = -α des n i , 0 < i < k, ∂ t n k = -α col N N0 n k -α des n k , (3)
where the rates α ion , α des , α col can be found in Sec. 2.7. The instantaneous number of neutrals is denoted by N . Let us also note that we have N + i n i = N 0 at every time. The evolution equation on N is given in Eq. ( 6). The factor 2 in the first equation of system (3) comes from the fact that through collisional ionization, one electron in the k-th energy level ionizes one electron in the valence band: these two electrons are then in the lowest level of the conduction band.

Remark 2 There is no diffusion between groups, the cascading effect is not considered yet. These transitions will be accounted for later to ensure the conservation of energy in Sec. 2.5.

To each energy group, we associate a mean velocity denoted by v i for the free electrons in the energy group E i . The evolutions of the different velocities of the free electrons are calculated using Drude's law [START_REF] Mao | Dynamics of femtosecond laser interactions with dielectrics[END_REF]:

∂ t v i = Q e m e E -ν (i) e-v i , 0 ≤ i ≤ k, (4)
where Q e is the charge of an electron, m e its mass and ν

(i)
e the collision rates containing electron-electron and neutral-electron collisions. Their expressions are obtained from

ν (i) e = σ eff 2E i m e n i , (5)
as in [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF].

Neutrals and bound electrons

We only consider one electron by atom. This electron is in the highest energy level. For each point x and time t, we denote by N (t, x) the number of neutrals. According to our assumptions, this also corresponds to the number of bound electrons in the highest level of the valence band.

The number of neutrals increases with the energy lost by radiation by the free electrons. This number decreases by multi-photon and collisional ionization. The total density N of atoms obeys

∂ t N = k i=0 α des n i -α ion N -α col n k N N 0 , (6)
N 0 being the number of neutrals at equilibrium without ionization. The transition rates α ion , α des , α col are the same as those found in Eq. ( 3).

Their expressions are given in Sec. 2.7. In this work, we do not consider that other species than neutrals could be produced by this recombination (e.g. negative ions).

To validate our rendering of the ionization process, we choose to use two different approaches for describing laser-matter interaction. The first one at the quantum-mechanical level uses the Bloch equation to describe the evolution of the material and the Maxwell equations for the light where the polarization term comes from the microscopic description. The second one is a model based on the Maxwell equations where the polarization term takes the dispersive effects into account.

Using the Bloch equation and the density matrix formulation [START_REF] Boyd | Nonlinear Optics[END_REF][START_REF] Newell | Nonlinear optics. Advanced Topics in the Interdisciplinary Mathematical Sciences[END_REF]: We describe the atoms according to their energy using the density matrix formalism. In this formalism, the material is statistically described at the quantum-mechanical level. Each molecule of the medium is considered as a quantum system with 2 discrete levels of energy.

The density matrix ρ represents the various populations in each of the energy levels of the free Hamiltonian in its diagonal terms. The off-diagonal terms of the density matrix represent the quantum coherences of a set of atomic states. It was shown in [START_REF] Besse | A Maxwell-Bloch model with discrete symmetries for wave propagation in nonlinear crystals: an application to KDP[END_REF], that this approach is adequate to describe a wide range of dispersive dielectrics.

Let us denote by ω jk = ω j -ω k the difference between the two frequencies associated to levels j and k. We shall also denote by µ the dipolar matrix of the material (depending on its microscopic structure) [START_REF] Boyd | Nonlinear Optics[END_REF][START_REF] Newell | Nonlinear optics. Advanced Topics in the Interdisciplinary Mathematical Sciences[END_REF]. The potential of interaction of the wave-field with the material is given by V = -µ • E.

The density matrix ρ satisfies

∂ t ρ ℓℓ = i [V, ρ] ℓℓ , (7)
(where we remind that for two operators A and B, [A, B] = AB -BA) for the population terms and the following equations for the coherences.

∂ t ρ jk = -ıω jk ρ jk + ı [V, ρ] jk , j = k. ( 8
)
The polarization is given by

P a = N tr (µρ),
hence its time derivative is computed through

∂ t P a = tr (µρ)∂ t N + N tr (µ∂ t ρ). (9)
This approach will be useful for introducing the influence of the lower level electrons.

Using a nonlinear polarization in the Maxwell equations: The polarization in ( 1) is obtained from

P = N N 0 (P L + P N L ), (10)
where the linear P L and nonlinear P N L terms are computed as in [START_REF] Bourgeade | Dynamic Load Balancing Computation of Pulses Propagating in a Nonlinear Medium[END_REF]. The number of neutrals without the laser pulse is denoted by N 0 . This approach is useful to account for a Kerr-type nonlinearity.

Interaction terms

The currents involved in Maxwell equations ( 1) are decomposed into J = J e + J mpi . The first one is due to the movement of free electrons and the ionization is responsible for the second one. The first current is obtained from the expression

J e = i0 i=0 Q e n i v i . (11)
Then we have to compute the current due to photonic ionization. Its expression is given by

J mpi = α ion W ion N E 2 E, (12)
where for silica W ion = 9|Q e |, Q e being the charge of an electron, E•J mpi = α ion W ion is the ionization energy [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF][START_REF] Stuart | Nanosecond-to-femtosecond laser-induced breakdown in dielectrics[END_REF]. If the electric field is vanishing, we take J mpi = 0. The expression of α ion is found in Sec. 2.7 and depends on the intensity of the incoming wave-field.

Diffusion between energy groups: ohmic heating

At this stage, we have not taken into account the effect of the acceleration of electrons by the laser and the transition it generates between energy groups. In [START_REF] Rethfeld | Free-electron generation in laser-irradiated dielectrics[END_REF], this is achieved directly in the rate equations. However, the corresponding probability or avalanche coefficient is difficult to recover experimentally and is generally intuitively estimated [START_REF] Arnold | Acoustic-phonon runaway and impact ionization by hot electrons in silicon dioxide[END_REF][START_REF] Stuart | Nanosecond-to-femtosecond laser-induced breakdown in dielectrics[END_REF]. We choose another approach to render these physical effects. Equations (3) do not contain cascading terms. Free electrons are heated by the laser pulse and their energy increases. We have assumed that there are n i electrons with mean energy E i (or with an energy between i ω L and (i + 1) ω L ). If these electrons are heated they may leave this group to a group of electrons with higher energy. We choose to compute the number of these electrons leaving the current group and add them to the appropriate group. Yet, if we add these electrons to a different group, the quantity of movement of the first group changes. Then, we also have to adjust the velocities in each group.

Initially, we consider that all electrons are in the lowest level n 0 . We will denote by i 0 (t, x), the highest populated energy level at the point x and time t. We have i 0 (0, x) = 0 and i 0 (t, x) ≤ k.

In the sequel, we use the notation ω c = Wion k , where W ion is the energy of ionization of the material.

Renormalization of densities: The laser heating makes electrons gain energy and thus they may change from one energy level to another. We will denote by E * i the energy obtained from E i . If E * i > E, there will be a flux of electrons from the energy group E i to E i+1 .

The principle of the projection is shown in Fig. 2. We consider that the electrons are found uniformly in each group. So the transition E 0 → E * 0 gives rise to a flux

δn 0 = E * 0 -E 0 E * 0 -0 n 0 , (13)
corresponding to the second gray rectangle in Fig. 2 and for the higher groups we have

δn j = E * j -E j E * j -E * j-1 n j . (14)
We have

E * j -E * j-1 = E * j -E j + E j -E j-1 + E j1 -E * j-1 , = ∆E j + ω c + ∆E j-1 ,
where ∆E j = E * j -E j . To compute the variation of energy, we use

∂ t E i = Q e v i • E, (15)
which gives the flux between two groups

∆E i ≡ E * i -E i = δtQ e E • v i + v i-1 2 , ( 16
)
Density Energy 0 !" 2!" n 0 n 1 ! 0 , v 0 ! 1 , v 1 ! 0 * , v 0 * Figure 2:
Principle of the computation of the new densities to take heating into account. The red and blue rectangles represents the two initial groups corresponding to energies E 0 and E 1 . After heating, these two groups are shifted and have to be projected to compute the new densities of each energy level. We assume that the electrons are uniformly distributed in each group.

where we have approximated the velocity of the interface between groups i -1 and i by vi+vi-1

2

. We compute

   ñ0 = n 0 -δn 0 , ñi = n i + δn i-1 -δn i , 1 ≤ i ≤ i 0 -1 ñi0 = n i0 + δn i0-1 . ( 17
)
Remark 3 It shall be noted that this process does not change the total number of free electrons.

Renormalization of velocities: The total momentum per group must also be computed as we have changed the electron densities. The momentum of the energy group E i is given by

p i = m e n i v i . ( 18
)
The quantity of movement leaving the first group at the interface is given by

δp 0 = -m e δn 0 v 0 + v 1 2 , (19)
and for the other levels we obtain

δp i = m e δn i-1 v i-1 + v i 2 -δn i v i + v i+1 2 . ( 20
)
Then we compute the adjusted velocities with

ṽi = m e n i v i + δp i m e ñi , 0 ≤ i ≤ i 0 . (21)
Renormalization of last group: As electrons are heated, their energy can exceed the threshold E i0 . In this case, we need to add an additional energy level (for instance, with i 0 = 1, we have to add an energy level to account for free electrons in the right-most gray rectangle of Fig. 2). For this matter, we use the conservation of energy.

Let us denote by U i the total energy of the group i. We have U i = (i + 1 2 ) ω c n i for 0 ≤ i < i 0 , the energy in the last group ensures the conservation of the total energy.

First, we compute the total energy in the material

S o = i0 i=0 U i , S = i0-1 i=0 (i + 1 2 ) ω c n i . (22)
The total energy created by ohmic heating during dt is given by

S c = i0 i=0 Q e n i (v i • E)dt, (23)
and the energy due to radiation is

S r = i0-1 i=0 (i + 1 2 ) ω c n e (1 -e -ν des dt ) + U (i 0 )(1 -e -ν des dt ), (24)
where we have used Eq. ( 3) to evaluate the energy lost by radiation during a time dt.

The energy of the last level is given by

U (i 0 ) = S o -S -S r + S c + δn k ω c - U (k) n k . ( 25
)
We shall now decide if there is a need to add another energy level (assuming that i 0 < 6). If U (i 0 ) > n i0 (i 0 + 1 2 ) ω c , we have to create a new level. To ensure that the last level is not almost empty, we will only add a level if

U (i 0 ) > n i0 (i 0 + 1) ω c .
The energy per electron is given by U (i0) ni 0

. Then assuming that the electrons are uniformly distributed we get

δn(i 0 ) = n i0 U (i0) ni 0 -(2i 0 + 1) ω c U (i0) ni 0 -i 0 ω c , (26)
and we set the new values of the densities n * (i 0 ), n * (i 0 + 1) and energies U * (i 0 ), U * (i 0 + 1):

n * i0 = n i0 -δn(i 0 ) ; n * i0+1 = δn(i 0 ) U * (i 0 + 1) = U (i 0 ) -(i 0 + 1 2 ω c )n * i0 ; U * (i 0 ) = (i 0 + 1 2 ω c )n * i0 . (27)
Finally, we have to compute the adequate velocities. For that matter, we use the same method as in the previous paragraph. We ensure that we have the conservation of momentum. The momentum of the group we have just added must correspond to the initial momentum of the group i 0 minus the momentum of the updated group i 0 .

For the sake of simplicity, let us assume that i 0 = 0. We denote by v the velocity of the interface between groups 0 and 1. As in Eq. ( 19), the momentum leaving the first group is given by

δp 0 = -m e δn 0 v, (28)
but this time we have to compute the velocity of the interface. We choose to use a linear approximation i.e. we take v = 2v 0 . We have conservation of the momentum, which (using the same notations as before) leads to

m e n * 0 v * 0 = m e n 0 v 0 + δp 0 . (29)
We replace δp 0 and v by their expressions and we obtain

v * 0 = 1 - δn 0 n * 0 v 0 . (30)
Finally, we take

v * 1 = v. ( 31
)
Remark 4 One can also use an approximation of higher order.

Adding a relaxation process toward thermal equilibrium

We now assume that the electronic densities are relaxing toward a Maxwellian distribution [START_REF] Stuart | Nanosecond-to-femtosecond laser-induced breakdown in dielectrics[END_REF]. That is that we assume a Maxwellian distribution for the electron velocity i.e.

F (v) = 1 √ πv th exp(- v 2 v th ). ( 32
)
To obtain the corresponding energy distribution function f e , we take

W = W th v 2 v 2 th ,
where W th = 1 2 mv 2 th = 3 2 k B T e (T e being the electronic temperature). This gives

f e (W ) = 1 2 √ πW th W exp(- W W th ). (33)
The equations on electron densities become

   ∂ t n 0 = α ion N -α des n 0 + 2α col N N0 n k -1 τm (n 0 -n eq 0 ), ∂ t n i = -α des n i -1 τm (n i -n eq i ), 0 < i < k, ∂ t n k = -α col N N0 n k -α des n k -1 τm (n k -n eq k ), (34)
where we have denoted by τ m the relaxation time toward equilibrium (identical for each level).

We shall now compute the densities at equilibrium. We have

n eq i = j n j kT (i+1) ω i ω exp(- W kT )dW , 0 ≤ i < k = ( j n j ) 1 -exp(- ω kT ) exp(- i ω kT ), 0 ≤ i < k, (35)
and

n eq k = ( j n j ) exp(- k ω kT ). ( 36
)
Let us note that we have j n eq j = j n j .

Remark 5 In our experiments, most of the time the highest level is not populated. Therefore this relaxation process toward equilibrium will accelerate the ionization.

Probabilities

We have to determine the probabilities associated to each of the three processes involved. We have taken the expression found in [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF]. For reader's convenience, we recall these expressions shortly.

For the probability of collision, we have used the expressions from [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF], namely

α col = 2α 0 col √ π N N 0 ( 15 4
β 3 2 - 1 2 β)e -1 β + π 2 ( 15 4
β 2 -3β + 1)erfc( 1 √ β ) (37) where β = W th
Wion , N 0 = 2×10 28 m -3 and α 0 col = 1.5×10 17 . We have denoted by erfc the function 1 -2 √ π x 0 exp(-t 2 )dt. The probability of radiation is also a constant value:

α des = 1 τ trap = (1.5 × 10 -13 ) -1 s -1 ,
see [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF] for details. The influence of this phenomenon is probably limited with ultrashort pulses.

The probability of ionization depends on the intensity of the incoming wave-field.

α ion = σ 6 N 0 10 -12 E × H 6 , (38)
as in [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF] we took σ 6 = 2 × 10 7 m -3 s -1 (TW/m 2 ) -6 and N 0 = 2 × 10 28 m -3 .

3 Numerical method

Discretization of the Maxwell equations

For discretizing Maxwell equations, we use a finite-difference scheme based on Yee's work [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF]. The physical variables are computed using the scheme shown in Figure 3. For a detailed presentation of this scheme, one can see [START_REF] Bourgeade | Numerical methods for the bidimensional Maxwell-Bloch equations in nonlinear crystals[END_REF]. The discretization of the model based on nonlinear Maxwell equations is described in [START_REF] Bourgeade | Dynamic Load Balancing Computation of Pulses Propagating in a Nonlinear Medium[END_REF].

z y t = nt t =(n + d r d 1 2 -) t D x D y D z H z H y H x ,, P

Computation of free electrons densities

The computation of the free electron densities given by ( 3) is decomposed in two parts using a splitting method. The first one is completely linear

         ∂ t n 0 = α ion N -α des n 0 -1 τm (n 0 -n eq 0 ), ∂ t n i = -α des n i -1 τm (n i -n eq i ), 0 < i < k, ∂ t n k = -α des n k -1 τm (n k -n eq k ), (39)
and the second one is nonlinear but only involves the lowest and highest energy level:

   ∂ t n 0 = 2α col N N0 n k , ∂ t n k = -α col N N0 n k , N = N 0 -n 0 - k-1 i=1 n i -n k .
(40)

Resolution of (39)

As the system corresponding to (39) is linear, we will use an exponential scheme. Hence we will compute

n n+1 j = n n j e -(α des (j)+ 1 τm )dt + n eq i τ m α des (j) + 1 1 -e -(α des (j)+ 1 τm )dt , 1 ≤ j ≤ k,
(41) and for the lowest level we have

n n+1 0 = n n 0 e -(α des (0)+ 1 τm )dt + α ion N n+ 1 2 + 1 τm n eq 0 α des (0) + 1 τm 1 -e -(α des (0)+ 1 τm )dt .
(42)

Resolution of (40)

We have to compute the densities from (40) (i.e. we are splitting the equations on n 0 and n k in time). Summing the two first equations of (40), we get

∂ t (n 0 + 2n k ) = 0, hence n n+1 0 + 2n n+1 k = n n 0 + 2n n k . We compute n n+1 k = e -α col N n+ 1 2 N 0 dt n n k , (43)
and

n n+1 0 = n n 0 + 2(n n k -n n+1 k ). (44)

Computation of the electrons velocities

One has to compute the velocities with equations (4). For this purpose, we use an exponential scheme yielding

v n+1 i = e -νedt v n i + Q e m e ν e (1 -e -νedt )E, (45)
where this equation is taken component by component.

Numerical results

The two approaches described in Sec. 2.3 yield identical results. Both will be used in this section.

Propagation of an ultrashort pulse in fused silica

In this experiment, we compute the propagation of an ultrashort pulse in fused silica. The intensity of the 50 fs pulse (λ = 800 nm) is 1.5×10 10 V/m. We have considered 3 µm of silica. The initial density of atoms is N 0 n = 2 × 10 28 m -3 , and we have n 0 e = 10 10 m -3 free electrons at each point of the dielectric.

We have taken 100 points per wavelength in the direction of propagation and 150 points in the transverse direction. The time step is computed to ensure the CFL condition.

The initial wave-field is given by

E(t, y) = E max sin( πt τ ) cos(ω 0 (t - τ 2
))e

-1 2 ( y wy ) 2
H(τ -t), where w y = 0.5 µm, E max = 1.5×10 10 V/m, τ = 50 fs. Thus the input laser intensity is roughly 8.8 × 10 13 W/cm 2 and the peak power 689 kW.

The amplitude of the electric field E y after 35 fs is shown in Fig. 4. The absorption at the center of the pulse is clearly visible. To underline the nonlinear effects, the spectrum of this field is shown in Fig. 5. Harmonic generation is clearly observed and only odd harmonics appear as expected. The ionization process is a source of nonlinearities, the nonlinear index of the material is also taken into account.

We now show the evolutions of the different populations of electrons at the front of the dielectric in the center of the pulse. These evolutions are shown in Fig. 6. The medium is almost completely ionized by the laser beam and a large number of electrons reach the highest energy level of the conduction band.

Influence of Maxwellian relaxation

In this section, we study the effect of the relaxation toward a Maxwellian equilibrium. This process is added in our model in section 2.6.

We use the same experimental setup than in the previous experiment with two different values for the relaxation time (τ m in Eq. ( 34) ). Let us note that the value τ m = 1 fs is shorter than one optical cycle. We choose such a small value to observe the effect of the relaxation process.

The result of the run is shown in Fig. 7.

As initially higher levels are not populated, the relaxation process increases ionization. In order to be able to observe the effect of this relaxation, we had to take a very short relaxation time (τ m = 1 fs) which is probably not realistic. With a more reasonable value, there is almost no difference with the evolution we observe without this Maxwellian relaxation (in Fig. 6).

The effect of this relaxation on the number of free electrons is quite small, we believe it can be safely neglected. The electrons are not very hot and the relaxation time is thus too long for this process to be influent.

Influence of the discretization of the conduction band

Contrary to other models like [START_REF] Rethfeld | Free-electron generation in laser-irradiated dielectrics[END_REF], the one-photon absorption by free electrons is not directly accounted in the equations (3). These energy levels do not have to be separated exactly by the energy of one photon. Thus we can discretize the conduction band as we wish as long as the difference between the first and the last level is conserved.

In this experiment, we wish to study the influence of this discretization on the ionization. For instance, we can see if a finer discretization is needed to better account for the cascading effect. We have run the same experiment with three different discretizations of the conduction band. One has 5 levels, the second one is the same we have used before (7 levels) and last one has 2 additional levels.

It shall also be noted that, of course, refining the discretization increases the computation time and that finer discretizations may require smaller time steps (however it is not the case in this test).

We have shown the evolution of the number of neutrals and the number of free electrons in the lowest energy level in Fig. 8.

There is almost no difference on the number of neutrals (or the total number of free electrons). The population of electrons in the level E 0 is a 9: Energy fluxes in J/m 2 (1D, first picture), J/m (2D, second picture) over width (m) for our model (i.e. quantities prefixed with "Mg") and the model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF] (quantities prefixed with "Pen") in 1D and 2D (in the center of the dielectric). We plot for each approach the absorbed, transmitted and total energy. Figure 10: Energy (J) over width (m) obtained from our model (quantities prefixed with "Mg") and the model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF] (quantities prefixed with "Pen") in 3D. We plot for each approach the absorbed, transmitted and total energy.

Experimental comparison

In this section, we will use the experimental results presented in [START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF] to study the accuracy of our model.

The detailed description of the experimental setup is given in [START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF]. The authors analyze the phase shift of probe laser pulse propagating in a sample of SiO 2 by observing the pattern resulting from interferences between this pulse and a reference beam. The probe wavelength is 618 nm, the sample temperature is 300 K. We use the fitting curve they have built from their experimental measurements,

We have plotted the phase shift predicted by the multigroup model and the temperature model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF] and the results measured in [START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF]. The results are shown in Fig. 14 for an intensity of 1.4 TW/cm 2 and in Fig. 15 for an intensity of 2.1 TW/cm 2 .

From these two tests, we may conclude that the multigroup approach appears to be more accurate than the temperature model. However the central part of the phase shift differs a little bit between the multigroup approach and the results of [START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF].

With the multigroup model, collisional ionization only occurs once the higher energy level is populated. Thus for weak ionization, collisional ionization happens later than with the temperature model. In [START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF], the ex- perimental results are fitted with a model neglecting collisional ionization. For the lower intensity 1.4 TW/cm 2 on Fig. 14, there is almost a perfect match between the results of [START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF] and our model.

Conclusion

We have presented a mathematical model for early ionization in dielectrics. It uses multiple-rate equations for unbound electrons and neutrals. The laser pulse is described by Maxwell equations. The cascading effect and ohmic heating are rendered by a procedure ensuring the conservation of energy.

Comparing to pre-existing works, we have a full time-space domain model (contrary to [START_REF] Rethfeld | Free-electron generation in laser-irradiated dielectrics[END_REF]) where we do not need an intuitive expression for the one-photon intraband absorption and using multiple rate equations we can render cascading effects (contrary to [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF]).

The model was kept (relatively) simple : several phenomena, which we believe do not play necessarily an important role with ultrashort pulses, were neglected. For instance, the trapped electrons (accounted for by the term in α des in Eq. (3)) are not set aside whereas they could be considered as excitons. We have also neglected the fact that when the highest energy 14: Phase shift (in rad) in SiO 2 for the three approaches for an intensity of 1.4 TW/cm 2 . "Pen" denotes the curve obtained with the temperature model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF], "Mg" the curve obtained with the multigroup model described in this paper and "Exp. results" the results from [START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF].

Figure 15: Phase shift (in rad) in SiO 2 for the three approaches for an intensity of 2.1 TW/cm 2 . "Pen" denotes the curve obtained with the temperature model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF], "Mg" the curve obtained with the multigroup model described in this paper and "Exp. results" the results from [START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF]. level of the valence band is empty, the pulse can interact with the lower layers of the valence band to repopulate this level. However we believe that for the pulse durations we consider, the role of this phenomenon is probably quite small. Yet, this is the subject of a forthcoming work, which will also introduce tunneling ionization we have not considered here.

Some parameters that we have used are difficult to determine precisely. In particular the ionization rate σ 6 in Eq. (38) comes from [START_REF] Keldysh | Ionization in the field of a strong electromagnetic wave[END_REF] and the collisional rate α ion in Eq. ( 37) is not rigorously justified. The results obtained also show that the observed phenomena depend on the model. The recovery of parameters has to be coupled to a specific model.

Even if our numerical results are qualitative, they show that with strong ionization, our multi-group approach yields different simulations that the temperature model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF]. This may prove that there is no electronic temperature. However, our model is much more computationally expensive than [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF]. Our numerical codes are adapted to run on parallel architectures.

Figure 3 :

 3 Figure 3: Scheme used to discretize the physical variables.

Figure 4 :

 4 Figure 4: Amplitude of the electric field E y after 35 fs of propagation. The pulse propagates from bottom left to top right.

Figure 5 :

 5 Figure 5: Fourier transform of the wave-field in the y polarization in the center of the material after 35 fs of propagation. The amplitude is shown in arbitrary units in a log scale. The frequency scale is normalized by the frequency of the incoming wave.

Figure 6 :

 6 Figure 6: Evolution of the populations of electrons at the center of the pulse at the front of the material over time. The plain lines denote the different populations of free electrons corresponding to the seven energy levels (0,. . . 6) and the dotted line represents the population of neutrals.

Figure 7 :

 7 Figure 7: Electron populations over time in the setup of Sec. 4.1 with τ m = 100 fs (first picture) and τ m = 1 fs (second picture) in the relaxation process added in Sec. 2.6.

Figure 8 :

 8 Figure 8: Electron populations over time in the setup of Sec. 4.1 with one finer (9 points), one standard (7 points) and one coarser discretization (5 points) of the conduction band.

Figure

  Figure 9: Energy fluxes in J/m2 (1D, first picture), J/m (2D, second picture) over width (m) for our model (i.e. quantities prefixed with "Mg") and the model[START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF] (quantities prefixed with "Pen") in 1D and 2D (in the center of the dielectric). We plot for each approach the absorbed, transmitted and total energy.

Figure 11 :

 11 Figure 11: Total number of free electrons over time for the two models for three different points in the dielectric over time in 1D (first picture) and 3D (second picture). The first point is located at the front of the dielectric, the second one after 1.5 µm and the last one at the back.

Figure 12 :

 12 Figure 12: Absorbed energy for the temperature model. The absorbed energy is represented in J. The vertical slice at the back of the picture is located at z = 10 -7 m of the front of the material. The two other slices are at the center of the x and y directions. The same color scale is used for the three slices.

Figure 13 :

 13 Figure 13: Absorbed energy for the multigroup model. The contours shown represents five uniformly distributed values ranging from 2.6 × 10 10 J (i.e. the maximum obtained with the temperature model) to 3.2 × 10 11 J (not included). The absorbed energy is represented in J. The vertical slice at the back of the picture is located at z = 10 -7 m of the front of the material. The two other slices are at the center of the x and y directions. The same color scale is used for the three slices.

Figure

  Figure14: Phase shift (in rad) in SiO 2 for the three approaches for an intensity of 1.4 TW/cm 2 . "Pen" denotes the curve obtained with the temperature model[START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF], "Mg" the curve obtained with the multigroup model described in this paper and "Exp. results" the results from[START_REF] Martin | Subpicosecond study of carrier trapping dynamics in wide-band-gap crystals[END_REF].

little smaller with the coarser discretization. As the second level is closer to the first one with a finer grid, it is probably easier to reach for an electron. Hence, the population in the lowest level is smaller.

4.4

Comparison with the model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF] In this experiment, we compare our model with the model described in [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF]. This model only uses one group for the free electrons and assumes that these electrons are at thermal equilibrium.

In Fig. 9, we have plotted the evolution of energy fluxes of our model and the model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF]. Computations were performed in 1D and 2D (the picture represents the evolution of the flux inside the dielectric as a function of the depth).

The same experiment is also run in 3D. The result is shown in Fig. 10. In Fig. 11, we have plotted the total number of free electrons for three different points in the dielectric over time for the two models. These points are uniformely distributed at the center of the material in the 3 µm of the width.

Finally, we have plotted, in Fig. 12 and 13, the absorbed energy. The pulse propagates from the back of the figure in the z-direction. The first picture represents the absorbed energy (in J) for the temperature model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF]. We can observe that the absorption occurs mainly at the surface of the dielectric. The width of the material is as before 3 µm.

In Fig. 13, we show the result of the same experiment with the multigroup model. As before, one can observe that the material is more ionized with this model than with the temperature model. The largest contour shown in this picture represents the maximal value of absorbed energy obtained with the temperature model (more precisely this contour corresponds to E abs = 2.6 × 10 10 J).

Even if the results look similar, with the model [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF], the material is less ionized than with our model. Without additional information from experiments, it is hard to determine which model is the most accurate. However, let us note that our model is obtained from less restrictive assumptions than [START_REF] Peñano | Transmission of intense femtosecond laser pulses into dielectrics[END_REF]. Indeed, the electronic distribution observed with the multi-group approach does not correspond to a thermic equilibrium.

It is nevertheless interesting to note that because the absorption occurs more in the bulk for the temperature model, the reflection of the pulse is less effective. This explains why the energy inside the dielectric is greater for the temperature model although the pulses are the same. This effect is more visible in 1D than in 2D or 3D. Actually, the reflection in 3D is quite negligible.