Alaeddine El Fawal
email: alaeddine.elfawal@epfl.ch

Jean-Yves Le Boudec
email: jean-yves.leboudec@epfl.ch

Kave Salamatian
email: kave.salamatian@epfl.ch

Vulnerabilities in Epidemic Forwarding

HAL is

Introduction

In this paper, we identify vulnerabilities in epidemic forwarding over ad-hoc networks. We are interested in broadcast applications such as chatting in a traffic jam or coupon advertisements [START_REF] Garyfalos | Coupons: Wide scale information distribution for wireless ad hoc networks[END_REF]. The principle of epidemic forwarding is that nodes repeat with some probability the information they hear from others, thus propagating fresh information.

Epidemic forwarding employs several mechanisms such as inhibition, which prevents a node from forwarding oversent or over-received packets in order to minimize redundancy. Each mechanism can be implemented using different alternative methods. Thus, the existence of vulnerabilities is highly dependent on the mechanisms employed and on the methods adopted to achieve them. We examine the links between these methods and the vulnerabilities.

We classify the vulnerabilities into two categories: malicious and rational. A malicious attacker harms other nodes and does not look for personal benefit. It aims at decreasing other nodes' throughput and/or spread. In contrast, a rational attacker aims at increasing its personal profit from the network. It tries to increase its throughput or save power.

We evaluate the vulnerabilities by simulations. We show that a malicious attacker does not have much effect in highly mobile networks, but it might be very harmful in static networks. In contrast, independently of mobility, the rational attacker increases dramatically its throughput. Attacks that are otherwise very harmful lose their efficiency in the presence of some epidemic forwarding mechanisms such as adaptive spread control and adaptive inhibition. Moreover, many elements such as the attacker position and the node density influence the malicious attacker's impact.

The simulations are carried out on networks with up to 1000 nodes using a JAVA implementation of the epidemic forwarding system in JIST-SWANS [START_REF]Java in simulation time / scalable wireless ad hoc network simulator[END_REF]. Beside static scenarios, we applied the epidemic forwarding system to the vehicular network using an extension of JIST-SWANS called STRAW [START_REF]Street random waypoint / vehicular mobility model for network simulations[END_REF], which provides a mobility model based on the operation of real vehicular traffic.

Epidemic Forwarding Mechanisms

In this section, we explain the different mechanisms used in epidemic forwarding, in order to understand their vulnerabilities.

Inhibition Mechanisms

Inhibition mechanisms aim at preventing nodes from forwarding over-sent or over-received packets in order to minimize redundancy. We classify them into two sets: rigid and adaptive. With the former set, the mechanisms cannot adapt themselves to different network settings: when the settings change, their parameters need to change. The adaptive mechanisms ensure a good performance in a wide range of settings without changing their parameters.

2.1.1. Rigid Inhibition. Within this set we find Gossipbased epidemic forwarding [START_REF] Modiano | Maximizing throughput in wireless networks via gossiping[END_REF] where a node decides to forward a packet with a fixed probability p and drop it with (1 -p).T h ev a l u eo fp depends on the setting but Gossip does not involve any mechanism to adapt p.

Adaptive Inhibition.

Within this set we distinguish between two methods.

2.1.2.a. Counter Based Inhibition:

This method is essentially the one proposed in [START_REF] Ni | The broadcast storm problem in a mobile ad hoc network[END_REF]. A packet stored in the epidemic buffer has a counter called "Receive Count" incremented by 1 when a duplicate of this packet is received. Initially, i.e. when the packet is created by the application or received for the first time, the counter is set to 0. When the counter reaches a maximum value, the packet is discarded from the epidemic buffer. When a packet is transmitted, the value of Receive Count is lost.

2.1.2.b. Virtual Rate Based Inhibition: This method is proposed in [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF]. With this method, a packet in the epidemic buffer is retransmitted with a probability that depends on its "virtual rate"; it is equal to c 0 a R b S where c 0 is a constant (inverse of a time), R [resp. S] is the number of times this packet or a duplicate was received [resp. sent] and a and b are unit-less constants less than 1. Thus the virtual rate of a packet decreases exponentially with any send/receive event of the same packet. A scheduler decides which packet is selected next for transmission by the MAC layer; it serves packets with a rates not exceeding their virtual rates. Hence, a packet in the epidemic buffer, which has seen many send/receive events, is scheduled at a very low rate and it is more likely that it will be dropped by the used buffer management mechanism before being transmitted [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF]. The constant c 0 is equal to the nominal packet rate of the MAC layer.

Spread Control Mechanisms

Spread control mechanisms are essential for epidemic forwarding, as the broadcast capacity does not scale with the population. Spread control can be implemented using one of the following methods.

2.2.1. Classic TTL. This is the method that comes by default with the Internet Protocol (IP). When a packet is created by a source and placed into the epidemic buffer, it receives a TTL value equal to some positive constant "MaxTTL". When the packet is accepted for transmission by the MAC layer, the TTL field of the transmitted packet is equal to the value of the TTL field in the packet in the epidemic buffer, minus 1. The TTL field in the packet stored in the epidemic buffer is unchanged.

When a packet created by some other node is received for the first time at this node, the packet is delivered to the application, and the value of the TTL is screened. If it is equal to 0, it cannot be retransmitted and the packet is discarded. Else (TTL≥ 1), the packet is stored in the epidemic buffer, with TTL equal to the value present in the received packet. When and if the packet is later accepted for transmission by the MAC layer, the transmitted TTL field is equal to the stored TTL minus 1, and the stored TTL is unchanged.

2.2.2. Aging. This method is proposed in [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF] in a different but essentially equivalent form. We give here a presentation that combines different options in one single frame-

The method uses the TTL field like Classic TTL, but the TTL of a packet may be decremented while it is stored in the epidemic buffer, depending on receive and send events. Formally, every packet in the epidemic buffer has an "age" field, which is a fixed decimal positive number less than 256. When a packet, created by some other node, is received by this node for the first time, its age is set to the complement to 255 of the received TTL: age= 255-TTL. When a packet is transmitted, its stored age is incremented by a fixed amount K 0 and then its TTL is set to 255age. When a duplicate packet is received, the received TTL is ignored but the stored age is incremented by K 1 : age = age+K 1 . When any packet is received, the stored age of all packets in the epidemic buffer is incremented by K 2 : age = age+K 2 . The node drops packets with age larger than 255.

Scheduler

Epidemic forwarding needs a scheduler for buffer management. To our knowledge, the only scheduler that is explicitly detailed in the literature is in [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF]. It is used with the virtual-rate based inhibition (Sect. 2.1.2.b). It decides which packet in the epidemic buffer is selected for transmission, i.e. to be passed to the MAC layer. In order to ensure some level of fairness, the scheduler serves packets per source Id, using a processor sharing approach. Moreover, every packet should be served at a rate not exceeding its virtual rate computed in Sect. 2.1.2.b.

Control of Injection Rate

The only explicitly defined method to achieve control of injection rate is the one in [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF]. It is used together with the aforementioned scheduler. The packets generated by the application at a given node are placed into the epidemic buffer, where they compete with the other packets for transmission (but with a larger virtual rate, having R = S =0 ). The application rate is controlled by a windowing system : The number of outstanding packets the application is allowed to have in the epidemic buffer at this node is limited to -at most-2 [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF]; a packet is deleted from the epidemic buffer when a duplicate is received, which serves as implicit acknowledgment (Ack).

Attacks

In this section, we describe the vulnerabilities that are specific to epidemic forwarding. We distinguish between two types of attackers: malicious and rational. The former does not look for a personal benefit but aims to harm other nodes. In contrast, the latter seeks to increase its personal profit from the network. Most of the attacks are described by drawing (Figs. 1 and 2) using a generic example where the attacker is M and the victim in malicious case is A.

Malicious Attacks

A malicious attacker aims at decreasing the spread and/or the injection rate of the victim by exploiting vulnerabilities in epidemic forwarding mechanisms. In the following we identify five attacks and map them to their corresponding epidemic forwarding mechanisms.

Send on behave

M is the malicious node and A is the victim. We refer by Self to packets generated at the node transmitting them, by Fwd to packets forwarded by the node but generated by others and by Fake to packets that are generated by M but carrying the victim identity (Id = A). We will explain only the "Inhibit by forwarding and TTL" attack, as other attacks have similar explanation. In "Inhibit by Forwarding and TTL", A sends a Self packet i with TTL = MaxTTL -1, that is received by M and B. M forwards the packet i (Fwd: i) 3 times with TTL = 0, the packet (Fwd: i) is received by B and C. Thus, the inhibition mechanism at B will inhibit packet i, as it is received 4 times, and C will drop the packet, as its TTL = 0.

