P Kumar
email: prkumar@stanford.edu

S Tsikata

K Hara
email: kenhara@stanford.edu

Effects of multiply charged ions on microturbulence-driven electron transport in partially magnetized plasmas

Nonlinear interaction between kinetic instabilities in partially magnetized plasmas in the presence of multiply-charged ion streams is investigated using kinetic simulations. It was observed by Hara and Tsikata [Phys. Rev. E 102, 023202 (2020)] that the axial ion-ion two-stream instability (IITSI) due to singly and doubly charged ion streams, coupled with the azimuthal electron cyclotron drift instability (ECDI), enhances cross-field electron transport. In the present study, it is observed that the addition of triply charged ions (as a third ion species) contributes to damping of the excited modes, leading to a reduction in the cross-field electron transport. The net instability-driven electron transport is shown to be a function not only of the azimuthal modes, such as the ECDI, but of the multiple ion species that dictate the development of additional plasma waves. It is found that trapping of the higher ion charge states within the plasma waves results in a broadening of the ion velocity distribution functions.

I. INTRODUCTION

Cross-field electron transport across magnetic field lines plays an important role in fusion, astrophysical, and E × B discharge plasmas [START_REF] Doyle | Chapter 2: Plasma confinement and transport[END_REF][START_REF] Boeuf | Tutorial: Physics and modeling of Hall thrusters[END_REF] . While cross-field discharges have been widely studied theoretically [START_REF] Morozov | The conceptual development of stationary plasma thrusters[END_REF][START_REF] Boeuf | Rotating instability in low-temperature magnetized plasmas[END_REF][START_REF] Cavalier | Hall thruster plasma fluctuations identified as the E × B electron drift instability: Modeling and fitting on experimental data[END_REF][START_REF] Smolyakov | Fluid theory and simulations of instabilities, turbulent transport and coherent structures in partially-magnetized plasmas of E × B discharges[END_REF] , experimentally [START_REF] Smith | Time and space-correlated plasma potential measurements in the near field of a coaxial Hall plasma discharge[END_REF][START_REF] Yoshikawa | Anomalous diffusion of a plasma across a magnetic field[END_REF][START_REF] Parker | Transition in electron transport in a cylindrical Hall thruster[END_REF][START_REF] Ellison | Cross-field electron transport induced by a rotating spoke in a cylindrical Hall thruster[END_REF] , and numerically [START_REF] Mikellides | Numerical simulations of Hall-effect plasma accelerators on a magnetic-field-aligned mesh[END_REF][START_REF] Hara | An overview of discharge plasma modeling for Hall effect thrusters[END_REF][START_REF] Charoy | 2D axialazimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF][START_REF] Boeuf | E × B electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF] , turbulent transport of electrons across the magnetic field is still not well understood. Experimentally measured electron transport has been shown to be significantly higher than that predicted from classical collisional transport processes.

Over the last few decades, a number of different models, using fluid [START_REF] Mikellides | Numerical simulations of Hall-effect plasma accelerators on a magnetic-field-aligned mesh[END_REF][START_REF] Barral | Low-frequency model of breathing oscillations in Hall discharges[END_REF][START_REF] Keidar | Plasma flow and plasma-wall transition in Hall thruster channel[END_REF] , kinetic [START_REF] Adam | Study of stationary plasma thrusters using two-dimensional fully kinetic simulations[END_REF][START_REF] Héron | Anomalous conductivity in Hall thrusters: Effects of the non-linear coupling of the electron-cyclotron drift instability with secondary electron emission of the walls[END_REF][START_REF] Yan | Two-dimensional modeling of the cathode sheath formation during the streamer-cathode interaction[END_REF][START_REF] Cho | Kinetic particle simulation of discharge and wall erosion of a Hall thruster[END_REF] , and hybrid approaches [START_REF] Fife | Hybrid-PIC modeling and electrostatic probe survey of Hall thrusters[END_REF][START_REF] Bareilles | Critical assessment of a two-dimensional hybrid Hall thruster model: Comparisons with experiments[END_REF][START_REF] Parra | A twodimensional hybrid model of the Hall thruster discharge[END_REF] , have been proposed to explain these phenomena. Fully kinetic particle-in-cell (PIC) models have been recently used to investigate the role of plasma turbulence induced by kinetic instabilities on the cross-field transport of electrons. While electron diffusion across the magnetic field lines can be caused by fluctuations in the electric field in the E × B direction alone, microturbulence may be driven by nonlinear coupling between different linear instabilities [START_REF] Muschietti | Microturbulence in the electron cyclotron frequency range at perpendicular supercritical shocks[END_REF][START_REF] Tsikata | An axially propagating two-stream instability in the Hall thruster plasma[END_REF][START_REF] Villafana | 2D radial-azimuthal particlein-cell benchmark for E×B discharges[END_REF] and can influence the cross-field electron mobility [START_REF] Boeuf | E × B electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF][START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF] . In developing reducedorder models for electron transport, accounting for the effects of such mode coupling is necessary.

The existence of a variety of charge states within such plasmas has additional wide-ranging implications. Depending on the discharge characteristics, such as the operating power, triply and even quadruply charged ions can be generated in addition to the singly and doubly charged ions in HETs [START_REF] Gulczinski | Near-field ion energy and species measurements of a 5-kw Hall thruster[END_REF] . The generation of higher charge states is also a known feature of other E × B devices such as planar magnetrons [START_REF] Bohlmark | The ion energy distributions and ion flux composition from a high power impulse magnetron sputtering dis-charge[END_REF] operating in high-current pulsed regimes. Ions of different species also result in complex kinetic phenomenon, such as the mixture of ion species near the plasma sheaths [START_REF] Baalrud | Instability-enhanced collisional friction can determine the bohm criterion in multiple-ion-species plasmas[END_REF] and light ions potentially stabilizing plasma turbulence [START_REF] Kasparek | Co2-laser scattering from thermal fluctuations in a plasma with two ion components[END_REF] .

Previous work by Hara and Tsikata in Ref. 32 investigated the coupling between the electron cyclotron drift instability (ECDI) and ion-ion two-stream instability (IITSI) in Hall effect thrusters (HETs) due to the mixture of singly and doubly charged ion streams using kinetic theory and simulations. The relative velocity between singly and doubly charged ion streams causes the excitation of the IITSI, first identified and studied in coherent Thomson scattering measurements [START_REF] Tsikata | An axially propagating two-stream instability in the Hall thruster plasma[END_REF] . In addition, it was shown in Ref. 32 that cross-field electron transport is enhanced due to the multidimensional plasma waves initiated by the coupling between ECDI in the azimuthal direction (in the direction of the E×B drift) and IITSI in the axial direction (in the direction of the applied E field).

In the present paper, we now analyze a more complex, but physical, case involving three ion populations using the axialazimuthal PIC simulation. We consider the case in which a large-amplitude IITSI due to singly and doubly charged xenon ions is present while the ECDI is initiated in the cross-field direction. A third species of triply charged xenon ions is then introduced in the simulations and their effects on the multidimensional plasma wave and cross-field plasma transport are investigated by varying the triply charged ion fraction. Charged species streaming with different velocities nonlinearly interact with the plasma wave, leading to particle trapping that influences the nonlinear saturation of these instabilities. The trapping mechanisms for the three species considered in the simulation and their impact on the corresponding velocity distribution function (VDF) is also analyzed. This paper is organized as follows. The simulation method and physical and numerical parameters are described in Section II. Section III A discusses the coupling of ECDI and IITSI. In Section III B, the wave damping due to triply charged ions is analyzed. The effect of the wave damping on crossfield electron transport is discussed in Section III C. Finally, broadening of ion VDFs and nonlinear trapping of ion species This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0067305

is presented in Section III D.

II. SIMULATION METHOD AND PARAMETERS

The 2D axial-azimuthal Particle-in-Cell (PIC) model originally proposed by Boeuf and Garrigues 14 and later utilized by Charoy et al. [START_REF] Charoy | 2D axialazimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] and Hara et al. [START_REF] Hara | Cross-field electron diffusion due to the coupling of drift-driven microinstabilities[END_REF] is adapted and further extended to incorporate a triply charged ion species in the simulations.

