Yuxiang Sun(孙宇翔) and Di Hu(胡地)

School of Physics,Beihang University,Beijing 100191,China

Keywords: tokamak,magnetic island,tearing mode,magneto-hydrodynamic(MHD)instability

1.Introduction

Due to the nonlinear growth of global magnetohydrodynamic (MHD) instabilities, tokamak plasmas would eventually experience disruptions which not only interrupt the energy output but also cause irreversible damage to the device.[1]Such damage during the thermal quench(TQ)phase of disruptions is characterized by a large amount of thermal load on the plasma facing components (PFCs),[1]while the most feared disruptive consequence in the current quench(CQ) phase is the localized deposition of a large amount of runaway electrons.[2]To avoid such damages and safely guard the continued operation of future high performance devices,various disruption avoidance systems (DAS) are designed,such as the electron cyclotron current drive(ECCD)or Heating (ECRH),[3–6]to suppress the instability of the dominant magnetic islands by modifying the current and pressure profile within.Facing this,disruption mitigation systems(DMS)such as the shattered pellet injection(SPI)[7–17]or other massive material injection(MMI)[18–25]techniques are used to deplete the majority of thermal energy via radiation loss during the TQ, and raise the density to help with runaway electron suppression during the CQ phase.

There have been a lot of recent numerical explorations[11,26–35]on how the DMS would interplay with global MHD modes.Such interplay could strongly impact the mitigation efficiency and these studies provided valuable insight into the disruption mitigation dynamics.However,current simulations often investigate the scenario of injection into initially axisymmetric plasmas,and all MHD response is provoked by the injection itself.In reality,it is most likely that by the time the DMS is triggered, there is already unstable long wavelength global MHD modes, such as tearing modes,existing in the target plasma.The interplay between such finite amplitude mode and the injected materials is of great interest,since it may have a significant impact on the density transport and the outgoing heat flux during the mitigation process.Specifically, the islands could be further destabilized by the helical profile modification within the O-point,[28,34]prompting premature triggering of the TQ before the fragments penetrate into the plasma core.Such a scenario is detrimental to the disruption mitigation efficiency.

While expensive 3D nonlinear MHD codes could be arranged to explore the aforementioned stability impact from profile modifications caused by either the DAS and the DMS,quasi-linear stability analysis using helical Grad–Shafranov equation[36]provides a more agile and cost-effective tool.In our simulation environment, which was run on a virtual machine on AMD R9 5900HX CPU, the typical convergence time ranges from 50 ms to 300 ms,depending on the strength of profile modification.With better hardware this could be further improved.Though it may still not be quick enough for real-time feedback control,it should be possible to obtain a surrogate model for quick equilibria solution,[37,38]which could potentially be used to satisfy the demand of the realtime control.This requires dedicated work and is part of our future work.

The justification for using such a perturbed equilibrium approach for stability analysis is that the growth time of resistive modes is usually on the order of millisecond timescale,much longer than the Alfv´enic timescale which is the perpendicular inertia time for magnetized plasma,so that the plasma can still be considered as a quasi-equilibrium.The advantage of this approach is its capability of avoiding the singularity in solutions close to the resonant surface,and accurately describing the impact of profile modification within the island on the stability criterion∆′.[36]

Considering a dominant finite sized island in a cylindrical analogue of the tokamak plasma with a toroidal effect multiplier 1−m2/n2,[39]the 2D temperature and pressure profile modification with helical symmetry are described and their impact on the island stability is obtained through the calculation of∆′up to the resonant surface.Several different combinations of the current density and pressure gradient modification are considered, and it is found that positive helical current perturbation on island O-point is always strongly stabilizing while the negative one does the opposite, consistent with previous analytical[40]and numerical results.[28,34]Furthermore, the pressure gradient modification is found to have a non-monotonic impact on the stability depending on the island size.It should be noted that,in the scope of this paper,we only analyse simple perturbation of current density and pressure gradient without consideration of specific physical mechanism.In reality,the current and pressure profile modification could be coupled together.For example, the bootstrap current modification means certain current profile modification accompanying any given pressure profile modification.Selfconsistent modelling of such profile modifications is needed for more comprehensive exploration of the island stability,and this is an immediate future development we will pursue.