Artificial High Density (AHD).

In this attack, we exploit the adaptability of the spread control to the conges-tion and node density. The attacker places itself close to the victim. It acts like any node: it has its self packets (packets that are generated at this node) to send and relays others packets. But it does not forward victim packets. By generating much traffic in the very close surrounding of the victim, the attacker incites the spread-control mechanism at the victim's good neighbors to react negatively and prevent the victim packets from going farther.

3.1.2. Inhibit by Forwarding (IbF) Attack. With IbF (Fig. 1), the attacker exploits the adaptive inhibition. It immediatly forwards the victim packets a number of times (this number is called Attack-Persistency) to inhibit its neighbrhood from forwarding the same packets (see Sect. 2.1.2). With the counter based inhibition (see Sect. 2.1.2.a), the Attack-Persistency is equal to the maximum value the counter can reach. With the virtual rate based inhibition, this Attack-Persistency should be large enough to make the corresponding virtual rate close to zero (in practice two times are enough).

3.1.3. Inhibit by TTL (IbTTL) Attack. This attack exploits the spread control using TTL. As the attacker receives a victim packet, it forwards it immediatly with a TTL = 0. In Fig. 1, B and M receives a packet from A with TTL = MaxTTL -1. M forwards it with TTL = 0 instead of (MaxTTL -2). Hence, the attacker decreases the chance the packet has to travel beyond C as B is inhibited and C drops the packet. Even if B succeeds in forwarding the packet after M, this will change nothing with C. This example considers the use of "classic TTL" (see Sect. 2.2.1). With "aging", the attack is exactly the same. Note that in Fig. 1, B applies the first strategy (see Sect. 2.2.1) upon receiving a duplicate of the packet and thus it keeps the old TTL.

3.1.4. Inhibit by Forwarding and TTL (IbFTTL) Attack. This is a combination of IbF and IbTTL (Fig. 1). In this case M forwards the victim packets Attack-Persistency times with TTL = 0 to insure that the victim packets at B are well inhibited and thus the packet loses any chance of travelling beyond C.

3.1.5. Send on Behalf of the Victim (SoB) Attack. The attacker exploits the scheduler and the aging mechanism. In Fig. 1, M sends fake packets with A's Id. As the scheduler serves packets per source Id, A's packets are delayed in the epidemic buffer and they will be dropped either by the aging mechanism (they become too old) or by buffer overflow.

Rational Attakcs

A rational attacker tries to increase its injection rate while maintaining large spread. In the following, we identify two rational attacks.

Do Not Cooperate (DNC) Attack.

When a new packet is injected by the application at a given node, it is placed in the epidemic buffer, where it competes with pack-ets received from other nodes. This competition prevents the application from injecting at the full rate allowed by the packet injection control mechanism because of the additional delay in the epidemic buffer. Thus, an attacker decides to not cooperate and to keep only its self packets (packets that are generated at this node) in the epidemic buffer. Note that, if the attacker tries to go beyond the allowed rate, its packets will be accumulated in other nodes, which are not able to serve them at the same rate. Thus, it risks killing its packets for the same reason as in Sect. 3.1.5.

Sybil Attack.

We refer to the Sybil attacker as the node that forges multiple identities [START_REF] Douceur | The sybil attack[END_REF]. This is a wellknown attack in networking, but the way it is exploited in this paper is new and very specific to epidemic forwarding. As the scheduler serves packets per source Id, the attacker sends its self packets with different Ids and thus it increases their share of the bandwidth. In Fig. 2, we present the scheduler as a process sharing approach where queues are per source Id. In this case, M's packets receive larger bandwidth share than A's packet at B.

Performance Evaluation

In this section, we evaluate the impact of the aforementioned attacks by simulation. We apply them in static scenarios, as well as in highly mobile networks. We consider vehicular mobility on the highway. Our metrics are based on the spread and the injection rate: a malicious attacker aims at reducing the victim spread and a rational one tries to increase its rate while maintaining large spread.

In our simulation we consider the epidemic forwarding system proposed in [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF], called SLEF. To our knowledge, SLEF is the only complete system proposed for a wide range of settings. Furthermore, SLEF implements all epidemic forwarding mechanisms already discussed in Sect. 2: The virtual rate based inhibition, spread control by TTL and aging, injection rate control and the scheduler discussed in Sect. 2.3. The parameter values of the virtual rate based inhibition are a = b =0 .15, c 0 corresponds to 802.11b basic rate (1Mbps) with a packet length of 1500 bytes. As for the aging, we use K 0 = K 1 =25and K 2 =0.5. The victim is in the middle, the attackers start filling the positions around the victim according to their index in an increasing order: if one attacker, it fills position 1. If 2 attackers, they fill positions 1 and 2, and so on. R can be either 25m or 100m. d = 25m.

Settings

Figure 3. Malicious attackers positions.

With the static scenarios, we simulate from 200 up to 600 nodes uniformly distributed over a square of 500x500 m 2 , but, in most cases, we show the results only for 400 nodes as the others are similar. The transmission range is around 50 m (PDA transmission range).

In the case of a malicious attack, the victim is in the middle of the square and attackers take place around it as it is indicated in Fig. 3. We want to evaluate the impact of the distance between attackers and the victim. Therefore, the radius R in Fig. 3 can have one of two values: 25m and 100m. With the former, the attackers of the corresponding circle are within the transmission range of the victim and they are outside it with the latter.

In the case of a rational attack, there exists only one attacker, which is in the middle.

The network can be either congested, where all nodes are sources sending at full rate (capacity allowed by the channel) or non-congested, where the victim is the only source in the network and it is sending at full rate. Beside the victim, only attackers can act as sources in the non-congested scenario, based on the attack they want to achieve.

In the following we will use the following notations: "close" [resp. "far"] to indicate that R is equal to 25m [resp. 100m] and "one" [resp. "all"] to indicate that the network is non-congested [resp. congested].

As for the mobile scenario, we simulate 1000 vehicles in an urban two-lane road. The speed limit is 80 km/h. The car density is 12.5 cars/km in each direction. The transmission range is 300m, which is typical for vehicular network.

Our simulations are carried out through JIST-SWANS [START_REF]Java in simulation time / scalable wireless ad hoc network simulator[END_REF], an open source simulator for ad hoc networks. The MAC layer is a very accurate implementation of 802.11b in DCF mode with the basic rate of 1 Mbps as we transmit in broadcast (pseudo-broadcast [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF]). As for the radio, we use the capture effect to approach the real WIFI cards that all implement it [START_REF] Kochut | Sniffing out the correct physical layer capture model in 802.11b, berlin, germany[END_REF]. We consider fading channels with free space path-loss. As for the mobile network, we use an extension of JIST-SWANS called STRAW [START_REF]Street random waypoint / vehicular mobility model for network simulations[END_REF], which simulates the vehicular traffic and provides a mobility model based on the operation of real vehicular traffic.

Static Scenarios

Malicious Attacks.

4.2.1.a. AHD:

The results are shown in Fig. 4. Let us begin with the "all" scenario (see Sect. 4.1) where the attackers are sources and act as any other node, except that they do not forward victim packets. In the "close" case, the impact of the attack is considerable in both scenarios, "all" and "one", and it increases with the number of attackers. In contrast, in the "far" + "all" scenario, the attackers do not have a major impact; the reason is twofold: (1) the attackers are far from the victim and thus they do not increase the density as much as in the "close" scenario; (2) the inhibition mechanism is adapted. Indeed, the attackers are numerous and they cooperate in forwarding all packets except the victims, hence they inhibit their neighbors from forwarding packets except those of the victim. Thus, the increase in the victim spread, which we notice for 8 attackers, is due to the fact that victim packets are less inhibited than others. If the inhibition were rigid, we expect that AHD would have more impact on the victim. In the "far" + "one" scenario, attackers are still injecting new packets in the network as before. They reduce considerably the victim spread.

4.2.1.b. IbF:

The attackers do not generate fresh packets: their role is merely to forward victim packets as it is explained in Sect. 3.1.2. The results are shown in Fig. 4. It is clear that IbF does not achieve its goal. This can be explained as follows: When an attacker receives a new victim packet, it immediately forwards it Attack-Persistency times (if the MAC layer allows). If the same attacker receives another victim packet, before it finishes forwarding the previous packet, it cancels the previous and it starts anew with the newest. Thus, let us consider a scenario that happens frequently. The victim sends a new packet. The attacker forwards it immediately. The victim receives a duplicate of its self packet and considers it as an implicit Ack. Hence, it injects a new self packet that will be received by the attacker before finishing the forwarding process and it will be received by other attackers even before beginning the forwarding process. Thus, all attackers cancel the previous packet, which explains why it is not inhibited.