A 2D Cartesian domain in the axial (x) and azimuthal (y) plane with axial length L x = 2.5 cm and azimuthal length L x = 1.28 cm is discretized uniformly using 500 and 256 cells, respectively. At the beginning of the simulations, electrons and ions are loaded uniformly with a density of 5.0 × 10 16 m -3 and a Maxwellian velocity distribution with electron temperature, T e = 10 eV, and ion temperature, T i = 0.5 eV. While the collisions are neglected in this setup, ionization events are modeled using an ionization rate function, S i (x), assumed constant in time, leading to a constant ion current density in steady state [START_REF] Boeuf | E × B electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF] . The steady-state continuity equation for the ions can be written as,

∇ • (j i1 + j i2 + j i3 ) = e(S i1 + 2S i2 + 3S i3 ),
where j ik is the ion current density, e is the elementary charge, and S ik is the ionization rate for k = 1, 2, 3 corresponding to singly, doubly, and triply charged ion species, respectively. We define α d and α t to be the fraction of doubly and triply charged ion current densities with respect to the total ion current density, such that the individual source terms are assigned as S i1 = (1α dα t ) S i (x), S i2 = α d S i (x)/2, and S i3 = α t S i (x)/3. Here, the total ionization rate S i (x) is defined as

S i (x) = S 0 cos π x -x M x 2 -x 1 ,
where x 1 = 0.25 cm, x 2 = 1.0 cm, x M = (x 1 + x 2 )/2, and S 0 is adjusted to impose j M = 400 A/m 2 to be the total ion current density.

The radial magnetic field is applied in the axial (x) direction, given by B(x) = B max + B 0 ξ (x), where

ξ (x) = 1 -exp 1 - x -x L 2σ b 2 and B 0 = (B a -B max )/ξ (0) if x < x L and B 0 = (B c - B max )/ξ (x L ) if x ≥ x L .
Here, x L = 0.75 cm, σ b = 0.625 cm, B max = 100 G, B a = 60 G, and B c = 10 G. Neutral atom dynamics and transport is not taken into account. Periodic boundary conditions are used in the azimuthal direction. The anode plane, located at the left axial boundary is kept at a fixed potential of 200 V. At the cathode plane, i.e., x = 2.4 cm, the azimuthally averaged potential is set to be zero [START_REF] Boeuf | E × B electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF] , and the electrons are injected from a random position in the azimuthal direction maintaining current conservation, i.e., Γ ec = Γ ea -Γ i1a -2Γ i2a -3Γ i3a , where Γ ec is the number of electrons re-injected from the cathode plane, and Γ ea , Γ i1a , Γ i2a , and Γ i3a are the number of electrons, singly charged ions, doubly charged ions, and triply charged ions absorbed at the anode plane, respectively. An explicit PIC code, developed within the Plasma Dynamics Modeling Laboratory at Stanford University and benchmarked against several different codes in the community, is utilized for the simulations presented in this paper [START_REF] Charoy | 2D axialazimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF][START_REF] Hara | Cross-field electron diffusion due to the coupling of drift-driven microinstabilities[END_REF] . The code is written in C++ and utilizes Message Passing Interface (MPI) for particle decomposition. The Poisson equation is solved using the scalable linear solver library, HYPRE. Simulations are performed with a grid size of 50 µm in both directions and the number of particles per cell reaches approximately 250 in the steady state. Numerical convergence studies have shown that these parameters are sufficient for resolving the physics processes of interest in this paper.

III. RESULTS AND DISCUSSION

The effect of increasing the fraction of triply charged ion species (α t ) in addition to singly and doubly charged ions is studied. The results presented in the paper use a condition in which the doubly charged ion fraction is fixed (α d = 20%) such that a large-amplitude IITSI coexists with ECDI [START_REF] Hara | Cross-field electron diffusion due to the coupling of drift-driven microinstabilities[END_REF] and the triply charged ion fraction is varied α t = 0 -10%. The plasma wave behavior reaches a steady state around 10 -15 µs and the simulations are run up to 30 µs to ensure that the plasma wave characteristics do not diverge.

A. Coupling between ECDI and IITSI

Figure 1 shows the instantaneous azimuthal (E y ) and axial (E x ) electric field due to coupling between ECDI and IITSI for α t = 0 in the top row and 10% in the bottom row, respectively, in steady state. The results without the triply charged ions (top row) are identical to the cases presented in Ref. 32.

Here, the addition of triply charged ions (bottom row) similarly shows the coexistence of ECDI (modulation along y) and IITSI (modulation along x). The azimuthal plasma wave observed in the simulation exhibits wavelength corresponding to the resonance at k y λ D ≈ 0.9, where the ECDI growth rate is maximum. The dominant wavelength of axial fluctuation is on the scale of 1 mm driven by the IITSI. As discussed in Ref. 32, the simulation results show that the ECDI is primarily driven in the region where E × B drift is largest (x ∈ [0.5 cm,1 cm] in this simulation) and the IITSI is driven in the region where the velocity difference between Xe + and Xe 2+ increases (x ≥ 1 cm in this simulation).

B. Wave damping due to addition of triply charged ions

As the fraction of triply charged xenon ions (Xe 3+ ) is increased, damping of the the plasma wave amplitude can be seen. To further quantitatively investigate such damping in the presence of triply charged ions, we perform fast Fourier transform (FFT) analysis of the two-dimensional electric field This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Amplitudes of the FFT of azimuthal field (E y ) and axial field (E x ) are shown in the (k x , k y ) plane for the 0% Xe 3+ case in Figs. 2(a1) and (b1), respectively. The peaks at k x = 0 in Fig. 2(a1) correspond to the resonance modes of the ECDI. Additionally, signatures of axial fluctuation of E y are evident in the bright spots at k x = 0, suggesting that the axial oscillations of E y are correlated with the axial oscillations of E x due to the IITSI, as shown in Fig. 2(b1). However, the ECDI resonances in the E y spectra are not found in the E x spectra, indicating that the axial plasma waves (E x ) are not influenced by the ECDI. This is because the ions are accelerated to velocities that can initiate the IITSI downstream of the region in which the ECDI occurs at k x = 0 due to the large E × B drift.

PLEASE CITE THIS ARTICLE AS

Figures 2(a2) and (b2) show the results when the fraction of triply charged ions is increased from 0% to 10%. While the FFT results show a structure similar to the case without triply charged ions, there is a slight change in the amplitude of the electric field fluctuations. Figure 2(a3) shows the FFT of E y as a function of k y , i.e., the results in Figs. 2(a1) and (a2) integrated in k x . It can be seen that the amplitude of the azimuthal field E y in the plume decreases approximately 30%, when 10% triply charged xenon ions are added (α t = 0.1). This is noteworthy since the electric field investigated here is within the plume. As shown later in Fig. 3(d), the timeaveraged, azimuthally-averaged axial electric field is not mod-ified due to the presence of triply charged ions, meaning that the maximum E × B drift around x ∈ [0.5 cm, 1 cm] remains unchanged. Hence, the azimuthal wave initiated by the E × B drift would not be locally affected, if there is no influence from other regions of the plasma flow. Additionally, using the 2D linear dispersion theory in Ref. 32, the ECDI growth rate was found to be unaffected by the presence of multiply charged ions. The decrease in the E y amplitude shown in Fig. 2(a3) is thus unexpected, as the results indicate that the axial IITSI in the presence of of Xe 3+ affects the azimuthal wave (the E y fluctuation) that is primarily driven by ECDI in the region where E × B drift is largest. The effects of mode damping in the presence of light and heavy ion species have been shown in different contexts, as was mentioned in the introduction [START_REF] Kasparek | Co2-laser scattering from thermal fluctuations in a plasma with two ion components[END_REF][START_REF] Du | Impurity effects on trapped electron mode in tokamak plasmas[END_REF] . It is interesting to note that in the present case, significant mode damping effects occur due to the presence of multiple ion species of the same mass but with different velocities. With differing charge states, these species reshape the distribution of energy in the fluctuations via their response to the applied accelerating field.

Finally, Fig. 3(b3) shows that the amplitude of the axial field FFT|(E x )|, integrated along k y and plotted against k x , decreases nearly 8% when triply charged ions are present. Since the linear dispersion theory indicates that the growth rate of the IITSI is unaffected by the presence of triply charged ions, the observed reduction of the plasma wave amplitude is due to the nonlinear dynamics between multiply charged ion streams. As will be shown later, multiply charged ions have a different bounce frequency due to particle trapping within the plasma wave. It is likely that such nonlinear phenomena lead to the waves having more temporal and spatial modes, leading to the damping of the plasma wave amplitude.

