David Galdeano

Ahmed Chemori
email: chemori@lirmm.fr

Sébastien Krut
email: krut@lirmm.fr

Philippe Fraisse
email: fraisse@lirmm.fr

Task-based whole-body control of humanoid robots with ZMP regulation, real-time application to a squat-like motion

In this paper, a task based whole-body control strategy is proposed for humanoid robots. Its basic idea lies in the control of (i) the CoM with a ZMP regulation, (ii) the relative pose of robot's feet and (iii) joint's limit avoidance. Through the proposed study, it is shown that these tasks allow to produce smooth whole-body motions. Real-time experimentation results are presented to show the effectiveness of the proposed control scheme.

I. INTRODUCTION

Since the early days of the humanoid robotics, most of studies has been focused on two field: (i) the study of biped walking control and stability issues targeting the lower-limbs considering the upper-limbs as a dummy mass [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF], (ii) the manipulation of items or other tasks involving upper-limbs motions only without making a step [START_REF] Asfour | Toward humanoid manipulation in human-centred environments[END_REF].

This separation has allowed a simpler study of humanoid control but trying to combine upper and lower-limbs motion together is not an easy process if these are designed separately since stability is not ensured [START_REF] Nakaoka | Generating whole body motions for a biped humanoid robot from captured human dances[END_REF]. To tackle this issue and improve the efficiency of control, the humanoid robot has to be considered as a whole system.

In the literature, several approaches have been proposed to deal with whole-body control; most of them can be classified into three classes: optimization based control, human-capture based control and task-based control.

Optimization based control relies on accurate model of the robot and his environment to compute off-line the best joints' trajectories. Several studies [START_REF] Lengagne | Guaranteed computation of constraints for safe path planning[END_REF], [START_REF] Salini | Lqp controller design for generic whole body motion[END_REF] use this methodology to produce whole-body motions. However, this method has several drawbacks: (i) some optimization can take a very long time to compute [START_REF] Lengagne | Planning and fast replanning safe motions for humanoid robots[END_REF], (ii) accurate model of the environment are difficult to produce and (iii) this approach is not reactive hence difficult to use in a dynamic environment.

Human capture based control relies in using sensors to record human motion data. Then these last ones are used to generate human-like motions such as balancing [START_REF] Montecillo-Puente | On realtime whole-body human to humanoid motion transfer[END_REF], walking [START_REF] Kim | Stable whole-body motion generation for humanoid robots to imitate human motions[END_REF], or dancing [START_REF] Montecillo-Puente | On realtime whole-body human to humanoid motion transfer[END_REF], [START_REF] Kim | Stable whole-body motion generation for humanoid robots to imitate human motions[END_REF]. However, despite evident similarities, the human locomotor system is much more complex than the humanoid's one. Therefore, a direct mapping from human captured motion to humanoid robot will fail in most cases and often lead to instable motions. This implies that data need to be adapted to take into account the specific structures of humanoid robots such as joint limits [START_REF] Pollard | Adapting human motion for the control of a humanoid robot[END_REF].

Task-based control uses the principle of task as defined by Nakamura [START_REF] Nakamura | Task-priority based redundancy control of robot manipulators[END_REF] and Siciliano [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF] to track different objectives in operational space using a joint space control. Applications to humanoid robot can be found in Mansard, Ramos and Sentis works [START_REF] Ramos | Dynamic motion capture and edition using a stack of tasks[END_REF], [START_REF] Mansard | A versatile generalized inverted kinematics implementation for collaborative working humanoid robots: The stack of tasks[END_REF], [START_REF] Sentis | A whole-body control framework for humanoids operating in human environments[END_REF]. It is worth to note that Montecillo [START_REF] Montecillo-Puente | On realtime whole-body human to humanoid motion transfer[END_REF] had based his work on both human capture data and tasks hierarchy controls.

In this paper, a whole-body control strategy is proposed. The basic idea of the proposed control strategy is to use the position of the Center of Mass (CoM) and the relative pose of the feet to produce stable motions. By using Siciliano's recursive formulation of task formalism [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF], we can track these two quantities. We added (i) a zero moment point (ZMP) stabilizer to improve the stability of the robot and (ii) a joints' limit avoidance to deal with singularities issue.

Every approach proposed in the literature make a decomposition of walking motion in different phases. The usual ones are simple support, double support and impact phases. Our approach allows to avoid this phenomenon. We can create a continuous whole-body control without having to switch from phases like double support to single support when a foot is lifted. The dynamic stability is ensured by the ZMP regulator.

This paper is organized as follows: in next section, our demonstrator HOAP-3 is introduced. Section III is devoted to the proposed control scheme, where its basic principle is presented and discussed. Real-time experimental results are presented in section IV, with a presentation and discussion of the obtained results. The paper ends with some concluding remarks and future work.

II. DESCRIPTION AND MODELING OF HOAP-3

In this study, we propose to use the humanoid robot HOAP-3 from Fujitsu company as a testbed for the illustration of the proposed control solution (cf. Fig. 1-(a)).

This robot is a small humanoid of 60 cm tall, 8.8 kg weight, and has 28 degrees of freedom (dof) distributed as follows: Each leg is composed of six dof, three dof at the hip, one dof at the knee and two dof at the ankle. The torso is composed of 1 dof and the neck is composed of 3 dof. Both arms are composed of 6 dof, 3 dof at the shoulder, one dof at the elbow and two dof at the hand. The spatial distribution of the different dofs is illustrated in Fig. 1-(b). The robot is equipped with optical incremental encoders at all joints, a 3-axis accelerometer and 3-axis gyroscope used for posture sensing, four force sensors per foot for ground contact forces' measurement and two cameras for visual feedback.

The actuators are brushless motors paired with microcontrollers, which are controlled by an embedded PC running with a real time kernel RT-Linux at a sample frequency of 1 kHz.

III. PROPOSED CONTROL SCHEME

As previously mentioned in the introduction, the challenge to perform a whole-body control law is still an open problem. The use of task formalism to control a robot for tracking several objectives in the operational space has been proposed by Nakamura in [START_REF] Nakamura | Task-priority based redundancy control of robot manipulators[END_REF] and Siciliano in [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]. The task formalism has been used on humanoid robots by Sentis and Khatib in [START_REF] Sentis | A whole-body control framework for humanoids operating in human environments[END_REF] for multi-contact dynamic motions. Recently, Mansard [START_REF] Mansard | A versatile generalized inverted kinematics implementation for collaborative working humanoid robots: The stack of tasks[END_REF], [START_REF] Mansard | A unified approach to integrate unilateral constraints in the stack of tasks[END_REF] has used this approach on humanoid robots and has generalized the formalism by using the addition and removal of tasks during the control execution. He has also provide a flexible C++ implementation of the proposed algorithm.

The actual study is also based on the task formalism. The main contributions of this paper is a whole-body control architecture including a ZMP regulation to ensure stability. A first study has been proposed in [START_REF] Galdeano | Partial human data in design of human-like walking control in humanoid robotics[END_REF] to deal with whole-body control, but without any experimental results.

In the literature, several tasks are needed to produce stable whole-body motions. However, in this work, the proposed architecture is focused on only three main tasks: the relative feet position and orientation tracking, CoM position tracking with ZMP regulation and the joints' limits avoidance.

It is worth to note that additional tasks such as manipulation or Human-Robot Interaction (HRI) tasks can easily be added within this framework.

A. First task: Relative feet pose

The relative feet pose task allows to place a foot relatively to the other one (cf. Fig. 2). It is based on the tracking error of the relative feet pose:

ε r = [E T pos E T ori ] T (1)
where E pos ∈ R 3×1 is the position error and E ori ∈ R 3×1 is the orientation error.

The position error is defined by:

E pos = P rd -P r (2)
where P rd = [x rd y rd z rd ] T expresses the desired feet relative position, designed to produce a stable motion and P r = P rf -P lf is the estimated feet relative position where P rf = [x rf y rf z rf ] T and P lf = [x lf y lf z lf ] T are respectively the positions of the right and left foot obtained using the forward kinematic model (F KM r ).

The orientation error is defined by:

E ori = R rf (ln(R -1 rf R lf R rd )) ∨ (3)
where R rd ∈ R 3×3 is the desired feet relative orientation matrix, R rf ∈ R 3×3 and R lf ∈ R 3×3 are respectively the right and left foot orientation matrix obtained using F KM r . The logarithmic map operator ω = (lnR) ∨ is defined by: (5) where I ∈ R 3×3 is the identity matrix. More details on this operator can be found in [START_REF] Kajita | Introduction à la commande des robots humanoides[END_REF] and [START_REF] Li | A mathematical introduction to robotic manipulation[END_REF].