4.2.1.c. IbTTL: Our implementation of IbTTL is similar to the one of IbF with the difference that it modifies the TTL before forwarding, as it is explained in Sect. 3.1.3. This attack is more harmful than IbF. The attacker needs to forward the packet only once with TTL = 0. Thus, nodes that receive the packet from the attacker for the first time are not able to forward it due to its TTL. This makes the difference with IbF, which needs to forward several times to inhibit the packet in its neighborhood.

4.2.1.d. IbFTTL:

This attack has approximately the same impact as IbTTL, which is to be expected as IbF has little effect on the victim.

4.2.1.e. SoB:

The attackers send only fake packets at full rate. Fig. 4 shows a significant decrease in spread and rate. The spread reduction is due to the fact that victim packets are killed in the epidemic buffers before being forwarded, which is due to the delay caused by the fake packets (for more explanation see Sect. 3.1.5). Moreover, the decrease in rate is due to the delay of the implicit Ack that controls the injection rate as explained in Sect. 2. [START_REF] Fawal | Self-Limiting Epidemic Forwarding[END_REF].

From what we have seen in this section, we can conclude that the attackers are not able to harm the victim in the presence of mobility for two reasons. The first is that the impact of the attackers is very position-dependent. The second is that, even with the most harmful attack, the attackers could reduce the spread of the victim, but its packets still reach a few tens of nodes. If these nodes are mobile, they will carry the victim packets beyond the barrier imposed by the attacker. This conclusion is well verified later in the vehicular network scenario.

Rational Attacks.

4.2.2.a. DNC:

We evaluate the impact of DNC only in the "all" scenario, where increasing the injection rate is a challenge. In the "one" scenario, the attacker is the only source in the network and he has the entire network capacity, thus it is meaningless to evaluate its impact in this case. In Fig. 6, the performance of a DNC attacker is compared with a well-behaved node that is very close to it and thus they both experience the same network conditions. We show the spread of both nodes and the rate ratio (DNC over well-behaved). The DNC rate is four times larger than a well-behaved node. But, surprisingly, the DNC spread is much larger when the network is very dense (600 nodes). The reason is as follows: The attacker does not forward others' packets. Thus, when it receives others' packets, it drops them without updating the age of its self packets in the epidemic buffer. Hence, the age of its self packets does not increase during their stay in its epidemic buffer by K 2 (see Sect. 2.2.2), which allows them to travel farther.

4.2.2.b. Sybil: We evaluate the impact of Sybil in only "all" scenario for the same reason as with DNC. The attacker uses five different identities. In addition, it does not forward others' packets. So, our implementation is in fact a combination of both attacks, Sybil and DNC, explained in Sect. 3. This implementation gives the attacker a much larger advantage than using DNC alone (up to 10 times larger than a well-behaved node and 2.5 larger than the DNC attacker), which explains the impact of Sybil alone. In this scenario, nodes are highly mobile and the position of the victim is not known. Thus, the attackers are chosen randomly. Beside the attackers, the network contains 1000 well-behaved nodes, all of them are sources ("all" scenario). All nodes cross the same urban road.

Vehicular Network Scenario

4.3.1. Malicious Attacks. Fig. 7(a) shows the impact of malicious attacks. The spread of the victim is drawn according to number of attackers. The Attack-Persistency of IbF and IbFTTL is 2. Other values give the same results. As we notice, the effect of the attackers is negligible even in the presence of 100 attackers. In the most harmful case, the IbFTTL attacker reduces the victim spread from 50 to 30 nodes, which is not significant. This can be explained by the presence of the spread control mechanism; the attacker can affect the victim only if their spreads interfere, i.e. there exist common nodes that receive the attacker and the victim packets. And the amount of harm is proportional to the amount of interference. As the spread is limited by the spread control mechanism, this interference is not considerable and does not happen frequently.

4.3.2. Rational Attacks. Contrary to malicious attacks, rational attacks are still powerful even in highly mobile network. The results are shown in Fig. 7(b). Sybil still ensures higher gain than DNC.

Conclusions

We identify vulnerabilities that are specific to epidemic forwarding over wireless ad-hoc networks. We classify these vulnerabilities into two categories: malicious and rational. We evaluate their impact according to the number of attackers and the different network settings. We find that the impact of malicious attacks depends on the position of the attacker relative to the victim, the network density, the traffic load and mobility. In static scenarios, we identify the attacks that reduce dramatically the victim spread, whereas the harm of other attacks is reduced due to the adaptive inhibition and the injection rate control. In highly mobile vehicular network, the impact of malicious attacks are minimized due to the spread control.

We have studied the rational case in presence of only one attacker in the network. The attacker could achieve considerable profit in all scenarios.

Our work can be extended in different directions. We plan to examine the impact of the presence of several rational attackers on the network. Another extension is to find solutions to recover from these vulnerabilities.

Figure 1 .Figure 2 .

 12 Figure 1. Malicious attacks.

Figure 5 .

 5 Figure 4. Malicious attacks in static scenario.

Figure 6 .

 6 Figure 6. Rational attacks in static scenario.

Figure 7 .

 7 Figure 7. Vehicular Network. Malicious (a) and Rational (b) attacks, "all" scenario.


Claire David

Problèmes inverses : Géométrie, Intégrabilité et formalisme de Lax

Introduction

Les travaux de C. Lacave et D. Iftimie ..., dans la lignée de ceux de J.Y. Chemin ..., permettent d'étudier l'influence d'obstacles immergés sur des champs de vitesses soit connus, soit solutions d'équations classiques régissant le mouvement d'un fluide incompressible (divergence nulle) : l'application de la loi de Biot et Savart permet en effet d'obtenir l'expression du champ de vitesses résultant : si l'on considère un domaine D de IR 2 , constitué de IR 2 privé d'un obstacle Ω et de sa frontière Ω, où se trouve un fluide incompressible dont le champ de vitesse U est régi par les équations d'Euler, alors, le champ résultant v est donné par :

v = -→ r ∧ curl U r 3 dV (1)
Dans ce type de problème, la forme "géométrique" de l'obstacle est connue.

Il serait intéressant de considérer le problème "à l'envers" : connaissant le champ de vitesses initial, et le champ de vitesses final, peut-on en déduire la forme de l'obstacle ? ? ? ? Il va sans dire que les applications pourraient être très intéressantes, ainsi, de "forcer la forme de l'obstacle" en fonction du champ de vitesses que l'on souhaite obtenir.

A cet effet, il serait intéressant de considérer, pour le champ initial, non pas un champ spécifiquement donné (ce qui, en général, n'est d'ailleurs pas le cas puisque celui est, le plus souvent, considéré en tant que solution d'un système d'équations, Euler ou Navier-Stokes), mais, bien en tant que solution d'un système donné, mais, non pas, classiquement, comme un système d'équations aux dérivées partielles, mais, plutôt, comme solution de "l'équation de Lax" correspondante, sous la forme

dL dt = [L, B] (2)
où L et B sont deux opérateurs différentiels, symétrique pour L, antisymétrique pour B.

Le système dont est solution le champ de départ, qui est un système intégrable, peut aussi se traduire par une condition de compatibilité entre deux problèmes linéaires.

De même : le champ que l'on souhaite obtenir, et qui est donc connu, peut être lui aussi considéré comme la solution d'un système intégrable, que l'on met aussi sous forme de Lax.

Pour ce type de problèmes, existent donc des surfaces remarquables (condition de courbure nulle notamment, tores invariants, etc ...), passer du champ initial au champ final via une approche "géométrico-différentielle" devrait donc donner des résulats intéressants.

Notamment, utiliser le fait que les trajectoires du système hamiltonien intégrable correspondant s'enroulent sur des tores.

Un début : couple de Lax pour les équations d'Euler 3D

Sous forme rotationnelle, les équations d'Euler tridimensionnelles s'écrivent :

où Ω désigne le rotationnel du champ u :

Un couple de Lax correspondant est alors donné par [START_REF] Li | Isospectral Theory of Euler Equations[END_REF] :

avec :


Carole Garnier

Wu Ke

Jean-Louis Dillenseger

K Wu

J.-L Dillenseger

Bladder Segmentation in MR Images using Active Region Growing Model

Prostate segmentation in MRI may be difficult at the interface with the bladder where the contrast is poor. Coupled-models that segment simultaneously both organs with non-overlapping constraints offer a good solution. As a presegmentation of the structures of interest is required, we propose in this paper a fast deformable model to segment the bladder. The combination of inflation and internal forces, locally adapted according to the gray levels, allow to deform the mesh toward the boundaries while overcoming the leakage issues that can occur at weak edges. The algorithm, evaluated on 33 MRI volumes from 5 different devices, has shown good performance providing a smooth and accurate surface.