The paper is arranged as follows.In Section 2,we introduce our quasi-linear perturbed equilibrium equations for the mode structure calculation considering the current and pressure gradient profile modifications.This will form the governing equation of our numerical solution.Then in Section 3,we describe our numerical method,set up our simulation and demonstrate the stability impact of the aforementioned profile modifications.Finally, conclusions and discussion are given in Section 4.

2.Quasi-linear perturbed equilibrium and profile modification for finite island

The fundamental approach of the perturbed equilibrium method is already detailed in Ref.[36],but we explicitly write out the essential parts here for the completeness of the paper.The helical Grad–Shafranov equation for the cylindrical plasma is[41]

where thek=n/(mR) and ¯F ≡Bz+krBθ,mandnare the wave numbers of cylindrical plasma, andBzorBθmeans the magnetic induction intensity inzorθdirection, respectively.Here the toroidal effect multiplier (1−m2/n2) is added to mimic the toroidal coupling effect.[39]The helical fluxΨ∗is the magnetic flux through a helical membrane, and we letΨ∗=(r)+ψwithψbeing the helical flux perturbation.We take the current density to be≡µ0j,µ0=4π×10−7H·m−1.¯p=µ0p(¯Ψ∗)is the plasma pressure in the helical Grad–Shafranov equation.The right hand side can reflect the helical current term asjh= ¯F¯F′(Ψ∗)+(1+k2r2)¯p′(Ψ∗).

With the long wavelength limitk2r2≪1 and use the 2D version of the current and pressure profiles,we can rewrite the Eq.(1)as[36]

We takeP(Ψ∗)≡dp/dΨ∗as pressure gradient,andT(Ψ∗)≡.Here the 2D proflie provides a more accurat e marginal stability boundary compared with results obtained from simple 1D perturbed profile, with the 1D profile predicting instability while the 2D island is actually marginal stable.[36]Combining Eqs.(2) and (3), the helical magnetic fluxΨ∗can be calculated for given boundary conditions once we know the current and pressure gradient profile as a function ofΨ∗.We can see from these equations that the current density perturbation will directly affect the derivative of the mode structure and thus change the instability criterion∆′.But the pressure gradient perturbation has to enter the right hand side together with the (r2−) term.In a symmetric island, the contribution of the pressure gradient should thus be symmetric on both sides of the resonant surface, so that such a symmetrical island takes no effect from this term in the calculation of∆′.And only with the in–out asymmetry of the finite island could this term have an impact on the stability.

Separating the zeroth and the principle harmonics, we could write down the perturbed equilibrium equations[36]

Here,〈f〉0represents the zeroth harmonic component of functionf,while〈f〉1represents the principle harmonicm,n.We are only considering a single helicitym/nhere,since our equilibrium approach requires a given symmetric direction which is determined by the helicity.Although the multi-helicity problem is important in the triggering of the TQ, the pre-TQ phase often sees a single helicity dominant mode which nonlinearly drives the others helicity modes before the TQ onset.[42]If we wish to retain the nonlinear coupling by these subdominant multi-helicity modes in the future,we should retain the principle harmonic of the nonlinear term involving the others helicity modes in the right hand side of Eq.(5), this is a potential future development of this work.The above equations could be re-written as

Here we defined

With a suitable description of the 2D helical modification of the helical current densityj∗and the pressure gradientP, all terms on the right hand side (RHS) of Eqs.(6) and (7) are known, we can then obtain the perturbed helical fluxψand the corresponding∆′from given boundary condition.

The aforementioned profile modifications can be expressed as follows.Since both the helical current density and the pressure gradient are still flux functions,we could first obtain the description of the helical flux near the resonant surface,then extract the corresponding profiles’zeroth and principle harmonics,which will be used in the RHS terms of Eqs.(6)and(7),respectively.The flux function near the resonant point can be represented by the following equation:[36]