(lnR) ∨ =        [0 0 0] T if R = I θ 2 sin θ   r 32 -r 23 r 13 -r 31 r 21 -r 12   if R = I (4
The orientation error is computed form the rotation matrix (R ∈ R 3×3 ) to have a well defined error without gimbal lock limitations which appears with Euler angles notation. These orientation matrix needs only to be defined in a common reference coordinates. Here, we take the reference foot orientation as a common reference coordinates.

The Jacobian of the relative feet pose task is then defined as follows:

ε r = J r ε q (6)
where ε q ∈ R 28×1 represents the vector of articular position errors of the robot, J r ∈ R 6×28 is the Jacobian matrix of the relative translation between the feet.

B. Second task: Center of mass position

The Center of Mass (CoM) position tracking allows to improve the stability in the static case. This is a common task in the task formalism applied to humanoid robots [START_REF] Montecillo-Puente | On realtime whole-body human to humanoid motion transfer[END_REF], [START_REF] Cotton | On the manipulability of the center of mass of humanoid robots: Application to design[END_REF], [START_REF] Sentis | Compliant control of multicontact and center-of-mass behaviors in humanoid robots[END_REF] since it is a usually used stability indicator. By tracking a CoM designed to achieve stability, we can produce stable motions.

The tracking error of the CoM is expressed by:

ε CoM = d -CoM (7)
where Both of these two positions are expressed in the reference foot coordinates. This allows to design a cyclic trajectory for walking without having to take care about the position of the robot in the global coordinate system. The position of the robot in the global coordinate system can be estimated using odometry method based on foot position.

CoM d = [x CoM d y CoM d z CoM d ] T is
The jacobian of the CoM task is then defined as follows:

ε CoM = J CoM ε q (8)
where J CoM ∈ R 3×28 is the CoM Jacobian matrix.

C. ZMP regulation

Now, to improve dynamic stability of the robot, we propose to consider a ZMP regulation. This regulation is obtained by injecting the ZMP error in the CoM tracking task using a Proportional Derivative (PD) controller.

Let us defined the tracking error on the ZMP position, ε ZM P , as follows:

ε ZMP = α dZMP lef t + (1 -α) dZMP right (9)
α = AZMP lef t AZMP lef t + AZMP right [START_REF] Nakamura | Task-priority based redundancy control of robot manipulators[END_REF] dZMP lef t = ZMP lef t desired -ZMP lef t measured [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF] dZMP right = ZMP P right desired -ZMP right measured [START_REF] Ramos | Dynamic motion capture and edition using a stack of tasks[END_REF] where dZMP lef t and dZMP right are illustrated in Fig. 3. AZMP lef t and AZMP right are the amplitudes of forces In equation ( 9), the ZMP position tracking error ε ZMP is a weighted distribution of the ZMP error on both feet. The weight for this distribution is based on the repartition of force amplitude on both feet as defined in equation [START_REF] Nakamura | Task-priority based redundancy control of robot manipulators[END_REF]. This allows to neglect the error on one foot when it is lifted off for walking scenarios (for instance) and reinforce the stabilization on the support foot.

The tracking error of the ZMP: ε ZM P , is then injected in the CoM task using a Proportional Derivative (PD) controller.

ε CoM &ZM P = ε CoM + K p ε ZMP + K d dε ZMP dt (13)
The effectiveness of this part of the controller is demonstrated in section IV-A.

D. Third task: Joints' limits avoidance

Using only the two defined previous tasks control can lead the robot to non desired behavior because of drift on joint positions while tracking perfectly these tasks. If the drift leads the robot joints to a singular pose or to joints' articular limits, the induced loss of dof can alter the efficiency of tasks tracking. To overcome this problem, it will be necessary to add some constraints on the control.

The proposed solution to tackle this problem is based on attractive potential fields who allow to defined a comfort position (cf. Fig. 4).

The attractive potential fields are defined as:

ε Joints = β [ε q1 ε q2 ε q3 ... ε q22 ] T ( 14
)
Fig. 4: Robot's comfort position is defined with joints' configuration far from their limits.

with:

ε qi = 2 (q i -q imed ) (q imax -q imin ) 2 ; q imed = q imax + q imin 2 ( 15
)
where β ∈ R * + is a convergence gain, q imax and q imin are respectively the upper and lower limit on the joint i, q i is the current joint i position and q imed is its comfort position.

This specific task must be added after all operational space tasks because it affect all the degrees of freedom. The null space used to add tasks is filled meaning no other task can be tracked.

E. Architecture of the control scheme

In order to control the robot, we have designed the control scheme summarized in Fig. 5. The task computation algorithm is based on Siciliano's recursive formulation detailed in [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]. The feet relative pose is the task with the highest priority since an error in foot placement often lead to fall of the robot. The control law ε q1 with this task only is given by the following:

ε q1 = J + r ε r (16)
The CoM task with ZMP regulation has the second priority since small error on the CoM position are tolerable. Indeed, these errors can often be corrected without a loss of equilibrium. The control with the two tasks becomes ε q2 given by:

ε q2 = ε q1 + ( J CoM ) + (ε CoM -J CoM ε q1 ) (17)
with:

J CoM = J CoM P r and P r = (I -J + r J r ) (18)
Other tasks like items manipulations and upper body position and/or orientation can be added now without compromising the motion's stability using Siciliano's recursive formulation [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF].

The last task to be considered is the joints' limits avoidance because it provides a good way to avoid the loss of dof due to singularity. The main inconvenient is that the control's nullspace is filled by this task, meaning that no other task can be added after this one. The final control law becomes the following:

ε q3 = ε q2 + (P Joints ) + (ε Joints -ε q2 ) (19)
with:

P Joints = (P r -J + CoM J CoM ) (20)
This final control law allows to create continuous wholebody motions while ensuring dynamic stability, without having to switch between phases like double support and single support in dynamic walking tasks.

IV. REAL-TIME EXPERIMENTAL RESULTS

The proposed control scheme (cf. Fig. 5) has been developed and implemented in the HOAP-3 robot (presented in section II) with reference signals obtained by numerical computation. These reference signals depends on the two scenarios to be performed.

In the first one, the reference signals given for the relative feet pose are constant, which means that we want the feet to be static, and the reference CoM position is set to the initial robot's CoM position and kept constant which means that the robot should not move excepted for ZMP regulation.

In the second one, the reference signals given for the relative feet pose are constant, which means that we want the feet to be also static, and the reference CoM position to be moved up and down using a sinus signal to generate a squat-like motion. The ZMP regulation ensure that the produced motions are stable.

A. First scenario: static standing under disturbing forces

The objective of this scenario is to prove the effectiveness of the ZMP regulation presented in section III-C.

We deliberately push the robot while controlling the wholebody to keep constant feet relative pose and a fixed CoM (cf. Fig. 6).

The only motion from the robot is due to the ZMP regulation inside CoM task.

The ZMP measurement resulting from this scenario are displayed in Fig. 7. The objective of this scenario is to produce a whole-body control of the HOAP-3 by moving the CoM position up and down according to a sinusoidal signal (cf. Fig. 8). In Fig. 9, the joints' trajectories are displayed. It is worth to note that the trajectories converge to a periodic motion since the desired tasks are periodic, some joints are not solicited since not all of them are useful for this whole-body motion. In Fig. 10, W lef t and W right are the amplitudes of the measured force for the left and right foot respectively, x and y frontal and lateral positions of the ZMP relative to the center of each foot.

[rad] q2 [rad] q3 [rad] q4 [rad] q5 [rad] q6 [rad] q7 [rad] q8 [rad] q9 [rad] q10 [rad] q11 [rad] q12 [rad] q13 [rad] q14 [rad] q15 [rad] q16 [rad] q17 [rad] q18 [rad] q19 [rad] q20 [rad] q21 [rad]
According to the obtained results, one can observe that the produced motions are continuous and smooth, without switching phases, and the robot produces the desired squatlike motion (cf. Fig. 8).

V. CONCLUSION AND FUTURE WORK

This paper deals with a control scheme for humanoid robots. The proposed control solution is based on three tasks controlling (i) the relative distance between feet of the robot, (ii) the trajectory of its CoM with a ZMP regulation and (iii) the avoidance of joints' limits. One of the main advantages of such proposed control scheme lies in that fact that it provides a continuous control framework for wholebody motions. Indeed, the proposed principle don't use any decomposition of the movement in different phases, which avoid eventual discontinuities. The obtained results are very promising.

In future work, we aim at using different scenarios such as walking or following several objectives including for instance carrying or manipulating objects while walking.

Fig. 1 :

 1 Fig. 1: The HOAP-3 humanoid robot used as demonstrator for real-time experiments.

Fig. 2 :

 2 Fig. 2: Illustration of the tracking error on the relative feet pose.

11 r

 11 12 r 13 r 21 r 22 r 23 r 31 r 32 r 33   and θ = cos -1 r 11 + r 22 + r 33 -1 2

  the desired position of the center of mass and CoM = [x CoM y CoM z CoM ] T is the real position of the center of mass obtained from the forward kinematic model of center of mass (F KM CoM ).

Fig. 3 :

 3 Fig. 3: Graphical illustration of tracking errors on ZMP position.

Fig. 5 :

 5 Fig. 5: Block diagram of the proposed control scheme with the control of (i) the relative pose of robot's feet, (ii) the CoM with a ZMP regulation and (iii) joint's limits avoidance.

Fig. 6 :

 6 Fig. 6: Static standing under disturbing forces.

  x right [mm] y right [mm] W left [N] x left [mm] y left [mm]

Fig. 7 :Fig. 8 :

 78 Fig. 7: ZMP measurements. In Fig. 7, W lef t and W right are the amplitudes of the measured forces for respectively the left and right foot. x and

Fig. 9 :

 9 Fig. 9: Joints trajectories.

Fig. 10 :

 10 Fig. 10: ZMP measurements.

ACKNOWLEDGMENT

This research was supported by the French National Research Agency, within the project R2A2 (ANR-09-SEGI-011).


J Rosato
email: joel.rosato@univ-amu.fr

H Capes

R Stamm

Influence of correlated collisions on Stark-broadened lines in plasmas

Keywords: numbers: 32.70.Jz, 32.60.+i, 52.25.Dg, 52.70.Kz

An investigation of spectral line broadening in plasmas is carried out within a kinetic-theory approach, based on the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy. The model employs a resummation procedure to account for correlated emitter-perturber collisions. Applications to hydrogen lines indicate that such collisions strongly affect the width and the shape in the core region. This argument is supported by comparisons to numerical simulations. It is also shown that the usual collision operator models, based on a binary description of emitter-perturber collisions, can be extremely inaccurate. The present model, in a better agreement with numerical simulations, is suggested as an extension suitable for the design of fast and accurate numerical routines for plasma diagnostics.

I. INTRODUCTION

The Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy provides a valuable approach to the study of dilute medium such as gases or plasmas. The typical scheme for obtaining a kinetic equation consists in writing a closure relation based on suitable assumptions for the two-particle correlation function, and inserting this relation in the first equation of the hierarchy. This approach has been successfully used in the past for describing transport effects in plasmas (e.g., [START_REF] Balescu | Statistical Dynamics -Matter Out of Equilibrium[END_REF] for a review). In this work, we consider an adaptation of the BBGKY formalism to the description of Stark line shapes for plasma spectroscopy purposes and we focus on an extension devoted to retain correlated collisions between an emitter and the surrounding perturbers. Plasma spectroscopy is of general interest for diagnostics, given that the light naturally emitted from the atomic species (including multi-charged ions) contains information on the plasma parameters (N , T e , etc.) [START_REF] Lister | [END_REF][3][4][5]. Detailed line shapes are also of interest for radiative transfer simulations, where they serve in transport codes as cross sections for the photon absorption and emission processes [6,7]. It has been a long-standing challenge for plasma spectroscopy to accurately describe Stark-broadened line shapes [START_REF] Griem | Principles of Plasma Spectroscopy[END_REF].

The problem consists in a correct description of the atomic dipole autocorrelation function in the presence of the plasma microscopic electric field. The latter is created by numerous charged particles, so that one is confronted with a statistical time-dependent quantum problem which, even today, has no general analytical solution. A basic description of a line shape, sometimes referred to as the "standard model", consists in assuming the ions motionless during the dipole correlation time (or time "of interest") τ dd and, at the same time, assuming the electrons moving rapidly. This short-time assumption for the electrons allows one to describe their contribution through a non-Hermitian part in the Hamiltonian, usually referred to as "collision operator". This quantity can be calculated using standard perturbation approaches, either semi-classical [START_REF] Griem | [END_REF] or fully quantum-mechanical [10], using a binary atom-electron interaction model (impact approximation). The ionic contribution is described for its part with a constant Stark effect term in the Hamiltonian, and a statistical average over the ionic electric field is then performed using an appropriate probability density function. As a rule, deviations to the standard model are expected if the ordering τ e ≪ τ dd ≪ τ i is not satisfied. Here, τ e and τ i denote the electron and ion collision times, estimated as r 0 /v e , r 0 /v i with v e , v i being the electron and ion thermal velocities, respec-tively, and r 0 = N -1/3 being the characteristic scale for the mean interparticle distance. Ion dynamics effects (that arise, for instance, if τ i ∼ τ dd ) have been observed [11] and have been extensively investigated, using numerical simulations [4,[12][13][14][START_REF] Stambulchik | High Energy Density Phys[END_REF] or ad hoc models based on the statistical properties of the electric field (such as the model microfield method "MMM" [START_REF] Brissaud | [END_REF]17] or the frequency fluctuation model "FFM" [18][19][20]). Incomplete electron collisions (that are expected when τ e ≪ τ dd is not satisfied) have also been investigated, in particular by using refined models for the collision operator, either based on kinetic theory (such as the "unified theory" [21,22]) or semiempirical procedures (e.g. using the Lewis cutoff [23]).

Such models reproduce the result of the impact approximation at the line center (i.e., with the frequency detuning ∆ω = 0) and yield an asymptotic behavior identical to that expected within the quasi-static approximation. Several works in the past have been devoted to generalizations of the collision operator models, able to account for non-impact effects, e.g. [24][25][26]. A major issue concerns the modeling of correlated collisions between the emitter and the charged particles surrounding it. An interpretation of the correlated collisions is that an emitter suffering a binary collision also "feels" the presence of the other particles.

They are important when the characteristic collision frequency becomes of the same order as or larger than the inverse correlation time of the emitter-perturber interaction potential [27,28]. The latter is of the order of the inverse plasma frequency ω -1 p , and the collision frequency can be estimated as N b 2 W v ln(λ D /b W ), where N is the particle density, v is the thermal velocity of the perturbers, λ D is the Debye length, and b W = hn 2 /m e v is the socalled Weisskopf radius (here n is the principal quantum number of the upper level). It determines an effective cross section which corresponds to collisions yielding coherence loss of the atomic wavefunction (strong collisions). The usual binary models assume that there are no simultaneous strong collisions, i.e. that the parameter g = N b 3 W is small compared to unity. In general, correlated collisions concern both strong and weak collisions and, hence, may be important even for small g. In this work, we address the problem of modeling correlated collision effects on the line broadening with a first-principles approach inspired from the unified theory. The theoretical development makes extensive use of kinetic equations of BBGKY-type and follows early works reported on in the seventies [27,28]. A preliminary attempt to retain correlated collisions has been performed by using a generalization of the Kirkwood truncature hypothesis to quantum operator [29]. Here, we follow a more rigorous approach based on a partial resummation of the correlation terms that arise in the BBGKY hierarchy. This procedure amounts to a renormalization of the atomic energy levels involved in the collision operator, accounting for their perturbation due to collisions. The method is designed for regimes where g < ∼ 1 and is suitable either for electron or ion collisions. We introduce the general formalism in Sec. II and, for the sake of simplicity, we specify it to one-component plasma in the weak coupling limit. This assumption allows one to consider the perturbers as non-interacting quasi-particles evolving on straight lines (ideal gas) and generating a Debye electric field. The resummation procedure is presented in Sec. III. We show that the collision operator is determined by a nonlinear equation and we examine an iterative scheme for solving it numerically. In order to examine the role of correlated collisions, we perform calculations in ideal cases and compare the results to the unified theory (Sec. IV). Comparisons with numerical simulations are also made. Finally, Sec. V reports on applications to realistic cases in the framework of magnetic fusion research.

II. FORMALISM

We give here a brief overview of the BBGKY approach to the line broadening problem, following the early works reported on in [21,27,28]. An atom immersed in a set of N charged particles is considered. For simplicity, we assume one-component plasma and we consider the particles evolving along straight lines and generating a Debye electric field. A generalization accounting for correlations between the perturbers in a rigorous way can be developed following [27]. A Stark line shape I(ω) is written as

I(ω) = 1 π Re ∑ ε ⟨⟨d • ε| Φ0 (-iω)|ρd • ε⟩⟩. (1)
Here, the double ket notation for Liouville space |...⟩⟩ has been used, ρ is the atomic density operator, d • ε is the dipole projected onto the polarization vector ε, and Φ0 (-iω) =

∫ +∞ 0 dte iωt Φ 0 (t)
Here, 1...N stand for the phase space coordinates (r 1 , v 1 ...r N , v N ) of the perturbers, d1...dN are the corresponding volume elements, f N is the N -particle phase space distribution, and

   ∂ ∂t + iL 0 + N ∑ j=1 [ v j • ∂ ∂r j + iV (j) ]    U = 0, (3)
with the initial condition U (t = 0) = 1. All the quantities present in Eq. ( 3) are operators acting in the atomic Liouville space formed by the dyadics |ab⟩⟩ ≡ |a⟩|b⟩. The term L 0 is the Liouvillian accounting for the atomic energy level structure and

V (j) = V (r j ) = -d • E(r j )
denotes the Stark term resulting from the electric field due to the j-th perturber.

Equation ( 2) provides generalizations of the reduced phase space distributions, which account for the presence of the atom. It is customary to introduce a cluster expansion (t is not written explicitly)

                         Φ 1 (1) = f 1 (1)Φ 0 + Γ 1 (1) Φ 2 (1, 2) = f 2 (1, 2)Φ 0 + f 1 (1)Γ 1 (2) + f 1 (2)Γ 1 (1) + Γ 2 (1, 2) Φ 3 (1, 2, 3) = f 3 (1, 2, 3)Φ 0 + f 2 (1, 2)Γ 1 (3) + f 2 (1, 3)Γ 1 (2) + f 2 (2, 3)Γ 1 (1) +f 1 (1)Γ 2 (2, 3) + f 1 (2)Γ 2 (1, 3) + f 1 (3)Γ 2 (1, 2) + Γ 3 (1, 2, 3) ... (4)
This decomposition allows one to write a hierarchy of equations for the correlations

   ∂ ∂t + iL 0 + p ∑ j=1 [ v j • ∂ ∂r j + iV (j) ]    Γ p (1...p) = -i p ∑ j=1 f 1 (j)V (j)Γ p-1 (1...j -1, j + 1...p) -iN ∫ d(p + 1)V (p + 1)Γ p+1 (1...p + 1), (5)
with the initial condition Γ p (1...p; t = 0) = 0 for p ≥ 1. It has been assumed that the phase space distributions are space independent and factorize as f p (1...p) = f 1 (1)...f 1 (p). Also, the thermodynamic limit (N → ∞, V → ∞, N /V = cst with V being the volume of the system) is assumed. In Eq. ( 5), by convention Γ 0 ≡ Φ 0 and Γ -1 ≡ 0.

Originally, the unified theory was developed for applications in cases where the collisions are uncorrelated. The treatment involves equations for Φ 0 and Γ 1 only, setting Γ 2 ≡ 0:

( ∂ ∂t + iL 0 ) Φ 0 = -iN ∫ d1V (1)Γ 1 (1), (6)
and

[ ∂ ∂t + iL 0 + v 1 • ∂ ∂r 1 + iV (1) ] Γ 1 (1) = -if 1 (1)V (1)Φ 0 . ( 7
)
Equation ( 7) can be solved formally by using the propagator of the atom under the influence of one perturber, Q(1, t). It obeys the following equation

[ ∂ ∂t + iL 0 + iV (r 1 + v 1 t) ] Q(1, t) = 0, (8)
with the initial condition Q(1, t = 0) = 1, and it is proportional to a time-ordered exponential (Dyson series)

Q(1, t) = e -iL 0 t T exp [ -i ∫ t 0 dτ e iL 0 τ V (r 1 + v 1 τ )e -iL 0 τ ] , (9)
T being the time-ordering operator. The correlation Γ 1 is obtained from a convolution

Γ 1 (1; t) = -if 1 (1) ∫ t 0 dτ Q(r 1 -v 1 τ, v 1 , τ )V (r 1 -v 1 τ )Φ 0 (t -τ ). ( 10
)
Inserting the solution Eq. ( 10) into the right-hand side of Eq. ( 6) provides a closed, integrodifferential equation for the averaged evolution operator Φ 0 (t). Its solution takes a simple form in the Laplace space:

Φ0 (s) = [s + iL 0 + K(s)] -1 . ( 11
)
Inserting this relation in Eq. ( 1) provides an analytical expression for the line shape. Decomposing the double bra and kets onto an appropriate base of the Liouville space indicates that the line shape can be written as a sum of generalized Lorentzian functions, whose widths are frequency-dependent and given by matrix elements of K(s = -iω). The latter quantity is a collision operator that accounts for incomplete collisions. It is completely determined in terms of the propagator Q and the interaction term V :

K(s) = N ∫ ∞ 0 dte -st ∫ d1V (r 1 + v 1 t)Q(1, t)V (1)f 1 (1). ( 12
)

III. CORRELATED COLLISIONS

The assumption Γ 2 ≡ 0 considered in the unified theory is not valid in regimes where correlated collisions are present. Such correlations occur if, during the characteristic time for a binary collision (estimated as λ D /v ∼ ω -1 p ), another perturber affects the atomic wavefunction. Such a perturbation occurs at a time scale of the order of the inverse matrix elements of the collision operator. A complete treatment of correlated collisions involving the infinite hierarchy Eq. ( 5) seems out of reach. To proceed further we propose a simplification that uses the singular role of the p + 1-th particle in the p + 1-correlation present in the integral term of Eq. ( 5). The approach is inspired from diagrammatic techniques used in kinetic theory ("ring approximation", e.g. [START_REF] Balescu | Equilibrium and Nonequilibrium Statistical Mechanics[END_REF]). For all p ≥ 1, we assume that p-correlations are created or destroyed due to the p-th particle only. Therefore, the sums involving V in Eq. ( 5) are replaced by their last term. The following evolution equation is obtained

  ∂ ∂t + iL 0 + p ∑ j=1 v j • ∂ ∂r j + iV (p)   Γ p (1...p) = -if 1 (p)V (p)Γ p-1 (1...p -1) -iN ∫ d(p + 1)V (p + 1)Γ p+1 (1...p + 1). ( 13
)
To solve this system (p ≥ 1) it is practical to use operator techniques in the Laplace space.

We define the resolvent

G p (z) =   z -L 0 + i p ∑ j=1 v j • ∂ ∂r j -V (p)   -1 , ( 14
)
the inversion being taken in the operator sense, and we introduce two operators A and P acting on the correlations

(A Γp-1 )(1...p) = G p (is)f 1 (p)V (p) Γp-1 (1...p -1), (15)
(P Γp+1 )(1...p) = G p (is)N ∫ d(p + 1)V (p + 1) Γp+1 (1...p + 1). ( 16
)
These operators can be interpreted as creating and destroying correlations, respectively. The evolution equation ( 13) can be written in the Laplace space in a compact form

Γp = A Γp-1 + P Γp+1 . ( 17
)
We solve this equation using the ansatz P Γp+1 = C Γp where C depends on A and P . This relation stems from the fact that lim p→+∞ Γ p = 0 is assumed. A recurrence relation for Γp

is obtained Γp = (1 -C) -1 A Γp-1 . ( 18
)
This leads to a formal solution of Eq. ( 17):

Γp = [(1 -C) -1 A] p Φ0 . ( 19
)
Here we have again used the convention Γ0 ≡ Φ0 . The latter quantity is obtained from Eq. ( 17) with p = 0, A Γ-1 ≡ iG 0 (is), using that

P Γ1 = C Φ0 Φ0 = i(1 -C) -1 G 0 (is). (20)
A closed relation for C is obtained from Eq. ( 18), multiplying on the left by P and again using the ansatz (P Γp = C Γp-1 ):

C = P (1 -C) -1 A. ( 21
)
This equation is transcendental because P is not invertible. A resolution should therefore involve an iterative scheme.

The collision operator is obtained from Eq. ( 20), performing algebraic manipulations and identifying with the terms in Eq. ( 11)

K(s) = iG -1 0 (is)C 0 . ( 22
)
Here, C 0 denotes a restriction of C such that CΦ 0 ≡ C 0 Φ 0 . This quantity is determined using the recurrence relation (p ≥ 0)

C p = P p (1 -C p+1 ) -1 A p+1 , ( 23
)
which stems from Eq. ( 21) and where A p , P p are restrictions of A and P such that AΓ p ≡ A p+1 Γ p and P Γ p ≡ P p-1 Γ p . We define a set of operators K 0 , K 1 , K 2 ... generalizing the collision operator, in such a way that the following property is satisfied

K p (1...p; s) = iG -1 p (is)C p . ( 24
)
These operators obey the following recurrence relation (the dependence on 1...p, p + 1 is not written explicitly)

K p (s) = iG -1 p (is)P p [1 + iG p+1 (is)K p+1 (s)] -1 A p+1 , (25)
or, using the explicit definition of A p+1 and P p and making elementary algebra

K p (s) = iN ∫ d(p + 1)V (p + 1) [ G -1 p+1 (is) + iK p+1 (s) ] -1 V (p + 1)f 1 (p + 1). ( 26
)
We simplify this relation by using the ansatz K p (1...p; s) = K(s

+ ∑ p j=1 v j • ∂/∂r j
), which stems from the presence of is + i ∑ p j=1 v j • ∂/∂r j in the denominator. Setting p ≡ 0 in Eq. ( 26) yields a closed equation for the collision operator

K(s) = iN ∫ d1V (1) Ḡ1 (is)V (1)f 1 (1), ( 27
)
where Ḡ1 is a modified resolvent

Ḡ1 (z) = [ z -L 0 + iv 1 • ∂ ∂r 1 -V (1) + iK ( -iz + v 1 • ∂ ∂r 1 )] -1 . ( 28
)
The latter can be interpreted as a propagator involving energy levels "dressed" by the presence of the other particles, through the presence of the collision operator as a non-Hermitian part in the Liouvillian L 0 . In the case where the correlated collisions are neglected, Ḡ1 ≡ G 1

and the collision operator reduces to that obtained within the binary approximation, Eq. ( 12). The set of equations ( 27) and ( 28) is nonlinear. It is practical to define an effective propagator Q eff (1, t), satisfying the dressed one-particle Schrödinger equation

[ ∂ ∂t + iL 0 + iV (r 1 + v 1 t) ] Q eff (1, t) + ∫ t 0 dτ M (τ )Q eff (1, t -τ ) = 0, (29)
with the initial condition Q eff (1, t = 0) = 1. Here M (t) is the inverse Laplace transform of the collision operator. Equation ( 29) is interpretable as describing the evolution of the atom under the influence of one collision represented by the interaction term V , given a set of collisions occurring in its past history. These collisions act accumulatively and are taken into account by the kernel M (t). The absence of this term in the unified theory [Eq. ( 8)] stems from the assumption Γ 2 ≡ 0. The latter is valid provided the characteristic evolution time for Q eff be much shorter than the time between two collisions. In the case of weakly correlated plasma, the evolution time is of the order of λ D /v, so that a validity criterion is provided by the relation λ D γ/v ≪ 1 where γ is a typical matrix element of the collision operator, e.g. estimated as γ = N b 2 W v ln(λ D /b W ). This result is in agreement with the analysis reported on in [27,28] and it indicates that the present extension of the unified theory allows one to explore regimes where correlated collisions are present. In terms of the adimensional parameter h = N b 2 W λ D ln(λ D /b W ), such regimes correspond to h > ∼ 1. The collision operator accounting for correlated collisions can be written in a fashion similar to that obtained within the unified theory:

K(s) = N ∫ ∞ 0 dte -st ∫ d1V (r 1 + v 1 t)Q eff (1, t)V (1)f 1 (1). ( 30
)
The substitution of Q by Q eff denotes a renormalization of the atomic energy levels induced by the correlated collisions. This stems from the structure of the modified resolvent Eq. ( 28). This is also illustrated by rewriting Eq. ( 29) in the Fourier and Laplace domain:

Qeff (k 1 , v 1 , s) = [s+iL 0 +K(s)] -1 [ (2π) 3 δ(k 1 ) -i ∫ d 3 k 2 (2π) 3 V (k 2 ) Qeff (k 1 -k 2 , v 1 , s -ik 2 • v 1 ) ] . ( 31
)
Here the convention F (k) = ∫ d 3 re -ik•r F (r) has been used for any function of space F (r). As can be seen, the collision operator appears in the denominator as a non-Hermitian contribution to the Liouvillian. Equation ( 31) presents similarities with the result of the resonance broadening theory used for plasma turbulence [START_REF] Dupree | [END_REF]32], where the quasi-linearized Vlasov equation plays a role similar to Eq. ( 29). In this theory, the coupling between the oneparticle distribution function and the plasma's electric field is described through a diffusion coefficient in the velocity space and the latter obeys a nonlinear equation as does our collision operator. Resonance broadening models have also been used in solid state physics for the calculation of dielectric functions (e.g., [33]).