I. INTRODUCTION

Prostate cancer is one of the most frequent cancer in men [START_REF]Cancer Facts and Figures[END_REF], [2]. Mainly treated by surgery or radiotherapy, other alternatives as brachytherapy or high intensity focused ultrasound (HIFU) are also proposed. Except for surgery, these treatments use preoperative imaging to target prostate and preserve sensitive organs in its surrounding as the bladder and the rectum. These organs of interest are delineated, most often manually, before or during the operation. This process is time consuming and is prone to intra-and inter-expert variabilities. Automation of the segmentation step would save significant time and improve the treatment accuracy.

Compared to CT-scan or ultrasound, MRI remains the imaging process that offers the best contrast between structures. As the boundary between organs may be weak, coupled-models techniques, proposed for CT device [START_REF] Costa | Automatic Segmentation of Bladder and Prostate Using Coupled 3D Deformable Models[END_REF], [START_REF] Song | Graph search with appearance and shape information for 3-D prostate and bladder segmentation[END_REF], that simultaneously segment the organs by applying non-overlapping constraints at the interface are promising. However, first approximations of the organs shapes and positions are needed. This paper focuses on the segmentation of the bladder. As the bladder appears with a high contrast compared to its surroundings, the region growing technique is used in [START_REF] Pasquier | Automatic Segmentation of Pelvic Structures From Magnetic Resonance Images for Prostate Cancer Radiotherapy[END_REF]. However, due in part to the partial volume effect, the bladder boundary may appear badly defined. The region growing method is not adapted to these cases leading to leakage issues. A level set approach is proposed in [START_REF] Duan | A Coupled Level Set Framework for Bladder Wall Segmentation with Application to MR Cystography[END_REF] to segment the bladder wall but this kind of method can be time consuming.

In this paper, we proposed a fast deformable model that solves the leakage issue. The algorithm combines an inflation force which acts like a region growing process and an internal force which constraints the shape of the surface. The method, with its initialization, its deformation and its resampling processes, is detailed in section II. In section III, the characteristics of the MRI volumes from different devices that have been used to evaluate the method are described. The results of the deformable model compared with those of the region growing method are then presented and discussed. Finally, we conclude in section IV.

II. METHOD

A. Initialization

Given an initial point C inside the bladder, an initial mesh centered at C is built. It is shaped as a small ellipsoid surface (Fig. 1) to roughly match the future bladder global morphology. Its radii, aligned along x, y and z directions, were initialized to 12, 6 and 6 millimeters respectively. Its vertices are organized along M a = 12 parallels and M b = 16 meridians. In the next sections, i = (a, b) will represent the index of a vertex located on the parallel a and the meridian b.

A bladder intensity prior I ref is defined to serve as a reference in the segmentation algorithm. This one is given as the average of the gray levels of the voxels included in the initial mesh.

B. Deformation

Once the initial mesh is built, the deformation acts like a Discrete Dynamic Contour [START_REF] Lobregt | A discrete dynamic contour model[END_REF], [START_REF] Hu | Prostate boundary segmentation from 3D ultrasound images[END_REF]. A total force f tot i is applied at each vertex i. Then, the dynamic equations derived from Newton's second law enable to calculate the acceleration a i , the velocity v i and the new position p i of each vertex i from the previous time step:

p i (t + △t) = p i (t) + v i (t)△t (1)
v i (t + △t) = v i (t) + a i (t)△t (2)
a i (t + △t) = 1 m i f tot i (t + △t) (3)
where m i is the mass at vertex i and △t the time step. This process is iterated while the deformation of the mesh is significant. Every 10 iterations, the vertex to vertex distances between the previous (at t -10△t) and the current mesh are computed. While vertices that have moved more than 2 mm exist, the deformation carries on.

In our study, f tot i corresponds to a weighted sum of inflation (f inf i ), internal (f int i ) and damping (f d i ) forces:

f tot i = w inf i f inf i + w int i f int i + f d i (4)
Inflation forces enlarge the mesh according to the intensity of the image at vertex positions. In order to constrain the mesh deformation, the inflation force at a vertex is not only based on the gray level of the voxel corresponding to its location but also on the gray levels of its neighboring vertices. Given a coefficient α j :

α j (x j , y j , z j ) = 1, if I(x j , y j , z j ) ≥ T surf 0, otherwise (5)
where I(x j , y j , z j ) is the intensity of the voxel (x j , y j , z j ) and T surf a threshold that characterizes the intensity of the bladder close to its wall, the inflation force is defined as:

f inf i =   1 N i + 1 N i j=0 α j   r i (6)
where r i is the unit radial vector at vertex i. j = 0 corresponds to the vertex i and j = 1, ..., N i to its N i neighbors.

Internal forces aim at smoothing the mesh during deformation. The shape of the mesh is thus constrained avoiding the segmented surface leakage at weak edges. These forces represent the surface tension at each vertex i:

f int i = - 1 N i Ni j=1 e ij e ij .r i .r i (7)
where e ij = p i -p j is the edge vector connecting the adjacent vextices with coordinates p i and p j . Damping forces ensure stability of dynamic behavior and prevent oscillations during deformation. They are proportional to the velocity:

f d i = w d i v i (8)
with w d i , a small negative constant. In order to better control the inflation, weights are also adjusted according to the image intensity.

w inf i are defined so that a vertex located inside the bladder moves quickly but slows down when it gets closer to the bladder wall:

w inf i =    w inf, in i , if I(x i , y i , z i ) ≥ T in w inf, surf i , if T in > I(x i , y i , z i ) ≥ T surf 0 otherwise. (9)
with T in a threshold that represents the gray levels of the inner part of the bladder and w inf, in i , w inf, surf i are positive weights that satisfy w inf, in i > w inf, surf i . Thus, considering this definition and those of the inflation force (6) the displacement of a vertex will be more important if the number of its neighbors that satisfy the threshold criteria is large and its corresponding gray level is intense. On the contrary, a vertex whose intensity is weak and/or whose neighbors are already stopped will not keep on moving. This process also contributes to solve leakage issues.

w int i and w d i are adapted to stop the deformation of the surface in the areas where the vertices enter the bladder wall:

w int i = w int, surf i , if I(x i , y i , z i ) ≥ T surf 0 otherwise. ( 10
)
w d i = w d, surf i , if I(x i , y i , z i ) ≥ T surf -m i otherwise. (11)
Indeed, the total force of a vertex whose intensity is less than T surf will be defined as:

f tot i = -m i v i , (12)
and, according to (2), the velocity at the next iteration will become nil.

Thresholds and weights values will be deduced from a set of various qualities images and given in section III.

C. Resampling

The initial resolution of the ellipsoidal mesh may become insufficient depending on the size of the bladder. Indeed, the resolution must be fine enough to capture the shape variations. But, according to the definition of the inflation force [START_REF] Duan | A Coupled Level Set Framework for Bladder Wall Segmentation with Application to MR Cystography[END_REF], the area covered by a vertex and its neighbors must be large enough so the neighboring vertices can stop the current one in case of a local bad definition of the boundary.

Let d a b,b+1 denotes the edge vector oriented from the vertex (a, b) to the vertex (a, b + 1). Given a parallel a, the average distance da between its vertices is given by:

da ← 1 M b M b -1 b=0 d a b,b+1 . (13)
The average distance db along a meridian b is computed in the same way. A good compromise in the resolution of the mesh was found by adapted regularly, when necessary, the number of parallels and/or meridians so that da=Ma/2 and db=0 keep equal on average to 5 mm.

Because of the deformation of the mesh, the distance between vertices will change. The resolution may become locally smaller or bigger depending on the areas (Fig. 2 (a)). As the bladder may present local variations, the regularity of the sampling must be preserved to better delineate the organ. A process was thus applied to keep the distances between vertices homogeneous successively along parallels and meridians.

First, the positions of the vertices are moved along the edge vectors of parallels a so that the distances between them are equal to da in the way described on Fig. 3. Then, the accelerations and the velocities at new positions are obtained using a bilinear interpolation from the previous vertices locations. The same process is then applied along the meridians from the new vertices positions. The resulting mesh after resampling is illustrated on Fig. 2 (b).

III. MATERIALS AND RESULTS

The algorithm was tested on 33 T2 axial MRI volumes from 5 different devices. The MRI systems and volumes characteristics are summarized in Tab. III. MRI 3D images were interpolated along z before segmenting the bladder.

This work comes within the context of a segmentation with a coupled-model technique. Before segmenting the bladder, the prostate shape and position is approximately known. Thus, in this study, the initial point C inside the bladder is located automatically using geometrical (location of the bladder with respect to those of the prostate) and image (high intensity the bladder) features.

