Linear NIMROD simulations of n = 0 modes for straight tokamak configuration and comparison with analytic results Open Access

Vertical, axisymmetric modes with toroidal mode number n = 0 are of concern for the safe operation of tokamak plasma discharges, as these modes, if growing uncontrolled, may lead to vertical displacement events (VDE) and plasma current disruptions. Therefore, a considerable amount of work has been dedicated to the study of this problem^1–8 (references quoted above are a non-exhaustive sample). In general, active feedback control systems have been developed, which under normal circumstances can prevent the occurrence of VDE.^9–14 

Recently, there has been renewed interest in the theory of n = 0 modes, prompted by the observation of saturated, high frequency n = 0 oscillations in the JET tokamak.^15,16 Analytic theory work, which was possible on the basis of simplified plasma equilibria and geometry, indicated that magnetic X-points, which are intrinsic in modern tokamak plasmas with divertor configurations, can have a significant impact on the dynamics of these modes.^17,18 In addition, a new type of n = 0 mode, dubbed¹⁹ vertical displacement oscillatory mode (VDOM), which could be driven unstable by resonant interaction with fast ion orbits,²⁰ was proposed. Clearly, the treatment of more realistic equilibria requires numerical work. It must be pointed out that treating X-point effects on n = 0 perturbations is challenging from the numerical point of view, since the toroidal field line going through a magnetic X-point is resonant, in the sense that $B eq · k = 0$ for n = 0 at magnetic X-points, with $k$ being the perturbation wavevector and $B e q$ being the equilibrium magnetic field.

The purpose of this article is to provide a numerical benchmark between analytic theory and linear, initial value simulations performed by means of the extended-MHD code NIMROD.²¹ Some preliminary tests of the numerical model for n = 0 vertical modes with NIMROD have already been presented in Ref. 22, with this manuscript providing a more complete and more detailed benchmark. The benchmark will include scenarios where n = 0 oscillations are either unstable according to ideal-MHD, or executing stable oscillations with a frequency close to the poloidal Alfvén frequency. In these numerical simulations, we will mostly stay away from the difficult numerical treatment of X-point effects, in the sense that the hot plasma density is assumed to vanish before the magnetic separatrix is actually reached. Nevertheless, evidence of X-point effects can be seen in the simulations presented in this article, as the X-points actually lie in a low-density, high-resistivity halo plasma, rather than in vacuum. Thus, X-point effects manifest themselves through the appearance of X-point currents carried by the halo plasma.

In order to better clarify the scope of the present article, we point out that four different scenarios come to mind for the occurrence of n = 0 modes in tokamak plasmas. These modes can appear as follows: (1) Unstable vertical displacement perturbations, potentially leading to VDE;^1–8 (2) Stable, vertical displacement oscillatory modes (VDOMs),^14,19 with a frequency of the order of the poloidal Alfvén frequency (a few hundred kHz in a tokamak such as JET), which could be driven unstable by fast ions;²⁰ (3) Global Alfvén eigenmodes²³ (GAE), also with a frequency of the order of the poloidal Alfvén frequency, typically somewhat higher than the VDOM frequency, depending on details of the q profile; and (4) Geodesic acoustic modes (GAM),^24,25 with a relatively low frequency of the order of the sound wave frequency (e.g., a few tens of kHz in a tokamak plasma such as JET). GAE and GAM may also be driven unstable by fast ions.^15,26 In this article, we will be dealing only with the first two types of n = 0 perturbations, i.e., unstable vertical displacements and VDOMs. Also, fast ions physics, which could lead to VDOM destabilization, is not included in this article.

This article is organized as follows. In Sec. II, we briefly review analytic theory results on unstable vertical displacements and VDOMs that are up for comparison with the numerical simulations presented in Sec. IV. In Sec. III, we recall the main features of the NIMROD code. Differences between the assumptions and approximations used in the analytic model and the physics underlining the NIMROD code are highlighted. The straight tokamak equilibrium adopted in analytic theory and in NIMROD simulations is presented. In all linear NIMROD simulations discussed in this article, the confining wall is taken as ideally conducting. Results are presented in Sec. IV, subdivided into five subsections. Subsection IV A considers the no-wall limit, where vertical displacements are ideal-MHD unstable, with a growth rate that depends on plasma ellipticity. In Subsection IV B, we present wall position scans at fixed values of plasma ellipticity. As the wall is moved closer to the plasma, the ideal-MHD vertical instability is suppressed (passive wall stabilization) and the n = 0 perturbation transitions into stable VDOMs. Subsection IV C presents the results obtained by scanning plasma elongation (plasma triangularity is not considered in this article) at fixed wall position, showing again the transition from ideally unstable vertical displacements to stable VDOMs. Subsection IV D focuses on oscillatory modes. Subsection IV E examines the question of X-point currents carried by the halo plasma. Conclusions are presented in Sec. V.