In practice, a calculation of the collision operator from Eq. ( 30) should be done by iterations. A simplification, practical for numerical applications, is provided by assuming K(s) ≃ K(-iω 0 ) ≡ K 0 (with ω 0 being the central frequency of the line under consideration) in Eq. ( 31), using that the collision operator is governed by the values of Qeff near the resonance. This amounts to setting M (t) ≡ K 0 δ(t) in Eq. ( 29), and it leads to a simple expression for Q eff , with a structure similar to that in the binary case Eq. ( 9):

Q eff (1, t) = e (-iL 0 -K 0 )t T exp [ -i ∫ t 0 dτ e iL 0 τ V (r 1 + v 1 τ )e -iL 0 τ ] . (32)
A simple expression of the collision operator for hydrogen Lyman lines, suitable for small r 0 /λ D , is derived in Appendix A. Figure 1 shows the convergence of the iterations in a specific case. The diagonal element of the collision operator for Ly-α is plotted in terms of the number of iterations for r 0 /λ D = 0.1, assuming b W /r 0 = 0.5. These conditions correspond to h ≃ 7.5, i.e. a regime such that correlated collisions are present. As can be seen, the calculation converges after only a few iterations.

IV. COMPARISONS

We have applied the collision operator formula Eq. ( 30) to calculations of hydrogen line shapes in ideal cases. The effective propagator has been estimated from Eq. ( 32) and the simplification presented in Appendix A [Eq. (A4)] has been used. respectively. A numerical result from an ab initio simulation code is also shown in the figure.