The parameters of the deformable model were adjusted from the 33 MRI volumes. Given the different qualities of the MRI systems, the bladder wall definition varies from one volume to another. Also, the position of the bladder at the beginning of the volume and the thickness of the axial slices lead to strong partial volume effects. Thus, the bladder intensity can not be considered as really homogeneous. The weights and thresholds values defined in section II were adapted to take into account this heterogeneity. The chosen parameters values give good results on the whole set of data. The intensity of the bladder decrease from its center to its boundary. The threshold T in (9) characterizes the intensity of the inner part of the bladder except the region close to its wall. It must be high while taking into account the inhomogeneity of the gray levels. It was thus set to 0.8I ref .

A value of T surf = 0.6I ref (Eqs. 5, 9-11) that represents the nearness of the bladder wall produces good segmentations.

A voxel whose gray level is higher than T in have a great likelihood to belong to the bladder and can move quickly. w inf,in i = 1 and w inf,surf i = 0.25 (9) allow to deform the mesh fast inside the bladder and to slow down progressively near its boundary. The internal forces must have a weight in the same level of magnitude but slightly larger than the one of the inflation force. Indeed, if an inflation force attempts to pass through the bladder wall because of a local bad definition of the boundary, the internal force must be strong enough to stop it. A value of 1.2 was thus chosen for w int,surf i (10). The weight w d,surf i (11) which is less sensitive [START_REF] Ladak | Prostate boundary segmentation from 2D ultrasound images[END_REF] was set to -0.6.

The mass and the time step of the dynamic process (Eqs. 1)-(3) were set to m i = 1, △t = 1 for simplicity.

Five examples of the resulting segmentations obtained with the proposed deformable model are displayed Fig. 4 with axial (X-Y), sagittal (Y-Z) and coronal (X-Z) sections (from left to right). White arrows indicate some areas where a region growing method would leak. Figs. 5 and6 illustrate this issue on the same volumes than Figs. 4(a) and 4(b) respectively. The three slices on the top of Fig. 5 were obtained with a threshold equal to I ref . We observe that this value is insufficient to segment the whole bladder due to its heterogeneity. Decreasing the threshold to T in (bottom of Fig. 5) helps to delineate the organ but also extracts a part of the surroundings and the inside of the prostate (indicated on Fig. 4(a) with a white arrow). Fig. 6 shows a result obtained with a threshold of I ref . Once again, only the half of the bladder is segmented and, even with this high threshold, a part of the surroundings is already included in the segmentation result.

The deformable model gives good results for all these critical cases. Even if in Fig. 4(a) a little part of the bladder is missing due to is thinness, the boundary close to the prostate is always well delineated. Besides, no leakage issue occurs in spite of the weakness of the bladder wall in some areas and of the similarity between the bladder and the background intensities. A good example is given on Fig. 4(d) where the is really thin and badly defined on a large area at the left bottom of the bladder. Even in this case, the surface stopped at the boundary.

Concerning the computational time, the segmentation process lasts 1 second on average and 3 second to the maximum on a 2.33 GHz Xeon PC.

IV. CONCLUSION

A fast and efficient 3D deformable model was proposed in this paper to segment the bladder in T2 MRI volumes. The algorithm was evaluated on 33 patients from different MRI devices and has shown good performance for all specimens. Without any pre-or post-processing, the method produces a smooth segmented surface. Besides, the adaptive definitions of forces and weights according to the intensity of the image help to solve the leakage issues that can occur with a region growing technique. The resulting surface can then be included in a coupled bladder-prostate segmentation method in order to constrain and improve the segmentation at the interface between the two organs.

Fig. 1 .

 1 Fig. 1. Initialization of the algorithm: ellipsoidal mesh built from an initial point.

Fig. 2 .

 2 Fig. 2. Resampling process at the 50 th iteration: (a) before resampling, (b) after resampling.

Fig. 3 .

 3 Fig. 3. Resampling of a parallel a. The vertices are moved along the edge vectors to be spaced of da.

Fig. 4 .

 4 Fig. 4. Five examples of segmentation results obtained with the proposed deformable model. From left to right: axial (X-Y), sagittal (Y-Z) and coronal (X-Z) sections. White arrows indicate some areas where a region growing method would leak.

Fig. 5 .

 5 Fig. 5. Segmentation obtained with a region growing method on the same volume than Fig. 4(a). The threshold is set to I ref at the top and T in at the bottom. From left to right: axial (X-Y), sagittal (Y-Z) and coronal (X-Z) sections. Leakage can be observed, notably in the area indicated with a white arrow in Fig. 4(a).

Fig. 6 .

 6 Fig. 6.Segmentation obtained with a region growing method with a threshold of I ref on the same volume than Fig.4(b). From left to right: axial (X-Y), and sagittal (Y-Z) sections.

TABLE I MRI

 I SYSTEMS AND VOLUMES CHARACTERISTICS. n: NUMBER OF VOLUMES, sa: AXIAL PIXEL SIZE, st: TRANSVERSE PIXEL SIZE.

	MRI System	n	Size	sa	st
			(voxel)	(mm)	(mm)
	3T	10	720 × 720 × 20	0.416	4.001
	Philips Achieva
	3T		320 × 320 × 20	0.625	3.600
	Siemsens Vero	5
	1.5T Siemens		256 × 256 × 24	0.781	3.000
	Magnetom Symphony 16
	3T	1	352 × 352 × 24	0.540	2.630
	Philips Achieva
	GE Signa HDxt	1	512 × 512 × 16	0.391	4.500

ACKNOWLEDGMENTS The authors would like to thank E. Blanc and C. Vurpillot from EDAP-TMS, Pr. O. Rouvière from the Edouard Herriot hospital of Lyon and R. Mathieu from the Pontchaillou hospital of Rennes for providing the datasets.

This work is part of the French MULTIP project supported by an ANR Grant (ANR-09-TECS-011-06).


Redamy Pérez-Ramos
email: redamy.perez@uv.es

Vincent Mathieu
email: vincent.mathieu@ific.uv.es

Miguel-Angel Sanchis-Lozano
email: miguel.angel.sanchis@ific.uv.es

Three-particle correlations in QCD parton showers

Three-particle correlations in quark and gluon jets are computed for the first time in perturbative QCD. We give results in the double logarithmic approximation and the modified leading logarithmic approximation. In both resummation schemes, we use the formalism of the generating functional and solve the evolution equations analytically from the steepest descent evaluation of the one-particle distribution. We thus provide a further test of the local parton hadron duality and make predictions for the LHC.

Introduction

The observation of quark and gluon jets has played a crucial role in establishing QCD as the theory of strong interaction within the standard model of particle physics. The jets, narrowly collimated bundles of hadrons, reflect configurations of quarks and gluons at short distances. Powerful schemes, like the double logarithmic approximation (DLA) and the modified leading logarithmic approximation (MLLA), which allow for the perturbative resummation of soft-collinear and hard-collinear gluons before the hadronization occurs, have been developed over the past 30 years (for a review see [START_REF] Yu | Basics of Perturbative QCD[END_REF]). One of the most striking predictions of perturbative QCD, which follows as a consequence of angular ordering (AO) within the MLLA and the local parton hadron duality (LPHD) hypothesis [START_REF] Ya | [END_REF], is the existence of the hump-backed shape [START_REF] Yu | Basics of Perturbative QCD[END_REF] of the inclusive energy distribution of hadrons, later confirmed by experiments at colliders. Indeed, the shape and normalization of single inclusive distributions are compared with an experiment; a constant K ch , which normalizes the number of soft gluons to the number of charged detected hadrons (mostly pions and kaons), turns out to be close to unity (K ch ∼ 1), giving support to the similarity between parton and hadron spectra [START_REF] Yu | Basics of Perturbative QCD[END_REF]. Thus, the study of inclusive observables like the inclusive energy distribution and the transverse momentum k ⊥ spectra of hadrons [3] has shown that the perturbative stage of the process, which evolves from the hard scale or leading parton virtuality Q ∼ E to the hadronization scale Q 0 , is dominant. In particular, these issues suggest that the hadronization stage of the QCD cascade plays a subleading role and, therefore, that the LPHD hypothesis is successful while treating one-particle inclusive observables.