wherexis the radial distance to the resonant surface, and the subscript s means the value on the resonant surface.The island semi-widthwand asymmetrysare

We usexSlandxSrdenote the left and right boundaries of the separatrix

Those equations show thatwis the island semi-width,whilesinduces a periodic asymmetry for the island along the poloidal direction.It is obvious that the asymmetry results in the outward shift of X-point and inward shift of O-point of the island,causing a corresponding asymmetry in the (r2−)Pterms discussed above.

With the island geometry known, we introduce the following flux-coordinate with the dimension of length squared:

So that the current density profile can be represented as[36]

wherewl≡w+s=−xSl(0),wr≡w −s=xSr(0), anddl≡wl+s,dr≡wr−s.Like-wise the pressure gradient could be described as the following:

Here we name the different zones of the MHD by the following way.The area near the resonant surfacexS0is a magnetic island as the magnetic lines are reconnected to form a shape of island.Then the zones at−dl

In Eq.(16), the parameters,,,are set to match the value and derivative of unperturbedjprofile at the edge of the buffer zone.Besides, thanks to the fact that conductivity in the plasma is quite high, current perturbations can be limited close to the island region.So we also choose to keep the total current in quasi-linear zone the same as the unperturbed results.And the parameters,,,in Eq.(17)are also iterated to satisfy the matching condition.Here, we consider the cooling or heating timescale to be fast compared with the global transport time,hence the pressure profiles outside of the buffer zone would not be affected significantly.A typical zeroth harmonic result of such profile modification is shown in Fig.1,represented by the green solid lines.

Fig.1.The current density and pressure of 2/1 mode with setting island semi-width to w=0.08 m and initial parameters to f1 =0.5, f2 =0.35.(a)Different j′′i can change the zeroth harmony of current density greatly.The green lines are the profile of the current density which is calculated by the ideal theory.The red dash lines and blue lines are calculated by the quasi-linear theory,while the former set to zero and the latter one set to 50.(b)As the former,shows the effect of different at zeroth harmony, when the general pressure is set as 0.239 MPa.Vertical lines indicate the different zone of the island.

In Ref.[36], the second-order gradientsandwere set to zero to represent the profile flattening within the island.As discussed in Section 1,many physical processes could result in a bump or hole on the O-point for the corresponding profiles.For example,the radiative cooling from impurity SPI into the island could result in both current holes and pressure holes on the O-point,[16,40]while the electron cyclotron wave injection(like ECRH)could result in bumps within the island for both profiles.[3–6]More bizarre cases could exist, for example during hydrogen isotope injections, where the dilution cooling resulted in current expulsion from the island interior but didn’t deplete the pressure due to the low radiation power.Indeed, there might even exist a minor pressure peak within the island due to the heat transfers into the low temperature,but high density,island region.

Some examples of proflie modifciations with selected values ofandare shown in Fig.1 for the zeroth harmonic which is calculated with the setting as semi-widthw=0.08 m and initial parametersf1=0.5,f2=0.35.In Fig.1(a),the unperturbed current density profile is shown in purple,a flattened profile is shown in green and a current hole case is shown in blue.Likewise, the unperturbed, flattened and pressure hole cases are represented by the same colours in Fig.1(b).The red vertical dashed lines represent the location of the resonant surface, the blue dashed lines represent the maximum island width, and the green dashed lines are the 1D buffer region away from the island which is impacted by the 2D island profile modification.

We wish to remind the readers that these profiles are merely some examples of possible profile modification caused by material injection or external drives.The exact shapes of these profile modifications are subjected to detailed analysis of heating,cooling,current drive and transport processes,which are beyond the scope of this paper.In the rest of the paper,we will perform scans on the strength of such profile modifications without delving too much into the feasibility of such modifications,as our focus is on stability analysis.

With the above description,we can finally solve the perturbed equilibrium equation from both the magnetic axis and the plasma boundary,up to both sides of the resonant surface and obtain the stability criterion