This code uses the "collision-time technique", which has been developed and designed in the past so as to provide correct statistics at regimes where b W /r 0 ≪ 1 [34] (in particular, it avoids the bias induced by some prescribed boundary conditions). The technique assumes particles moving along straight lines, in agreement with the quasi-particle model used in the collision operator model. As can be seen, the binary model overestimates the width and predicts a different shape structure, with a dip at the center that increases with b W /r 0 .

This dip is a consequence of the inadequacy of the binary model. Correlated collisions are important in this region because they govern the average atomic evolution operator at long times, hence, by virtue of the Laplace transform, at small frequency detuning. The typical range for the dip corresponds to frequencies smaller than the matrix elements of K (see discussion in Appendix B). The renormalized model gives a much better result, with no dip and with an overestimate of the width no larger than 10% at b W /r 0 = 0.5. Practically, this better result stems from the presence of K into the denominator in Eq. [START_REF] Dupree | [END_REF]. The small discrepancy remaining at the center is due to the part of the correlations that is not retained in the resummation procedure. This discrepancy increases with the ratio b W /r 0 . The binary model overestimates the width and predicts a different shape, with a dip at the center that increases with b W /r 0 . The renormalized model gives a much better result, with no dip and with an overestimate of the width no larger than 10% at b W /r 0 = 0.5.