The study of particle correlations in intrajet cascades, which are less inclusive observables, focuses on providing tests of the partonic dynamics and the LPHD. In [4], this observable was computed for the first time at small x (energy fraction of the jet virtuality taken away by one parton) in MLLA for particles staying close to the maximum of the one-particle distribution. In [5], the previous solutions were extended, at MLLA, to all possible values of x by exactly solving the QCD evolution equations. This observable was measured by the OPAL Collaboration in e + e -annihilation at the Z 0 peak, that is, for √ s = 91.2 GeV at LEP [6]. Though the agreement with predictions presented in [5] was improved, a discrepancy still subsists pointing out a possible failure of the LPHD for less inclusive observables. However, these measurements were redone by the CDF Collaboration in pp collisions at the Tevatron for mixed samples of quark and gluon jets [7]. The agreement with predictions presented in [4] turned out to be rather good, especially for particles having very close energy fractions (x 1 ≈ x 2 ). A discrepancy between the OPAL and CDF analysis showed up and still stays unclear. Therefore, the measurement of the two-particle correlations at higher energies at the LHC becomes crucial. Furthermore, going one step beyond, in this article we give predictions for the three-particle correlations inside quark and gluon jets. This observable and the two-particle correlations can be measured on equal footing at the LHC so as to provide further verifications of the LPHD for less inclusive observables.

Kinematics and evolution equations

A generating functional Z(E, Θ; {u}) can be constructed [START_REF] Yu | Basics of Perturbative QCD[END_REF] that describes the azimuth averaged parton content of a jet of energy E with a given opening half-angle Θ; by virtue of the exact AO (MLLA), which satisfies an integro-differential system of evolution equations. In order to obtain exclusive n- 1: Three-particle yield and angular ordering inside a high energy jet.

x x 1 x 2 3 Θ 0 Θ 1 Θ 2 Θ3 E (1-z)E A zE C B Figure

particle distributions D (n)

A (k i , E) one takes n variational derivatives of Z A over u(k i ) with appropriate particle momenta, i = 1 . . . n, and sets u ≡ 0 afterwards; inclusive distributions are generated by taking variational derivatives around u ≡ 1. Let us introduce the n-particle differential correlations for

A = G, Q, Q jets as, A (n) 1...n (z) ≡ x 1 z . . . x n z D (n) A x 1 z . . . x n z , ln zQ Q 0 , (1)
together with

A (n) 1...n ≡ A (n) 1...n (1
) for later use; x i corresponds to the Feynman energy fraction of the jet taken away by one particle "i" and z is the energy fraction of the intermediate parton. For instance, for three-particle correlations n = 3, the observable to be measured reads C

(3)

A 123 = A (3) 123 A 1 A 2 A 3 .
The production of three hadrons is displayed in Fig. 1 after a quark or a gluon (A) jet of energy E with half opening angle Θ 0 and virtuality Q = EΘ 0 has been produced in a high energy collision. The kinematical variable characterizing the process is given by the transverse momentum

k ⊥ = zEΘ 1 ≥ Q 0 [or (1 - z)EΘ 1 ≥ Q 0 ]
of the first splitting A → BC. The parton C fragments into three offspring such that three hadrons of energy fractions x 1 , x 2 , and x 3 can be triggered from the same cascade following the condition Θ 0 ≥ Θ 1 ≥ Θ 2 ≥ Θ 3 , which arises from exact AO in MLLA [START_REF] Yu | Basics of Perturbative QCD[END_REF]. We make use of variables, ℓ = ln z x 1 , y = ln

x 3 EΘ 1 Q 0 , ℓ i = ln 1 x i , y j = ln x j EΘ 0 Q 0 , η ij = ln x i x j , Y = ℓ i + y j + η ij = ln(Q/Q 0 ) and λ = ln(Q 0 /Λ QCD ).
The two variables entering the evolution equations are z and Θ 1 , such that

x 1 ≤ z ≤ 1 ⇒ 0 ≤ ℓ ≤ ℓ 1 .
Accordingly, the anomalous dimension related to the coupling constant can be parametrized as

γ 2 0 (Q 2 ) = 2N c α s (Q 2 ) π , γ 2 0 (ℓ + y) = 1 β 0 (ℓ + y + η ij + λ)
,

where β 0 = 1 4Nc 11 3 N c -4 3 n f T R , with T R = 1/2
and n f the number of light quark flavors. From AO and the initial condition at threshold

x 3 EΘ 0 ≥ x 3 EΘ 1 ≥ x 3 EΘ 3 ≥ Q 0 , one has the bounds Q 0 x 3 E ≤ Θ 1 ≤ Θ 0 ⇒ 0 ≤ y ≤ y 3
for the integrated evolution equations. The evolution equations satisfied by [START_REF] Yu | Basics of Perturbative QCD[END_REF] are derived from the MLLA master equation for the generating functional Z A (E, Θ; u(k i )). For three-particle correlations, one takes the first δZ A δu(k 1 ) , second

δ 2 Z A δu(k 1 )δu(k 2 )
, and finally third

δ 3 Z A δu(k 1 )...δu(k 3 )
functional derivatives of Z A (E, Θ; u(k i )) over the probing functions u(k i ) so as to obtain the differential system of evolution equations:

Q(3) ℓy = C F N c γ 2 0 G (3) - 3 4 C F N c γ 2 0 G (3) ℓ -β 0 γ 2 0 G (3) , (2)
Ĝ(3) ℓy = γ 2 0 G (3) -aγ 2 0 G (3) ℓ -β 0 γ 2 0 G (3) + (a -b)γ 2 0 Ĝ(2) 12 G 3 + Ĝ(2) 13 G 2 + Ĝ(2) 23 G 1 ℓ (3) -β 0 γ 2 0 Ĝ(2) 12 G 3 + Ĝ(2) 13 G 2 + Ĝ(2) 23 G 1 + (a-c)γ 2 0 (G 1 G 2 G 3 ) ℓ -β 0 γ 2 0 G 1 G 2 G 3 , where Â(2) ij = A (2) ij -A i A j and Â(3) = A (3) -A 1 A 2 A 3 - Â(2) 12 A 3 - Â(2) 13 A 2 - Â(2) 23 A 1 .
The subscripts ℓ and y in Eqs. ( 2) and (3) denote ∂/∂ℓ and ∂/∂y, respectively. The first terms of Eqs. ( 2) and ( 3) are of classical origin and, therefore, universal. Corrections ∝ - 3 4 , a, (a -b), and (a -c), which are O( √ α s ) suppressed, better account for energy conservation at each vertex of the splitting process, as compared with the DLA O(1). The hard constants are obtained after integration over the regular part of the DGLAP splitting functions [START_REF] Yu | Basics of Perturbative QCD[END_REF] as performed in [4,5]. In the equation for the gluon initiated jet (3), the first and second constants a(n f = 3) = 0.935 and b(n f = 3) = 0.915 were obtained in the frame of the single inclusive distribution and two-particle correlations respectively [4]. The third constant c(n f ) appearing for the first time in this frame reads

c(n f ) = 1 4N c 11 3 N c + 4 3 n f T R 1 -2 C F N c 3 n f =3
= 0.917.

MLLA and DLA solutions of the evolution equations

Equation ( 3) is self-contained and can be solved iteratively by setting

G (3) = C (3)
G 123 G 1 G 2 G 3 and G (2) ij = C (2)
G ij G i G j in the left-and right-hand sides of (3). Accordingly, the solution of ( 2) is also obtained by setting

Q (3) = C (3) Q 123 Q 1 Q 2 Q 3 and Q (2) ij = C (2) Q ij Q i Q j in the left-hand side of (2) and G (3) = C (3) G 123 G 1 G 2 G 3 in
the right-hand side of the same equation such that the iterative solutions can be written in the compact form

C (3) A 123 = C (2) A 12 -1 F (2) A 12 + C (2) A 13 -1 F (2) A 13 + C (2) A 23 -1 F (2) A 23 + N 2 c C 2 A F (3) A 123 . (4)
The MLLA two-particle correlators C

(2)

A 12 will be taken from [5] for the computation of C

(3)

A 123 . Moreover, F (2)
G ij = 1 + 1 -bΨ ℓ + ξ ij 1 -ǫ 1 2 + ∆ 12 + ∆ 13 + ∆ 23 + ǫ 1 , (5)
F (3) G 123 = 1 -cΨ ℓ + ξ 12 1 + ξ 13 1 + ξ 23 1 -ǫ 1 2 + ∆ 12 + ∆ 13 + ∆ 23 + ǫ 1 (6)
and for the quark jet

F ( 2
)
Q ij = 1 + ξij 1 -ǫ1 3 + ∆ 12 + ∆ 13 + ∆ 23 -aΨ ℓ + ǫ1 , (7)

F

(3)

Q 123 = C (3) G 123 (1 -aΨ ℓ ) + ξ12 1 + ξ13 1 + ξ23 1 -ǫ1 3 + ∆ 12 + ∆ 13 + ∆ 23 -aΨ ℓ + ǫ1 , (8)
where

