A Review of Fusion and Tokamak Research Towards Steady-State Operation: A JAEA Contribution

Providing a historical overview of 50 years of fusion research, a review of the fundamentals and concepts of fusion and research efforts towards the implementation of a steady state tokamak reactor is presented. In 1990, a steady-state tokamak reactor (SSTR) best utilizing the bootstrap current was developed. Since then, significant efforts have been made in major tokamaks, including JT-60U, exploring advanced regimes relevant to the steady state operation of tokamaks. In this paper, the fundamentals of fusion and plasma confinement, and the concepts and research on current drive and MHD stability of advanced tokamaks towards realization of a steady-state tokamak reactor are reviewed, with an emphasis on the contributions of the JAEA. Finally, a view of fusion energy utilization in the 21st century is introduced.


Introduction
Fusion research to provide a scientific basis for fusion energy production has been carried out for 50 years [1]. World fusion research is ultimately directed to the construction and demonstration of a fusion power station in DEMO and to play a meaningful role in the energy supply at the end of this century [2,3]. Fusion research is a combination of: (1) the development of a scientific basis for plasma confinement and enabling technologies for fusion; (2) demonstrations of fusion-relevant plasma cores and technologies for fusion energy and (3) the integrated demonstration of fusion power production. The continued interaction among these three areas of research is fundamental to the success of fusion research.
To achieve such a goal, a fundamental understanding of fusion and plasma confinement is essential [3]. Over the past 50 years of fusion research, significant progress has been made in the area of (1) and we have reached the stage of the demonstration of the scientific and technological feasibility of fusion energy by ITER [4]. Research area (1) can be developed independently if the issues are clearly understood.
The development of (2) requires an integrated understanding of the system concept since the system imposes many constraints and interlinkages among the various physics. This leads to the need for the development of the reactor concept. In 1990, we developed a steady-state tokamak (from the transliteration of the Russian sentence toroidal'naya kamera s aksial'nym magnitnym polem or toroidal chamber with an axial magnetic field) reactor concept called SSTR best utilizing the bootstrap current [5] and its conceptual design [6,7]. The conceptual development of the SSTR (Steady State Tokamak Reactor) was done for this purpose and to implement this philosophy into large tokamak experiments [8].
This paper starts from a review of fusion energy and the principles of magnetic confinement in Section 2, followed by the principles of the steady state tokamak reactor in Section 3, parallel collisional transport physics for the operation of a steady state tokamak in Section 4, MHD stability of the advanced tokamak in Section 5, the role of fusion energy in the 21st century in Section 6 and Summary and Acknowledgements in Section 7.

Fusion Energy
The question of energy development using fusion reactions appeared in the early 20th century when the theory of relativity and quantum mechanics was proposed. German physicist W. Heisenberg , who initiated quantum mechanics in early 20th century, recorded discussions with Danish physicist Niels Bohr  and Lord Rutherford  in Chapter 13-Atomic Power and Elementary Particles (1935)(1936)(1937) of his book "Physics and Beyond". In 1942, the inventor of the fission reactor, Italian physicist Enrico Fermi , while having a lunch at the Columbia Faculty Club, suggested to Edward Teller  the possibility of burning deuterium to develop a large source of energy. Based on his suggestion, Teller made various calculations and found that fusion between deuterium (D) and tritium (T) is a possibility. Tritium and deuterium react at relatively low energy, creating Helium and a neutron. Since helium has a higher binding energy, we can generate a huge amount of energy.
The deuterium nucleus consists of one proton and one neutron. Among combinations of the two nuclei, p-p, n-n and p-n, a bound state is only possible for p-n, which is deuterium. Deuterium was discovered in 1932 by American chemist H.C. Urey 1934 Nobel Prize for Chemistry), who showed that one out of 7,000 hydrogen atoms is deuterium. The bound state of a proton and a neutron to form deuterium can be treated as a two-body problem of nuclear force between neutron and proton originated from the meson exchange forces, as predicted by Japanese Nobel Prize Winner Hideki Yukawa . He explained this strong force as an exchange force of a finite mass meson to get short-range force in marked difference with Coulomb force explained by the mass-less photon exchange force. The wave function inside the box type nuclear potential may be approximated by a sine function AsinK(r-c) (A and K are constants, c is minimum distance). The wave function outside the nuclear potential may be Be −kr (B, k are constants). Decay rate k of the wave function is related with the binding energy as k~E b 0.5 , which implies that radial decay o the wave function is weak and two-nucleon probability stays outside of nuclear potential, which becomes larger (~/ru/ 2~e−2kr from Born relation) if E b is small. In other words, the two nucleons are not strongly bonded and are easy to separate. Conversely, if the binding energy E b is large, the wave function decays rapidly outside the nuclear potential, leaving the separation probability of nucleons as small. The binding energy of deuterium (E b = 2.225MeV) is small and the wave function decays slowly in the r direction. The small binding energy of deuterium resulted in the dissociation of deuterium into a proton and a neutron at the high radiation temperature of the early Universe, by which formation of Helium is delayed to three minutes after the Big Bang and the reason we have so much hydrogen remaining in the Universe. Hydrogen with mass number 3 is called tritium. The word "tritium" comes from the Greek word meaning "third". The nucleus of tritium consists of one proton and two neutrons. Tritium is an unstable isotope and decays into helium-3 by emitting a high-energy electron and a neutrino (T→ 3 He + e − + ). This is called beta decay and has a half-life of 12.26 years. Tritium was first produced in the laboratory by Australian physicist M. Oliphant  in 1934 by colliding deuterium.
Tritium as a fuel for DT fusion is generated by the nuclear reaction of neutrons with lithium. Lithium has two isotopes ( 6 Li and 7 Li) and the abundance of 6 Li and 7 Li in natural Li is 7.4% and 92.4%, respectively. The 6 Li reaction 6 Li + n 3 T + 4 He + 4.8MeV is an exothermic reaction, while the 7 Li reaction 7 Li + n 3 T + 4 He + n'−2.5MeV is an endothermic reaction. The cross section of the 6 Li reaction is of 1/v type and easy to react at low energy. Meanwhile, the 7 Li reaction is called a threshold reaction whose cross section becomes nonzero above a critical energy. The reaction rate for 6 Li is much larger than that for 7 Li. The amount of lithium in seawater is about 233 billion tons and it can be considered an infinite resource if low-cost technologies for the recovery of lithium from seawater were established.
Neutrons were discovered by the British physicist J. Chadwick (1891-1974) in 1932. Neutrons have no net charge, but in the center there is slightly positive charge distribution, which is cancelled by a slightly negative charge distribution in the periphery. The mass is distributed within a radius of about 0.8 fermi. The neutron is slightly heavier than the proton and the difference is about twice the electron mass (1.29MeV). Neutrons alone can't exist stably, and decay to a proton emitting the electron and a neutrino with a half-life of about 12 minutes (n→p + e − + ). The mass of the neutron is greater than the sum of the mass of a proton and an electron, and the mass difference leads to the energy release. This reaction was confirmed in 1948 by the observation of the electrical bending of protons and electrons from the beta decay of a strong neutron beam in a large cylindrical tank. Natural decay of isolated particles always seems to end with a decrease in mass. Since the mass of a neutron is larger than that of proton by about two times the electron mass, a neutron easily decays to a proton while it is difficult for a proton to decay to a neutron. Indeed, the question of whether proton decay occurs has been an important research subject in physics for many years.
Helium is an element with two protons and two neutrons and its mass number is four. The origin of the name of helium is Greek word "helios" meaning the Sun. During a solar eclipse observed in India on August 18 in 1868, British astronomer J.N. Lockyer (1836-1920), who launched the prominent scientific journal Nature, observed the solar corona and discovered a new emission spectrum. He thought that the emission came from an unknown element, which he called as "helium". The binding energy of 4 He is extremely large compared with those of hydrogen and lithium. Such a large binding energy for this particular nucleus is explained by the nuclear "shell model".
The potential shape for a neutron is different from that of a proton since only the nuclear potential operates on neutrons while the Coulomb potential is superimposed for protons. In Helium, two nucleons with + 1/2 and −1/2 spins can sit in the ground state for a neutron and proton independently. 4 He is the first element in the closed-shell state and "2" is the smallest "magic number". This leads to the large binding energy of 4 He. Thus, the 4 He nucleus is particularly stable, and abundant in the Universe created by the Big Bang.
The encounter of deuterium and tritium results in the formation of the "compound nucleus" 5 He ( D T 2 5 H e * 2 4 He n), as shown in Figure 1, by the tunnel effect at a fractional energy of 500 keV Coulomb barrier potential. The compound nucleus has a high reaction probability near 80 keV due to the resonance phenomenon. In this way, Nature gives human beings a chance to use this reaction for practical purposes. Strong nuclear force can operate beyond the Coulomb barrier when the distance r is less than 3 fermi if we use the nuclear radius formula R c = 1.1A 1/3 fermi (1fermi = 10 −15 m). The kinetic energy of the incident nuclei is distributed to the nuclei in the compound nucleus. A neutron and helium that achieve a large energy by chance will escape from the compound nucleus. Fusion cross section considering tunnel effect and resonance is given as: where, 2 = 2 (2ME) −1 is de Broglie wave length/2 of incident nuclei. P(E/E c ) is a Coulomb barrier penetration probability, and the last factor is the Breit-Wigner nuclear resonance cross section. An analytical form of the Coulomb barrier penetration probability was given by Gamov as The point for the derivation of above formula is the long-range nature of the Coulomb force. According to conventional wisdom, we assume the plane wave solution exp[ikz] as a boundary condition of the incident wave at infinity. Due to the long-range nature of Coulomb force, the incident wave is distorted even at infinity. The wave front should be perpendicular to the classical hyperbolic orbit and the incident wave is modified as exp[ik{z + b 0 lnk(r − z)}], where b 0 = e i e j /(4 0 m r u 2 ) = 7.2 × 10 −10 Z i Z j /E r (eV), where (m) is the Landau parameter. In fact, as seen in Figure 2, the measured fusion cross-section is fitted well by the truncated summation of equation (1)   A recent paper [58] revisited this fusion reaction questioning the discussion above and gave a simpler three-parameter fitting based on the approximate analytical collision cross section using the optical potential (= U r + iU i ) for D + T reaction. This simpler formula gives better agreement with the fusion cross-section data at low energy.