V. APPLICATION TO REALISTIC CASES

We consider an application of the model to realistic cases in the framework of magnetic fusion research. Observations of spectral lines with a high upper principal quantum number (typically n > 7, and up to n = 20 in some cases) have been reported on in experiments involving recombining plasmas, in Alcator C-Mod [35][36][37][38], in ASDEX Upgrade [39], in NAGDIS-II [40][41][42], in JET [43], and in NSTX [44,45]. For such lines, the ratio b W /r 0 for the electrons can be as high as 0.5 and the h parameter can be of the order of unity. This makes the binary approximation questionable. We have calculated the profile of the Ly 20 transition (n = 20 → 1) at N e = 4 × 10 13 cm -3 and T e = T i = 0.2 eV (typical conditions of divertor recombining plasmas), using the binary approximation and the renormalized model. The ionic contribution has been retained through statistical average over the microscopic electric field (quasi-static model). Practically, this amounts to making the formal substitution

L 0 → L 0 -d • E i ≡ L 0 -d z E i , (33)
I(ω) → I(ω, E i ), (34)
and to calculate the total line shape from the following integral

I tot (ω) = ∫ ∞ 0 dE i W (E i )I(ω, E i ), (35)
where W (E i ) is the probability density function of the ionic electric field. For the latter, we have used the Hooper model [46], which is suitable for moderately coupled plasmas such as those of magnetic fusion experiments. The profile obtained within the binary approximation again differs from that obtained within the renormalized model. The dip present at the center stems from the neglect of correlated collisions. The typical frequency range is given by the matrix elements of -iK -d z E 0 (with E 0 being the Holtsmark microfield). In a diagnostic context, this means that the use of the binary unified theory could lead to significant errors in the interpretation of spectra. We have also applied the model to the ion broadening of Ly-α in the atom's frame of reference (i.e. without Doppler broadening), motivated by opacity calculations for ITER modeling (e.g. [47,48]). The plasma conditions expected in the ITER divertor are such that the ions are in a "near impact" regime, with b W /r 0 < ∼ 0.5 [49]. The typical magnetic field is so strong that the Zeeman effect splits the line into well separated components. Figure 4 shows a plot of the central component obtained with the two models at N e = 2 × 10 15 cm -3 and T e = T i = 1 eV. Here again, the use of the binary approximation yields an incorrect result with a dip at the center. The influence of correlated collisions on Stark broadening is not specific to Lyman lines. We have applied the model to the H-α transition (n = 3 → 2), a line in the visible range. Figure 5 shows a plot of the spectral profile in the atom's frame of reference at N e = 10 15 cm -3 and T e = T i = 1 eV.

Only the transition 3d → 2p is presented in the figure, and no magnetic field is retained.

The profile has been compared to numerical simulations. As can be seen, the binary theory strongly overestimates the width and again provides a dip. The renormalized model is in a much better agreement, with a discrepancy not larger than 10%.

VI. CONCLUSION

We have investigated the role of correlated collisions on Stark-broadened line shapes in plasmas using a kinetic-theory approach. Correlated collisions occur if the collision frequency becomes of the same order as or larger than the plasma frequency. In the model, the average atomic evolution operator is obtained from a hierarchy of equations similar to BBGKY.

The presence of correlated collisions results in a renormalization of the atomic energy levels considered in a binary model. The results of the present model are in a very good agreement with numerical simulations. In contrast, the use of a binary model at conditions such that correlated collisions are present yields an erroneous line shape structure, with a dip present at the center. The method presented in this work provides an important correction to the usual models based on collision operators. The renormalization gives a structure of the new collision operator convenient for physical interpretation and suitable for numerical calculations, in particular for an implementation in line shape codes that use a binary collision operator (e.g. [20,50]). The model can be applied either to ions or electrons, and it can be generalized to isolated lines. Also, the kinetic theory treatment is transposable to quantum plasma [51]. provided by

Q(1, t) ≃      0 if r 1 < r m e -iL 0 t otherwise , ( A1
)
where r m is a characteristic radius that remains to be determined. A possible approach consists in estimating the matrix elements of the integral ∫ t 0 dτ e iL 0 τ V (r 1 + v 1 τ )e -iL 0 τ that appears in the Dyson series and comparing them to unity [52]. Such a procedure yields the following relation

r m =      b W 1 if b W 1 < v 1 t √ b W 1 v 1 t otherwise , ( A2
)
where b W 1 = hn 2 /m e v 1 is the Weisskopf radius specific to velocity v 1 . This relation exhibits the separation between the short and long time regimes, which correspond to incomplete and complete collisions, respectively. The interaction term has been estimated as

V (r 1 + v 1 τ ) ∼ H(r m /v 1 -τ )b W 1 v 1 /r 2 m
, H being the Heaviside function. From the approximation Eq. (A1), the space integral in Eq. ( 30) can be calculated analytically, in a fashion similar to what is done for the electric field autocorrelation function (e.g. [53]). We assume a Coulomb field and retain Debye screening through a cut-off at λ D . Hence,

∫ d 3 r 1 V (r 1 + v 1 t)Q eff (1, t)V (1) = 4π(b W 1 v 1 ) 2 3(n 2 ea 0 /h) 2 d • e (-iL 0 -K 0 )t d ×            (1/r m -1/λ D ) if v 1 t ≤ r m ≤ λ D (1/v 1 t -1/λ D ) if r m ≤ v 1 t ≤ λ D 0 otherwise . (A3)
The time integral in Eq. ( 30) can be evaluated analytically. The matrix element γ(s) = ⟨n10|K(s)|n10⟩, which provides the line width and shift of Lyman lines, is given by

γ(s) = 4πN (b W v) 2 3(n 2 ea 0 ) 2 ∑ lm |⟨n10|d|nlm⟩| 2 ×    [ erf ( b W λ D ) - 2 √ π b W λ D e -(b W /λ D ) 2 ]   √ π s ′ b W v erf   √ λ 2 D s ′ b W v   - 1 -e -λ 2 D s ′ /b W v λ D s ′   + ∫ ∞ b W v/λ D dv 1 F (v 1 ) [√ π s ′ b W v erf ( √ s ′ b W v v 1 ) + 1 v 1 E 1 ( s ′ b W v v 2 1 ) - 1 v 1 E 1 ( λ D s ′ v 1 ) - 1 -e -λ D s ′ /v 1 λ D s ′ ]} , ( A4
)
where b W is the Weisskopf radius evaluated at the thermal velocity v = √ 2k B T /µ (µ is the reduced mass of the emitter-perturber system), s ′ = -i∆ω + ⟨nlm|K 0 |nlm⟩ with ∆ω being the frequency detuning,

F (v 1 ) = 4(v 1 /v) 2 exp[-(v 1 /v) 2 ]/(v √ π)
is the Maxwellian velocity distribution function, and erf, E 1 are the error and first exponential integral functions, respectively. The Hilbert space formalism is used here, i.e. the bras, kets, and operators refer to atomic states instead of dyadics in the Liouville space. Equation (A4) is suitable for numerical calculations. The binary approximation ("usual" unified theory) corresponds to setting s ′ ≡ -i∆ω, and the complete collision approximation ("impact" theory) corresponds to setting s ′ → 0. An explicit calculation for b W ≪ λ D yields

γ(s) ≃ 4 √ πN b 2 W v 3(n 2 ea 0 ) 2 ∑ lm |⟨n10|d|nlm⟩| 2    2 + E 1   ( b W λ D ) 2      , ( A5
)
which corresponds to the standard result for hydrogen [START_REF] Griem | [END_REF]. In the case where a static electric field is present (Sec. V), the substitution Eq. ( 33) leads to the presence of a static Stark effect term in the Liouvillian in Eq. (A3). This term can be accounted for in the formula Eq.