We can solve the mode structure up to the resonant surface in this approach since we avoid the singularity in the mode structure solution in the perturbed equilibrium method.[36]We wish to emphasize here that we are calculating∆′on both sides of the reconnecting resonant surface,instead of both sides of the island boundary as conventional wisdom held.[40,43]Such a∆′can be used as an indicator of the island stability.But to obtain the growth rate of the island,one may have to further carry out the flux-average of Ohm’s law carefully, which is beyond the scope of this paper.

As a further note, in toroidal geometry, one should also subtract the Pfirsh–Schluter current contribution[44]to obtain the accurate island stability.Since we are considering∆′on both side of the resonant surface, the flux-average of the Pfirsh–Schluter current on the separatrix must be used instead of the form given in Ref.[44].In this study, however, due to the large major radius we chose, this additional current contribution is small.[45]Hence we simply neglect this additional contribution here, and we are satisfied with merely obtaining∆′as the indication of island stability.We will take this term into calculation in future researches to analyse the growth rate for higher inverse aspect-ration tokamak configuration.

3.Stability impact from the current density and pressure profile modifications

3.1.Numerical setup

We consider a simple unperturbed equilibrium profile with the following form:

whereais the minor radius.The parametersf1andf2are given manually to produce an unperturbed profile,upon which the finite island profile modification would occur.Besides those initial equilibrium profiles,there are some variables like the island semi-widthwand perturbation parametersand, which we will vary in our numerical solution to find out their impact on the stability.In our numerical solution,we also choose some constant parameter values as shown in Table 1.

We are considering cases where there are dominant single helicity modes in the study.As an example,we investigate the 2/1 mode and the 3/1 mode each with a set of current profiles respectively.Thef1andf2will be given for each case.The safety factor for each current profile can be obtained by

The current density, pressure andqprofiles for each case are shown in Fig.2 for the two cases respectively.The initial parameters of 2/1 mode in Fig.2(a)aref1=0.5,f2=0.35,edge safety factorq(a)=3 and the initial parameters of 3/1 mode in Fig.2(b) aref1=−0.30,f2=−0.20, edge safety factorq(a)=3.95.

Table 1.The initialized parameters for the quasi-linear simulation in this study.

Since Eqs.(16) and (17) are 2D in nature, we have to extract their zeroth and principle harmonics to be used in the RHS of Eqs.(6) and (7).We begin with the current density profile.For regions outside of the island but within the buffer region as is marked by the green vertical dashed lines in Fig.1,we have

so the current density outside the island becomes

Fig.2.Use Eq.(19) to calculate the initial profile used for our simulation, and use Eq.(21) to calculate the safety factor in island.Panels (a)and(b)are structures of m=2,n=1,initial parameters in Eq.(19)are f1=0.5, f2=0.35,edge safety factor q(a)=3.Panels(c)and(d)are structures of m=3,n=1,initial parameters in Eq.(19)are f1=−0.30, f2=−0.20,edge safety factor q(a)=3.95.

Here we have defined

whereE(x)is the complete elliptic integral of the second kind.As those equations are functions ofx,we can calculate the zeroth harmonic ofjprofile for any given radius.The perturbed current density in this buffer region is then calculated by taking the derivative of the zeroth harmonic current density just like in linear theory.

The profile modifications within the island region as is marked by the blue dashed lines in Fig.1 are more complicated due to its 2D nature.We have to take into account both the outside island region and the inside island region expressions in Eq.(16).We definet ≡|x−s|/dl,then Eq.(16)could be re-written as

It is obvious that the Eqs.(28)and(29)are functions of both radiusxand poloidal angleθas mentioned above,to use them in Eqs.(6)and(7),one have to extract the zeroth and principle Fourier harmonics,respectively.

The pressure gradient profile can be modified like-wise.As we have discussed above,we consider the pressure increase or decrease within the island to have a faster timescale compared with that of the global transport, so that the pressure bump or hole is located within the island and its immediate vicinity.Once the pressure gradient profileP(χ) is extracted from Eq.(17),the perturbed flux could be solved.The pressure profile itself can be obtained by simply taking the integration from the boundary conditionp(a)=0