Ψ ℓ = ψ 1,ℓ + ψ 2,ℓ + ψ 3,ℓ = O(γ 0 ) and ψ = ln[G(ℓ, y)].
Higher order corrections arising from the solution of the system of Eqs. 2 and 3 have been neglected in (4). In this case, G(ℓ, y) is the inclusive energy distribution, which will be inserted from the steepest descent method presented in [5]. The other functions appearing in ( 5) and ( 6) are ∆ ij = γ -2 0 (ψ i,ℓ ψ j,y + ψ i,y ψ j,ℓ ) = O(1) and

ζ ℓ = Ċ(3) G 123 ,ℓ Ċ(3) G 123 = O(γ 2 0 ), ζ y = Ċ(3) G 123 ,y Ċ(3) G 123 = O(γ 2 0 ), χ ij ℓ = Ċ(2) G ij ,ℓ Ċ (2)
G ij = O(γ 2 0 ), χ ij y = Ċ(2) G ij ,y Ċ (2)
G ij = O(γ 2 0 ), ξ ij 1 = γ -2 0 χ ij ℓ Ψ y + χ ij y Ψ ℓ = O(γ 0 ), ǫ 1 = γ -2 0 (ζ ℓ Ψ y + ζ y Ψ ℓ ) = O(γ 0 ), with ζ = ln Ċ(3) G 123 and χ = ln Ċ(2)
G . The set of functions appearing in ( 7) and ( 8) is obtained from the previous by replacing

ζ → ζ, χ → χ, ξ → ξ, Ċ(2) G ij → Ċ(2) Q ij and Ċ(3) G ij → Ċ(3) Q ij where the dotted Ċ(2)
A ij and Ċ(3) A ij are the DLA solutions of the two-and three-particle correlators; that is why this solution is said to be iterative. Moreover, corrections ǫ 1 , ǫ1 and ξ ij 1 , ξij 1 are very small and do not play a significant role in the shape and normalization of the three-particle correlations.

The DLA two-particle correlators are taken from [8] and the DLA expression for Ċ(3) A ij can be obtained from ( 4) by setting all MLLA O(γ 0 ) corrections to zero:

Ċ(2) A ij -1 = N c C A 1 1 + ∆ ij ; (9)
Ċ(3) A 123 -1 - Ċ(2) A 12 -1 - Ċ(2) A 13 -1 - Ċ(2) A 23 -1 (10)
= N c C A Ċ(2) A 12 -1 + Ċ(2) A 13 -1 + Ċ(2) A 23 -1 2 + ∆ 12 + ∆ 13 + ∆ 23 + N 2 c C 2 A 1 2 + ∆ 12 + ∆ 13 + ∆ 23 .
The solutions have the following simple physical interpretation: the first term (= -1) in the left-hand side translates the independent or decorrelated emission of three hadrons in the shower. After inserting the two-particle correlator with color factor ∝ Nc C A (9) in the left-hand side of (10), terms ∝ Nc C A correspond to the case where two partons are correlated inside the same subjet, while the other one is emitted independently from the rest. Next, replacing (9) in the right-hand side of (10), one obtains a contribution

∝ N 2 c C 2 A
describing the independent emission of two partons inside the same subjet. The last term ∝ N 2 c C 2 A involves three particles strongly correlated inside the same partonic shower as depicted in Fig. 1. This term is indeed the cumulants of genuine correlations, first obtained in this article for this observable.

The evaluation of (4), which is expressed in terms of the logarithmic derivatives of the single inclusive distribution ln[G(ℓ, y)], will be performed using the steepest descent method to determine G(ℓ, y) [5,8]. Thus, the MLLA logarithmic derivatives were written in [5] in the form:

ψ i,ℓ (µ i , ν i ) = γ 0 e µ i + 1 2 aγ 2 0 e µ i Q(µ i , ν i ) -tanh ν i -tanh ν i coth µ i 1 + e µ i Q(µ i , ν i ) (11) - 1 2 β 0 γ 2 0 1 + tanh ν i 1 + K(µ i , ν i ) + C(µ i , ν i ) 1 + e µ i Q(µ i , ν i ) + O(γ 2 0 ), ψ i,y (µ, ν) = γ 0 e -µ i - 1 2 aγ 2 0 2 + e -µ i Q(µ i , ν i ) + tanh ν i -tanh ν i coth µ i 1 + e -µ i Q(µ i , ν i ) - 1 2 β 0 γ 2 0 1 + tanh ν i 1 + K(µ i , ν i ) -C(µ i , ν i ) 1 + e -µ i Q(µ i , ν i ) + O(γ 2 0 ), (12)
where the functions Q(µ i , ν i ), C(µ i , ν i ) and K(µ i , ν i ) are defined in [5] and (µ i , ν i ) are expressed as functions of the original variables (ℓ i , y j ) by inverting the nonlinear system of equations [8]:

y i -ℓ i ℓ i + y i = (sinh 2µ i -2µ i ) -(sinh 2ν i -2ν i ) 2(sinh 2 µ i -sinh 2 ν i ) , sinh ν i √ λ = sinh µ i √ ℓ i + y i + λ .
In particular, this method allows for the estimation of the observable for particles with energies near the maximum or hump (ℓ max = Y /2) of the one-particle distribution | ℓ -Y /2 |≪ σ ∝ Y 3/2 , which applied to the three-particle correlations will appear in a forthcoming paper. For instance, at DLA one has ∆ ij = 2 cosh(µ i -µ j ) with such a parametrization of the logarithmic derivatives of the inclusive spectrum. Close to the hump one has ∆ ij ≃ (ℓ i -ℓ j ) 2 ; thus the correlations are expected to be quadratic as a function of (ℓ i -ℓ j ) and to have a maximum for particles with the same energy x i = x j . In this frame, the role of MLLA corrections should be expected to be larger than for the two-particle correlations. Indeed, higher order corrections increase with the rank of the correlator, which is known from the Koba-Nielsen-Olesen problem for intrajet multiplicity fluctuations [9]. For the two-particle correlations, for instance, one has ∝ -b(ψ 1,ℓ + ψ 2,ℓ ) and for the three-particle correlator one has the larger correction ∝ -c(ψ 1,ℓ + ψ 2,ℓ + ψ 3,ℓ ).

Phenomenology and comparison with existing e + e -and p p data

The study of n-particle correlations is very important because, being defined as the n-particle cross section normalized by the product of the single inclusive distribution of each parton

C (n) A 1...n = A (n) 1...n A 1 . . . A n ,
the resulting observable becomes independent of the constant K ch , thus providing a refined test of QCD dynamics at the parton level. Since our study of three-particle correlations depends on previous results for two-particle correlations, we briefly review recent results about this observable. The MLLA evolution equations for two-particle correlations, quite similar to those leading to the hump-backed plateau, were solved iteratively in terms of the logarithmic derivatives of G(ℓ, y) [5]. That is how, the result previously obtained by Fong and Webber in [4], only valid in the vicinity of the maximum ℓ max of the distribution, was extended to all possible values of x. Consequently, as displayed in Fig. 2, the normalization of the more accurate solution of the evolution equations is lower and reproduces some features of the OPAL data at the Z 0 peak Q = 91.2 GeV of the e + e -annihilation, like the flattening of the slopes towards smaller values of x [5]. Qualitatively, our MLLA expectations agree better with available OPAL data than the Fong-Webber predictions [5]. There remains however a significant discrepancy, markedly at very small x. In this region nonperturbative effects are likely to be more pronounced. They may undermine the applicability to particle correlations of the LPHD considerations that were successful in translating parton level predictions to hadronic observations in the case of more inclusive single particle energy spectra [START_REF] Yu | Basics of Perturbative QCD[END_REF].

These measurements were redone by the CDF Collaboration for pp collisions at center of mass energy √ s = 1.96 TeV for mixed samples of quark and gluon jets [7]. For comparison with CDF data, the two-particle correlator was normalized by the corresponding multiplicity correlator of the second rank, (2) Q [5] in the process e + e -→ q q as a function of

ℓ 1 + ℓ 2 = | ln(x 1 x 2 )| for ℓ 1 -ℓ 2 = ln(x 2 /x 1 ) = 1.0 (left) and ℓ 1 -ℓ 2 = | ln(x 2 /x 1 )| for ℓ 1 + ℓ 2 = ln(x 1 x 2 ) = 6.0 (right).
which defines the dispersion of the mean average multiplicity inside the jet. In this case, the MLLA solution by Fong and Webber [4], the more accurate MLLA solution [5], and the NMLLA solution [3] were compared with the CDF data. The Fong-Webber predictions turned out to be in good agreement with CDF data in a range from large to small x, also covering the region of the phase space where MLLA predictions should normally not be reliable, that is, for x > 0.1 (see Fig. 3). As these figures were taken from [7], different notations have been used in this case, for instance,