Topology
In the natural fusion reactor, the Sun, a dense and hot plasma, is confined with a gravitational field. The characteristic of this force is that it is a central force field and the force acts in the direction of the field lines. For this reason, the confinement bottle has a "Sphere" topology. In the man-made fusion reactor, a high temperature plasma is confined by trapping charged particles with the Lorentz force in a magnetic field to sustain the reaction within a small dimension a 100 millionth that of the Sun. The characteristic of this force is that it acts in the direction perpendicular to the field line. For this reason, the confinement bottle has a "torus" topology. In the MHD framework, plasma equilibrium is governed by a balance between the pressure gradient force and the Lorentz force. This requires "symmetry" for the field line structure. The tokamak configuration has symmetry in a toroidal direction, while "hidden" symmetry is required for a helical configuration.
Considering the magnetic confinement of a hot plasma in a region of three-dimensional space, the boundary must be a closed surface. The sphere is a typical closed surface but it can't be covered with a non-zero vector field and it is always associated with a fixed point (or null point). In the torus, however, the surface can be covered with a non-zero vector field. In our case, we consider the magnetic field as a vector field. Mathematically speaking, all surfaces homeomorphic to a sphere will have a fixed point. This means that a sphere and a torus have different topologies. This surface property of a sphere and a torus shown in Figure 3 does not change, even if they are bent or stretched. A geometrical property, which does not change by continuous deformations, is called "topology". French mathematician Henri Poincare (1854Poincare ( -1912 proved the theorem "A closed surface that can be covered with a vector field without a fixed point is restricted to a torus." This is called the "Poincare theorem" [59]. The meaning of Poincare theorem is important for high temperature plasma confinement. Considering the boundary surface of the magnetic confinement, the plasma will leak from the zero point of a magnetic field vector. To confine the hot plasma, the surface must be covered by non-zero magnetic field. This is why we use toroidal geometry for magnetic confinement.

Integrability and Symmetry in Plasma Equilibrium
The magnetic field is characterized by its incompressibility ( B = 0). This leads to the existence of a vector potential A( × A = B) given by A = − ζ + G ( and ζ are arbitrary poloidal and toroidal angles, G is gauge transformation part). This leads to the Hamilton structure for the magnetic field evolution in the direction of toroidal angle ζ. The magnetic field line trajectory is given by d /dζ = / , d /dζ = − / and can be regarded as the Hamilton equation if we regard as the Hamiltonian, θ as the canonical coordinate, f as the canonical angular momentum, ζ as time. Variational principle of a field line is given by the analogy to the Hamilton action integral, S = Ldt = [p dq/dt − H]dt by substituting the relationship p , q , H , t ζ as so variational principle to give a field line is given by the following formula: In plasma equilibrium, the plasma's expansion force (− P) is balanced with the Lorentz force (J × B). Here, J is the current flowing in the plasma, B is the magnetic field, P is the plasma pressure. This is the basic principle of the magnetic confinement fusion: In this case, the magnetic field lies on the constant pressure surface ( ), called "integrable" and the surface is called the "magnetic surface". This means that B is a linear combination of tangent vectors x/ and x/ ζ on the flux surface. The incompressibility condition of B leads to the existence of the flow function and the coordinate transformation of ( m ) by which B becomes a straight line gives the Clebsch form for the magnetic field B = × , where = m −ζ/q. Then, and becomes 1/2 of toroidal and poloidal fluxes inside the constant P surface and P = P( ). The coordinates ( , m , ζ) are called flux coordinates. The Lagrangian of the magnetic field line L = d /dζ − becomes L = d m /dζ − ( ). Therefore, the Lagrangian has no explicit dependence on "canonical coordinate" m nor "time" ζ and they ( m and ζ) are ignorable coordinates. The action integral of a field line in flux coordinates is given by: The existence of such ignorable coordinates is essential for the existence of plasma equilibrium. If there is apparent geometrical symmetry such as axisymmetry in a tokamak, it is easy to show the existence of a flux function. But, hidden symmetry can be found in general 3D toroidal geometry, called 3D equilibrium. For general 3D equilibrium, the action integral S is given by [60] as follows: The equivalence of variational principle S = 0 to equilibrium equation (4) is apparent from the following equation as shown by [60]: This variational principle is implemented as numerical VMEC code for general toroidal equilibrium [61]. Plasma equilibrium with apparent symmetry in a torus is the axisymmetric magnetic configuration such as the tokamak configuration shown in Figure 4, which is a major object of present fusion research. In the cylindrical coordinate system (R, ζ, Z), the ζ is a cyclic coordinate and / ζ = 0. Hamiltonian (or poloidal flux/2 ) is given by using the ζ component of the vector potential A as = RA ζ (R, Z). The variational principle of plasma equilibrium in axisymmetric geometry is given by S = 0, where S L dRdZ R(B p 2 /2 0 B 2 /2 0 P) dRdZ [62]. The Euler-Lagrange equation for is given by L/ ( L/ R ) / R ( L/ Z ) / Z 0 , which gives the following Grad-Shafranov equation [63,64]: This equation is simply the component of 2 A = − 0 J and can be solved numerically by giving functional form for P and F, which are determined by the transport processes. This plasma current becomes one of free energy sources in addition to plasma pressure to drive various MHD instabilities, which are described in Section 5.