where〈P〉∗1 marks the conjugation of the principle harmonic of the pressure gradient.The contribution from the〈Pψ/dxterm is small compared with the average field term for our parameters considered, however.An example of the final pressure profile thus obtained has been shown in Fig.1.

3.2.Impact from current density profile modification within the island

Fig.3.An example of the 2/1 mode structure for an island semi-wide w=0.0801 m,initial parameters f1=0.5, f2=0.35 and without pressure ingredient.(a) The unperturbed current density calculated by a quasilinear way with five vertical lines indicating the different zones of the island.Panel(b)is the helical flux ψ and(c)its radial gradient ∂rψ.

Fig.4.The current density of 2/1 mode with setting island semi-width to w=0.08 m and initial parameters f1=0.5, f2=0.35.The five vertical lines still indicate the five zones of the island,and the green line in each figure is the current density for perturbation is zero and the purple line is for the proflie calculated with(a)set the perturbation to 50;(c)set the perturbation to −50.Panels(b)and(d)are the ∂rψ for the structures(a)and(c).In panels(a)and(c),the red lines mark the resonant surface,the blue lines divide the magnetic island and the green lines are the edge of the buffer zones.In panels(b)and(d)the black lines mark the resonant surface.

Fig.5.Value of ∆′ varies by semi-width w and perturbations Panels (a) and (b) are results of 2/1 mode with edge safety factor q(a)=3 and initial parameters f1 =0.5, f2 =0.35.Panels(c)and(d)are results of 3/1 mode with edge safety factor q(a)=4.5 and initial parameters f1=0.5, f2=0.35.

In this subsection, we first demonstrate the stability impact from only the current profile modification,and show that our perturbed equilibrium result is consistent with previous finite island analyses using constantψassumptions on island boundary.Figure 3 shows an example profile we got by our code which has the island semi-widthw=0.0801 m and no pressure.Figure 3(a)presents the current density which is flattened within the magnetic island.As the figure indicated,the blue dashed lines mark the edge of the magnetic island,and the green dashed lines divide the area where we need to consider the 2D perturbations,buffer zones are the areas which between the blue lines and green lines.While the red dashed line in the middle shows the place of the resonant surface.Those vertical dashed lines have the same meaning in the following figures.What’s more,Figs.3(b)and 3(c)show the calculated magnetic flux and its first derivative.A notable feature is that the perturbed flux is not constant-ψacross the island region.The jump in thecorresponds to the stability criterion∆′.

To demonstrate the effect of the current bump or hole,we show two further examples in Fig.4.As we can see, if we give a negative helical current perturbation to the island,the∆′increases, which means the instability grows.On the other hand, if we give a positive perturbation, the∆′will decrease to negative and indicates the island becomes more stable.This is consistent with previous findings in both the ECCD island stabilizing studies and the cooling island destabilizing studies.[3,4,16,40]A positive helical current perturbation on the island O-point tends to stabilize the island, while a negative helical current perturbation on the O-point tends to do the opposite.

A more comprehensive look at the helical current perturbations impact can be obtained from a parameter scan of island semi-width and current perturbation amplitude, the result of which is shown in Fig.5, assuming no pressure contribution.In Figs.5(a)and 5(c),we show a scan of current perturbation strength with different fixed island widths.First,as the figure indicates, the stronger the perturbation, the bigger the impact to∆′.Consistent with the results in Fig.4,the helical current bump (negative) results in a decrease of∆′and stabilizes the island, while the helical current hole does the opposite.Second,the effect of the current bump or hole depends on the island width,so that for the samethe wider the island,the stronger the effect.This is understandable as with a samethe more current are perturbed in the larger island, thus the stronger impact.Such a trend is common for both modes.In Figs.5(b) and 5(d), these trends are further confrimed by the contour of∆′fordifferent islandsemi-widthwandIn these fgiures,themarginalstableboundary∆′=0 is marked by a solid line,while the others contours are in dashed lines.