C. Effect of multiply charged ions on cross-field electron transport

Electron transport in the cross-field direction can be characterized by the axial electron bulk velocity [START_REF] Hara | An overview of discharge plasma modeling for Hall effect thrusters[END_REF][START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] :

u e,⊥ = -µ ⊥ E ⊥ + 1 en e ∇ ⊥ p e , (1)
considering the drift-diffusion approximation, i.e., in steady state and when the inertia terms are negligible. Here, u e,⊥ is the cross-field electron bulk velocity, µ ⊥ is the cross-field mobility, E ⊥ = E x is the axial electric field, and p e is the electron pressure. Note that when collisional processes are dominant, the cross-field mobility can be derived as

µ ⊥ = µ 0 (1 + Ω 2 ) -1
, where µ 0 = e/(m e ν m ) is the unmagnetized electron mobility, m e is the electron mass, ν m is the electron momentum transfer collision frequency, Ω = ω B /ν m is the Hall parameter, and ω B = eB/m e is the electron gyrofrequency.

Figure 3 shows the plasma properties averaged in the y direction and over a 5 µs period in the steady state. The changes in the Xe + density as α t is varied while keeping α d = 20% can be seen in Fig. 3(a). There is no notable change in the doubly charged ion density as α d = 20 % is used for all cases, as shown in Fig. 3 ion density increases proportional to α t . While the axial electric field, shown in Fig. 3(d), is slightly damped, the overall structure is not modified and the difference at the peak is within 4% with increasing fraction of triply charged ions. Since the contribution from the pressure gradient, i.e., diffusion flux, is smaller than the drift flux due to the electric field [START_REF] Hara | One-dimensional hybrid-direct kinetic simulation of the discharge plasma in a Hall thruster[END_REF][START_REF] Hara | Mode transition of a Hall thruster discharge plasma[END_REF] within the acceleration and plume regions, the results shown in Fig. 3 indicate that the cross-field electron mobility is proportional to |u e,⊥ |.

Figure 4 shows the time-averaged, azimuthally-averaged axial electron velocity, u e,⊥ , for different triply charged ion fractions. The electron bulk velocities shown here are obtained by averaging in the y direction and over a 5 µs period in the steady state between 15 and 20 µs. Note that the values are negative as the electrons travel toward the anode (to the left of the computational domain). The baseline case, α d = 20% and α t = 0%, is identical to the results in Ref. 32. As the fraction of triply charged ions is increased by varying α t to 5% and 10%, a consistent reduction in the absolute value of axial electron velocity is observed. While the electron bulk veloc-ity in the upstream, e.g., x ∈ [0 cm, 0.3 cm], is not affected by the triply charged ions, the decrease in the electron bulk speed is approximately 20% when α t = 10% from the case when α t = 0.

Since the bulk plasma properties remain virtually unchanged to this incremental change in the triply charged ion fraction, as discussed in Fig. 3, it can be concluded that the decreased axial electron bulk velocity is caused by the reduction in the anomalous electron mobility. It is to be noted that the net ion current density in the system is constant, as discussed in Section II. While it was shown in Ref. 32 that the generation of IITSI due to Xe 2+ and Xe + in addition to the ECDI enhances the cross-field electron mobility, Fig. 4 shows that the addition of Xe 3+ results in a slight reduction of the enhanced cross-field electron mobility. It was discussed in Ref. 32 that the cross-field electron flux density in the presence of plasma turbulence can be written as, Γ e,x = n e E y /B z , where n e and E y are the fluctuations of the electron density and electric field, respectively, and the bracket denotes temporal and spatial average. The reduced electron transport across magnetic field is This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. therefore likely to be attributed to the dampened electric field amplitudes, shown in Fig. 2. The results indicate that the amplitude of the microturbulence developed due to the coupling between ECDI and IITSI affects the cross-field electron transport in the presence of ion species streaming with different velocities.

PLEASE CITE THIS ARTICLE AS

D. Nonlinear ion wave trapping

The charged particles can be trapped within the axial plasma wave initiated by the IITSI and the azimuthal plasma wave initiated by the ECDI, resulting in a non-Maxwellian ion velocity distribution (IVDF). Here, the axial and azimuthal IVDFs of Xe + , Xe 2+ , and Xe 3+ are presented.

Axial IVDF

The bulk velocity of ions that accelerate across a constant discharge voltage V d can be calculated as

U Z+ x = 2ZeV d m i 1/2 , ( 2
)
where Z = 1, 2, and 3 for Xe + , Xe 2+ , and Xe 3+ respectively, e is electron charge, m i is ion mass. For V d = 200 V, the ion bulk velocities U + x ≈ 1.7 × 10 4 m/s, U 2+ x ≈ 2.4 × 10 4 m/s, and U 3+

x ≈ 2.96 × 10 4 m/s. The maximum velocity attained by the ions is limited by the applied voltage in the absence of any axial oscillations. However, some ions gain energy from FIG. 4: Cross-field electron transport is reduced when increasing triply charged ion species fraction. Plasma properties are averaged in y direction and over 5 µs from t = 15 µs -20 µs. For all cases, α d = 20% is used and the results for α t = 0, 5%, and 10% are shown.

the particle-wave interaction and have energy exceeding the limit imposed by the applied voltage.

Figure 5 shows IVDFs in x-v x plane in the presence of 20% Xe 2+ and 10% Xe 3+ . Ion VDFs for the x-component of the velocity (v x ) is shown in Figs. 5(a), (b), and(c) for Xe + , Xe 2+ , and Xe 3+ , respectively. The particles are binned into a discretized phase space with ∆x = 50 µm and ∆v = 100 m/s. In addition, Fig. 5(d) shows the IVDFs in the downstream, averaged from x = 2 cm to x = 2.3 cm.

Despite the presence of triply charged ions, the linear dispersion relation suggests that the growth rate and real frequency are not significantly modified. This is primarily because the doubly charged ion density is larger than triply charged ion density and therefore the resonant condition is mostly dictated by the interaction between singly and doubly charged ion streams. Similar to the observations in Ref. 32, the phase velocity, v φ , of axial plasma wave can be estimated as v φ = ω/k x ≈ U + x + c s , where ω and k are angular frequency and wave number of the plasma wave and c s is ion acoustic speed. Figures 5(a)-(c) show the trapping of charged species around the phase velocity, v φ ≈ 2.15 × 10 4 m/s. As U +

x ≈ 1.7×10 4 m/s, this means that c s ≈ 4.5×10 3 m/s, which is consistent with the ion acoustic speed evaluated using 30 eV electron temperature observed around x ≈ 1 cm, where IITSI starts to grow [START_REF] Hara | Cross-field electron diffusion due to the coupling of drift-driven microinstabilities[END_REF] . The plasma wave results in Xe + ions that are higher energy than U +

x and a significant portion of Xe 2+ to be decelerated, which are similarly observed in Ref. 32. The most notable result in the present study is that the trapping of Xe 3+ can be seen from Fig. 5(d). It is interesting to note that the axial plasma wave results in (i) Xe 3+ particles that are trapped around the phase velocity, forming a plateau This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. in the Xe 3+ IVDF, and (ii) Xe 3+ particles with higher velocities than U 3+ x , i.e., the ballistic velocity predicted in Eq. ( 2), shifting the peak of the IVDF to a higher velocity.

PLEASE CITE THIS ARTICLE AS

To quantitatively estimate the efficiency of the ion acceleration, one can calculate the bulk velocity of the ions evaluated from the first moment of the IVDFs, i.e., u = ∞ 0 v f (v)dv, where v is the particle velocity and f is the normalized VDF. The bulk velocities of Xe + , Xe 2+ , and Xe 3+ can be obtained from the IVDFs, shown in Fig. 5(d),as u + x = 1.61 × 10 4 m/s, u 2+

x = 2.06 × 10 4 m/s, and u 3+ x = 2.62 × 10 4 m/s, respectively. Using the ballistic velocities as shown in Eq. ( 2), if the beam acceleration efficiency is defined to be η b = (u x /U x ) 2 , one obtains η b = 0.94, 0.86, and 0.88, for Xe + , Xe 2+ , and Xe 3+ , respectively. The lower beam acceleration efficiencies for Xe 2+ and Xe 3+ ions are consistent with the axial deceleration and collisionless thermalization due to the plasma wave, as shown in Fig. 1.

Azimuthal IVDF