Tokamak Confinement and Inductive Operation
Tokamak has geometrical symmetry in the toroidal direction and this symmetry provides robustness in maintaining a nested flux surface against various parametrical changes leading this configuration to be a front-runner in fusion research. Tokamak achieved equivalent break-even conditions in large tokamaks such as JT-60U [8] and JET or produced significant fusion power (>10 MW) in TFTR [65] and JET [66], while other magnetic confinement fusion experiments remained much lower, as shown in Figure 5. Geometrical symmetry provides good confinement of energetic charged particles as well as thermal plasmas. This is a reason why the tokamak concept was selected for ITER.
However, this configuration requires a net toroidal plasma current [the right hand side of (8) is proportional to toroidal plasma current], which is driven mainly by inductive means. This method is quite effective since the electrical conductivity of 10 keV plasma is 20 times higher than that of Cu at room temperature. But induction of toroidal electric field is limited to a finite pulse length (300-500 seconds in ITER) due to current limits in the transformer (the CS coil current in ITER). This means that a tokamak fusion power station may be pulsed as shown in Figure 6 or it requires huge energy reservoir. The tokamak reactor design based on inductive operation was first made in UWMAK studies by R. Conn [67]. . Tokamak confinement showed significant progresses by decades (1970s, 1980s, 1990s) to reach break-even conditions. The key is to achieve good confinement with collisionless high temperature plasmas. Figure 6. Schematic view of the inductive operation of a tokamak. Change in the primary current induces a toroidal electric field to drive and sustain the plasma current.

Tokamak Continuous Operation
Since present power sources such as oil/coal/natural gas fired plants and fission plants operate continuously, it is highly desirable for a tokamak reactor to be a steady-state power station. To achieve continuous operation in a tokamak, a non-inductive current drive is essential. After the theoretical development of a current drive using lower hybrid waves by N.J. Fisch [68] and subsequent experimental demonstration in the JFT-2 [6], the STARFIRE design [70] was made to realize continuous fusion power production, which in turn showed that recirculation power is larger than expected. Nature blesses human being by providing an interesting physical process, called "bootstrap current" to make tokamaks operate in steady-state [71,72]. The bootstrap current is driven by the collisional relaxation of a distorted velocity distribution function in a rare collision regime (called collisionless plasma), which is a kind of thermo EMF that drives the plasma current in a toroidal direction and 80% of the plasma current is thus driven by the bootstrap current, as shown in Figure 7 [73]. The physics of non-inductive current drive are governed by the collisional transport along the magnetic field. Since the power required for current drive reduces the net electricity from a fusion power station, efficient and also reliable non-inductive current drive methods have to be developed for steady state tokamak fusion plants, as described in Section 4.

The Steady State Tokamak Reactor
Observation of a high bootstrap current fraction, up to 80% in the JT-60 high-β p discharges [73] stimulated the design development of a SSTR, consistent with updated scientific and technological knowledge at that time [5]. The SSTR concept was originally developed in 1989 as a DEMO concept (aiming to demonstrate sustained electric power generation from fusion) with minimum extrapolation from the knowledge in those days [5,6]. A power reactor concept with similar core plasma assumptions but more aggressive technical provisions, ARIES-I [74], was developed independently. Since then various concepts of a tokamak fusion power system have been developed on the basis of advanced tokamak scenarios with high bootstrap current and high confinement performance; namely A-SSTR [75], ARIES-RS [76], CREST [77]. The SSTR design was made with strong involvement of industries [78]. The plant layout is shown in Figure 8 and a bird's eye view of the tokamak core is shown in Figure 9 [79].

Reactor Power Balance
The SSTR concept as well as ARIES-I have provided good scientific and technical guidelines for fusion research and development in the World. Reactor power balance is an important aspect in SSTR [80]. The energy flow diagram of SSTR is shown in Figure 10. Here, P f is the plasma fusion power, P CD is the current drive power, Q is the energy gain of the confined plasma Q = P f /P CD , BD is the energy multiplication factor in the blanket-divertor system, P Ge is the gross electric power, P Net is the net electric power to the grid, rP Ge is the re-circulating power (r is the re-circulating power fraction), CD is the overall system efficiency of the CD system, and aux P th is the power required for auxiliary equipment. From the relations in Figure 10, the plant efficiency η Plant (= P Net /P th ), the ratio of net electric power output to thermal power output, is given by: p l a n t t h a u x Here, the thermal conversion efficiency th is 0.345 for water cooling in a fission light water reactor, while it is 0.49 for an advanced high temperature He cooling system. The second term of rhs is the reduction of plant efficiency due to auxiliary equipment, which ranges from 0.02-0.04. Power for auxiliary systems in SSTR is assumed to be 80 MW, similar to the STARFIRE design but an accurate evaluation is require,d based on ITER experience. The third term of rhs is specific for a tokamak power system and represents the current drive and heating requirements. The system efficiency of N-NBI is estimated to be 0.50, which is much higher than that of the laser driver in ICF, but the electric power of 120 MW is not small and represents one of drawbacks of helical system, that still requires basic research to improve plasma confinement at a reactor-relevant high temperature. The energy multiplication factor in a Blanket-divertor system BD has to be evaluated considering various processes in the blanket and divertor. The nominal BD value for SSTR is 1.21 (= 3710 MW/3060 MW) [78]. The η plant -Q diagram is shown in Figure 11 for a pressurized-water cooling and a high temperature helium gas cooling system. It must be noted that η plant depends weakly on Q at Q = 30 to 50, η plant~0 .3 for pressurized water and η plant~0 .4 for high temperature helium. Efficiency of the current drive by non-inductive means (NBCD or RFCD) is expressed by the current drive efficiency η CD , defined by η CD = I pCD R p < n e >/P CD and has a certain limit η CD~5 × 10 19 A/Wm 2 which is much less than the efficiency of an inductive current drive. In this sense, it is difficult to achieve the required Q level of Q = 30 to 50 by only using non-inductive current drive by external means. This is a fundamental reason why we have to utilize the bootstrap current to achieve the efficient steady state operation of a tokamak reactor. Operation at some Q has two meanings since the CD system has two functions, current drive and heating. The first one is that current drive efficiency must be compatible with this Q value, to support the plasma current Ip with this current drive power. The second one is that the energy confinement time of the plasma must be compatible with this Q value.