(A4) by using a decomposition onto the parabolic base. This procedure yields the following substitutions where n 1 , n 2 are the parabolic quantum numbers. A similar procedure can be applied to the other matrix elements of the collision operator, which are required if a static Stark effect term is present.

APPENDIX B

An analysis of the failure of the binary model when correlated collisions are present can be performed from an estimate using a simplification of the collision operator formula Eq.

(A4). We ignore the velocity average, so that F (v 1 ) can be replaced by δ(v 1 -v), and we assume b ′ W ≪ λ D . We also ignore the imaginary part and we focus on frequencies much smaller than v/b W and larger than v/λ D , in such a way to investigate a region dominated by incomplete weak collisions. With these orderings the diagonal matrix elements of the binary collision operator can be evaluated as follows

γ(∆ω) ≃ γ 0 ln ( v b W ∆ω ) , ( B1
)
γ 0 = 4πN b 2 W v 3(n 2 ea 0 ) 2 ∑ lm |⟨n10|d|nlm⟩| 2 . ( B2
)
This expression is similar to the Lewis model [23]. The line shape is given by

I(∆ω) = C γ(∆ω) π 1 ∆ω 2 + γ(∆ω) 2 , ( B3
)
where C is a normalization constant. A dip occurs if the line shape has extrema at a finite

Figure 2 FIG. 1 :

 21 FIG.1: Diagonal matrix element of K 0 for Ly-α in terms of the number of iterations for r 0 /λ D = 0.1, assuming b W /r 0 = 0.5. The calculation converges after only a few iterations.

FIG. 2 :

 2 FIG. 2: Profile of Ly-α at b W /r 0 = 0.2 (a) and 0.5 (b), i.e. when correlated collisions are present.

FIG. 3 :

 3 FIG.3: Profile of Ly 20 in recombining divertor plasma conditions. The dip present at the center stems from the neglect of correlated collisions. No dip is present if the renormalized model is used.

FIG. 4 :FIG. 5 :

 45 FIG. 4: Profile of Ly-α, without Doppler broadening, in plasma conditions relevant to the ITER divertor. Only the central Zeeman component is shown here. The binary model again overestimates the width and predicts a different shape.

1 |⟨nlm|n 1 n 2

 12 m⟩| 2 × ... (A6) s ′ → s ′ -i⟨n 1 n 2 m|d z E|n 1 n 2 m⟩/h, (A7)

( 1 +

 1 value for ∆ω, larger than v/λ D . The resolution of the equation dI(∆ω)/d∆ω = 0 shows that this occurs for ∆ω = γ(∆ω)/ √ (1 + 2 ln(v/b W ∆ω)). The right-hand side of this expression can be estimated taking a constant value Λ of the order of unity for the logarithm. This provides a position for the local maximum of the line shape and a criterion for the presence of a dip. If Λ is estimated such that Λ/ √ 2Λ) = ln(λ D /b W ), this criterion reduces to h > ∼ 1 where h is the parameter introduced in the main text.

  is the Laplace transform of the evolution operator averaged over the per-

	turber trajectories (classical path assumption). The latter is obtained from the following
	definition with p = 0
	∫
	Φ s (1...p; t) =	d(p + 1)...dN f N (1...N )U (1...N ; t).

Acknowledgments

This work was carried out within the framework of the European Fusion Development Agreement and the French Research Federation for Fusion Studies. It is supported by the European Communities under the contract of Association between Euratom and CEA. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

APPENDIX A

We derive here a simplification of the collision operator formula applicable to the hydrogen Lyman lines in the limit of small r 0 /λ D . We assume that collisions with a small impact parameter (strong collisions) give no significant contribution to the integrals in Eq. [START_REF] Balescu | Equilibrium and Nonequilibrium Statistical Mechanics[END_REF].

On the other hand, we consider that the atomic wavefunction is not sensitive to perturbers located at large distance. In this framework, an approximation for the evolution operator is


J Rosato

H Capes

L Godbert-Mouret

M Koubiti

Y Marandet

R Stamm

Accuracy of impact broadening models in low-density magnetized hydrogen plasmas

The impact approximation used in the modelling of Stark profiles is examined when a magnetic field is present. Motivated by tokamak plasma spectroscopy, we calculate line shapes and S-matrix elements for the first Lyman lines of hydrogen with two models proposed for retaining simultaneously Stark and Zeeman effects in the impact limit. An evaluation of the accuracy of the two approaches is made with the help of a numerical simulation.

Introduction

Detailed line shapes are of interest for accurate plasma diagnostics in magnetic fusion [1]. A renewal of interest in the modelling of Stark broadening has been provided recently in the context of radiation transport simulations for tokamak divertor physics. In divertors operating at high density regimes (up to 10 15 cm -3 and higher: Alcator C-Mod [2] and soon ITER [3]), the plasma close to the wall contains a sufficiently large amount of neutrals so as to be optically thick to the first Lyman lines of the hydrogen isotopes. Coupled plasma-radiation-neutrals simulations have indicated the possibility for a shift of the ionization-recombination balance induced by the photoabsorption, in particular due to the Ly-α transition [4]. The need for reliability in these simulations has led to the development of accurate models for the spectral profiles which are involved in the photon emission and absorption rates appearing in the equation of radiative transfer. Efforts have been carried out in order to properly describe the ion dynamics on the Stark shape in the atoms' frame of reference [5]. The latter reference discusses in details the formation of Ly-α and includes the combined action of the Zeeman effect and the fine structure, which can be of the same order of magnitude (e.g. [6,7] for early works).

This paper deals with Stark broadening at plasma conditions such that the impact approximation can be used for ions. In a recent work [8], a collision operator model accounting for the Zeeman degeneracy removal of the energy levels has been developed in the framework of the standard impact theory (Griem-Baranger-Kolb-Oertel or GBKOmodel [9]), hereafter referred to as ST. An application to the Lorentz triplet of Ly-α (strong field regime) has indicated the possibility for having the lateral components much narrower than the central component. In a work published last year ( [10], note added in proof), it has been mentioned that the use of the ST leads to an overestimate of the "primary, adiabatic" and "secondary, nonadiabatic" contributions. Details have been published at the same time in [11]. We review in this paper the main features of the ST and of the alternative impact theory referred to as "generalized theory" (GT) in Refs. [10,11]. The GT uses a different interaction representation than the ST and avoids the use of a cut-off for small impact parameters. In this work, we perform a detailed comparison between the two approaches, and discuss their respective accuracy in the treatment of strong collisions. We propose a way for a fair comparison between these approaches, and calculate the two resulting line shapes. Other comparisons are made for the S-matrix and line broadening cross section, for which the two models are also confronted to a numerical simulation enabling an evaluation of their accuracy. Since the GT is designed for parabolic states, our calculations neglect fine structure. Section 2 gives an introduction to the impact formalism, Sec. 3 and Sec. 4 present the ST and the GT, respectively, and Secs. 5, 6 report on comparisons with simulations. Details and clarifications for a meaningful comparison of the two approaches are given in the Appendix.

Ion impact broadening: setting up the formalism

The spectral profile of an atomic line, I(ω), is given by the Fourier transform of the atomic dipole autocorrelation function C(t):

I(ω) = 1 π Re ∫ +∞ 0 dtC(t)e iωt , (1)
C(t) = { Tr(ρ u ⃗ d ⊥ • U † l (t) ⃗ d ⊥ U u (t)) } . (2)
Here, ⃗ d ⊥ is the dipole operator projected into the polarization plane; U (t) is the evolution operator; ρ is the density operator; the brackets {...} stand for a statistical average over the perturbers' states and the trace Tr(...) is performed over the atomic states. The subscripts u and l refer to the restriction of operators in the subspace corresponding to the upper and lower energy levels of the transition, respectively. For the sake of simplicity when examining the ST and the GT we will neglect the broadening of the lower levels (strictly speaking this is applicable to Lyman lines in the absence of fine structure and neglecting quenching interactions), so that U l (t) can be replaced by 1 l , the identity matrix in the lower levels' subspace. A more general treatment, accounting for the lower levels, would consist in replacing all operators by super-operators acting in the Liouville "line" space [12]. Since all of the development done hereafter will exclusively concern the upper levels' subspace, we will not write the subscript u anymore. For notational simplicity we will use SI-atomic units (m e = e = h = a 0 = 1, ε 0 µ 0 c 2 = 1, [START_REF] Bransden | Physics of Atoms and Molecules[END_REF]).

Within the impact approximation for ions, the evolution operator is given by U (t) = exp[(-iH 0 -Φ)t] where H 0 is the atomic Hamiltonian accounting for the Zeeman effect, i.e. H 0 = H at -⃗ µ • ⃗ B with ⃗ µ, ⃗ B being respectively the magnetic moment and the magnetic field; and Φ is the ion collision operator {the electron contribution is neglected since weaker by a factor (2m e /m i ) 1/2 at equal temperatures, m e , m i being the electron and ion masses, respectively [8]}. Its matrix elements denote typical collision frequencies and provide the width and the shift of the line under consideration, of Lorentzian shape. This operator is defined by

Φ = N ∫ +∞ 0 dvf (v)v ∫ λ D 0 db2πb { 1 -S( ⃗ b, ⃗ v) } angle , (3)
where the integrals are done over the velocity (module) v and the impact parameter (module) b, N and f stand for the perturbers' density and velocity distribution function, the brackets denote an angular average and S is the scattering matrix corresponding to a binary collision. The emitter motion is accounted for through a reduced mass model. The upper bound λ D is the Debye length, which accounts phenomenologically for the screening of the Coulomb field due to the electrons. In the following we consider the cross section operator

σ(v) = ∫ λ D 0 db2πb { 1 -S( ⃗ b, ⃗ v) } angle
and its matrix elements