In this section, we briefly summarize the analytic results on vertical modes and VDOMs obtained in Refs. 17–19, 22, and 27, limiting our considerations to those results that are relevant for the numerical benchmarking presented in this article. The magnetic field is represented as $B = e ϕ × ∇ ψ + B ϕ e ϕ$⁠, where $e ϕ$ is the unit vector along the ignorable toroidal direction, and $B ϕ$ is constant. The plasma flow is represented as $v = e ϕ × ∇ φ + v ϕ e ϕ$⁠. Neglecting toroidal curvature effects, an idealized, low-beta, straight tokamak equilibrium is assumed, with an equilibrium, toroidal current density of the type $J e q = J 0 H ( ψ b − ψ e q )$⁠, where H(x) is the Heaviside (unit step) function, $ψ e q ( x , y )$ is the equilibrium magnetic flux function, and $ψ e q = ψ b = const$ is an elliptical flux surface with minor semiaxis a and major semiaxis b representing the plasma boundary. Beyond the plasma boundary, the density also drops to zero, and, therefore, in analytic work, the magnetic X-points associated with this elongated equilibrium are assumed to lie in vacuum. The plasma and vacuum regions are bounded by a resistive wall, which is also modeled as an ellipse, with major and minor semiaxes b_w and a_w, respectively, confocal with the elliptical plasma boundary, i.e., $b w 2 − a w 2 = b 2 − a 2$⁠. It is also assumed that equilibrium flows are absent.

The normal mode analysis for n = 0 modes is carried out on the basis of the linearized, reduced ideal-MHD model.²⁸ In the low-beta limit that is standard for tokamak plasmas, velocity flows and magnetic field perturbations along the toroidal direction decouple from the perpendicular dynamics. Thus, the model equations further reduce to a two-field model for the variables ψ and $φ$⁠. It is convenient to introduce elliptical coordinates $( μ , θ )$⁠, which are related to Cartesian coordinates through the transformation relations $x = A sinh ( μ ) cos ( θ )$ and $y = A cosh ( μ ) sin ( θ )$⁠, with $A = b 2 − a 2$⁠. The elliptical plasma boundary corresponds to $μ = μ b$⁠, such that $a = A sinh μ b$ and $b = A cosh μ b$⁠.

Plots of $ω 2 τ A 2$ as function of κ and $b / b w$ are shown in Sec. IV here below.

Considering finite wall resistivity in the relevant limits D > 1 and $ω 0 τ η w ≫ 1$⁠, one finds that the two oscillatory roots, $ω ≈ ± ω 0$⁠, are weakly damped by wall resistivity, while the zero-frequency root is purely growing on the resistive wall time scale. In tokamak experiments, active feedback stabilization by means of external currents is used to suppress the n = 0 resistive wall mode. However, in this article, we focus our attention on the two other roots of the n = 0 dispersion relation in the limit of an ideal wall.

So far, the analytic results described above assume that the X-points of the equilibrium flux function lie in vacuum. If the hot plasma density extends to the magnetic separatrix, axisymmetric X-point currents can be driven. It was shown in Ref. 17 that these currents can suppress the growth of vertical displacements on the ideal-MHD Alfvén time scale, even if the ideal wall is moved to infinity (D = 0). From the numerical point of view, this requires a very accurate treatment of the plasma response in the vicinity of the magnetic X-points, which will be the subject of a future publication. Nevertheless, in Sec. IV of this article, where the X-points are modeled to lie in a low-density, high-resistivity halo plasma, we will show evidence of the existence of X-point currents carried by the halo plasma.