High Bootstrap and High Poloidal Beta Operation
The major feature of the efficient steady state operation of a tokamak is the maximum utilization of the bootstrap current. Since the bootstrap current fraction is proportional to the poloidal beta β p* = (4/(μ 0 I p 2 R p ) PdV (or f boot~( a/R) 0.5 β p* ,), the reactor should operate in a high β p* regime. This is a marked difference with high current and high toroidal beta research made before the initiation of the advanced tokamak research. The constraint on the plasma poloidal beta comes from the ideal and resistive MHD stability. The most important constraint is the so-called Troyon scaling described by Prof. F. Troyon in 1984 [81]: where β N is a constant named "normalized beta", β t = < P > /(B t 2 /2μ 0 ) is the volume averaged toroidal beta, I p is the plasma current, a is plasma minor radius. The combination of Troyon scaling with the definition of poloidal beta β p* = (4/(μ 0 I p 2 R p ) PdV gives the following relation: where is vertical plasma elongation. Figure 12 shows the β p* -β t diagram in which the solid curve indicates the β p* -β t relation for = 1.8 and β N = 3.5 for reference. From equation (11), the toroidal beta β t is inversely proportional to β p* at a fixed β N . This scaling is confirmed most beautifully in DIII-D as shown in Figure 12 [82]. The regime corresponding to lower Ip and high q (typically q 95~5 ) regime is called the advanced tokamak regime. . Experimentally confirmed operating regime in ( p , t ) to confirm Troyon scaling in Doublet III (right) [82].
It is found in this figure that the steady-state fusion power concepts SSTR and ARIES-1 (β p* = 2 to 2.1) as well as the current ITER-ss design (β p* ∼1.5) adopt high β p* operation to increase the bootstrap current fraction [8]. While a comprehensive theory of the bootstrap current is given later, a simple expression of the bootstrap current fraction is given by Cordey [83] as follows: With current profile control, it is possible to achieve f bs~7 5% at β p*~2 . This enables current drive of the remaining 25% of the plasma current by a high-energy beam or RF. It is important to notice that before 1990, world tokamak research was focused on increasing plasma current to improve energy confinement, typically represented by the design change from INTOR to ITER-CDA. Also, frontier research was directed towards achieving a high toroidal beta close to 10% with a highly normalized-current I p /aB t in DIII-D [84].

Current Profile Control of High Bootstrap Current Fraction Plasma
Since the bootstrap current has a hollow current profile, control of current profile is important for MHD stability and confinement improvement. Since the current profile is frozen and difficult to change when the plasma temperature becomes high, current profile control before intensive auxiliary heating is important. Access scenarios to advanced confinement regimes such as positive shear, weak shear and negative shear are shown in Figure 13. Figure 13. Access to steady state operation regime with high bootstrap current fraction from OH regime. Since current profile is frozen at high temperature, current profile control in the target OH regime is important, according to Kishimoto [8].
One of important operational diagrams of a tokamak is the Cheng Diagram [85] or (q,li) diagram. Figure 14 shows the (q,li) diagram for JT-60U [86]. The low li boundary is limited by surface kink modes or locked mode while the upper li boundary is limited by tearing mode activities for low regime. Advanced tokamak operation modes such as weak positive shear regime (or high p regime) and reversed shear regimes are below the sawtooth boundary [87,88].  [86]. Here, q eff ~0.8q 95 . High N regime is characterized by relatively higher li and high p regime is characterized by a relatively lower li below the sawtooth boundary [87].

Weak Positive Shear Regime
To eliminate the hollow current profile, an active central current drive with N-NBI was first proposed for the steady state operation with elevated q(0) to stabilize ballooning modes [5], now called weak positive shear regime, as shown in Figure 15. In this regime, improved core confinement is observed in high power neutral-beam heated large tokamaks such as JT-60U [87][88][89] and TFTR [90,91]. Ideal MHD stability of a weak shear regime was studied by Ozeki [78], showing that a SSTR relevant regime can be stable for ideal MHD modes with q(0) > 2 and wall stabilization.

Negative Shear and Current Hole Regimes
After the SSTR proposal, Ozeki [92] found for the first time that a hollow current profile with peaked pressure profile and reduced pressure gradient near q min can be stable for ideal MHD modes and this is called reversed shear operation or negative shear operation. He proposed to use the off-axis NBCD to realize a reversed shear profile as shown in Figure 16. Since the bootstrap current is a hollow current profile and has 1/B p dependence, it is easier to obtain a higher bootstrap current fraction with a hollow current profile. Since then, many works have been done for optimization of the reversed shear scenario in both theory [93,94] and experiments in TFTR [95], DIII-D [96], and JT-60U [97]. To achieve a hollow current profile, a step-up scenario of P NB power during the current ramp-up (dIp/dt~0.5MA/s) is the key aspect to obtain a high electron temperature and to enhance the skin current effect. A key issue for stable evolution of an RS plasma is its stability when q min passes through a low m/n rational surface. When the plasma current is ramped-up, q min also decreases with time and may pass through a low m/n rational surface such as q min = 4 and 3 and tends to disrupt at q min = 2 due to beta collapse [98,99]. Figure 16. Negative shear operation scenario using off-axis current drive by Ozeki [92].
As an extreme situation in a negative shear configuration, equilibrium with zero plasma current in the central regime called "Current Hole" was formed in JT-60U [100] and JET [101], as shown in Figure 17. This current-hole regime can be stably sustained for a few seconds. This regime is interesting from the control viewpoint in that it has low li and is easier to get an elongated plasma and also is easier to get a high bootstrap current fraction. On the other hand, the Current Hole regime is subject to higher ripple loss and sets severe constraints on maximum toroidal field ripple as well as the low no-wall beta limit seen in negative shear regime.

Advanced Tokamak Research
In 1993, the first comprehensive review on the prospects of the steady state tokamak reactor was given, which covered the physics requirements of a high bootstrap current fraction, confinement enhancement factors, non-inductive current drive, MHD stability including disruption probability, power and particle control and the need for new research directions was stressed, in addition to some engineering features of the magnet, neutral beam, coolant and material selection [102]. In 1994, Goldston [103] gave a talk on advanced tokamak physics in TPX design activity to establish the physics basis for the steady state tokamak. Since then, steady state tokamak research is called by the name of "Advanced Tokamak".
In 1997, Taylor [104] gave an EPS invited talk on the physics of advanced tokamaks, which is typically shown by the upper right regime in Figure18 (left). This regime corresponds to a lower Ip and high q (typically q 95~5 ) regime, which was called advanced tokamak regime. He surveyed improved confinement based on ExB shearing of microscopic turbulence and improved MHD stability to achieve higher β N using shaping and pressure profile control, as shown in Figure 18 (right). Ozeki also addressed the physics issues of high bootstrap current tokamaks, including TAE stability in EPS [105]. Recent DIII-D advanced tokamak experiments show 100% noninductive plasma with β t = 3.6%, β N = 3.5, H 89 = 2.4 [106].

Parallel Collisional Transport Physics for Steady State Tokamak Operation
Since plasma current plays an essential role in tokamak confinement, it is quite important to understand the parallel transport physics, especially generalized Ohm's law. Fortunately, most of the parallel transport in a tokamak is governed by the collisional transport and we have developed Matrix Inversion (MI) Method [73] based on the Hirshman-Sigmar neoclassical transport theory [107].

Moment Equation
As Hirshman and Sigmar [107] showed, we obtain the following moment equations for momentum and heat flux by taking v and v 2 v moments of Vlasov-Fokker Planck equation of species a after subtracting convective heat flux from the v 2 v moment: m a n a du a d t e a n a (E u a B Here, n a , u a , q a , P a , a , a , F a1 , F a2 , M a , Q a are density, velocity, conduction heat flux, the average plasma pressure, viscosity tensor (anisotropic component of the pressure) and viscous heat tensor, friction force, heat friction force, momentum source, and heat momentum source, respectively. Velocity distribution function in strong magnetic fields shows anisotropy parallel and perpendicular to the magnetic field, as shown by Chew, Goldberg, and Low [108], and a and a can be expressed as: a ( // a a )(bb Here, b = B/B is the unit vector parallel to B and = a /L is a smallness parameter as a ratio of Larmor radius and the macroscopic length scale L of the plasma, a = v Ta / a is the Larmor radius, v Ta = (2T a /m a ) 1/2 is the thermal velocity, a = e a B/m a is the cyclotron angular frequency. Taking the cross product B with (10) and (14), and neglecting the time derivative for the drift time scale O(( 2 ) -1 ) much longer than the Alfven time scale O(( ) −1 ), and neglecting other O( 2 ) terms smaller than P, , T terms, the major components of particle and heat flows of particle species a perpendicular to B are given as: Here, first term of rhs of (17) is the ExB drift flow, the second term is the diamagnetic drift flow caused by the pressure gradient. Equation (18) Here B F a1 and B F a2 are the frictional forces on species a by the parallel flow on the magnetic surface, and B a and B a are viscous forces parallel to B, which originate from the relaxation of velocity space anisotropy between parallel and perpendicular to B. Substitution of (15) and (16) In 1995, a comprehensive review on the experimental evidence for the bootstrap current was given [109]. In this paper, the origin of this velocity space anisotropy is pictorially explained in case of the electron as in Figure 19. The magnetic moment is conserved in a high temperature plasma when the electron moves along the magnetic field. So, the orbit of the electrons satisfying B max ≥ E/ is trapped in the weak magnetic field regime reflected by the magnetic mirror (trapped particle orbit: or called -Banana Orbit‖ from its shape). Consider the case density is lowering towards the outside (dn/dr < 0). Consider the velocity distribution function in a magnetic surface. There are less electrons for the trapped electrons with v // > 0 since it comes radially from outside, while there are more electrons for the trapped electrons with v // < 0 since it comes radially from inside. Meanwhile, the orbit of un-trapped electrons stays much closer to the magnetic surface and the number of electrons for v // > 0 is roughly equal to that for v // < 0. Then, there appears a discontinuity in the trapped/un-trapped boundary of the velocity distribution function. Small Coulomb collision smoothes this gap and causes the particle diffusion in the velocity space. This collisional diffusion in velocity space acts as a viscous force in the magnetic field direction.

Flux surface Averaged Momentum and Heat Flow Balance
The friction and heat friction forces in (19) and (20) are given by: Here, l ij ab is called the friction coefficient, which has the symmetry l ij ab l ji ba due to self-adjoint property of the Coulomb collision term: Here, the friction coefficient has the following symmetry due to the self-adjointness of the Coulomb collision operator: Since viscous force operates when the particle moves poloidally to feel the variation of the toroidal field, the viscous force is proportional to the poloidal flows: Here, a1 , a2 , a3 are called parallel viscosity coefficients. Collisional transport regimes in tokamak are divided into three regimes: (1) Banana regime where the collision time is longer than the bounce time of trapped particle orbit ( c < b ; c : collision frequency, b : bounce frequency), (2) Pfirsh-Schlüter regime where the collision time is shorter than the transit time of the un-trapped particle ( c > t ; c : collision frequency, t : transit frequency ~v Ta /Rq), (3) plateau region between the two. Expression of the viscosity coefficient is derived for each regime and the velocity partitioned approximate viscosity coefficient valid for all velocity region is derived. The viscosity coefficient is obtained by the integration in the velocity space as follows: Here, first term of the denominator of (33) is the correction term to connect the banana and plateau regimes, while the second term is the Pfirsh-Schlüter correction term. T Also, a * is the collisionality defined by the ratio of collision frequency 1/ aa , the transit frequency Rq v T a a a (38) Here, f t is a trapped particle fraction and is related to un-trapped particle fraction f c through f t + f c = 1. The f c is given by the following equation: Substituting these formulas for viscosity and friction coefficients into (19) and (20) Here, ai , l ij ab , u //a , q //a , M a// , Q //a , V ia (BV 1a = −F( )(d /d + (dP a /d )/e a n a ), BV 2a = −F( )(dT a /d )/e a n a , and is electrostatic potential) are viscosity coefficient, friction coefficient, parallel flow, parallel heat flow, parallel momentum source, parallel heat source, thermodynamic forces, respectively.

Generalized Ohm's Law
The following system of linear equations is obtained by writing down (39) for electrons, ions, and impurities [109]: Here, M, L, U // , V , E*, S // are viscosity matrix, friction matrix, parallel flow vector, thermodynamic force vector, electric field acceleration vector, parallel source vector, and given as follows: Then, plasma current density parallel to the magnetic field is expressed in the following form and is called generalized Ohm's law: Here terms on the right hand side are called bootstrap current, ohmic current, beam and RF driven currents, respectively and are given as follows:

Electrical Conductivity
The inductive operation is quite effective to drive the plasma current since the electrical conductivity of high temperature plasma (e.g., = 10 9 /ohm-m for Te = 10 keV) is 20 times larger than that of copper (e.g., = 5 × 10 7 /ohm-m at room temperature). From equation (45), we obtain the expression of electrical conductivity (fast ion contribution is neglected): // NC e a n a e b n b (M L) 1 ab b e,i,I a e,i,I (48) Here, L represents the collisional friction forces among various species, and M represents effect of trapped particle. If there is no trapped particles, the viscosity matrix M = 0 and conductivity σ is given in this case as follows: which is called -Spitzer conductivity‖. The electrical conductivity is reduced by viscosity M. The trapped particle is trapped in the banana orbit as shown in Figure 19 and does not contribute to the current. It creates a frictional force by the velocity relative to un-trapped electrons, which are drifting by the electric field. The electrical conductivity given by L. Spitzer Jr. [110] Here, Z eff is the effective charge. Hirshman [111] also gave an approximate analytic expression for electrical conductivity as follows: Here, the trapped particle fraction f t (54) is not accurate enough for a non-circular cross-section plasma and in such a case (38) must be used.

Bootstrap Current
The generalized Ohm's law in (42) includes current driven by the thermodynamic forces V 1a and V 2a (BV 1a = −F( )(d /d + (dP a /d )/e a n a ), BV 2a = −F( )(dT a /d )/e a n a ) as follows: Although V 1a includes an electrostatic potential term, this term ( J~ F( Z b n b ba )d /d ) vanishes for the axisymmetric plasma due to charge neutrality. Distortion of the velocity distribution function occurs in collision-less plasma. The electron distribution function is drifting in the direction of v // < 0 while the ion velocity distribution function is drifting in the direction of v // > 0, as seen in Figure 20. This produces the plasma current and is named bootstrap current.  Figure 21 shows a comparison of plasma surface voltages between measurement and 1.5 dimensional transport simulation (called so by the coupling of two dimensional equilibrium and one dimensional current diffusion transport simulation) using the measured plasma parameters. If we do not include bootstrap current in the transport simulation, the simulation does not reproduce the measurement and the existence of the bootstrap current is confirmed [109].  (42) of JT-60 discharge [109]. Up to 80% of the plasma current is carried by the bootstrap current from the simulation. Without including bootstrap current, the measured surface voltage can't be reproduced. [112] 4.6.1. Neutral Beam Current Drive Theory When a fast neutral beam is injected tangentially to the torus, circulating fast ions produce a fast ion current by multiple circulations around the torus. Collision with bulk electrons produces a shielding current by induced drift in the same direction as the fast ions. This shielding is not perfect due to the existence of trapped electrons and impurities. The sum of the fast ion and shielding currents is called beam-driven current J bd , which is a current in response to the external momentum source S //a given in (40) but momentum balance of fast ion has to be included as follows [109]:  (63) where, e b , n b S b are electrical charge, density and momentum source density of beam ion, respectively. Also, F is called stacking factor. Parametric dependences of F on Z eff and e in arbitrary aspect ratio (0 1) are calculated by (18) and shown in comparison with Start-Cordey calculation [113] in Figure 22. Figure 22. Stacking factor F for neutral beam current drive as a function of r/R and Z eff [109] in comparison with Calculation by Start-Cordey [113].

Neutral Beam Current Drive
To evaluate the fast ion current < B J > fast , we have to solve the fast ion Fokker-Planck equation for the velocity distribution function of fast ions f b valid for v Ti ≤ v b ≤ v Te . Fast ion Fokker-Planck equation in non-uniform magnetic field is give by Connor [114]:  (64) is given in J. Cordey [115] as follows: where, a n (v) is the analytical solution for uniform magnetic field by Gaffey [116]: Here, S(v, ) = S 0 (v − v b )k( ) = S 0 (v − v b ) k n c n ( ) and fast ion distribution function above beam energy (v > v b ) comes from energy diffusion in velocity space. The following equation determines the eigen-value l n and eigen-function c n ( ): This equation is solved numerically by the Rayleigh-Ritz method in the ACCOME code [117]. Using above formulas, the flux surface averaged parallel fast ion current multiplied by B, < B J > fast can be calculated as: Bds v The dependence of the current drive efficiency h = J bd /SE b0 on T e and n e are given as follows by taking se~Te 1.5 /n e and E b0 /v b0~vb0~( v b0 /v c )v c , v c~Te 0.5 : where f(x) = f 0 (x)/x 1/2 .

Demonstration of Current Drive with N-NBI
The SSTR design utilizes high energy neutral beam injection for the non-inductive current drive. Therefore it is necessary to demonstrate the effectiveness of heating and current drive by high energy neutral beam for the ITER and the steady state tokamak reactor. Consequently, the first world project to install a 500 keV negative ion based neutral beam injection system (N-NBI) was started in JT-60U to demonstrate its feasibility [118]. Key advantages of N-NBI are: (1) the ability to place the injector far from the reactor for easy maintenance, (2) no accessibility issues across the plasma boundary, (3) robustness against plasma conditions, especially no dependence on the toroidal field, (4) high energy conversion efficiency using a gas neutralizer or plasma neutralized for the reactor, (5) acceptable current drive efficiency = < n e > RI p CD /P CD = 5 × 10 19 A/W/m 2 . Figure 23 shows a bird's eye view of the N-NBI system in JT-60U, which started operation in 1996 [15]. Table 1 shows the major specifications of the N-NBI system in JT-60U [118]. Since then, progress has been made to increase beam voltage by conditioning [119], to improve ion source uniformity [120,121] and to reduce heat load to the electrodes [122,123]. Table 1. Design specification of JT-60U N-NBI system [118]. Current drive efficiency was extensively investigated for a range of beam energies (≤350 keV) and electron temperatures (Te ≤ 14 keV) much wider than in previous experiments [124]. A method to calculate plasma equilibrium with inductive and non-inductive (beam-driven and bootstrap) currents is established in numerical codes such as ACCOME [117]. A method of experimental determination of non-inductively driven current was established by Forest [125]. Ohm's law in a general toroidal geometry without bootstrap current and RFCD is given by:

Item
where the electric field is related to the time variation of poloidal flux at constant F surface as follows: Time evolution of the poloidal flux (r, t) and < B J > can be measured by the MSE diagnostics with equilibrium magnetic fitting code such as EFIT [126]. By calculating the Ohmic current using measured / t, density, temperature and Zeff profiles, non-inductively driven current < B J> NBCD can be -measured‖ as difference < B J > NBm = < B J > − //NC < B E >. This current density profile can be compared with the numerical calculation of < B J > NBc using measured density, temperature and Zeff profiles. A comparison is shown in Figure 24. It is worth noting that calculated and -measured‖ non-inductively driven current profiles agree with each other when multi-step ionization effect is taken into account in the calculation of < B J > NB .
Theoretical NB current drive efficiency increases with electron temperature T e in Equation (70). This dependence is also confirmed experimentally, as shown in Figure 24. Maximum NBCD efficiency h CD = 1.55 × 10 19 A/W/m 2 is achieved at T e (0) = 14 keV with the beam energy of 360 keV. Projected NBCD efficiency for E b = 1 MeV in ITER is (2 to 3) × 10 19 A/W/m 2 for T e (0) = 10 to 20 keV. For the DEMO, higher central temperature T e (0)~30 keV and higher beam energy E b~2 MeV might be necessary to have high NBCD efficiency of 5 × 10 19 A/W/m 2 . These experimental results are encouraging for steady-state operation of ITER, DEMO and beyond. Redistribution of beam driven current and/or reduced driven current have been observed in various tokamaks [127] if the discharge is associated with MHD activities, such as toroidicity-induced Alfven eigenmodes (TAE), sawteeth, fishbones and tearing modes, so it is important to control MHD activity so that NBCD is not deteriorated. Figure 24. Comparison of calculated and measured beam driven current (left) and NBCD current drive efficiency as a function of central electron temperature T e (0) [124,127].

MHD Stability of Advanced Tokamak [128]
Plasma current in a tokamak acts as a free energy source to drive MHD modes such as kink modes. The bootstrap current is linked to the pressure gradient and this linkage produces new MHD modes. Theoretical and experimental progresses to understand and control MHD modes in advanced tokamak are reviewed in this section.

Energy Principle and 2D Newcomb Equation
The stability of plasma equilibrium (8) has been an important subject in fusion research and can be studied by the variation of action integral (6) subject to equilibrium constraint (7). Then variation S becomes second order in displacement and given as follows: where the force operator F known as Hermite (self-adjoint) operator. The W = −(1/2) F( )dv is a change of potential energy and is given by Furth [129] as follows: Here, W SA is the bending energy of the magnetic field and is a source of shear Alfven wave. W MS is the compressing energy of the magnetic field and is a source of magnetosonic waves. W SW is the compressing energy of the plasma and a source of the sound wave. All these terms are positive and stabilizing. Meanwhile, W IC is interchange energy of plasma pressure in the curved magnetic field and can take a positive or a negative value. W kI is the kinking energy by the current and can take a positive or negative value. Here, the curvature vector is given by = b b. If P < 0, the interchange energy is the source of instability. The MHD instability in tokamaks comes from the kink energy term W kI at low beta, called current driven kink modes and interchange/ballooning term W IC is added at high beta, called ballooning modes, while MHD instability in helical system comes mainly from interchange/ballooning term W IC .
In the case of axisymmetric torus, energy integral is minimized under the incompressibility condition = 0 as in the case of cylindrical symmetry. The energy integral W under = 0 can be expressed in the following form by using X = r and V = r ( − /q) in the flux coordinates (r, , ) with r = [2R 0 0 (q/F)d ] 1/2 as formulated by Tokuda [130]: Absence of V/ r term leads to simpler Euler-Lagrange equation for V and its solvability condition on leads to the following two-dimensional Newcomb equation for X: Here, X = (--, X -2 , X -1 , X 0 , X 1 , X 2 ,--) t (t: transposed) where X m is the Fourier component of X and f, g, h are constant matrices. Two solutions of this 2M + 1 equation are singular at a rational surface and others are analytical. MARG2D [130] solves this 2D Newcomb equation for the analysis of peeling modes with high n numbers. Here, peeling mode is an external mode localized near the plasma edge driven by the finite edge current. This mode can be coupled to pressure driven ballooning mode and thought to be a cause of ELM (Edge Localized Modes) in a tokamak.

Tearing and Neoclassical Tearing Modes
In tokamaks, most unstable kink modes with poloidal mode number m = 1,2,3 can be stabilized if the safety factor at 95% flux surface q 95 is above 3 (q 95 > 3) which is the ITER standard operation condition. Operation of q(0) < 1 gives rise to the m = 1 internal kink instability. Stability of these linear ideal MHD modes can be analyzed using 2D Newcomb equation.
Finite resistivity in a high temperature plasma gives rise to resistive instabilities. This finite resistivity enables kink-like deformation with its resonant surface inside the plasma by changing the magnetic field topology at rational surface through magnetic reconnection to create magnetic island and is called -tearing instability‖. This instability with poloidal and toroidal mode numbers m/n = 2/1, 3/1 3/2 is particularly important in a tokamak. The resistive instability with mode number m/n = 1/1 is somewhat different from others due to the breakdown of constant approximation at the rational surface and is called resistive kink mode as cause of sawtooth oscillation through reconnection. When constant approximation at resonant surface is valid, perturbed flux = irB r /m (B r : perturbed radial magnetic field) is approximately governed by following diffusion equation: t 0 2 r 2 (78) This gives rise to the evolution of magnetic island width w as dw/dt = '(w)/2 0 [or more accurately dw/dt = 1.66 ( '(w) − w)/ 0 ), where '(w) = [d /dr(r s + w/2) − d /dr(r s − w/2)]/ (r s )]. As seen from Figure 25, the perturbed current inside the magnetic island is antiparallel to the equilibrium plasma current to form counter clockwise field lines around the island for the case of positive magnetic shear s > 0 (s = (r/q)dq/dr). The formation of magnetic islands reduces the pressure gradient and the reduction of the bootstrap current occurs and accelerates the growth of magnetic islands. This mode is called neoclassical tearing mode (NTM). On the other hand, the perturbed current is parallel to the equilibrium plasma current and reduction of the bootstrap current reduces the magnetic island for s < 0.

Double Tearing Modes in Negative Shear Plasma
Negative shear or reversed shear shown in Figure 26 (ii) is stable to NTM but is subject to double tearing mode (DTM) since there are two rational surfaces. DTM can grow explosively if mode coupling between two rational surfaces is strong. An explosive growth of DTM can occur after quiescent Rutherford regime by nonlinear destabilization of high m/n modes for intermediate separation of two rational surfaces [131]. If the separation between two rational surfaces is large enough, modes are decoupled and island will not grow explosively. Important implication for the plasma control is to pass through low m/n rational q min as quick as possible under reduced pressure gradient and keep wider separation in quasi steady state [132].

Resistive Wall Modes
When the SSTR concept was presented at the 1990 IAEA conference with q(0) > 2 for ballooning mode stability, this calls for Ramos' idea on ideal MHD stability that the free boundary beta limit is inversely proportional to q(0) due to low m, n modes at this conference, later published in 1991 [133]. He proposed a modification to Troyon scaling as follows: This is true for free boundary modes, but the situation with a stabilizing wall is quite different. The MHD stability of SSTR for low m, n kink-ballooning modes are analyzed and found to be stable when q(0) > 2 with wall stabilization [78]. At that time, most people believed that wall stabilization may not work for realistic tokamak circumstances since a wall necessarily has finite resistance and penetration of the magnetic field will nullify wall stabilization. Even if plasma is rotating in the toroidal direction, it was believed that the mode [called resistive wall mode (RWM)] is attached to the wall and the mode will slip with respect to the plasma rotation [134,135].
It was unfortunate that it was not well recognized at that time that it was already shown shear Alfven wave has a continuous spectrum in the inhomogeneous plasma as shown in Figure 27 and RWM mode may damp by the phase mixing process or by the Landau damping through the mode conversion of shear Alfven wave to the kinetic Alfven wave, as already shown by Chen in 1974-1976 [136-139]. We should also note that mode conversion to KAW was already observed experimentally in 1989 [140]. The same mechanism works for TAE modes called continuum damping [141,142], which was also confirmed experimentally in 1995 [143]. In 2007, both DIII-D [144] and JT-60U [145] showed that RWM is stabilized with small toroidal rotation close to the expectation by the continuous damping of Alfven wave, as shown in Figure 28.  [145] and in DIII-D (right) [144].
Effect of wall stabilization in ideal MHD stability was studied by Manickam in 1994 for both positive and negative shear regimes, as shown in Figure 29 [146]. Beta limit is limited by n = 1 mode without wall stabilization and medium n for wall stabilized plasma. Changes in the beta limit between with and without wall stabilization is quite large for negative shear case. The reason for this is quite simple: wall stabilization is easy if the current is closer to the wall but such a surface current is unstable if wall is not effective.
One of practical issues of wall stabilization for tokamak system is the installation of an active feedback control coil, which is now actively discussed in ITER. Another important issue is the relative proximity of the resistive wall in the tokamak reactor [147]. Since the replaceable blanket should be replaced every few years, they must be segmented to reduce electromagnetic force and will not have good shell effect against kink like modes. So, a reactor will have r wall /a = 1.3 to 1.4. Hence, the integrated optimization of blanket and RWM stabilization is important. Figure 29. Comparison of beta limit with and without wall stabilization for positive shear regime q(0)~1 (left) and negative shear regime (right) by Manickam [146]. positive shear case negative shear case

Ballooning and Peeling Modes
The stability of the ballooning mode is determined by the balance between the bending energy W SA and the interchange energy W IC . Ballooning mode has a unique structure with a double periodic condition in poloidal and toroidal directions and is given as summation of quasi-modes in covering space (− , ) and is regarded as infinite radial overlapping of resonant modes.
The finite edge current can drive external modes localized near the plasma edge. This mode is called peeling mode. Peeling mode becomes most unstable when a rational surface is located just outside the plasma surface. This mode can be coupled to pressure driven ballooning mode and is thought to be a cause of ELM (Edge Localized Modes) in tokamaks.
One important attractive feature of a high β p regime is stable access to the second stability regime of ideal ballooning by increasing q(0) [148][149][150][151]. Second stability access near the plasma edge is possible with high edge current density, but the existence of peeling mode near the edge makes it complicated. This mode is destabilized by the edge current, especially edge bootstrap current by the steep pressure gradient.
Local confinement improvement in the edges and inside the plasma results in the edge transport barrier (ETB) (called H-mode) and internal transport barrier (ITB), respectively. These lead to increases of the local pressure gradient destabilizing ballooning mode through an increase in W IC and peeling mode as a combination of W IC and W kI near the confinement barriers. These localized instabilities at ETB and ITB are called edge localized mode (ELM) and barrier localized mode (BLM) [152], respectively.
In the collisionless edge plasmas, a steep pressure gradient near the edge produces a large edge bootstrap current. This edge bootstrap current destabilizes the peeling mode and the radial extent of peeling mode becomes larger with the magnitude of thye edge bootstrap current. This situation is particularly crucial for power control of ELM energy loss in ITER and beyond.

Infernal Modes
If the magnetic shear is finite, radial coupling of various resonant MHD modes (m, m 1, m 2,--) becomes strong since radial separation between modes is small. However, if the magnetic shear is very low s = rdq/dr/q~0, radial mode separation becomes larger and radial mode coupling becomes weaker and standard ballooning mode theory base on dense radial mode coupling breaks down [153] (Hastie and Taylor derived the applicable conditions of ballooning mode theory as n >> ( q'( )) −2 >> 1). Mode growth rate becomes oscillatory as a function of n (or toroidal wave number k z ) treated as a continuous parameter. Under such circumstances, an intermediate integer n mode may become unstable even if lower n modes are stable. This low n internal pressure-drive mode was named -infernal mode‖ by Manickam [154] and is responsible for the p collapses in JT-60U high p regime. Figure 30. From left, assumed pressure and q profiles for ideal MHD stability calculation (A'(li = 1.2), B'(li = 1), (C':li = 0.8)), calculated stable-unstable boundaries of normalized beta g (= N ) for a relatively peaked pressure profile C, and infernal mode stability regime in (n, q 0 ) space [94].
Ideal MHD stability of a high p plasma was investigated by T. Ozeki to explain p collapses as shown in Figure 30 [94]. A peaked pressure profile decreases the beta limit N to low values ( N = 1.5 to 2) in a range of current profiles (li = 0.8 to 1.2), as observed in high performance high p discharges [89]. In the case of high li~1.2 and peaked pressure profiles, beta limits observed as p collapses are identified as internal kink modes, especially by the infernal mode in the low q(0) regime. With increasing q(0) (lowering li), infernal mode is stabilized because of the increase in the magnetic well depth [94].

Alfven Eigenmodes
There are two Alfven waves in a uniform plasma immersed in a static magnetic field. One is the shear Alfven wave (displacement vector is perpendicular to B) and other is the compressional Alfven wave ( parallel to B). The shear Alfven wave at some location can have a resonance independent from neighboring field lines and the resonant condition (k ) can be expressed by = k // v A , where v A = B/( 0 n i m i ) 1/2 is the Alfven velocity. Near the shear Alfven resonance, the shear Alfven wave is mode-converted to a kinetic Alfven wave and strong damping of the wave occurs as long as resonance exists in the plasma, as shown in Figure 31 (left), but for a frequency slightly lower than min (k // v A ) (no shear Alfven resonance) in the plasma but cut-off still exists in the plasma, a shear Alfven wave can exist without strong damping and is called Global Alfven Eigenmodes (GAE).
In toroidal geometry where in-out inhomogeneity B~B 0 /(1 + cos ) exists, the Alfven resonance condition is given by a coupling of m (k //m = (n + m/q)/R) and m + 1 modes as: This gives a forbidden band of for Alfven wave resonance as shown in Figure 31. As we know, sin(m ) + sin((m + 1) ) = sin((m + 0.5) )cos(0.5 ), which implies periodicity after two circulation similar to Mobius band. The Mobius band is not orientable and there is no distinction between front and back surfaces. In this situation, it is difficult to have Alfven resonance, which is an explanation of the gap in the shear Alfven resonance. Although resonance does not occur, shear Alfven mode can propagate for a frequency range in the gap as point (or discrete) spectrum and can be destabilized by the interaction with fast ions as shown by Cheng [155,156]. This TAE was first observed in neutral beam heated plasma by Wong [157]. A review of Alfven Eigenmodes was given by Wong [158]. Alfven eigenmodes due to a particles in ITER are an important control issue. In JT-60U, TAE are observed during ICRF heating [159,160], and new Alfven Eigenmode Non-circular Alfven Eigenmode was observed [161]. TAE modes including chirping modes and reversed shear Alven eigenmodes are also observed during N-NBI heating [162,163].

Role of Fusion Energy in the 21st Century [3]
To prevent global warming, a low carbon society free from fossil energy must be achieved during the 21st century. Although attention has been paid to electricity and hydrogen power as clean energies, low CO 2 emission energy sources must be used to create electricity and hydrogen. CO 2 emission per electric power kwh is an important figure of merit and is called specific CO 2 emission. As shown in Figure 32, fusion is one of energy sources with low specific CO 2 emission next to hydro and light-water reactor. In light water reactors and fast reactors, iodine 131 ( 131 I) produced by the reaction tend to accumulate in certain organs of human body, and the concentration limit in the air (tolerance) is set to very low levels. The hazard potential of 131 I in a 1 GW fission power station is about 1,000 times larger than that of tritium (T) in a 1 GW fusion power station ( Figure 33). Therefore, fusion energy has favorable characteristics in terms of radiation hazards. Realization of a low carbon society by suppressing the use of fossil fuels is necessary to prevent global warming and emissions of greenhouse gases like carbon dioxide and methane. According to the 2100 nuclear vision proposed by the office of Strategy Research, Japan Atomic Energy Agency [164], an option to utilize previous R&D results and atomic energy technologies under development to reduce national CO 2 emission to 10% of present value while maintaining stable energy supply in 2100 as shown in Figure 34. In that scenario, use of fossil energy currently sharing 85% of primary energy is reduced to 30% by the end of this century and the other 70% will be shared by non-fossil energy with a dominant contribution from Atomic Energy. Figure 34. CO 2 reduction scenario (left) and Electricity supply scenario from various sources (right) in Atomic Energy vision 2100 [164].
In this scenario, the use of electricity and hydrogen is promoted as an energy source. Significant reduction in energy consumption is realized by improved energy efficiency through the expanded use of electric vehicles and fuel cell vehicles by way of hybrid vehicles in the transport sector. Coal and oil use in industrial sector is eliminated by substituting coke as a reducing agent in the steel industry and naphtha in the chemical industry with hydrogen. Energy use in civilian areas may become electric except fpr solar heat.
In this way, 60% of the final energy demand in the year 2100 becomes electricity. It is difficult to respond to such a huge power demand with only renewable energy and a large-scale stable supply from Atomic Energy is most promising. In this vision, it is assumed that ~30 fusion plants will be operational in Japan by the end of this century assuming the scientific and technical feasibility of fusion energy is demonstrated in ITER and the construction and operation of DEMO progresses in a timely fashion to start construction of the first commercial fusion plants in the 2050s. It should be noted that it is possible that fusion energy will not be on the energy market by the end of 21st century if fusion R&D does not go well or its economic efficiency and reliability of operation are not good enough.
The tokamak system adopted in ITER shows the best performance in high temperature plasma confinement. The operation of a tokamak stops when an inductive electric field can't be supplied since the confinement field is created by the inductive plasma current. As a power generation device it becomes pulsed. As a method to overcome this drawback and to produce energy continuously, the use of bootstrap current is considered. Since 80% of the plasma current is shed by the bootstrap current in JT-60 ( Figure 21), a steady-state tokamak reactor (SSTR) is designed as a practical way in cooperation with industries [78].
To achieve a continuous operation using the bootstrap current, the majority of the plasma current has to be driven by the bootstrap current and rest by a beam or in some other way, as shown in Figure 35. The heat including the generation inside the blanket is removed and converted to electricity in a steam turbine and a fraction is used as recirculating power for beam generation and other plant needs.
Achievement of continuous operation with Q > 5 given in the ITER technical objectives is important. In this case, about half of the plasma current should be driven by the bootstrap current while the rest should be driven by the beam-driven current. Figure 35. Principle to achieve high net generated power by achieving efficient continuous operation using bootstrap current.

Concluding Remarks
Fusion research has progressed significantly during the 50 years elapsed since its inception in 1958, especially thanks to tokamak research. In this review, a brief introduction of fusion reaction, topology of magnetic confinement and hidden symmetry in force equilibrium is provided leading to the importance of the apparent geometrical symmetry of the tokamak. This symmetry in a tokamak is created by the toroidal plasma current, leading to the question of continuous operation. The bootstrap current driven by the plasma itself provided an opportunity for continuous operation and the concept of a steady state tokamak reactor was formed in 1990, clarifying the necessary plasma regime for the operation of a steady-state tokamak reactor. In this review, the theoretical foundation of the generalized Ohm's law and its experimental confirmation, creation of advanced operational regimes for steady state tokamaks, and progress in our understanding of MHD stability issues of a steady-state tokamak are highlighted. Finally, a view of fusion energy utilization in 21st century is introduced by assuming smooth advances in ITER and DEMO development. This review is a partial review of research towards realization of a steady state tokamak fusion system with emphasis on the principles, physics of plasma currents in tokamaks and a description of important MHDs, and the role of fusion in the 21st century. A more comprehensive physics review will be given elsewhere.