⟨α|σ(v)|α ′ ⟩ ≡ σ αα ′ .

Standard impact theory

In the standard theory (ST) [9], the S-matrix is calculated in the interaction picture defined by the unitary transformation

|ψ(t)⟩ → | ψ(t)⟩ = exp(iH 0 t)|ψ(t)⟩.
This leads to a time-ordered exponential

S = T exp [ -i ∫ +∞ -∞ dtV (t) ] = 1 -i ∫ +∞ -∞ dtV (t) - ∫ +∞ -∞ dtV (t) ∫ t -∞ dt ′ V (t ′ ) + ... (4)
where

V (t) = -⃗ d(t) • ⃗ E(t)
is the Stark term in the interaction picture corresponding to a binary interaction with an ion

[i.e. ⃗ E(t) = -( ⃗ b + ⃗ vt)/(b 2 + v 2 t 2 ) 3/2 ]
. Following [9], we separate the b-integral in Eq. (3) into two parts, which correspond to two regimes. (i) Weak collisions: The time-integrated perturbation is sufficiently small so that a second order approximation can be used. A condition on ( ⃗ b, ⃗ v), which ensures the validity of this approximation, may be formulated as |{S (2) αα } angle | ≤ 1, where S (2) stands for the second order term of the S-matrix (with the convention S = S (0) + S (1) + S (2) + ...) and α refers to the state |α⟩. (ii) Strong collisions: When the second order approximation cannot be applied, the S-matrix becomes a strongly oscillating function of b, so that its integral over b can be neglected compared to the integral of the identity matrix (such a treatment is appropriate in particular when the contribution due to weak collisions is dominant, e.g. [START_REF] Alexiou | [END_REF] for a detailed discussion). For each state |α⟩, the two regimes (i) and (ii) correspond to impact parameters large and small compared to the strong collision radius b st α , respectively. The latter is defined by the following transcendental equation

