Direct preemptive stabilization of m , n = 2 , 1 neoclassical tearing modes by electron cyclotron current drive in the DIII-D low-torque ITER baseline scenario

We report the first preemptive direct stabilization of m , n = 2 , 1 neoclassical tearing modes (NTMs) in DIII-D low-torque ITER baseline scenario plasmas by electron cyclotron waves. These experiments strongly support that the observed stabilization is achieved by replacing the missing bootstrap current in the island by electron cyclotron current drive (ECCD). These plasmas (with and without ECCD) are stable to 2,1 NTMs when the differential rotation between the q = 1 and q = 2 surfaces ( ∆ f 1 , 2 ) is sustained above a critical value ( ∆ f CRIT ≈ 1 kHz), but those evolving to low differential rotation ( ∆ f 1 , 2 < ∆ f CRIT ) routinely develop a disruptive 2,1 island at the time and frequency of a sawtooth precursor. Preemptive, local and direct stabilization by ECCD is tested by scanning the current drive amplitude and location inside and outside of the q = 2 rational surface shot by shot. Analysis of only low differential rotation time windows, i.e. when stabilization from ∆ f 1 , 2 is absent, shows focused ECCD at q = 2 prevents the onset of 2,1 NTMs. These results are important, as they give the first demonstration of disruptive NTM control by ECCD in the plasma scenario planned for ITER.