For simplicity of NIMROD implementation, the analytic straight tokamak is represented by a Cartesian rectangular domain in the poloidal plane and a single toroidal n = 0 mode in the axisymmetric z-direction. This differs from the confocal elliptical wall in the analytic model. Instead, the four sides of the rectangular simulation domain are constrained by imposing the confocal condition, $b w 2 − a w 2 = b 2 − a 2$⁠, so that the half-height, b_w, and the half width, a_w, of the rectangular boundary match the major and minor axes of a confocal elliptical boundary. The rectangular simulation domain circumscribes the analytic confocal ellipse and as such, for equal values of the parameter $b / b w$⁠, the simulation domain is larger than the analytic one. As we will discuss in Sec. IV, this difference in domain geometry is a likely source of discrepancy between the numeric and analytic results, particularly as the wall is brought closer and the difference in area more pronounced.

NIMROD does not support a vacuum model and instead utilizes a “halo” plasma model to represent the analytic vacuum. The “halo” plasma is a cold region with low density and high resistivity, η_halo; a few orders of magnitude difference from hot core plasma. This “halo” plasma fills the entire domain outside of the hot plasma (see Fig. 1), both between the hot plasma and separatrix and throughout the open field line region. A measure of the halo resistivity is represented by the dimensionless Lundquist number, $S halo = τ R / τ A ∝ η halo − 1$⁠, where the resistive diffusion time $τ R = a b / η halo$ with a being the semi-minor axis and b the semi-major axis, and the Alfvén time $τ A = μ 0 ρ m 0 a b / B 0$ with B₀ and $ρ m 0$ being the magnetic field and mass density at the magnetic axis.

In this section, we present the benchmark comparison between the analytic theory results summarized in Sec. II and NIMROD numerical simulations. In numerical work, the wall is assumed to be ideal. In NIMROD, ideal conducting walls are implemented by enforcing vanishing normal components of the perturbed magnetic field at the wall.

First, we will compare unstable vertical displacements and scans of wall position and ellipticity. Then, we will consider vertical displacement oscillatory modes (VDOMs). In analytic work, where vacuum fills the space between the plasma boundary and the wall, and in the absence of dissipation and fast ions, VDOMs are undamped, purely oscillatory modes. In NIMROD, where finite dissipation is usually required for numerical stability, and a low-density, high-resistivity halo plasma is used instead of vacuum, some degree of damping is always found. This damping is shown to saturate to small values when the halo plasma approach vacuum parameters.

In most simulations presented, the grid resolution is (mx,my) = (360,360) with polynomial degree = 3 and the time step $Δ t = 0.5 τ A$⁠. The density diffusion $τ D / τ A ∼ O ( 10 2 )$ (⁠ $τ D = a b / D$⁠); viscosity $τ ν / τ A ∼ O ( 10 4 )$ (⁠ $τ ν = a b / ν$⁠); resistivity $S ∼ O ( 10 6 )$ inside the elliptical plasma; and $∇ · B$ diffusion $τ κ divb / τ A ∼ O ( 10 − 1 )$ (⁠ $τ κ divb = a b / κ divb$⁠) are found sufficient for numerical convergence. However, in some cases, as indicated below, a more refined mesh was necessary to obtain convergence.

This first scan examines the growth rate of the plasma unbound by any conducting wall, i.e., the no-wall limit, as the plasma ellipticity is varied. NIMROD with conducting wall boundaries does not simulate the analytic no-wall limit directly. Instead, for each value of ellipticity, NIMROD estimates the no-wall growth rate by performing successive simulations moving the domain boundary further away until a saturated value for the growth rate is obtained. As shown in Fig. 2(b), the no-wall growth rate for $e 0 = 0.3$ is reached in simulation for values of $b w ≥ 6 b$⁠. For $e 0 = [ 0.2 − 0.6 ]$⁠, the wall parameter $b / b w$ necessary for convergence are in the range (0.13, 0.19).

Figure 2(a) shows good agreement between NIMROD linear simulations and analytic results. However, agreement is better at larger values of ellipticity. For $e 0 = 0.1$⁠, the growth rate is not quite converged. The highest value $b / b w = 0.13$ was used. Convergence is constrained by the need to include and resolve the X-points in the ever increasing domain. As ellipticity is decreased and approaches 0, the X-points retreat away from the plasma and approach infinity, requiring ever larger simulation domains, the majority of which is then composed of halo plasma and open field lines. Values of $e 0 < 0.1$ become computationally prohibitive; the e₀ = 0 circular limit is numerically out of reach.