ℓ ≡ ξ = ln(1/x), ∆ξ = ξ -ξ max (ξ max ≡ ℓ max = 1 2 ln(Q/Q 0 )) such that ∆ξ 1 + ∆ξ 2 = ℓ 1 + ℓ 2 -ln(Q/Q 0 ) and ∆ξ 1 -∆ξ 2 = ℓ 1 -ℓ 2 . 2 ξ ∆ + 1 ξ ∆ -2 -1.5 -1 -0.5 0 0.5 1 1. 5 2 ) 2 ξ ∆ , 1 ξ ∆ C( 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 2 ξ
+ ∆ξ 2 = | ln(x 1 x 2 )| -ln(Q/Q 0 ) for ∆ξ 1 = ∆ξ 2 (left) and ∆ξ 1 -∆ξ 2 = | ln(x 2 /x 1 )| for ∆ξ 1 = -∆ξ 2 (right).
As observed in Fig. 3 (left), the data are well described by the three cases in the interval ∆ξ 1 + ∆ξ 2 > -0.5, that is, at very small x. However, the Fong and Webber's solution also describes the data for ∆ξ 1 + ∆ξ 2 < -0.5, that is, for larger values of x where the MLLA is no longer valid. QCD color coherence for Fig. 3 (left, the peak at ∆ξ 1 + ∆ξ 2 = -1.5 is due to numerical uncertainties) should be observed if the analysis is extended to ∆ξ 1 + ∆ξ 2 > 2.5. Moreover, the NMLLA solution [3] extends, like for the k ⊥ spectra, the region of applicability of such predictions for larger values of x. In [7], it was concluded that despite the disagreement with the OPAL data in Fig. 2, the LPHD stays successful for the description of less inclusive energy-momentum correlations. Therefore, in this paper we encourage the analysis of these observables by other collaborations like ALICE, ATLAS, and CMS in order to clarify this mismatch.

Predictions for three-particle correlations and phenomenology

Finally, in order to extend the applicability of the LPHD to a larger domain of observables, we perform theoretical predictions for three-particle correlations in the limiting spectrum approximation (Q 0 ≈ Λ QCD ). This observable and two-particle correlations can be measured on equal footing at the LHC.

We display the MLLA solutions (4) of the evolution equations ( 2) and (3). The correlators are functions of the variables ℓ i , y i and the virtuality of the jet Q = EΘ 0 . After setting y i = Y -ℓ i with fixed Y = ln(Q/Q 0 ) in the arguments of the solutions (4), the dependence can be reduced to the following: (virtuality Q = 450 GeV and Λ QCD = 250 MeV), which is realistic for LHC phenomenology [5]. The representative values ℓ 3 = ln(1/x 3 ) = 4.5, 5.5 (x 3 = 0.011, x 3 = 0.004) have been chosen according to the range of the energy fraction x i ≪ 0.1, where the MLLA scheme can only be applied.

C Q 1 Q 2 Q 3 are displayed,
In Fig. 5, the DLA (10) and MLLA (4) three-particle correlators for A = G and A = Q, Q are depicted, in this case as a function of the sum (ℓ 1 + ℓ 2 ) = | ln(x 1 x 2 )| for the same values of ℓ 3 = ln(1/x 3 ) = 4.5, 5.5, for x 1 = x 2 and Y = 7.5. As expected in both cases, the DLA and MLLA three-particle correlators are larger inside a quark than in a gluon jet. Of course, these plots will be the same and the interpretation will apply to all possible permutations of three particles (123). As remarked above, the difference between the DLA and MLLA results is quite important in pointing out that overall corrections in O( √ α s ) are large. Indeed, the last behavior is not surprising as it was already observed in the treatment of multiplicity fluctuations of the third kind given by R 3 = 4.52 1 -(2.280 -0.018n f ) √ α s [10].

For instance, for one quark jet produced at the Z 0 peak of the e + e -annihilation (Q = 45.6 GeV), one has α s = 0.134. Replacing this value into the previous formula for the quark jet multiplicity correlator, one obtains a variation from 4.52 (DLA) to 0.83 (MLLA). Because of this, DLA has been known to provide unreliable predictions which should not be compared with experiments. From Fig. 4, the correlation is observed to be the strongest when particles have the same energy and to decrease when one parton is harder than the others. Indeed, in this region of the phase space two competing constraints should be satisfied: as a consequence of gluon coherence and AO, gluon emission angles should decrease and on the other hand, the convergence of the perturbative series k ⊥ = x i EΘ i ≥ Q 0 should be guaranteed. That is why, as the collinear cutoff parameter Q 0 is reached, gluons are emitted at larger angles and destructive interferences with previous emissions occur. This effect is clearly observed in Fig. 4; the steep fall of the distribution is more pronounced in the quark jet than in the gluon jet. Moreover, the observable increases for softer partons with x 3 decreasing, which is for partons less sensitive to the energy balance. In Fig. 5 the MLLA correlations increase for softer partons, then flatten and decrease as a consequence of soft gluon coherence, reproducing for three-particle correlations the hump-backed shape of the one-particle distribution. Because of the limitation of phase space, one has C (3) ≤ 1 for harder partons.

Summary

In this paper we provide the first full perturbative QCD treatment of three-particle correlations in parton showers, provide a further test of the LPHD within the limiting spectrum approximation, and briefly revise the comparison of two-particle correlations with OPAL and CDF data. The correlations have been shown to be strongest for the softest hadrons having the same energy x 1 = x 2 = x 3 in both quark and gluon jets, increasing as a function of ln(x i /x j ) and | ln(x i x j )| when x k softens, that is for partons less sensitive to the energy balance. This result becomes therefore universal for n-particle correlations.

Coherence effects appear when one or two of the partons involved in the process are harder than the others, thus reproducing for this observable the hump-backed shape of the one-particle distribution. Also, the two-and three-particle correlations vanish (C (2) → 1) when one of the partons becomes very soft, thus describing the hump-backed shape of the one-particle distribution. The reason for that is dynamical rather than kinematical: radiation of a soft gluon occurs at large angles, which makes the radiation coherent and thus insensitive to the internal parton structure of the jet ensemble.

Figure 2 :

 2 Figure 2: Two-particle correlations in two quark jets R = 1 2 + 1 2 C

Figure 3 :

 3 Figure 3: Two-particle correlations in a mixed sample of gluon and quark jets in pp collisions as a function of ∆ξ 1 + ∆ξ 2 = | ln(x 1 x 2 )| -ln(Q/Q 0 ) for ∆ξ 1 = ∆ξ 2 (left) and ∆ξ 1 -∆ξ 2 = | ln(x 2 /x 1 )| for ∆ξ 1 = -∆ξ 2 (right).

Figure 4 :

 4 Figure 4: Gluon jet 3-particle correlator as a function of | ln(x 1 x 2 )| for x 1 = x 2 and ln(1/x 3 ) (left) and as a function of ln(x 2 /x 1 ) for fixed | ln(x 1 x 2 )| and ln(1/x 3 ) (right).

  respectively, as a function of the difference (ℓ 1 -ℓ 2 ) = ln(x 2 /x 1 ) for two fixed values of ℓ 3 = ln(1/x 3 ) = 4.5, 5.5, fixed sum (ℓ 1 + ℓ 2 ) = | ln(x 1 x 2 )| = 10, and, finally, fixed Y = 7.5

5 Figure 5 :

 55 Figure 5: Quark jet 3-particle correlator as a function of | ln(x 1 x 2 )| for x 1 = x 2 and ln(1/x 3 ) (left) and as a function of ln(x 2 /x 1 ) for fixed | ln(x 1 x 2 )| and ln(1/x 3 ) (right).

We give the first analytical predictions of this observable in view of forthcoming measurements by AT-LAS, CMS, and ALICE at the LHC. Further information from the comparison with forthcoming data may also help to improve Monte Carlo event generators in the soft region of the phase space in intrajet cascades, where PYTHIA, ARIADNE and HERWIG face difficulties while reproducing the data [11].

We gratefully acknowledge enlightening discussions with W. Ochs and E. Sarkisyan-Grinbaum as well as support from Generalitat Valenciana under Grant No. PROMETEO/2008/004. M.A.S. acknowledges support from FPA2008-02878 and GVPROMETEO2010-056. V.M acklowledges support from the grant HadronPhysics2, a FP7-Integrating Activities and Infrastructure Program of the European Commission under Grant No. 227431, by UE (Feder), and by the MICINN (Spain) Grant No. FPA 2010-21750-C02-01.