Figure 6 shows the IVDFs in y-v y plane of the three ion species considered in the simulation at two axial locations. The left column, Fig. 6(a), corresponds to the particles in the source region, i.e., 0.5 ≤ x ≤ 1 cm, and the right column, Fig. 6(b), corresponds to the particles in the plume region, i.e., 1.4 ≤ x ≤ 2.5 cm. The condition for nonlinear particle trapping within a sinusoidal potential wave can be written as,

|v φ -v Z+ x | ≥ ∆v = 2Zeφ 0 m i 1/2 , (3)
where v Z+ x is the particle velocity for the different charge state Z = 1, 2, and 3, and φ 0 is the amplitude of the plasma wave. Thus, Eq. (3) shows that the higher charge states are more likely to be trapped within the potential waves.

Figure 6(a) shows that most ions are not fully trapped by the azimuthal plasma wave (E y ) within the source region. The trapping frequency [START_REF] Hara | Quantitative study of the trapped particle bunching instability in langmuir waves[END_REF] is given by

ω b = k 2Zeφ 0 m i = k 2ZeE 0 m i , ( 4
)
where k is the wavenumber and E 0 = φ 0 /k is the amplitude of the electric field fluctuation assuming a sinusoidal potential wave. From Fig. 1, E 0 ≈ 2×10 4 V/m and k = 2π/λ y = 9.55× 10 3 rad/m in the source region, where λ y is the wavelength of the plasma wave. This yields ω b ≈ 1.7 × 10 7 rad/s. Hence, the characteristic time for azimuthal particle trapping, τ b = 2π/ω b is on the order of 0.4 -1 µs. On the other hand, the characteristic time for ion transit, τ t = L/v i , is on the order of 0.1 µs. Hence, as τ t < τ b , there is not sufficient time to be fully trapped within the source region. It is to be noted that the Xe 3+ show the formation of high velocity components, This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. which can be seen in Fig. 6(a3), because of the fastest trapping frequency (Z = 3), as shown in Eq. ( 3). However, such triply charged ions are not fully trapped either in the source region.

PLEASE CITE THIS ARTICLE AS

Figure 6(b) shows the azimuthal IVDFs in the plume region. The most interesting result is that the ion trapping feature is not observable for Xe + , as shown in Fig. 6(b1), while Xe 2+ and Xe 3+ are fully trapped within the plasma waves, as shown in Figs. 6(b2) and (b3), respectively. It is to be noted that this feature is consistently found in other conditions despite the amount of the Xe 3+ ions. It can be seen from Figs. 6(b2) and (b3) that the phase velocity of the azimuthal plasma wave is v φ ≈ 6 × 10 3 m/s, which is consistent with the ion acoustic speed evaluated using the electron temperature of 50 eV observed in x ∈ [0.5 cm, 1 cm]. Assuming a potential fluctuation amplitude of 10 V and using Eq. (3), the particles that are trapped within the potential wells must have a trapping velocity of ∆v = 3.5 × 10 3 m/s, 5.4 × 10 3 m/s, and 6.6 × 10 3 m/s for Xe + , Xe 2+ , and Xe 3+ , respectively. It can therefore be concluded that v φ > ∆v for singly charged ions, leading to virtually no ion trapping. On the other hand, doubly and triply charged ions are trapped in the azimuthal direction as v φ ≤ ∆v is satisfied.