The growth rate in every run is calculated only after the perturbed fields have reached a pure exponentially growing phase after an initial transient. A typical example for the perturbed x-component of the magnetic field, $B ̃ x$⁠, from an unstable run is plotted in Fig. 3.

Figure 4(a) shows perturbed pressure contours and Fig. 4(b) 2D arrow plots indicating the direction of the perturbed velocity projected onto a poloidal cross section. The pressure contours show a clear m = 1 signature caused by the upward movement of the plasma column. The arrow plot indicates the direction of flow in the poloidal plane. Within the hot plasma, all arrows point upward, consistent with a nearly rigid vertical shift of the entire plasma and in very good agreement with the mode structure obtained analytically. Outside the hot plasma, in the halo region, the arrows show the direction of the return flow. In analytic theory, where vacuum is assumed instead of halo plasma, the return flow is concentrated on the plasma boundary, giving rise to localized vorticity sheets.

We next examine the impact of moving the wall closer to the plasma. The results of a confocal wall scan is presented for two plasma elongations, $κ = 2.0$ and 1.4, corresponding to ellipticity values $e 0 = 0.6$ and $e 0 = 0.32$⁠, respectively. For these two scans, we preserve the confocal constraints (⁠ $b w 2 − a w 2 = b 2 − a 2$⁠) applied to the rectangular domain boundary. Figure 5, the confocal wall scan, plots the square of the normalized frequency vs wall parameter, $b / b w$⁠, and shows agreement between NIMROD and analytic theory. Positive values (shown in circles) indicate oscillating modes. Negative values (shown in triangles) indicate growing modes. The no-wall limit corresponds to values of $b / b w → 0$⁠, while the wall touches the plasma when $b / b w = 1$⁠. The zero crossing of the analytic theory curve (green line) indicates transition between purely oscillating and purely growing mode, which occurs at b_w = b_X, with b_X the vertical distance of the X-points from the magnetic axis. In other words, marginal stability (⁠ $ω 2 = 0$⁠) occurs when the confocal elliptical wall intercepts the X-points. Starting from small values of $b / b w$⁠, the growth rate, $γ = − i ω$⁠, decreases as the wall is brought closer to the plasma (increasing $b / b w$⁠), until the wall approaches the X-point. With the proximity of the wall to the X-point, the mode becomes marginal, and then transitions to a purely oscillating mode for values of $b / b w > b / b X$⁠, as the wall is brought within the X-points, i.e., as the X-points lie outside the simulation domain.

In the simulations, the transition from the unstable growing mode to stable oscillatory mode does not occur exactly when the rectangular wall intercepts the X-points, as in analytic theory. We also note that numerical points lie slightly below the analytic curve, and that this discrepancy is more pronounced at larger values of $b / b w$⁠, where the wall is closer to the plasma boundary. We believe that this is due to the use of the rectangular simulation domain, which has a larger area than the analytic confocal elliptical boundary for equal values of the parameter $b / b w$⁠. Therefore, the rectangular wall is on average a bit further away from the plasma boundary than the equivalent confocal elliptical wall. Hence, the stabilizing effect exerted by the rectangular wall is a little weaker and the growth rate found numerically is indeed a bit larger (i.e., smaller negative $ω 2$⁠) than the analytic result.

To test this idea, in Fig. 6 we present another wall scan for $κ = 2.0$⁠, where we adopted a self-similar wall, as shown in Fig. 1 as green dashed-dotted rectangle, i.e., $b w / b = a w / a$⁠, instead of the confocal constraint. Again, we see agreement between simulation and analytic theory. As expected, the numerical results for the self-similar wall lie slightly above the analytic curve, whereas in the prior two wall scans with confocal walls, the simulation results lie below the analytic curve. Again, the discrepancy is larger at larger values of $b / b w$⁠. We believe that this is due to the fact that, in this case, the self-similar wall is on average closer to the plasma than the confocal wall for equal values of $b / b w$⁠, and so the stabilizing influence of the wall is larger in this case.