Introduction
Tearing modes are resistive magnetohydrodynamic instabilitie tokamak plasmas [1]. Due to the fast parallel transport along * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the field lines, the pressure profile flattens within the magnetic islands [2][3][4][5][6]. This leads to a helical perturbation of the pressure gradient driven bootstrap current (j BS ), which reinforces the island growth [7][8][9]. Neoclassical tearing modes (NTMs) are tearing modes that are dominantly j BS driven. m, n = 2, 1 NTMs degrade the confinement [10] and are the most common single root cause of tokamak disruptions [11] (m and n are the poloidal and toroidal mode numbers), which makes their avoidance and active stabilization critically important tasks for scenario development of present day experiments and for operational regimes of future fusion reactors. Local electron cyclotron current drive (ECCD) may be used to replace j BS , which can therefore stabilize NTMs [12][13][14][15][16][17][18]. In addition, due to the nested flux surface topology and reduced cross-field diffusivity within the island, EC heating alone produces a temperature perturbation in the island [19][20][21][22], which may also modify the local resistivity and thus the ohmic current to partially compensate for the missing bootstrap current. On the other hand, preemptively driven ECCD may prevent 2,1 seed islands from surpassing the threshold for neoclassical growth and thereby avoid the growth of large 2,1 islands and concomitant disruptions. However, the effect of the EC heating on the resistivity alone is not expected to prevent the onset of NTMs as it only becomes effective once the island has grown to a large enough size to produce a peaked temperature perturbation in response to the localized heating [23].
The ECCD stabilization approach has become a well documented experimental method in many present day tokamak experiments demonstrating the direct and local effect of ECCD on tearing modes in different variations of high-performancemode (H-mode) scenarios [19,[24][25][26][27][28][29][30][31][32][33]. In the past two decades, a number of developmental tasks have been undertaken at DIII-D to advance the art of tearing stabilization by ECCD. Control algorithms in the DIII-D plasma control system have demonstrated successful real-time alignment of 3,2 and 2,1 islands separately with the ECCD current deposition location, allowing full suppression of tearing modes [34,35]. Suppression of the 2,1 NTM was demonstrated by ramping the toroidal field to move the ECCD location until it overlapped with the q = 2 surface in a standard H-mode plasma [36]. Preemptive stabilization of the 3,2 island was demonstrated allowing stable long pulse operation in plasmas prone to 3,2 islands [37]. Next, the ECCD was modulated in synchronization with the O-point which allowed tearing control with less net EC power [38]. The development continued by real time mirror steering as an actuator for alignment control [29,39]. In more recent years, locked island stabilization was also reported [30]. Despite the long history and broad development of the ECCD technique in DIII-D, the stabilization of 2,1 islands in the low-torque ITER baseline scenario [40,41] has never been reported. This is an important next step since ECCD is the leading candidate to stabilize NTMs in reactors, including ITER [42,43], which will be equipped with ECCD NTM suppression capability [44,45].
In this paper we report the first preemptive, direct avoidance of rotating m, n = 2, 1 NTMs by ECCD in DIII-D plasmas [46] using the ITER baseline scenario normalized parameters and shape at low-torque (T inj < 1 Nm). With and without the application of ECCD, these plasmas are stable to 2,1 NTMs (i.e. passively stable) when the differential rotation between the q = 1 and q = 2 surfaces (∆f 1,2 ) is sustained above a critical value of ∆f CRIT ≈ 1 kHz (q is the safety factor). Those plasmas evolving to ∆f 1,2 < ∆f CRIT routinely develop a disruptive 2,1 island at the time and frequency of an m, n = 1, 1 sawtooth precursor when the ECCD is not aligned with the q = 2 surface. In agreement, the rotation has been long known to be important for NTM stability [47][48][49][50]. Therefore, plasmas that reach ∆f 1,2 < ∆f CRIT provide a testbed for evaluating the efficacy of ECCD in NTM avoidance in low-torque DIII-D ITER baseline scenario plasmas. In this paper we report tests of preemptive, local and direct stabilization of m, n = 2, 1 NTMs by ECCD by scanning the current drive amplitude and location inside and outside of the q = 2 rational surface shot by shot. Analyzing time-points of ∆f 1,2 < ∆f CRIT , i.e. when the instability could occur, the ECCD is shown to prevent the onset of rotating 2,1 NTMs provided that the electron cyclotron beam overlaps with the q = 2 rational surface. These results are important, as they give the first demonstration of disruptive NTM control by ECCD in the plasma scenario planned for ITER.
The rest of the paper is organized as follows. Section 2 describes the discharge parameters and the ECCD aiming scheme. This is followed by the analysis of 2, 1 tearing onset conditions in section 3. The ECCD stabilization analysis is presented in section 4 and modeling in section 5. Finally, the discussion and summary is given in section 6.

Discharge parameters and aiming
We analyze a group of 17 DIII-D ITER baseline scenario discharges. These plasmas are selected from a dedicated experiment testing direct preemptive stabilization of the 2,1 NTM (referred to as 'first set'). A few additional shots are complementing this set by a follow-up experiment (referred to as 'second set'). Two examples are shown in figure 1, one with and one without ECCD [figure 1(e)]. In the shot with ECCD in figure 1, the ECCD density is about twice the local bootstrap current density and it is driven at q = 2. The relevant parameters include I P = 1. The q = 2 rational surface is at ρ ≈ 0.75-0.8 (ρ is the square root of the normalized toroidal flux surface label). The edge safety factor is q 95 ≈ 3. With the goal to test the efficiency of direct preemptive stabilization of 2,1 NTMs by ECCD, the driven current (j ECCD ) and its location (ρ ECCD ) were varied shot by shot. The current drive location, aiming geometry and deposited power were held constant within each shot. The ECCD amplitude and location is calculated with the GA-Toray ray tracing code [53] using magnetic equilibrium reconstructions constrained by internal Motional Stark Effect data [54], and density and temperature measurements (Thomson scattering [55]) without taking into account the helical structure of the islands.
Examples of deposited heat and current drive profiles are shown in figures 2(b) and (c), respectively. In the first set of experiments, the wall conditions lead to high density, which could have prevented access to q = 2 by the ECCD. In order to enable access one gyrotron was aimed at ρ ≈ 0.6 (q ≈ 1.5) starting prior to H-mode. This kept the density at a safe level throughout the discharge by leveraging the density pump-out effect of this beam and thereby enabling the rest of the beams to be deposited at the desired location. In the example shown in figure 2 the beam at ρ ≈ 0.6 drove about 6.55 kA ECCD and deposited about 0.5 MW heat. The rest of the 5 EC beams were focused at ρ q = 2 (where the safety factor is q = 2), they drove about 12.9 kA ECCD and deposited about 2.5 MW heat. The maximum ECCD is about 10 A cm 2 , which is about twice the local bootstrap current density (j BS , calculated by the AutoOneTwo transport code based on measured equilibrium and kinetic profiles entered into the Sauter bootstrap current model [56]). The current drive position is scanned in the radial range defined by where FWHM ECCD is the full width half maximum of the Gaussian EC beam (sum of the individual beams). This provides measurements with various degree of overlap between the q = 2 surface and the EC beam between the maximum and zero, including cases with partial overlap on both sides of the surface.
In the second set of experiments, the wall conditions led to favorable density allowing access of q = 2 and beyond without using a fixed EC beam at ρ = 0.6. In this case, 4 gyrotrons were available and they were moved together shot-by-shot delivering 2.3 MW power. The driven current was about 9.7 kA when all the EC beams were centered at q = 2, corresponding to about half the local j BS . Finally, to test the effect of the EC heating alone without current drive, one additional discharge is included in the analysis with near perpendicular EC beam trajectory to the q = 2 surface (173 643). In this case 1.1 MW EC heating was used which is about half the lowest power at which stabilization by ECCD is observed in this study.
3. Key parameter selection: the role of differential rotation in 2,1 tearing onset in low-torque ITER baseline scenario plasmas These plasmas commonly develop one or more non-disruptive core islands with n > 1 and m = n + 1 mode numbers which rotationally couple, phase-lock and flatten the rotation profile preceding the onset of the m, n = 2, 1 island as first reported in experiments [57] and explained theoretically [58]. Subsequent experimental work lead to the observation of disruptive NTM seeding by 3-wave coupling in the low-torque ITER baseline scenario in DIII-D [59]. In most of these plasmas the sawtooth instability is nearly periodic and is characterized by nearly uniform amplitude. Examples of steep and flat rotation profiles before and at the time of onset, respectively, are shown in figure 3(a) (charge exchange recombination, CER [60]). The m, n = 2, 1 tearing onset occurs at the time and frequency of a sawtooth precursor when ∆f 1,2 approaches zero. This suggest that ∆f 1,2 has a stabilizing effect on the driven 2,1 seed islands [59], possibly through ion polarization currents [61][62][63][64] and screening [65] the sawtooth magnetic perturbation at the 2,1 rational surface. The 2,1 islands, once destabilized, rapidly grow to 25 G or more and terminate the H-mode confinement [figure 1(h)]. In plasmas with intermittent sawtooth activity (not shown) the n = 1 perturbations can grow to significantly larger amplitude which may seed 2,1 islands at ∆f 1,2 ≈ 1.5 kHz. Figure 3(b) depicts the growth of the 2, 1 islands vs. ∆f 1,2 in all unstable discharges of our database. This includes shots with and without ECCD. The open (filled) circles represent time-slices prior to (after) the instability onset. The points with ∆f 1,2 > ∆f CRIT are all open, i.e. the plasmas are stable at sufficiently high differential rotation. When ∆f 1,2 drops below ∆f CRIT , the plasma is metastable until the next sawtooth occurs. This allows the plasma to stay at ∆f 1,2 < ∆f CRIT for a limited time before the 2,1 NTM is destabilized, which is comparable to the sawtooth period. The average time these plasmas spent at ∆f 1,2 < ∆f CRIT before the 2,1 onset w/o well aligned ECCD is t CRIT = 240 ms and the sawtooth period is τ ST ≈ 176 ms.
To further illustrate the importance of ∆f 1,2 in the tearing stability of these plasmas, three example discharges are shown in figure 4. The first discharge is stable without ECCD (left column). The contour shown in figure 4(a) depicts the evolution of the plasma toroidal rotation profile. These profiles are spline fits (ZipFit) to the carbon rotation measured via CER. The radial location of the q = 1 and q = 2 surfaces are extracted from internally constrained equilibrium reconstructions and are over plotted in figure 4(a) with white lines. The toroidal rotation is then interpolated at the q = 1 and q = 2 surfaces in figure 4(b) which shows that the rotation shear survives between these rational surfaces throughout the entire discharge. The time trace of the differential rotation in figure 4(c) shows that ∆f 1,2 > 1 kHz and figure 4(d) shows that 2,1 tearing modes do not develop despite the presence of regular sawtooth activity. The second example (middle column), also without ECCD, becomes unstable and disrupts due to a 2,1 tearing mode when ∆f 1,2 approaches zero. The third example (right column) depicts a case using j ECCD /j BS ≈ 2 at q = 2. In this case the plasma remained stable to rotating 2,1 tearing modes despite ∆f 1,2 falling below 1 kHz.
As the 2,1 tearing modes can only destabilize with ∆f 1,2 < ∆f CRIT in this database, one must only use the ∆f 1,2 < ∆f CRIT time points in order to test the effect of ECCD on 2,1 NTMs. In other words, points with ∆f 1,2 > ∆f CRIT should be removed, as they are stabilized by the differential rotation, regardless of the use of ECCD. Therefore the points with ∆f 1,2 > ∆f CRIT provide no information about the efficacy of ECCD. On the other hand, if one assumed that all time points are equivalent for 2,1 tearing destabilization (i.e. ignoring the stabilization of ∆f 1,2 ) then one could populate the stability map with stable points at any ECCD amplitude and location (when ∆f 1,2 > ∆f CRIT ), thereby masking possible correlations between ECCD and the NTM stability.

ECCD stabilization map of 2,1 NTMs
The ECCD stabilization map is constructed by plotting each stable and unstable time-slice in the plane of j ECCD /j BS and (ρ ECCD − ρ q = 2 )/FWHM ECCD , see figure 5. Recall that j ECCD /j BS is the ECCD amplitude relative to local bootstrap current and (ρ ECCD − ρ q = 2 )/FWHM ECCD is the location of the center of the Gaussian ECCD beam relative to the q = 2 surface, normalized to the full width half maximum of the EC beam. As mentioned in section 2, j ECCD , ρ ECCD and FWHM ECCD are calculated via the GA-Toray ray tracing code, j BS is calculated via the AutoOneTwo transport code and ρ q = 2 comes from internally constrained EFITs.
In figure 5 the third dimension, ∆f 1,2 /∆f CRIT , is represented by the color of the points. As mentioned in section 3, ∆f 1,2 is calculated from ZipFits using carbon impurity CER data. Open circles represent time-slices prior to the onset of the 2,1 tearing mode with 20 ms time resolution (stable points), while filled circles represent the time-slice of the 2,1 island onset (unstable points).
The stability map in figure 5 shows that the central region bounded by the dashed boundary is populated mostly with stable points which suggests that stabilization might be occurring in that region. One instability landed within the expected stabilized region near the boundary. This could be caused either by effects not captured by the model, or it could be explained by the noise or a systematic uncertainty of the q = 2 surface location by ρ ≈ 0.02 in the equilibrium reconstruction. We label this region as stable region, whose boundary is estimated by modeling in section 5 and is shown by a thick dashed line. However, outside of the stable region the stable points do not isolate from the unstable points. As such, the  and stable plasma w ECCD ((i)-(l)). Evolution of plasma toroidal rotation profile ((a), (e), (i)), toroidal rotation at q = 1 and q = 2 ((b), ( f ), ( j)) and differential rotation ((c), (g), (k)). n = 1 magnetic activity and the input EC power ((d), (h), (l)). data points in the j ECCD /j BS and (ρ ECCD − ρ q = 2 )/FWHM ECCD plane (without accounting for ∆f 1,2 /∆f CRIT ) do not show an obvious correlation of the stability with the driven EC current.
The ∆f 1,2 /∆f CRIT data provides additional insights. On one hand, unstable points occur in the vicinity of low ∆f 1,2 /∆f CRIT points (i.e. instabilities only occur when the rotation profile is flat). On the other hand, stable, low ∆f 1,2 /∆f CRIT points occur (i) either near the unstable points or (ii) they lie in the stable region where the ECCD is well aligned and unstable points are absent (i.e. plasmas with flat rotation profiles cannot survive w/o ECCD). (i) Implies that the stable low ∆f 1,2 /∆f CRIT points in the vicinity of unstable points belong to the same discharges as they approach the instability onset and therefore these points are transients, not operationally stable equilibria. This can be elucidated by excluding all points with ∆f 1,2 > ∆f CRIT from the stability map, as these points are stabilized by the rotation, i.e. they are stable with or without ECCD and therefore they do not carry information about the efficacy of ECCD stabilization. Figure 6 shows the stability map retaining only those points below the critical differential rotation (for clarity, these points correspond to the ones in the unstable shaded region in figure 3). In addition, points that belong to the same discharge are connected with straight lines in figure 6. Note that the unstable discharges spend an average of t CRIT = 240 ms at ∆f 1,2 < ∆f CRIT before the 2,1 onset. The color of the stable points reflects their time relative to t CRIT . The time-slices after t CRIT are shown with solid black. Figure 6 shows that without well aligned ECCD every stable point is connected to the onset of a 2,1 tearing instability; hence they are not operationally stable, as mentioned earlier. In other words, despite being labeled as 'stable', one cannot operate the plasma in those points, as the plasma quickly evolves into a state where the 2,1 island is destabilized. Hence, the region outside of the stable region is labeled as metastable. In contrast, it is also clear that the plasmas can survive ∆f 1,2 /∆f CRIT < 1 without the onset of a 2,1 tearing mode when the ECCD is well aligned with the q = 2 rational surface in the stable region.
For clarity, we note that one discharge (190 335) with well aligned ECCD was terminated by an m, n = 3, 2 island. Here all 4 available gyrotrons were focused at q = 2 (i.e. no j ECCD at q = 3/2). The plasma differential rotation approached zero while multiple sawtooth crashes occurred but 2,1 islands did not grow, provided that the ECCD was well aligned with q = 2. This suggests that stabilization is possible even at j ECCD /j BS < 1. However the plasma is prone to tearing modes at low differential rotation: at ∆f 1,2 /∆f CRIT < 1, the m, n = 2, 1 NTM is the most unstable tearing mode when ECCD is not used but when the ECCD is focused at q = 2, the second most unstable mode grows and causes disruption (the m, n = 3/2 tearing mode in this case). When j ECCD is driven both at q = 2 and q = 3/2 (169 546 and 169 539), the plasma survives the pre-programmed flattop.

Modeling the transition between stable and metastable states
In this section we use the modified Rutherfordequation (MRE) to estimate the boundary of the stable domain in the parameter space of y = j ECCD /j BS and x = (ρ ECCD − ρ q = 2 )/FWHM ECCD as in a former study [37], but now in the low-torque IBS, at low differential rotation. The MRE reads: Here τ r is the resistive time-scale, r s is the minor radius coordinate of the q = 2 rational surface, W is the full island width, t is the time, ∆ ′ 0 is the equilibrium classical stability index (here ∆ ′ 0 ≈ 0), δ∆ ′ r s = −(5π 3/2 /32)(3L q /FWHM ECCD )G(x)j ECCD /j tot is the perturbation of the classical tearing stability index caused by j ECCD [13,17], j BS /j tot = 0.16 is the bootstrap current density relative to the total current density, L q = 25 cm is the equilibrium magnetic field shear length-scale at q = 2, W 2 χ ≈ 3 cm is the transport threshold island width [2,6,21,[66][67][68], W marg = 8.5 cm is the marginal island width at which the neoclassical growth becomes dominant (this equals a few times the ion banana orbit width) and K 1 is the effectiveness parameter [14,23]. These parameters are determined in a previous work [62]. Setting dW/dt = 0 at W = W marg , equation (1) gives the boundary for just removing the metastable parameter space, y(x) = 0.96/(1.41K 1 (x) + 2.54G(x)). The effectiveness of replacing the missing bootstrap current depends on both the relative size of the island (W marg /FWHM ECCD ) and the alignment between the electron cyclotron beam and the q = 2 surface, K 1 (x) = −0.021x 3 − 0.254x 2 − 0.432x + 0.599 (at W marg /FWHM ECCD = 0.78). The effectiveness of changing ∆ ′ also depends on the alignment, G(x) = 0.903x 3 − 0.426x 2 − 1.765x + 0.991. Both factors can be negative, i.e. destabilizing, if the misalignment is large. The resulting curve y(x) is shown in figures 5 and 6 with a thick dashed line. This boundary gives a fairly good, although not perfect, isolation between most stable and unstable points.

Summary, implications and open questions
This paper reports the first preemptive stabilization of rotating m, n = 2, 1 NTMs in the low-torque ITER baseline scenario in DIII-D. Scans of ECCD amplitude and location were conducted and the onset condition of the 2,1 tearing modes was investigated both in shots with and without ECCD.
It was shown that in all unstable discharges of our database the onset of the 2,1 island occurs at low differential rotation, regardless of the use of ECCD. Explanation of the unfavorable property of low differential rotation in 2,1 tearing stability is offered as follows. The 2,1 seed islands are triggered by the sawtooth and hence they rotate at the 1,1 frequency. When ∆f 1,2 > 0, the seed island drifts with respect to the ions at q = 2 giving rise to stabilizing ion polarization currents, but this stabilization is weak when ∆f 1,2 ⪅ ∆f CRIT . As high differential rotation naturally stabilizes the plasmas, the effect of ECCD can only be tested when ∆f 1,2 /∆f CRIT < 1.
Next, the ECCD stabilization map was constructed in terms of the j ECCD /j BS , (ρ ECCD − ρ q = 2 )/FWHM ECCD and ∆f 1,2 /∆f CRIT variables, which shows that (i) there are no 2,1 onsets with well aligned ECCD and (ii) there are no stable equilibria at low differential rotation without well aligned ECCD. The boundary between the stable and metastable parameter space is estimated by the MRE, which reasonably, although not perfectly, agrees with the data.
We conclude that the 2,1 island avoidance at low differential rotation was achieved by the use of electron cyclotron waves in these experiments. The presented data are highly supportive that the observed stabilization is achieved by the electron cyclotron current replacing the missing bootstrap current in the island. However, given the limited size of our database, the use of additional ECCD at q = 3/2 in some discharges and the limited power in the ECH only shot, more data are needed to improve the confidence of the latter.
Further, we have reported signatures of stabilization at j ECCD /j BS < 1, although the amount of time this shot spent at low differential rotation is less than the average time the unstable shots spent. This is a potentially important observation, as it suggests that narrow islands can be suppress by lower electron cyclotron current density as long as the pressure profile has not yet completely flattened at the O-point. This would be favorable, as the impact of continuous steady ECCD operation on the fusion energy gain factor (Q) would be less than previously thought [44]. However, further experiments will be needed to determine the minimum j ECCD /j BS required to maintain stability and thus to precisely address its effect on Q in this scenario. In the discharge that achieved stability with 3.5 MW EC power at j ECCD /j BS ≈ 1.5, the H-mode confinement enhancement factor (H 98 ) drops by 14% (from 1.05 to 0.9). Stabilization at j ECCD /j BS < 1 could occur due to the facts that (i) the ECCD modifies ∆ ′ and/or (ii) that narrow islands are not entirely flat due to the finite cross-field transport in the O-point region. At finite cross-field transport the fast parallel transport along the field lines of internal nested flux surfaces of the island do not completely remove the bootstrap current in the O-point region, as opposed to islands much larger than the transport boundary layer within the island separatrices [2,6,20]. On the other hand, operating the ECCD only when ∆f 1,2 ⩽ ∆f CRIT (still preemptively) would further improve the time average of Q.
These results are consistent with and support the recently developed understanding of the 2,1 tearing stability in the DIII-D ITER baseline scenario [59,69], that the 2,1 onset occurs at the time and frequency of a sawtooth precursor when the rotation profile flattens at various times in the discharge regardless of the current profile relaxation phase. This study also suggests that preserving the plasma differential rotation is key to maintaining passive stability in the ITER baseline scenario in DIII-D. Experiments testing the actuators to prevent the differential rotation collapse will be discussed in a future publication.
Finally, how much ∆f CRIT is required for passive stability and if ITER will be able to maintain it remains an open question. This calls for simultaneous modeling of the sawtooth instability growth, the rotation profile evolution and seed-island stabilizing mechanisms arising from the differential rotation in ITER.

Acknowledgments
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, using the DIII-D National Fusion Facility, a DOE Office of Science user facility, under Awards DE-FC02-04ER54698, DE-AC52-07NA27344 and DE-FG02-86ER53218. The data that support the findings of this study are available from the corresponding author upon reasonable request.

Disclaimer
This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.