The azimuthal IVDFs for the cross-field plasma in the presence of multiply charged ions are notably different from previous kinetic simulations where only singly charged ions are considered [START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF][START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. i. insights from particle-in-cell simulations[END_REF] because the only trappable ions are singly charged ions if multiply charged ions are not present in the simulation. The present results indicate that the IVDF measurements of Xe + , e.g., using laser induced fluorescence, may not be able to detect such ion trapping phenomena. Based on the present kinetic simulations, Xe 2+ and Xe 3+ are more likely to exhibit such non-Maxwellian behavior. However, it is also to be noted that the ionization mechanism is simplified and is not self-consistently modeled in the present kinetic simulations. Thus, the degree of ion trapping may change depending on the location of the ion production within the discharge plasma. Additionally, as the present results show that the amplitude of the multidimensional plasma waves affects the electron transport, three-dimensional kinetic simulations are required to understand the electron transport in real sys-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0067305 tems.

IV. CONCLUSIONS

It was shown in Ref. 32 that singly and doubly charged ion streams excite the IITSI, which, through its nonlinear interaction with the ECDI, increases electron transport in the crossfield direction. In the present work, microturbulence developed due to these interactions between the ECDI and IITSI is now investigated accounting for a third (triply ionized) ion population, using a multidimensional kinetic simulation. Simulations were performed for several fractions of triply charged xenon ions while fixing the fraction of Xe 2+ and other plasma characteristics, which sets up the ECDI and a large-amplitude IITSI.

The simulation results show that the IITSI driven by the relative velocities between singly, doubly, and triply charged xenon ions interacts nonlinearly with the ECDI. The mode coupling between ECDI and IITSI is shown to decrease the amplitudes of both axial and azimuthal plasma waves as the fraction of triply charged ions is increased. The 2D dispersion relation accounting for singly, doubly, and triply charged ion species and magnetized electrons, predicts a negligible effect of adding ≤ 10% Xe 3+ on the linear growth rate. However, the nonlinear saturation of different instabilities plays an important role in how large the amplitudes of plasma waves become. In particular, these results indicate approximately 30% damping in azimuthal amplitude (in the direction of E × B drift) and nearly 8% in axial amplitude (in the direction of applied E field and across the B field) when 10% triply charged xenon ions are added to the simulation. With an incremental increase in the fraction of triply charged ions from 0 to 10%, the crossfield electron transport is shown to reduce accordingly.

Ion species streaming with different velocities interact nonlinearly with the plasma wave and undergo trapping and heavily influence the saturation of these instabilities. Trapping of the three ion species near the source and in the plume shows that when amplitude of plasma wave driven by the IITSI is large enough, resonant particle-wave interactions take place and higher charge state ions are deeply trapped. This leads to broadening of the axial ion VDFs of all three species considered, where the Xe 3+ VDF is shown to undergo the largest broadening due to the separation between the particle velocity and the phase velocity of the axial wave. In addition, the present simulations indicate that the singly charged ions may not be fully trapped within the azimuthal plasma wave driven by the ECDI. The multiply charged ions are more trapped because the trapping velocity scales as the charge state Z, which may play an important role in experimentally detecting such ion trapping features. The results of this paper are relevant to partially magnetized devices in which a diversity of ion charge states are generated. These findings show the extent to which the presence of multiple charge states influences the nonlinear behavior of key instabilities and their contribution to anomalous transport.

FIG. 1 :

 1 FIG. 1: Effect of increasing triply charged ion concentration on instantaneous azimuthal (E y ) and axial (E x ) electric field due to coupling between ECDI and IITSI. Top row: α d = 20% and α t = 0; bottom row: α d = 20% and α t = 10%. Results are shown for t = 18 µs, i.e., after the steady state is reached. The colorbar is saturated for better visualization of axial fluctuations in the downstream region, where min and max on the color scale correspond to -2.0 × 10 4 V/m and 2.0 × 10 4 V/m, respectively.

FIG. 2 :

 2 FIG. 2: Damping of ECDI and IITSI in the presence of singly, doubly, and triply charged ions. α d = 20% is fixed while α t is varied. Color plots show 2D FFT of time-averaged azimuthal (E y ) and axial (E x ) fields in steady state in the (k x , k y ) plane: (a1) FFT(E y ) with α t = 0%, (a2) FFT(E y ) with α t = 10%. (b1) FFT(E x ) with α t = 0%, and (b2) FFT(E x ) with α t = 10%. Min and max values on the colorbar correspond to 100 and 1500 for (a1) and (a2), respectively, and 100 and 4000 for (b1) and (b2), respectively. The line plots (a3) and (b3) are obtained from the 2D data by integrating with respect to k x and k y respectively: (a3) shows the variation in the azimuthal electric field amplitude as a function of k y ; (b3) shows the variation of the axial electric field amplitude as a function of k x .

FIG. 3 :

 3 FIG. 3: Spatially and temporally averaged plasma properties. All quantities are averaged in y direction and over 5 µs period in steady state. (a) Xe + density, (b) Xe 2+ density, (c) Xe 3+ density, and (d) axial electric field. All the physical and numerical parameters are fixed while α d = 20 % and α t is varied. Three cases corresponding to α t = 0, 5 and 10% are shown.

FIG. 5 :

 5 FIG. 5: Instantaneous ion velocity distribution function (IVDF) for the x-component (axial) of the velocity, averaged in the y direction is shown for (a) Xe + (b) Xe 2+ , and (c) Xe 3+ . The horizontal dashed lines indicate the ion bulk velocities assuming that the ions are accelerated across a constant voltage drop of 200 V. The colormap shows the VDFs normalized with the reference VDF values, which is chosen as f + max for Xe + , 0.1 f + max for Xe 2+ , and 0.01 f + max for Xe 3+ where f + max is the maximum value of Xe + VDF. Colorbar uses logarithmic scale, where min and max correspond to 1.0 and 1.0 × 10 -4 respectively. Shown in (d) are the axial VDFs in the plume (2.0 ≤ x ≤ 2.3 cm) for the three species.

FIG. 6 :

 6 FIG. 6: Instantaneous ion velocity distribution function (IVDF) for the y-component (azimuthal) of the velocity, averaged in the x direction is shown for Xe + in the top row, Xe 2+ in the middle row, and Xe 3+ in the bottom row. Two regions in the axial domain corresponding to (a) source (0.5 ≤ x ≤ 1 cm) and (b) plume (1.4 ≤ x ≤ 2.5 cm) are shown. Colormap is identical to Fig. 5.

ACKNOWLEDGMENTS

This material is based on work supported by the Air Force Office of Scientific Research under Award No. FA9550-18-1-0090 and by the US Department of Energy, Office of Science, Office of Fusion Energy Sciences, under Award No. DE-SC0020623. Simulations are performed using the high performance computing resources at the Stanford Research Computing Center.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.


Rexhep Durmo

Benedetta Donati

Louis Rebaud

Anne Segolene Cottereau

Alessia Ruffini

Maria Elena Nizzoli

Sabino Ciavarella

Maria Carmela Vegliante

Christophe Nioche

Michel Meignan

Francesco Merli

Annibale Versari

Alessia Ciarrocchi

Irene Buvat

Stefano Luminari

and dacarbazine-treated, interimPET-negative classical Hodgkin Lymphoma patients: A radio-genomic study

Introduction

Classical Hodgkin lymphoma (cHL) is a lymphoid neoplasm characterized by varied disease presentation and is accounted as one of the most curable cancers in humans. [START_REF] Thomas | Epidemiology and etiology of Hodgkin's lymphoma[END_REF] The current approach to cHL patients requires an accurate assessment of the trade-off between the intensity of the initial treatment and its safety. ABVD (doxorubicin, bleomycin, vinblastine, and dacarbazine) combination with or without radiotherapy is considered a standard regimen for HL due to its excellent tolerability. However, approximately 30% of patients relapse or are refractory to ABVD therapy requiring intensive and toxic salvage therapies. [START_REF] Viviani | ABVD versus BEACOPP for Hodgkin's lymphoma when high-dose salvage is planned[END_REF] The early identification of patients who are a high risk of treatment failure currently represent an open research question. [START_REF] Merli | Long-term results of the HD2000 trial comparing ABVD versus BEACOPP versus COPP-EBV-CAD in untreated patients with advanced hodgkin lymphoma: a study by Fondazione Italiana Linfomi[END_REF] Fluorine-18-fluorodeoxyglucose positron emission tomography/computed tomography (18F-FDG PET/CT) is a non-invasive tool useful for staging and treatment response purposes in HL. [START_REF] Cheson | Recommendations for initial evaluation, staging, and response assessment of hodgkin and non-hodgkin lymphoma: the lugano classification[END_REF] Interim FDG PET/CT (iPET) performed after 2 cycles of chemotherapy is currently used for early assessment of response and to tailor treatment intensity. [START_REF] Borchmann | Progression-free survival of early interim PETpositive patients with advanced stage Hodgkin's lymphoma treated with BEACOPPescalated alone or in combination with rituximab (HD18): an open-label, international, randomised phase 3 study by the German Hodgk[END_REF][START_REF] Kobe | Outcome-based interpretation of early interim PET in advanced-stage Hodgkin lymphoma[END_REF][START_REF] Straus | CALGB 50604: risk-adapted treatment of nonbulky early-stage Hodgkin lymphoma based on interim PET[END_REF][START_REF] Radford | Results of a trial of PET-directed therapy for earlystage Hodgkin's lymphoma[END_REF] Although highly prognostic, the accuracy of iPET as predictive factor is still suboptimal mainly due to the non-negligible risk of relapse observed among iPET negative (iPET-) patients. [START_REF] Stephens | Five-year follow-up of SWOG S0816: limitations and values of a PET-adapted approach with stage III/IV Hodgkin lymphoma[END_REF] Thus, clinical research is currently aimed at identifying additional prognostic factors among those available at time of diagnosis, that may help to increase the accuracy of individual risk prediction. Radiomic features measured at baseline by PET like metabolic tumor volume (MTV), total lesion glycolysis (TLG) and maximum standardized uptake value (SUV-max) have been proposed as promising prognostic biomarker in lymphoma. [START_REF] Mettler | Metabolic tumor volume for response prediction in advanced-stage hodgkin lymphoma[END_REF] However, a definitive consensus has not been achieved. [START_REF] Schöder | Metabolic tumor volume in lymphoma: hype or hope?[END_REF][START_REF] Guo | Prognostic value of baseline metabolic tumor volume and total lesion glycolysis in patients with lymphoma: a meta-analysis[END_REF] Among them, the distance between the two lesions that are the furthest apart (Dmax), that reflects lesion dissemination, has gained momentum based on new promising results. Indeed a recent study reported that in Diffuse Large B Cell Lymphoma (DLBCL) the combination of Dmax and MTV identified a group of high-risk patients before treatment. [START_REF] Cottereau | 18 F-FDG PET dissemination features in diffuse large B-cell lymphoma are predictive of outcome[END_REF][START_REF] Cottereau | Risk stratification in diffuse large B-cell lymphoma using lesion dissemination and metabolic tumor burden calculated from baseline PET/CT[END_REF] In order to assess the prognostic impact of Dmax in cHL, we conducted a retrospective study in a series of consecutive newly diagnosed patients with confirmed cHL who were treated at our institution with ABVD chemotherapy. We also explored the molecular bases underlying tumor dissemination defined by Dmax through a gene expression analysis on diagnostic biopsies of cHL patients.

Materials and methods

Patient selection

This study was conducted as a monocentric retrospective analysis. We included patients with confirmed cHL diagnosed in our institution between 2007 and 2020, who were initially treated with ABVD, who had available baseline FDG-PET (bPET) and a minimum follow up of 6 months from diagnosis. Upon centralized review patients with a single lesion at baseline PET were excluded. The main study endpoint was progression free survival (PFS).

This study was conducted in accordance with Declaration of Helsinki and approved by the Internal Review Board. Written consent was obtained from all living patients.

FDG PET/CT imaging and interpretation

bPET and iPET for all patients were revised by two expert nuclear medicine physicians (RD, AV) who were blinded to patient outcome. PET scanning was performed according to the Society of Nuclear Medicine and Molecular Imaging (SNMMI) and European Association of Nuclear Medicine (EANM) guidelines. [START_REF] Delbeke | Procedure guideline for tumor imaging with 18F-FDG PET/CT 1.0[END_REF][START_REF] Boellaard | FDG PET/CT: EANM procedure guidelines for tumour imaging: version 2.0[END_REF] MTV was calculated for bPET using Beth Israel PET/CT viewer plugin for FIJI (http://petctviewer.org) and defined using the 41% SUVmax threshold. Total lesion glycolysis was calculated mathematically as the sum of MTV x SUVmean of each lesion. Bone marrow involvement was included in the volume measurement only if there was focal uptake as previously described. [START_REF] Adams | Systematic review and meta-analysis on the diagnostic performance of FDG-PET/CT in detecting bone marrow involvement in newly diagnosed Hodgkin lymphoma: is bone marrow biopsy still necessary?[END_REF] The centroid was obtained from the three-dimensional coordinates of MTV of each lymphoma lesion. The distances between all pairs of lesions (including both nodal and extra nodal) were calculated by Euclidean formula with LIFEx software. [START_REF] Nioche | LIFEx: a freeware for radiomic feature calculation in multimodality imaging to accelerate advances in the characterization of tumor heterogeneity[END_REF] The largest lesion distance was deduced in each patient (Dmax).

iPET was classified according to Deauville Score (DS) as originally reported [START_REF] Meignan | Report on the first international workshop on interim-PET scan in lymphoma[END_REF] and interpreted according to Lugano criteria.4 Body Surface Area (BSA) was calculated according to the Du Bois method. [START_REF] Shuter | Body surface area: Du Bois and Du Bois revisited[END_REF] We evaluated the gene expression profiles by NanoString using the PanCancer Immune Profiling Panel (NanoString Technologies) as previously described [START_REF] Luminari | A gene expression-based model to predict metabolic response after two courses of ABVD in hodgkin lymphoma patients[END_REF][START_REF] Donati | Gene expression profile unveils diverse biological effect of serum vitamin D in Hodgkin's and diffuse large B-cell lymphoma[END_REF] starting from N = 5, 5 µm FFPE slides including whole tissue.

We compared the expression profiles of patients with low-and high-Dmax. pValue (as twotailed Student's t test) and the false discovery rate (FDR) obtained by the Benjamini-Hochberg method were calculated. Bioinformatic analyses were conducted by R Software v4.0.4 packages.

Gene expression profiles data are available at the Gene Expression Omnibus (GEO) repository (accession number: GSE184662).

CIBERSORT analysis

CIBERSORTx algorithm (https://cibersortx.stanford.edu/) [START_REF] Newman | Robust enumeration of cell subsets from tissue expression profiles[END_REF] was used to perform deconvolution of bulk expression data derived from 126 cHL. A 1016-gene customized signature matrix was built as previously described [START_REF] Ciavarella | Dissection of DLBCL microenvironment provides a gene expression-based predictor of survival applicable to formalin-fixed paraffinembedded tissue[END_REF] (Supplementary Table S1). Bulk-mode batch correction was applied to mixture samples before imputing cell fractions and 1000 permutations were set for significance analysis. Patients were stratified based on Dmax median value (20 cm).

Significantly different infiltrate percentage between the two groups was evaluated for each cell type by Mann-Whitney U test.

Statistical analysis

Associations between clinical variables and radiomic features and Dmax values were evaluated using Fisher's exact test for univariate analysis and by generalized linear model for multivariate analysis.

OS and PFS were estimated with Kaplan-Meier curves and log-rank test was used to assess statistical significance. Cox proportional-hazards model was used for univariate and multivariate survival analysis. For internal validation of Dmax cutoff we first applied Leave One Out cross validation and we calculated the best cutoff of Dmax, TMTV, TLG and SUVmax to predict PFS with two alternative methods: (1) by using ROC-AUC at 36 months of follow-up and selecting the cut-off by means of Youden index 25 ; (2) by using maximally selected log-rank test [START_REF] Lausen | Maximally selected rank statistics[END_REF] . Both the results were confirmed by 1000 bootstrap resamples. Statistical analysis was performed by R Software 4.0.4.

Results

Patients characteristics

A retrospective consecutive cohort of 183 cHLs treated with ABVD regimen was initially identified. bPET scan was available for 180 of the selected patients. Furthermore, 25 patients had a single lesion and were excluded. Thus 155 patients were fully eligible (Figure 1A). Baseline characteristics of eligible patients are reported in Supplementary Table S2.

All patients were evaluated with iPET after 2 cycles of ABVD treatment. iPET was positive for 27 patients (17.4%) and negative for the remaining 128 cases (82.6%). Among the iPET positive patients, treatment was intensified in 10 cases (37.0%), of whom 8 received BEACOPP and 2 IGEV followed by ASCT. A total of 50/155 patients (32.3%) received consolidation radiotherapy (RT) after systemic chemotherapy (Supplementary Table 3). Response to the induction treatment defined by final PET (fPET) was available in 153 patients and was confirmed as CMR in 133 patients (86.9%) and non-CMR in 20 (13.1%).

Radiomic features based on bPET were evaluated and are reported in Supplementary Table S2.

Median SUVmax was 14 (IQR: 10-17). Median baseline MTV was 100 cm3 (IQR: 47-184).

Median baseline TLG was 583 (IQR: 267-1265). The agreement between two operators who revised PET images was excellent (Pearson r = 0.92, p < 0.001).

Analysis of Dmax and correlation with baseline patient features

Median Dmax was 20 cm (IQR: 11-38). Patients were classified as low-Dmax (Dmax ≤20 cm) and high-Dmax (Dmax >20 cm). Dmax value was not influenced by BSA, as demonstrated by correlation plot (R = 0.99, p < 0.001) (Supplementary Figure 1A-B). In univariate analysis features correlated to high Dmax values were male gender, low serum albumin, low LDH, high ESR and stage. We also observed a direct correlation of Dmax with MTV and TLG. By contrast, no association between Dmax and SUVmax or iPET response was observed (Table 1).

Multivariate logistic analysis confirmed the correlation with high Dmax values only for MTV (p = 0.003). Different median Dmax values were observed for patients with early or advanced stage defined according to EORTC criteria: 12 cm (IQR: 8-18) and 37 cm (IQR: 22-54), respectively. Distribution of Dmax according to stage is reported in Supplementary Figure 1C.

Dmax correlation with survival probability

Median follow-up of the cohort was 63 months (range, 6.5-139 months). 39 patients experienced an event for PFS including 29 disease progressions and 10 deaths. The 4-year PFS rate was 77.7% (95%CI 71.1-84.9). Considering Overall Survival (OS), 13 patients died of whom 4 due to lymphoma progression and 9 due to other causes resulting in a 4-year OS rate of 92.9% (95% CI 88.5-97.5).

To establish the prognostic value of radiomic features in cHL, we evaluated the risk of disease progression (as PFS) associated to the median value of each parameter by Cox regression hazard model. In univariate analysis, no significant correlation was observed for MTV, TLG and SUVmax for the whole series. A significant correlation between MTV and PFS was observed when analysis was limited to patients younger than 60 years and adjusted for the main clinical risk factors (HR 2.66, 95%CI 1.11-6.38). By contrast, high Dmax was significantly associated with increased risk of disease progression (HR = 2.59, 95%CI 1.31-5.12) (Table 2). Indeed, patients with low Dmax had a 4 years PFS of 84.3% (95%CI 76.0-93.4), significantly better than patients with high Dmax (4 years PFS: 71.1%, 95%CI to 82.5%) (Figure 1C). Dmax confirmed its prognostic role also when considered as a continues variable (HR 1.54, 95%CI 1.03-2.32) and also when corrected for BSA (Supplementary Figure 1B). Among clinical variables, advanced stage (AA III-IV) (HR = 1.93, 95% CI 1-3.73) and iPET positivity (HR = 1.95, 95% CI 0.97-3.93) showed a significant association with reduced PFS (Table 2). Multivariate analysis for PFS showed that Dmax was the only variable remaining associated with disease progression, confirming the independent prognostic role of this metabolic parameter (HR = 2.70, 95%CI 1.1-6.63). These results were confirmed when the multivariate model was corrected for RT consolidation and for treatment intensification in iPET + patients (HR = 2.50, 95%CI 1.20-5.21) (Table 2).

Combined analysis of Dmax and iPET

We evaluated the prognostic role of Dmax for patients who achieved a negative iPET and continued their treatment with ABVD chemotherapy. In this group high Dmax significantly increased the risk of progression (HR = 4.53, 95%CI 1.81-11.3) and this association was confirmed in a multivariate analysis incuding Stage (III-IV) as covariate (HR = 4.98, 95%CI 1.64-15.16). These results were also confirmed when the multivariate model was corrected for consolidation radiotherapy (HR = 3.96, 95%CI 1.52-10.33) (Table 3).

Combining Dmax and iPET we observed that Dmax was able to identify patients at different risk of progression only among iPET-patients, while no additional prognostic role was seen for the small group of iPET + cases. Thus, combining iPET and Dmax we were able to identify three risk groups (Figure 1D). Sixty-three (41%) patients had iPET-and low Dmax and a 4 years PFS of 90% (95%CI 82.0-98.9) significantly better compared to the 65 (42%) patients with iPETand high Dmax (4 years PFS 72.4% (95%CI 61.9-84.6)). The 27 (17%) patients with a positive iPET had a 4 years PFS of 62.3% (95%CI 46.3-83.9).

Internal validation of Dmax prognostic cut-off

To avoid data overfitting, due to the limited sample size and a low number of events, we used the median value of our cohort to analyze the prognostic role of Dmax. However, we applied different approaches to provide an internal validation of our observations and to strength the reliability of this Dmax cut-off.

First, Leave One Out (LOO) internal validation showed a repeatable performance of this cut-off, indeed original Dmax >20 cm and LOO Dmax>20 cm had HR of 2.59 (95% CI 1.31-5.11) and 2.38 (95% CI 1.18-4.78) with a Gohen/Heller Concordance of 0.61 (SE 0.04) and 0.59 (SE 0.04), respectively. Then, with the Dmax as continuous value, we evaluated the ROC-AUC at 36 months of follow-up (3 years) in the overall cohort (AUC = 0.57) and we selected the cutpoint of 19.4 by means of Youden index (Figure 2A). After 1000 bootstrap resamples the 2.5-97.5 percentiles were 16.6-20.7 with a bias of 0.40. Furthermore, we took advantages of a second approach to define Dmax cutoff on the overall cohort based on maximally selected logrank test. Even in this case we obtained an estimated cut point of 19.4 cm (p = 0.029), that was confirmed also after 1000 bootstrapped resamples (Bias -0.0036, 2.5-97.5 percentiles were 10.4-21.8) (Figure 2B andC). All together these analyses indicate that the Dmax >20 cm is an accurate threshold within our cohort. A similar approach was used to validate cut-offs for TMTV, TLG and SUVmax (Supplementary Figure 2).

Gene expression profiling and CIBERSOT analysis

We employed digital profiling for investigating the biological bases of Dmax. 126 of the 155 cHLs analyzed samples had biological material suitable for the analysis. Of these, 61 (48.4%) were low-Dmax and 65 (51.6%) were high-Dmax (Figure 1A). 556 of the evaluated genes were expressed in at least one of the two groups. Differential analysis of the gene expression profiles between the high-and low-Dmax groups resulted in 154 differentially expressed genes, of which 46 (29.8%) upregulated and 108 (70.2%) downregulated in high-Dmax (Figure 3A). After p-value adjustment, 61 genes remained differentially distributed between the two groups. Principal Component Analysis (PCA) showed that gene expression profiles partially segregated high and low-Dmax samples (Figure 3B). 21 genes were up-regulated and 40 down-regulated in high-Dmax versus low-Dmax samples (Figure 3C).

Gene Ontology analysis showed that down-regulated genes in high Dmax samples were particularly enriched in processes linked to T-and B-cells homeostasis, activation and mediated cell response (Figure 4A). Of note, among the others, expression of B-cell specific markers like CD20 (MS4A1) and CD79A-B and T-cell specific factors like FOXP3 resulted significantly downregulated in high-Dmax patients. By contrast, among the up-regulated genes we observed the macrophage-specific marker CD68 (log2(FC) = 0.41, adj.pValue = 0.05) and several macrophage-associated factors including CXCL9, CCL24, CD163, CD14 and MRC1 (Figure 3C).

A deconvolution approach was used to infer the cell populations associated to Dmax on the basis of gene expression profiling.23 HL cells, B cells, CD8+ and CD4+ T cells were the predominant cell populations in our samples. Also, vascular components including endothelial cells, pericytes and myofibroblasts were well represented. Differential analysis of the distribution of these cell populations showed that naïve B and T cells were strongly enriched in the low-Dmax subgroup.

By contrast CD8+ T cells expressing high level of PDL1, dendritic cells, monocytes and eosinophil were enriched in the high-Dmax samples. Even if poorly represented in the tumor microenvironment, Cancer Associated Fibroblasts (CAFs) resulted also to be preferentially associated with low Dmax (Figure 4B and Supplementary Table 4). No differences in the distribution of HL cells was observed between high and low-Dmax cases. We also explored how the 7 Dmax-associated cytotypes correlate with other PET parameters (Supplementary Figure 3A-D). The median value of TMTV and TLG was significantly associated with 5 and 4 cell subsets, respectively, whereas SUVmax only correlated with monocytes. Moreover, the whole panel of cytotypes showed a prognostic trend in term of PFS, which was statistically significant for naïve B cells and CD8+/highPDL1 T cells, in a way consistent with their association to Dmax (Supplementary Figure 3E-F).

Discussion

In this study we evaluated the prognostic role of tumor spread defined by Dmax, a simple measurement of the largest distance between lymphoma sites, in a retrospective cohort of consecutive newly diagnosed patients with cHL treated with ABVD. We were able to demonstrate the independent correlation of Dmax with the risk of progression for all patients included in the study and for the subgroup of patients who achieved a CMR early during the treatment course with a five-fold increase in the risk of disease progression associated with high Dmax compared to low Dmax.

Our study should be considered as a proof of concept about the prognostic role of Dmax in cHL that may have some implications related to the initial evaluation of patients and for the future development of risk adapted therapies in cHL. In spite of the relevant progresses gained in the recent years in the understanding and evaluation of disease features, still staging and consequential management of cHL are relying on historical qualitative criteria, that poorly account for the complex biology of this tumor.

Compared to standard Ann Arbor staging that failed to predict outcome in our study, as well as in other studies, [START_REF] Cottereau | 18 F-FDG PET dissemination features in diffuse large B-cell lymphoma are predictive of outcome[END_REF] Dmax appeared more effective for risk stratification. Indeed, according to our analysis, Dmax contributes to a more precise evaluation of tumor spreading compared to Ann Arbor index, in particular for patients at intermediate risk.

Diagnostic imaging features collected at diagnosis by PET are regarded as potential source of relevant prognosis indicators. We showed that, differently from other radiomic parameters, including MTV, TLG and SUVmax, whose prognostic value was not confirmed by our analysis, Dmax had an independent association with the risk of progression and was confirmed in the subgroup of early responding patients. Noticeably, Dmax is a very simple feature to calculate, that is measured as the distance between centroids of two lesions, that is not impacted by the scanner or the operator, differently from the other PET-derived features.

Regarding the contribution of Dmax to response adapted therapy in cHL patients, a major finding of our analysis was the evidence that Dmax is an independent prognostic factor for PFS among the iPET-patients. Indeed, an unsolved need is currently identified in the non-negligible proportion of early responding patients who still experience relapse after completion of front-line therapy, regardless of initial stage. The risk of disease progression reflects some innate features of the tumor that are not captured by standard prognostic factors but can be partly revealed by the analysis of Dmax.

Being an index of tumor spreading, Dmax may be considered as a surrogate of tumor aggressiveness. To map this connection, we employed digital gene expression profiling to investigate the biologic features of Dmax at diagnosis. We chose to not confining the analysis to Reed-Sternberg (RS) cells but to consider whole tumor tissue since, in cHL tumor microenvironment accounts for a significant part of the biology of the disease and its behavior. [START_REF] Liu | The microenvironment in classical Hodgkin lymphoma: an actively shaped and essential tumor component[END_REF] We used the same approach in our previous work demonstrating that microenvironmental dynamics concur to define innate chemorefractoriness to ABVD. [START_REF] Luminari | A gene expression-based model to predict metabolic response after two courses of ABVD in hodgkin lymphoma patients[END_REF] Therefore, our explorative study represents an attempt to link peculiar aspects of cHL microenvironment to Dmax as a measure of the proneness of tumor to spread during progression.

While we recognize the use of a limited panel of genes as a limitation of our analysis, the obtained profiles were able to partly segregate high and low Dmax lesions. Genes associate to low-Dmax, were representative of pathways related to B-cells and T-cells presence and activation (Figure 3B-C). In particular, B-cells specific markers like CD20 (encoded by MS4A1 gene), CD79 A and B were significantly overexpressed in low-Dmax lesions. Also, lineage committing transcription factor like PAX5 [START_REF] Revilla-I-Domingo | The B-cell identity factor Pax5 regulates distinct transcriptional programmes in early and late B lymphopoiesis[END_REF] and FOXP3, [START_REF] Marson | Foxp3 occupancy and regulation of key target genes during T-cell stimulation[END_REF][START_REF] Josefowicz | Regulatory T cells: mechanisms of differentiation and function[END_REF] were significantly overexpressed in low-Dmax cHLs. By contrast, expression of macrophages/monocytes markers including CD68, CD163 resulted associated to high-Dmax lesions, coherent with a general role of these cells in tumor progression. [START_REF] Condeelis | Macrophages: obligate partners for tumor cell migration, invasion, and metastasis[END_REF] Of note the expression of FOXP3, CD68 and CD20 in the tumor microenvironment has been already reported to have a prognostic value in cHLs. [START_REF] Greaves | Expression of FOXP3, CD68, and CD20 at diagnosis in the microenvironment of classical hodgkin lymphoma is predictive of outcome[END_REF] In particular, high expression of FOXP3 and CD20 at diagnosis were reported to be associated to superior OS, while increased levels of CD68 were predictive of increased first-line treatment failure and reduce OS. Our results are coherent with these observations and add a new layer of information suggesting that this peculiar asset of the tumor microenvironment may contribute to cHL progression by creating favorable conditions to tumor dissemination.

Deconvolution analysis by CIBERSORTx confirmed this evidence. Indeed, higher Dmax characterizes those cHL with an under-represented immune microenvironment, particularly in naïve B cells and CD8+ T cells compartments. Interestingly, a subset of CD8+ T effector expressing higher PD-L1 level directly correlated with Dmax, featuring more aggressive disease.

As a critical immune checkpoint, PD-L1 promotes immune-escape [START_REF] Diskin | PD-L1 engagement on T cells promotes self-tolerance and suppression of neighboring macrophages and effector T cells in cancer[END_REF] and, in keeping with higher infiltration of CD68-positive cells, is negatively associated with cHL patients' survival. [START_REF] Agostinelli | The combined role of biomarkers and interim PET scan in prediction of treatment outcome in classical Hodgkin's lymphoma: a retrospective, European, multicentre cohort study[END_REF] This is also in line with the concept that weaker immune response correlates with proneness of distant spreading in many cancer types including B-cell lymphoma. [START_REF] Hiam-Galvez | Systemic immunity in cancer[END_REF][START_REF] De Charette | Hide or defend, the two strategies of lymphoma immune evasion: potential implications for immunotherapy[END_REF] Our results were generated from and should be referred to an unselected population of cHL patients homogeneously treated with ABVD that is the first-choice therapy for the majority of patients around world. The only difference in the treatment protocol was represented by the number of ABVD cycles, by the adoption of consolidation radiotherapy according to available guidelines and by the intensification of treatment for some iPET + patients that was defined on a patient basis and upon multidisciplinary discussion. This choice allowed to include in the study most of the patients seen in our institution and to significantly reduce interferences caused by treatment dependent biases. Indeed, our observation is referred to the majority of patients with the only exclusion of stage I disease who are not eligible for Dmax calculation. However, considering that other options are currently prescribed mostly for patients with advanced disease, larger studies, and a full validation of Dmax in patients treated with regimens other than ABVD are required before our results can be generalized to all cHL patients. Kanoun et al. recently provided confirmatory data about the prognostic role of Dmax among patients initially treated with BEACOPP within the AHL2011 randomized trial. [START_REF] Kanoun | Total metabolic tumor volume and tumor dissemination are independent prognostic factors in advanced hodgkin lymphoma[END_REF] Also, regarding the cut-off value of Dmax that was defined at 20 cm in our study further validation is needed. Indeed, our data should be considered as a proof of concept about the prognostic role of Dmax in cHL patients and to identify subgroups of patients that could have the greater utility from the use of Dmax, that is, iPET negative patients. We acknowledge all the limitations of our study population in term of size and representativeness. However, we also believe that our conclusions are well supported by data. Larger trials with well-defined validation cohort will be needed to define the best cut off for Dmax that will likely be different for different risk subgroups.

Among study results the lack of prognostic value for MTV or TLG was in contrast with available data. [START_REF] Cottereau | Prognostic value of baseline metabolic tumor volume in early-stage Hodgkin lymphoma in the standard arm of the H10 trial[END_REF] Indeed, neither MTV or TLG were confirmed as prognostic factors for the whole patients series in the univariate analysis of PFS and were not analyzed in the multivariate model. The small sample size and the few observed events reported in our study might be a reasonable explanation to this discordancy. Indeed an association between MTV with patients' outcome was found when the analysis was limited to patients younger than 60 years suggesting the presence of biases associated with older patients.

In conclusion, our study suggests Dmax as a relevant feature allowing to better predict patients' outcome among subjects otherwise defined at low risk based on their early response to ABVD chemotherapy. Pending further validation Dmax may be suggested as a new simple assessment that, added to the initial workup of patients, might contribute to a better definition of patient risk profile and to a more precise personalization of therapy in cHL. LMR (lymphocytes/ monocytes ratio), LDH (Lactate dehydrogenase), ULN (upper limit normal), ERS (eritrocytes sedimentation rate), MTV (total tumor metabolic volume), SUV (standardized uptake value), TLG (total lesion glycolysis), Dmax (largest distance between the two most distant lesions). Univariate pValue was calculated as twotailed Fisher exact test. Multivariate pValue, including variables significantly associated with Dmax at univariate analysis, was calculated by generalized linear model (glm). The pValues shown in boldface correspond to p<0.05.

Figure 1 Association

 1 Figure 1

Figure 2 Validation

 2 Figure 2

Figure 3 Gene

 3 Figure 3

Figure 4 Cell

 4 Figure 4

Table 1 :

 1 Association of Dmax with clinical-pathological variables in classical Hodgkin Lymphoma (cHL) patients (n = 155)

		Dmax ≤ 20cm	Dmax > 20cm
		(N=78)	(N=77)	NA	pValue	Multiv. pValue
	Gender			-	0.006	0.077
	Female	47 (60.3%)	29 (37.7%)
	Male	31 (39.7%)	48 (62.3%)
	Age			-	0.640	-
	≤65	69 (88.5%)	66 (85.7%)
	>65	9 (11.5%)	11 (14.3%)
	Age			-	0.624	-
	≤45	49 (62.8%)	45 (58.4%)
	>45	29 (37.2%)	32 (41.6%)
	LMR			1	0.333	-
	≤2.1	33 (42.3%)	39 (51.3%)
	>2.1	45 (57.7%)	37 (48.7%)
	Albumin			19	0.017	0.720
	≤4 g/dl	40 (58.0%)	52 (77.6%)
	>4 g/dl	29 (42.0%)	15 (22.4%)
	Hemoglobin			2	0.102	-
	≤10.5 g/dl	7 (9.0%)	14 (18.2%)
	>10.5 g/dl	71 (91.0%)	61 (81.3%)
	Leukocytes			2	0.059	-
	≤15x10 3 cells/mm 3	74 (94.9%)	64 (85.3%)
	>15x10 3 cells/mm 3	4 (5.1%)	11 (14.7%)
	LDH/ULN			6	0.013	0.233
	≤1	52 (68.4%)	35 (47.9%)
	>1	24 (31.6%)	38 (52.1%)
	ESR			9	0.004	0.677
	≤50	40 (54.8%)	22 (30.1%)
	>50	33 (45.2%)	51 (69.9%)
	STAGE			-	< 0.001	0.114
	I-II	62 (79.5%)	15 (19.5%)
	III-IV	16 (20.5%)	62 (80.5%)
	Bulky			-	0.383	-
	no	63 (80.8%)	67 (87.0%)
	yes	15 (19.2%)	10 (13.0%)
	Risk Group			1	< 0.001	0.446
	Local	59 (76.6%)	13 (16.9%)
	Advanced	18 (23.4%)	64 (83.1%)
	HASENCLEVER			-	< 0.001	0.208
	0-2	59 (75.6%)	32 (41.6%)
	>2	19 (24.4%)	45 (58.4%)
	MTV			-	< 0.001	0.003
	≤ 100 cm 3	57 (73.1%)	20 (26.0%)
	> 100 cm 3	21 (26.9%)	57 (74.0%)
	TLG			-	< 0.001	0.143
	≤ 583	52 (66.7%)	24 (31.2%)
	> 583	26 (33.3%)	53 (68.8%)
	SUVmax			-	0.200	-
	≤ 14	45 (57.7%)	36 (46.8%)
	> 14	33 (42.3%)	41 (53.2%)
	iPET			-	0.673	-

Acknowledgements

We are grateful to Dr. Luigi Marcheselli for his valuable support to the statistical analysis. We also wish to thank Dr. Savonarola for inspiring our writing. This work was supported by Fondazione Italiana per la Ricerca sul Cancro (AIRC-IG2021-25802) of Prof Stefano Luminari, and by Fondazione GRADE ONLUS Reggio Emilia to Dr. Alessia Ruffini.

Open Access Funding provided by Universita degli Studi di Modena e Reggio Emilia within the CRUI-CARE Agreement.

Disclosure statement

No potential conflict of interest was declared by the author(s).