A scan of the plasma elongation is performed for values of κ between κ = 1 and κ = 2 (e₀ between 0 and 0.6), with wall parameter fixed at $b / b w = 0.25$⁠. Figure 7 shows the squared normalized frequency plotted against κ; oscillatory cases are denoted by circles and unstable cases by triangles. The zero crossing of the analytic curve occurs for the value of $κ = κ X$ corresponding to the elliptical wall intersecting the X-points. This value can be obtained from Eq. (11), where specifically we set $b / b X = b / b w = 0.25$⁠, yielding $κ X ≈ 1.12$⁠. Agreement between analytic theory and numerical results is excellent, as the wall at $b w = 4 b$ is relatively far from the plasma, and so the effect due to the difference between elliptical and rectangular wall is negligible in this case.

As the conducting wall is brought closer to the plasma, so that the X-points lie beyond the wall and outside the integration domain, the mode transitions from unstable to a stable oscillatory mode. These stable oscillations were dubbed VDOMs, which stands for vertical displacement oscillatory modes.¹⁹ 

Here, we examine one such stable oscillatory case. Figure 8(a) shows the time history of perturbed fields: velocity (⁠ $v ̃ y$⁠), magnetic field (⁠ $B ̃ x$⁠), density (⁠ $n ̃$⁠), temperature (⁠ $T ̃$⁠), and Fig. 8(b) midplane profiles of vertical momentum per unit mass, $n 0 v ̃ y$⁠, at several times during a single period, for κ = 1.4 and $b / b w$ = 0.55. The single oscillation period is denoted by the yellow bar in the time history. The normalized oscillation frequency is $ω τ A = 0.31$⁠.

The time slices of the midplane profile show that the n = 0 oscillation is mostly confined to within the hot plasma, with a perturbed mass flow that well approximates the rigid-shift vertical oscillation found in analytic work. Note that the plasma extends to $x / b = a / b ≈ 0.7$⁠.

Recall that NIMROD does not support a true vacuum, but instead employs a low-density, high-resistivity halo plasma in the region surrounding the hot plasma. As a result, consistent with the flow patterns shown in Fig. 4(b), a return flow is carried by the halo plasma, which guarantees near-incompressibility, $∇ · v = 0$⁠. However, owing to the low halo density, the return mass flow carried by the halo plasma is very small compared with the mass flow in the hot plasma.

Simulations show that the oscillation is damped due to the finite halo density and resistivity. Figure 9(a) is a scan of halo plasma density n_h and Fig. 9(b) Lundquist number S_halo vs damping. For the halo density scan, the Lundquist number is fixed at S_halo = 0.4. A converged damping rate is approached for $n h < O ( 10 17 ) m − 3$⁠. For the Lundquist number scan, we fix the density at $n h = 1 × 10 17 m − 3$⁠. A converged damping rate is approached for $S halo < O ( 10 − 1 )$⁠.

These simulations required a more refined mesh with (mx,my) = (480,480) and polynomial degree = 3.

Here, we discuss an ideal-MHD unstable n = 0 vertical displacement mode with κ = 2.0 and $b / b w$ = 0.4 and focus on the currents that form in response to the growing vertical motion of the plasma. Figure 10(a) shows perturbed surface currents at the plasma boundary, as predicted by analytic work. The surface current density is very narrow, a consequence of having chosen a step-function equilibrium current density with a sharp drop at the plasma boundary.

Of greater interest to us are the X-point currents, as we believe that this is an interesting result with potential consequences on X-point topology and plasma edge MHD activity. Figure 10(a) shows a 2D color plot and Fig. 10(b) 1D profiles along the lower vertical midplane at x = 0 of the perturbed axial current density, $J ̃ z$ normalized by the value of perturbed $v ̃ y$ at the origin. The color plot shows both surface and X-point currents, having opposite polarity, as expected. The 1D profiles of $J ̃ z$ are drawn for different values of the halo Lundquist number. These profiles become sharper and narrower as halo resistivity is reduced, which is in qualitative agreement with analytic predictions.¹⁷ More specifically, as shown in Ref. 17, the perturbed X-point current density becomes a delta function peaked at the resonant X-point in the limit of vanishing resistivity.

There is, however, an important difference between our NIMROD simulations and the analytic work. In these NIMROD simulations, the X-points lie in a halo plasma and the X-point current is sensitive to halo resistivity and halo density. By contrast, in analytic work,¹⁷ the X-points lie in a vacuum.