Adding fniite pressure results in the stabilization of small islands due to the Glasser effect[46]as shown in Fig.6,where a simple flattened pressure proflie within the island is considered and we usem/n=2/1, island semi-widew=0.0801 m, the initial parametersf1=0.5,f2=0.35 to set the mode.Comparing Figs.6(a)and 5(b),it is apparent that a strong stabilizing effect exists for the small island limit,which is caused by the mode structure modifciation due to fniite pressure gradient near the resonant surface in toroidal geometry.[36]Such pressure contribution vanishes in the large island limit, as shown in Figs.6(a) and 6(b), resulting in a linearly stable but nonlinearly unstable scenario.As the pressure increases the island becomes more and more stable,even with a finite island there comes a threshold inp(0) beyond which there is no unstable island, as shown in Fig.6(c).Note that here the current perturbation from the neo-classical effect is not included,otherwise the island would still be unstable from the decrease of the bootstrap current.

Fig.6.(a) Value of ∆′ varies by semi-width and perturbations , with given general pressure p0 =0.239 MPa, for 2/1 mode and edge safety factor q(a)=3 and initial parameters f1 =0.5, f2 =0.35.(b)A line of=0 shows that with p(0) the ∆′ is not monotonous.(c) ∆′ changes with the general pressure and semi-width of the island at no current density or pressure gradient perturbation.

3.3.Additional contribution from the pressure profile modification

Fig.7.The pressure of the 2/1 mode with setting island semi-width to 0.0801 m and initial parameters f1=0.5, f2=0.35.The five vertical lines still indicate the five zones of the island,and the green line in each figure is the pressure gradient for perturbation is zero and the purple line is for the profile calculated with(a)set the set the perturbation to 50000;(b)set the perturbation to −50000.

A similar scan of the 3/1 mode is plotted in Fig.10.Since the 3/1 mode is stable enough, we do not consider applying negative, which only decreases∆′.With small island width, the 3/1 mode shares the same trend of destabilizing the island with increasingwhile the same non-monotonic feature appears in the large island limit,as can be seen by comparing Figs.10(a)and 10(b)with Figs.9(c)and 9(d).Similar ridge-like feature in∆′exists for the small current hole case(=10), as can be seen by comparing Figs.10(a) and 9(c).With a stronger current hole(=20),the largeis found to stabilize the island in the large island limit for the 3/1 mode,consistent with the 2/1 mode trend.We thus anticipate this to be a common feature that,for a large island,a significant pressure hole would tend to mitigate the destabilizing effect of a current hole within the island,while it contributes to instability in the small island cases.

The above mentioned trend could be more directly shown by considering a scan of differentandwith constant island width.We choose two different island semi-widthsw=0.0235 m andw=0.0801 m, the results are plotted in Fig.11.For the large island case shown in Fig.11(a), thestability impact is once again seen as non-monotonic,as large positiveis seen to mitigate the destabilizing effect of a helical current hole (positive), and large negativeis also found to stabilize the island.For the smaller island case shown in Fig.11(b), theP′′icontribution to island stability is mostly monotonic,with pressure bump consistently stabilizing the island and pressure hole consistently destabilizing the island.These different behaviours in different island width regimes shown in Fig.11 are consistent with what we discussed above regarding the pressure profile modification stability influence.

Fig.8.Value of ∆′ varies by semi-width w and perturbations .Panels(a)and(b)are results of 2/1 mode with edge safety factor q(a)=3 and initial parameters f1=0.5, f2=0.35.Panels(c)and(d)are results of 2/1 mode with edge safety factor q(a)=3 and initial parameters f1=−0.30,f2=−0.20.Panels(e)and(f)are results of 3/1 mode with edge safety factor q(a)=3.95 and initial parameters f1=−0.30, f2=−0.20.