[8] b st α = b W α    [ K // α + K ⊥ α A ST ( b st α b m )] 2 + [ K ′⊥ α B ST ( b st α b m )] 2    1/4 , (5)
where b W α = (2/3) 1/2 n 2 α /v is the Weisskopf radius (n α is the principal quantum number); b m = 2v/B is the magnetic cut-off; K // α , K ⊥ α , K ′⊥ α are atomic constants defined by

K // α = 1 n 4 α ∑ α ′ |z α ′ α | 2 , (6)
K ⊥ α = 1 n 4 α ∑ α ′ [ |x + α ′ α | 2 + |x - α ′ α | 2 ] , (7)
K ′⊥ α = 1 n 4 α ∑ α ′ [ |x + α ′ α | 2 -|x - α ′ α | 2 ] , (8)
with the convention x ± = (x ± iy)/ √ 2; and

A ST (s) = s 2 [K 2 0 (s) + K 2 1 (s)] and B ST (s) = s 2 [K 0 (s)I 0 (s) -K 1 (s)I 1 (s)
] with I 0,1 and K 0,1 being the modified Bessel functions.

In the following b st α ≤ λ D will be assumed. We can write the diagonal matrix elements of the cross section operator as follows:

σ αα = ∫ b st α 0 db2πb + ∫ λ D b st α db2πb { 1 -S (0) αα ( ⃗ b, ⃗ v) -S (1) αα ( ⃗ b, ⃗ v) -S (2) αα ( ⃗ b, ⃗ v) } angle = π ( b st α ) 2 - ∫ λ D b st α db2πb { S (2) αα ( ⃗ b, ⃗ v) } angle . (9)
Here we have used S (0) = 1 and {S (1) } angle = 0 (isotropic medium). The angle average and the integral over b present in the second term of the right-hand side can be performed analytically. An explicit calculation leads to [8] (10) where a ST (s) = sK 0 (s)K 1 (s) and b ST (s) = π/2 -πsK 0 (s)I 1 (s) are the GBKO functions (s > 0). In accordance with the model, there is a separation between the strong collisions (first term of the right-hand side) and the weak collisions (second term). The strong collisions are described by a Lorentz-Weisskopf term of characteristic radius b st α . The latter is calculated from Eq. ( 5). It accounts for the Zeeman degeneracy removal and reduces to b W α (up to a factor dependent on the quantum numbers) when B → 0. The weak collision term is split into two parts, which correspond to the electric field component parallel (//) and perpendicular (⊥) to ⃗ B. The magnetic field leads to a reduction of the perpendicular contribution. In the particular case where

σ αα = π ( b st α ) 2 +2πb 2 W α { K // α ln ( λ D b st α ) + [ K ⊥ α a ST (s) -iK ′⊥ α b ST (s) ] s=b st α /bm s=λ D /bm } ,
b st α ≪ b m ≪ λ D , a ST (s)| s=b st α /bm
s=λ D /bm behaves as ln(b m /b st ), indicating that the effective range of the perpendicular electric field is reduced from λ D to b m (see [8] for an interpretation with classical mechanics). In practice, such a reduction leads to a narrowing of the Zeeman components. As a particular example consider the lateral components of the Zeeman-Lorentz triplet of hydrogen Ly-α, which are broadened due to the perpendicular contribution only. Figure 1 shows the blue component calculated assuming degenerate 5), (10)] and assuming finite b m , at N e = 10 13 cm -3 , T e = T i = 3 eV, and B = 4 T. As can be seen, the model predicts a component much narrower (by more than a factor of two) than expected when neglecting the Zeeman degeneracy removal. The component is also shifted to larger frequencies [this stems from the presence of an imaginary part in Eq. ( 10)].

Generalized theory

The GT is based on an interaction picture different from that introduced above. We follow the synthesis presented in [15] (details on the derivation of formulas can be found in [16]). All matrix elements refer to the parabolic base. The unitary transformation is done by including the contribution of the electric field parallel to ⃗ B in the unperturbed Hamiltonian, namely

|ψ(t)⟩ → | ψ(t)⟩ = exp[iH 0 t -i ∫ t -∞ dt ′ d // E // (t ′ )]|ψ(t)⟩ ≡ Q † (t)|ψ(t)⟩.
This leads to the following expression for the S-matrix

S = S a T exp [ i ∫ +∞ -∞ dtQ † (t) ⃗ d ⊥ • ⃗ E ⊥ (t)Q(t) ] , (11)
S a = exp [ i ∫ +∞ -∞ dtd // E // (t) ] . (12)
The quantity S a denotes the solution that would be obtained in the B = 0-case if timeordering was neglected (adiabatic approximation). If B is finite, S a corresponds to the solution that would be obtained if the virtual transitions between the Zeeman energy sublevels (i.e., the inelastic collisions) were neglected. Using Eqs. ( 11) and ( 12), S can be decomposed into S a + S na , where the non-adiabatic contribution S na is given by

S na = S a { T exp [ i ∫ +∞ -∞ dtQ † (t) ⃗ d ⊥ • ⃗ E ⊥ (t)Q(t) ] -1 } . (13)
Accordingly, the cross section operator also splits into two parts σ = σ a + σ na . In the GT, S a is calculated exactly and a second-order approximation is done concerning S na .

No distinction is made between the weak and strong collisions, i.e., the integral over b is carried out from 0 to λ D even though the second-order expansion is not valid for small values of b. In [15,16], this procedure is justified by the fact that the integral over b converges without the need for a cut-off. Physically, this convergence is due to the dressing of the energy levels by E // (t). Collisions with small impact parameters lead to a strong degeneracy removal of the energy levels, so that the effect of ⃗ E ⊥ is reduced in such a way that the cross section remains finite. However, this does not guarantee that the second-order approximation provides a reliable estimate of S na , since the contribution of higher order terms is in principle not negligible. We will see in the next section that ignoring these higher order terms leads to overestimate the strong collision contribution and, eventually, the line width.

The adiabatic cross section operator is given by

σ a = ∫ λ D 0 db2πb{1 -S a } angle . ( 14
)
It is diagonal in the parabolic base. The corresponding cross section reads

σ a = 2π ( 2z αα v ) 2 I ( vλ D 2z αα ) , (15)
I(s) = 1 6 [ 3s 2 -2s 3 sin ( 1 s ) -s 2 cos ( 1
s ) + s sin ( 1 s ) -Ci ( 1 s )] , ( 16
)
where Ci is the cosine integral function. I(s) has the following asymptotic behavior for large s:

I(s) ∼ 11 36 - γ 6 + 1 6 ln s, (17)
where γ ≃ 0.577 is Euler's constant. In this expression the logarithm corresponds to the weak collisions. The non-adiabatic cross section operator is given by

σ na = - ∫ λ D 0 db2πb{S na } angle , (18)
and the corresponding cross section reads Figure 2. The non-adiabatic contribution is mitigated by the Zeeman degeneracy removal, which leads to line narrowing. Here, plot of the blue lateral component of Lyman-α obtained using the generalized theory neglecting (dashed line) and retaining the degeneracy removal (solid line). As in the previous section, a shift is also visible.

σ na αα ′ = 4π 3v 2 ∑ α ′′ [ x + αα ′′ x - α ′′ α ′ + x - αα ′′ x + α ′′ α ′ ] ∫ λ D /bm 0 ds s C ± (χ, Y, s), (19)
exponential is evaluated using the scaling-squaring method, i.e. we use the identity exp[-i∆tV (t)] ≡ {exp[-i∆tV (t)/m]} m with m integer and we choose m sufficiently large so that the matrix elements of ∆tV (t)/m are small in absolute value (details on the scaling-squaring method can be found in [17] and Refs. therein). Then, the exponential is approximated by a low-order expansion. The angle average is performed by Monte-Carlo integration. We have examined the matrix element ⟨211|{S( ⃗ b, ⃗ v)} angle |211⟩, which is affected by the non-adiabatic contribution exclusively (here a notation in terms of the spherical base quantum numbers is used). Figure 3 shows a plot of the real and imaginary part in terms of b (atomic units), assuming v = √ 2T i /m i , at N e = 10 13 cm -3 , T e = T i = 3 eV and B = 4 T, obtained from the numerical calculation and compared to the ST and the GT. The broadening function involved in the GT has been evaluated following Eq. (20). Here, χ = 0. The double integral has been calculated by Monte-Carlo method, interpreting [w(x)w(x ′ )] 3 as a joint probability density function for x, x ′ . As can be seen, both the ST and the GT agree well with the numerical result at impact parameters larger than ∼ 200, which roughly corresponds to the strong collision radius b st as defined in Eq. (5). Below this radius the numerical solution oscillates around zero, with an increasing frequency as b decreases. This reinforces the choice of a cut-off made by the ST as described in Sec. 2. The GT yields a solution much different from the Figure 3. Real and imaginary parts of ⟨211|{S( ⃗ b, ⃗ v)} angle |211⟩, obtained within the ST and the GT, and compared to a simulation. For impact parameters larger than the strong collision radius b st (of the order of 200), the ST and the GT agree well with the numerical calculation. The presence of oscillations at impact parameters smaller than b st corresponds to strong collisions which, within the ST, are accounted for by a cut-off procedure. In this domain the GT overestimates the real and imaginary parts in absolute value. numerical one, with a strong overestimate of the real and imaginary parts in absolute value (in particular the real part obtained with the GT diverges when b → 0). This may lead to an overestimate of the Stark width and shift. This point is illustrated in Fig. 4, where plots of the real and imaginary parts of the cross section are presented. The result of the GT is always larger than that of the numerical calculation. A calculation of the integral Eq. (3) shows that the GT overestimates the line width and shift by 60% and 20%, respectively, whereas the ST never deviates from the simulation by more than 10%. . Real part of ⟨311|{S( ⃗ b, ⃗ v)} angle |311⟩ and the corresponding cross section, obtained within the ST and the GT, and compared to a simulation. The S-matrix obtained within the GT is not unitary at small impact parameters. The cross section is not overestimated here because of the presence of oscillations.

States with permanent dipoles

The particular case χ = 0 considered above denotes broadening of a parabolic state with spherical symmetry (i.e., z αα = 0 in the parabolic base). For hydrogen lines other than Ly-α, the non-adiabatic contribution also involves states with a permanent dipole moment. We have examined such states by confronting the ST and the GT to simulations. Figure 5 shows a plot of the real part of ⟨311|{S( ⃗ b, ⃗ v)} angle |311⟩ and the corresponding cross section at the same conditions as in the previous section. This case involves the parabolic kets |011⟩ and |101⟩, i.e. states with a permanent dipole moment (e.g. [START_REF] Bethe | Quantum Mechanics of One-and Two-Electron Atoms[END_REF] for details on the parabolic states). As can be seen, the GT yields oscillations at small impact parameters, with increasing amplitude as b decreases. Although the S-matrix resulting from the GT is no longer unitary in this regime, the presence of oscillations prevents from the overestimation discussed above and eventually leads to a cross section similar to that obtained either within the ST or the numerical calculation.

In general, the spherical states contribute to the broadening in addition to those with a permanent dipole, so that serious misestimates of the line width should still be expected a priori when using the GT. For example, an application to the blue lateral component of H-α (which involves the spherical state |322⟩) shows that the GT already overestimates the width obtained by simulation by 70% when the broadening of the lower level is neglected, whereas deviations no larger than 10% are obtained with the ST.

Conclusion

The work presented in this paper provides a complement to previous investigations in Stark profile modelling for magnetic fusion plasmas. Two models for ion impact broadening -namely, the standard theory (ST) and the generalized theory (GT) -have been examined and compared to simulations. In general, the ST provides accurate results whereas the GT systematically overestimates the width and the shift of the Zeeman components because of an inconsistent use of the perturbation theory involved in the non-adiabatic part of the collision operator. In particular, a second-order expansion is used for the strong collisions, i.e., precisely, where a perturbative treatment does not apply. Although such collisions are not frequent (since b st ≪ r 0 is assumed in the impact theory, r 0 being the mean interparticle distance), the failure of the perturbation theory yields a sufficiently large error so as to globally affect the line shape. This suggests a subsequent refinement of this theory, either based on a cut-off procedure or on an analytical extension of the Dyson series expansion up to higher orders. Our conclusion is in agreement with previous studies regarding the GT at zero magnetic field [START_REF] Alexiou | [END_REF]. The broadening function used in the GT involves a multiple integral of strongly oscillating functions, not suitable for fast numerical evaluation. A further refinement of the GT should consist in establishing a simple, CPU-convenient approximation for this integral.

In magnetic fusion plasmas, the typical densities are sufficiently low so that the fine structure can be in competition both with the Stark and Zeeman effects. Another refinement of the GT should also address this issue. For opacity calculation purposes, accurate and analytical expressions are required for line shapes, so that the ST (and derived approaches accounting for non-impact effects, e.g., the unified theory [20,21]) presently seems to be the best approach. m = 2, so that K // α = 9/64 ≃ 0.14. A straightforward use of the formulas given by Eqs. (10) and (15) ,

so that, at N e = 10 13 cm -3 and T e = T i = 4 eV (plasma conditions assumed in [11]), one has σ ST αα /σ GT,a αα ≃ 0.65. The magnitude of this deviation falls in the typical uncertainty range expected by the use of the GT and presented in the present paper.

Figure 1 .

 1 Figure1. Blue lateral component of the Zeeman-Lorentz triplet of the Lyman-α line, obtained using the standard impact theory neglecting (dashed line) and retaining the degeneracy removal (solid line). A strong reduction of the line width, by more than a factor of two, is present. A shift is also visible.

Figure 4 .

 4 Figure 4. Real and imaginary parts of ⟨211|σ|211⟩, obtained within the ST and the GT, and compared to the simulation. v T = √ 4T i /m i stands for the relative thermal velocity. The GT systematically overestimates the cross section, which results in an overestimate of the Stark width and shift.

Figure 5

 5 Figure 5. Real part of ⟨311|{S( ⃗ b, ⃗ v)} angle |311⟩ and the corresponding cross section, obtained within the ST and the GT, and compared to a simulation. The S-matrix obtained within the GT is not unitary at small impact parameters. The cross section is not overestimated here because of the presence of oscillations.

  shows that

	σ ST αα =	32π v 2	[	3 16	+ ln	(	vλ D 2 √ 6	)]	,	(25)
	σ GT,a αα	=	48π v 2	[	11 6	-γ + ln	(	12 vλ D	)]

Acknowledgments

This work is supported by the collaboration LRC DSM99-14 between the PIIM laboratory and the CEA Cadarache research center, within the French federation of research on magnetic confinement fusion (FRFCM).

where χ = z αα /δz with δz = (z αα -z α ′′ α ′′ ), and Y = δz/vb m . In the integral, C ± is the so-called broadening function, given by [15] C ± (χ, Y, s) = 3 4

Here j 0 , j 1 , j 2 denote the first spherical Bessel functions. The upper signs correspond to the nondiagonal matrix elements (α ′ ̸ = α) and the lower signs correspond to the diagonal matrix elements (α ′ = α). In the limiting case B → 0, C ± reduces to A ST + iB ST , the integral Eq. ( 19) diverges at the lower bound s = 0, and the upper bound tends to 0. In this case, a cut-off denoting separation between the weak and strong collisions (as in the ST) is used. At finite B, the virtual transitions between the Zeeman sublevels mitigate the role of the non-adiabatic contribution. This leads to a narrowing of the line components. This point is illustrated in Fig. 2. The blue lateral component of Ly-α obtained within the GT is superimposed to the result obtained at B = 0, at the same conditions as in the previous section. In this particular case, the line broadening is given by the non-adiabatic contribution only. As can be seen, the model predicts a component narrower than expected when neglecting the Zeeman degeneracy removal and the component is shifted to larger frequencies. The width reduction is weaker than that predicted with the ST, with only 40% instead of a factor larger than two.

Comparison of the two models to a simulation

We have compared the GT and ST to the result of a fully numerical calculation of the collision operator. In this method, the S-matrix is calculated for a set of impact parameters ⃗ b and velocities ⃗ v chosen in such a way that the perturber follows a straight line in the Debye sphere surrounding the atom, in agreement with the straight line trajectory assumption done both within the GT and the ST. The algorithm used for the integration of the Schrödinger equation uses the infinitesimal evolution operator U (t + ∆t, t), i.e. we evaluate S ≃ U (T /2, -T /2) for large T by iterating the relation U (t + ∆t, -T /2) = U (t + ∆t, t)U (t, -T /2). The infinitesimal evolution operator is approximated as

Appendix

In [10,11], it is mentioned that the ST overestimates the adiabatic contribution "by up to two orders of magnitude". We have examined the relevance of this statement. The largest deviation claimed in [11] concerns the (102) -(101) component of the Paschenα line (parabolic numbers), hence here we focus our analysis on this transition. The calculation of the ratio σ ST αα /σ GT,a αα in Ref. [11] is inconclusive for two reasons: (i) Both the broadening of the lower level and interference effects have been retained in σ GT,a αα {cf. the definition of X αβ in [11], Eq. (A.1)} whereas not in σ ST αα . These effects, absent on Lyman lines, become important as the lower principal quantum number of the transition increases. A meaningful comparison between the two theories should involve, or not, the broadening of the lower level and the interference effects, at the same time. In the presentation of the ST done above as well as in that reported in [8], we have not retained the interference effects, but they can be easily treated within the Liouville line space formalism. (ii) K // α ≡ 1 has been used for calculating the ST formulas. This is a very drastic approximation. For the upper level of the transition, n = 4, n 1 = 1, n 2 = 0,