Figure 11 shows the behavior of the growth rate as a function of halo resistivity. In a relatively narrow range, we find that the growth rate decreases as the Lundquist number increases. This suggests a stabilizing effect from the X-point currents on vertical displacements, but this effect is found to be modest in these simulations. We may expect a more pronounced effect in situations where the hot plasma extends to the X-points and the magnetic separatrix.¹⁷ The growth rate saturates both at the lower end and at the upper end of the resistivity scan. The saturated growth rate at low values of the Lundquist number (high resistivity), where the halo resistivity approximates vacuum, is in good agreement with analytic prediction.^18,19 The saturated growth rate at the higher end has no analytic comparison, as halo effects, which were not considered in analytic work, are important in determining that saturated values.

The main result of this article is that we have successfully reproduced recent analytic results on n = 0 modes using the extended-MHD code NIMROD. The agreement between analytic theory and numerical results is very satisfactory, providing a successful benchmark and a useful starting point for future numerical investigations of n = 0 modes using more realistic tokamak geometry and plasma equilibria. Differences between analytic and numerical results have been interpreted as due mostly to the different shape of the wall used in analytic work vs NIMROD simulation, and the fact that NIMROD adopts a low-density, high-resistivity halo plasma in lieu of vacuum assumed in analytic work.

When ideal-MHD vertical displacements are passively stabilized due to the proximity of the ideal wall, vertical displacement oscillatory modes, dubbed VDOMs in Ref. 19, are found. In the present simulations, VDOMs are weakly damped, due to low values of plasma and halo resistivity. Viscosity and other dissipative terms also affect the damping rate, but at the low values of these parameters used in the simulations (see Sec. IV, before Sec. IV A), their effects on damping is negligibly small. The VDOM frequency is Alfvénic, and for realistic values of the wall parameter D defined in Eq. (8), it falls slightly below the poloidal Alfvén frequency, which should make these modes immune from continuum damping when an actual toroidal geometry with realistic equilibrium current density profiles is considered. VDOMs are of great interest as they can be driven unstable by their resonant interaction with fast ions. Indeed, we have proposed²⁰ that recent observations^15,16 of n = 0 modes on JET tokamak experiments with a significant energetic ion population can be interpreted as fast-ion-driven vertical modes. Work is in progress studying recent JET discharges to test this hypothesis. Ongoing linear toroidal simulations, using the NIMROD code, aim to extend our analysis to more realistic plasma geometries and profiles, starting from reconstructed equilibria of the relevant JET discharges. We point out that these JET discharges are single-null configurations, and yet, especially as far as VDOMs are concerned, preliminary NIMROD simulations are in good qualitative agreement with what is reported in the present article. These results will be presented in a future publication.

This article also discusses for the first time numerical results related to the occurrence of axisymmetric X-point currents. In this article, X-point currents are supported by the halo plasma. In future work, where we plan to extend the hot plasma to the magnetic separatrix, X-point currents may become more pronounced, leading to a more significant impact on the topology of the X-point region, as the analytic work of Refs. 17 and 18 suggests. We point out that treating X-point effects on n = 0 perturbation is numerically challenging, since the toroidal field line going through a magnetic X-point is resonant (in the sense that $B eq · k = 0$ at magnetic X-points), which is the reason why axisymmetric current sheets localized near X-points are to be expected.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom Research and Training Programme (Grant Agreement No. 101052200–EUROfusion). A part of this work was carried out using the EUROfusion High Performance Computer (Marconi-Fusion) under the Project No. $FUA 37 _ NSAM$⁠. The views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

This study was carried out within the Ministerial Decree No. 1062/2021 and received funding from the FSE REACT-EU - PON Ricerca e Innovazione 2014–2020. This manuscript reflects only the authors' views and opinions, neither the European Union nor the European Commission can be considered responsible for them.

We would also like to gratefully acknowledge Eni S.p.A. for the use of their High Performance Computing facility (HPC4), where NIMROD simulations were carried out, and Claudio Carati of Eni S.p.A. for his interest in our work. The author D.B. expresses his gratitude to the partial support from Eni S.p.A for his current researcher position.

The authors have no conflicts to disclose.

Debabrata Banerjee: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Funding acquisition (supporting); Investigation (lead); Methodology (equal); Software (equal); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Charlson C. Kim: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Tommaso Barberis: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Investigation (equal); Methodology (equal); Software (supporting); Validation (supporting); Visualization (supporting); Writing – review & editing (equal). Francesco Porcelli: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

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