Fig.9.Change of ∆′ figure at given current density perturbation in the meantime for 2/1 mode and edge safety factor q(a)=3, initial parameters f1=0.5, f2=0.35,(a) =−5,(b) =−20,(c) =+10,and(d) =+20.

Fig.10.Change of ∆′ figure at given current density perturbation in the meantime for 3/1 mode and edge safety factor q(a)=4.5, initial parameters f1=−0.30, f2=−0.20,(a) =+10 and(b) =+20.

Fig.11.Change of ∆′ figure at given semi-width in the meantime for 2/1 mode and edge safety factor q(a)=3, initial parameters f1 =0.5,f2=0.35,(a)semi-width w=0.0801 m and(b)semi-width w=0.0253 m.

4.Discussion and conclusion

The stability analysis of finite size island with current and pressure profile modifications near the O-point is carried out in this study using a quasi-linear perturbed equilibrium approach.The helical Grad–Shafranov equation is considered and the mode structure of the principle harmonic of the perturbed helical flux is solved up to the resonant surface to obtain the island stability criterion∆′,where the PS current term from the curvature and pressure is neglected, because we use a small aspect-ratio.The profile modification within the island is introduced by considering the 2D profiles and extracting their zeroth and principle harmonic to feed into the RHS of the helical Grad–Shafranov equation.Among the RHS terms,the helical current perturbation contributes consistently to the island’s stability,with positive helical current perturbation always stabilizing the island while the negative one does the opposite.These findings are in agreement with the previous constantψanalysis.The pressure contribution is more complicated due to the cancellation of its odd parity contribution in the second derivative of the perturbed helical flux.In the small island regime, the pressure modification contribution is monotonic,and the pressure bump is found to stabilize the island while the pressure hole does the opposite,but their overall stability impact is weak due to the aforementioned cancellation.In the large island regime,the pressure impact becomes non-monotonic,especially when we combine the pressure and current modifications.Both large pressure bumps and holes are found to stabilize the island in this regime.

Based on our conclusions above,for example,one would expect that the local heating by electron cyclotron wave within the island would result in an additional stabilization effect apart from the current drive caused by ECCD.Another example is the injection of deuterium pellets into the island O-point,which results in dilution cooling inside the island region.The cooling itself would result in resistive helical current expulsion from the island O-point, thus destabilizing.But the low temperature island region would also cause heat flux to flow into the island O-point,increasing its thermal density and pressure.Based on our result, one would anticipate the pressure bump caused by such an injection to mitigate the destabilization caused by the current expulsion, thus slowing down the island’s growth.These are some examples derived from our results above.

We would like to emphasize again that the current and pressure perturbation considered in this study are arbitrary,so long as they can be expressed as a function of the helical flux.Some immediate future development would be pursued to obtain a more self-consistent description of the island’s stability.For example, the bootstrap current modification could be automatically specified for a given pressure modification, so that the bootstrap contribution can be taken into account selfconsistently.[47]The same is true for the vortices mode contribution within the island which impacts the thermal transport,resulting in pressure profile modification thus bootstrap current modification.[48]Implementing consistent description of the aforementioned profile modifications will be part of our future works.Another future development of the current work is to retain the principal harmonic of the nonlinear terms involving the other helicity modes in the RHS of Eq.(5)to take into the nonlinear coupling with other islands.The non-circular shaping effect of the plasma cross-section can be like-wise added.Furthermore,the convective term could be added into the force balance equation to take the finite flow effect into account in the equilibrium solution.

Overall, the study carried out here provides insight into the island’s response to various disruption avoidance or mitigation scheme such as current drive and local heating caused by cyclotron waves, or local cooling and current expulsion caused by massive material injections.Understandings on the above dynamics would help to achieve more efficient disruption avoidance and mitigation in future machines.

Acknowledgments

The authors thank L.E.Zakharov for fruitful discussion and permission to carry out code development upon his code infrastructure.Project supported by the National MCF Energy Research and Development Program of China (Grant No.2019YFE03010001) and the National Natural Science Foundation of China(Grant No.11905004).