On the power and size of tokamak fusion pilot plants and reactors

It is generally accepted that the route to fusion power involves large devices of ITER scale or larger. However, we show, contrary to expectations, that for steady state tokamaks operating at fixed fractions of the density and beta limits, the fusion gain, Qfus, depends mainly on the absolute level of the fusion power and the energy confinement, and only weakly on the device size. Our investigations are carried out using a system code and also by analytical means. Further, we show that for the two qualitatively different global scalings that have been developed to fit the data contained in the ITER ELMy H-mode database, i.e. the normally used beta-dependent IPB98y2 scaling and the alternative beta-independent scalings, the power needed for high fusion performance differs substantially, typically by factors of three to four. Taken together, these two findings imply that lower power, smaller, and hence potentially lower cost, pilot plants and reactors than currently envisaged may be possible. The main parameters of a candidate low power (∼180 MW), high Qfus (∼5), relatively small (∼1.35 m major radius) device are given.


Introduction
With ITER under construction, the design of candidate pilot plants (Q fus = P fus /P aux ∼ 5-10) and demonstration reactors (Q fus ∼ 20-30) based on the tokamak magnetic configuration is an active area of fusion research [1][2][3]. For a successful design, many parameters, such as the plasma size and shape, current, and toroidal field, have to be taken into account simultaneously. The initial scoping is usually carried out using system codes that capture the main elements of the physics and engineering. Central to the investigations are predictions of the energy confinement time (τ E ) that are typically calculated using the IPB98y2 scaling law, derived from a free-fit to the experimental ELMy H-mode database [4]. One common conclusion of many system code investigations carried out to date is that devices of high fusion gain will have to be large and powerful, typically with major radius in the range 6-9 m, volume 1000 m 3 and fusion power 1 GW [1,5].
We have developed a system code based on a wellestablished physics model to explore possible steady state, high gain fusion devices. We discovered a dependence of some of the key parameters that is contrary to expectation and is likely to be of general applicability in pilot plant and reactor design: that Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. for steady state tokamaks operating at fixed fractions of the density and beta limits, the fusion gain depends mainly on the absolute level of the fusion power and the energy confinement, and only weakly on the device size. Appropriate treatment of the equations underlying the physics model reveals and explains this dependence.
When expressed in dimensionless variables, the IPB98y2 scaling has a negative dependence of τ E on the plasma beta (β −0.9 ). However, it is recognized that this scaling is inconsistent with the results of several dedicated experiments on individual devices in which the beta dependence of τ E has been probed directly. These have shown that τ E has a very weak dependence on beta [6][7][8]. Global scalings in which the beta dependence is constrained to be zero have been found to fit the experimental data almost as well as the free-fit, beta-dependent, scalings (for example, [7]). Using the beta-independent scalings in combination with our finding on the absolute power and energy confinement dependence of the fusion gain potentially leads to devices with high fusion performance at relatively low power and small size with many accompanying advantages. This is not an intuitive result but it is clear from the trials made with the system code. The outline of our code, our analysis and principal findings are presented in this paper. The full details of the code and results of benchmarking comparisons with the results of other codes are available in the accompanying online supplementary material (stacks.iop.org/NF/55/033001/mmedia). The paper has five main sections. In section 2 we outline the system code and in section 3 we present the results of our investigations of the dependence of Q fus on plasma and device parameters. We consider the impact of the beta-independent scalings in section 4, and in section 5 we consider the implications of our results for the design of pilot plants and reactors. A summary is given and conclusions are drawn in section 6.

System code
Our system code is based mainly on the physics model described in the papers by Stambaugh et al [9] and Petty et al [10]. In the description that follows we refer to the equations in those papers as Sn and Pn respectively, where n is the equation number in the papers [9] and [10] respectively, and give ours where we use different equations. The code operates in seven distinct calculation steps followed by several physics and engineering checks. In what follows the units are generally m, T, MW, MA, keV. Densities are in units of 10 20 m −3 .
First, the geometry of the device is defined by the radius of the central mechanical structure, R c , plasma aspect ratio A = R 0 /a, elongation, κ, triangularity δ and plasma-wall gap, g. We calculate the major radius R 0 , minor radius, a, plasma surface area, S p , first wall area, S w , and plasma volume, V p , using: We set the wall load due to the neutron power, n w , and determine the neutron power, P n , fusion power, P fus and alpha power, P α , using: P n = n w S w , P fus = P α + P n , P α = P n /4 = P fus /5 so P fus = 5n w S w /4, P α = n w S w /4.
The required Q fus is set and the needed auxiliary power P aux = P fus /Q fus is calculated. We assume that all the auxiliary power is used to drive the non-bootstrap part of the current, i.e. P aux = P cd . In effect we are assuming operation where both current and power balance apply. Luce has described this as the optimum operating point [11]. This is a common assumption in studies of fusion pilot plants and reactors operating in steady state.
The density and temperature profiles are assumed to be parabolic (S7) and we select the exponents of the density and temperature profiles, S n , S T . The central temperature T 0 is selected. The code calculates the fusion reactivity integral using the approximation log 10 σ v = y + C(log 10 T − x) 2 where C = −1.703 2425, x = 1.727 894 07 and y = −15.056 179.
The central DT fuel ion density, n dt0 , is calculated using P26. We assume 50/50 DT mix, T i = T e = T and average mass number M = 2.5. The helium ash fraction, f He , impurity ion charge, Z imp , and impurity ion fraction, f imp , are prescribed. The code calculates the central electron density using P15 and Z eff using charge balance from P30.
We select the value of the plasma toroidal field, B T0 and dimensionless current drive efficiency [12], ζ cd = 32.7nI cd R/P cd T typically taken as 0.5, which is a value that has been achieved in experiments and computation [13,14]. The code calculates the bootstrap current fraction, f bs , using equation (19) in [15] and the plasma current, I p , is determined self-consistently. We take the internal inductance l i = 0.5 and the ratio of the major radius to the magnetic axis as 0.8. The safety factor is calculated from q eng = 5B T0 a 2 κ/R 0 I p and the toroidal plasma beta, β T , is calculated using S12. The power loss due to bremsstrahlung radiation (P brem ) is calculated using S63. Finally, the wall reflectivity, R w , is prescribed for the calculation of the power loss due to cyclotron radiation (P cycl ) using equation (19) in [16]; typically we set R w = 0.6.
Through these steps all the main parameters of the plasma are determined. However, the plasma may exceed known limits and so may not be viable. Therefore a number of physics checks are included: for example, the line average electron density, n, and the Greenwald density limit, n GW = I p /π a 2 are calculated and the ratio is kept <1; the normalized beta, β N = β T /(I p /aB T0 ), is determined and to ensure stability we limit β N < 9/A = β N (max), which is a conservative limit based on other work [17]; similarly we limit κ to 0.9κ max where κ max = 2.4 + 65 exp(−A/0.376) [3]; and we check that q eng > 2. Although at this stage this is primarily a physics investigation, some engineering checks are also included to give an early indication of the feasibility of the devices examined: for example, an indicator of the power load in the divertor is the ratio P div /R 0 , where P div = P transp − P imp , P transp = P α + P cd − P brem − P cyl and P imp is the power radiated in the plasma edge due to impurities, which we have taken as 0.3P transp for the results presented in this paper. The radial build allows for some shielding on the high field side, which attenuates the neutron flux and protects the central core components. The effectiveness of the shielding is calculated using a parametric fit to the results of dedicated MCNP calculations on candidate shield materials [18]. Further engineering details are planned and will be incorporated in later versions of the code. The code has been benchmarked against the results of several independent codes using published data and good agreement has been obtained. Full details of the calculation steps, the incorporated physics and engineering checks, and the results of the benchmarking comparisons are in the online supplementary material (stacks.iop.org/NF/55/033001/mmedia).

Dependence of Q fus on plasma and device parameters
From the point of view of optimizing the design of fusion pilot plants and power reactors it is beneficial to know which plasma and device parameters have the most influence on the fusion performance and especially on the fusion gain, Q fus . The system code is ideally set up to carry out such an investigation.
Using the code we have explored the sensitivity of Q fus to a wide range of plasma and device parameters including, A, R 0 , B T0 , T 0 and P fus . In general, high density leads to high fusion performance and so it is common in pilot plant and reactor studies to assume operation at a fixed, high fraction of the Greenwald density, typically 0.8. High magnetic fields are technically demanding and expensive and so it is also common to assume operation at a fixed, high fraction of the β N limit-typically 0.9 of the maximum allowed for stabilityto maximize the use of the available field. We made the same assumptions in our investigation and in all cases assume steady state conditions. Contrary to expectations, we found that, under the conditions investigated, the dependence of Q fus on device and plasma parameters could be reduced mainly to a dependence on two parameters-the absolute value of P fus and the confinement enhancement factor, H. Notably Q fus is only weakly dependent on the device size. An example of the results is shown in figure 1 where we use the IPB98y2 scaling of energy confinement time. In the upper graph we plot the derived value of Q fus versus R 0 for fixed A = 3.2 and H = 1.5 and for three different levels of P fus (250, 500 and 1000 MW); in the middle graph we plot Q fus versus R 0 for A = 3.2, P fus = 500 MW and for three different values of H (1.1, 1.5 and 1.7). At each point in the scans we have adjusted the values of the plasma temperature and toroidal field so that the operating point corresponds to the density of 0.8n GW and β N of 0.9β N (max). The values of the remaining key plasma and device parameters, such as the plasma current, safety factor and bootstrap fraction, are derived. The values of the main plasma and device parameters for the P fus = 500 MW, H = 1.1, 1.5 and 1.7 cases, are shown in figure 2. Over the wide range of R 0 the volume changes by almost three orders of magnitude demonstrating the very small dependence of Q fus on device size under these operating conditions. Of course, at small volumes the wall load, divertor load, characterized by P div /R 0 , and the magnetic field are unrealistically high, but at moderate R 0 , say ∼3 m, these parameters approach possible achievable values.
There are, of course, many inter-parameter dependences in the equations within the code but the main dependences found with the scans can be readily demonstrated by an analysis of the simplified underlying equations that determine the fusion performance [19]. In the following the symbols have their usual meaning in tokamak physics.
We make the cylindrical approximation so that V ∝Ra 2 ∝ R 3 /A 2 and begin by noting that the fusion power P fus ∝ n 2 T 2 R 3 /A 2 the loss power P L ∝ nT R 3 /A 2 (τ E ) stored energy , the Greenwald density n ∝ I p A 2 /R 2 , β ∝ nT /B 2 ∝ β N I p A/RB, and safety factor q ∝ BR/A 2 I p . From the beta and safety factor expressions, B ∝ nT R/β N I p A ∝ nT Aq/β N B Hence We assume operation at the nominal optimum H factor, i.e. when P aux = P cd . In this situation, P 1/2 L H ∝ n 1/2 T R/A 3/2 I p . We take P L = P cd + P α = P fus (1/Q fus + 1/5) = P fus (Q fus + 5)/ 5Q fus so (Q fus + 5)/5Q fus ∝ P L /P fus ∝ (nT 2 Immediately we see similar dependences to those found with the system code, that is Q fus is strongly dependent on H and P fus and only weakly dependent on device size. It is notable that Q fus does not depend on A. Since q depends on B and β N depends on B −1 , there is no direct dependence on B either. Of course, if there is a strong dependence on B in (τ E ) scaling then the situation will be different but for the IPB98y2 scaling at least that is not the case. A comparison between the predictions of the equation with the predictions of the system code for one particular set of conditions (P fus = 500 MW, A = 3.2, H (IPB98y2) = 1.5), is shown in figure 3. In order to make this comparison it is necessary to normalize the results of the equation to the results of the code, and that is done at R 0 = 2.46 m. The good agreement obtained confirms that equation (1) accurately captures the main dependences. It is clear that the origin of the main dependences is the chosen operation at fixed fractions of the density and beta limits: equation (1) follows from those limits. Although these conditions will not apply precisely for all operating modes both the density and beta will be limited at some level and so these findings have general qualitative applicability for steady state tokamaks.

Impact of beta-independent scaling
Operation at high beta is clearly desirable. The fusion triple product nT τ E ∝ βτ E B 2 T ; the bootstrap fraction is ∝ β and clearly needs to be high to minimize P cd and hence increase Q fus , and, looking further ahead, studies have shown that the cost of electricity is ∼ 1/β [20]. There are, therefore, at least three potential benefits from high beta operation. However, if τ E ∝ β −α as it is in the case of the IPB98y2 scaling, there is potentially a competition in requirements. Without the density and beta limits this could be overcome by using the size and field dependence of other variables in the τ E scaling, mainly the normalized Larmor radius, ρ * ∝ (RB) −1 ; but, as shown above, for steady state tokamaks, increasing the size and field have limited direct impact on Q fus . To achieve high Q fus it is necessary to increase P fus and that requires an increase in the density, which in turn requires an increase in B T and I P . The negative beta dependence in the ITER scaling therefore leads to powerful devices operating at high/moderate field and current, and to be able to handle the high power the devices have to be large. The dependence of τ E on β is therefore a critical matter with a high leverage on the device parameters.
Dedicated experiments on individual tokamaks designed and executed specifically to probe the dependence of τ E on β have, in general, not confirmed the negative dependence: many have shown that τ E is independent of β [7]. In some cases a negative dependence has been seen, and it has been suggested that this may be due to specific plasma shaping and fuelling effects present in those experiments [21]. The inconsistency between the results of the dedicated experiments and the global IPB98y2 scaling has led to many comments in the literature (for example in [22]), and this is recognized as an unresolved problem in the field. Since the beta dependence in the scaling has a substantial effect on the device parameters, the absence of this dependence is significant. The system code is ideally set up to investigate this impact.
In figure 4 we plot P fus versus R 0 for fixed Q fus = 30 for the same conditions as figure 1 but with the addition of the results of the beta-independent scalings. The beta-independent scalings differ but a factor of three to four reduction in P fus relative to the IPB98y2 scaling is typical. We have carried out many such scans and found this reduction to be a general result. For example, in figure 5 we show the results of another scan but in this case with A = 1.8. The dependences on R of the key plasma parameters for this latter case are shown in figure 6. These findings have significant implications for the design of pilot plants and reactors.

Implications for the design of pilot plants and reactors
Although large and powerful devices may ultimately be required for efficient net power production, it is clearly desirable through the development phase if relatively small, high Q fus devices, can be realized since they would make possible multiple, improving, relatively inexpensive development cycles as usually employed in the development of new technological devices. The results shown in figures 5 and 6 offer an attractive possibility. In the region of R 0 ∼ 1.5 m there is a class of devices that, according to the beta-independent scalings, have a high Q fus ∼ 30 but at relatively modest power level. Fortuitously these devices have a low aspect ratio and apparently have an energy confinement that is resiliently independent of beta. A series of investigations coordinated by ITPA Topical Group on Transport and Confinement thoroughly probed the beta dependence of the energy confinement to variations in several key parameters including triangularity, normalized pedestal height, and ion collisionality. The energy    confinement in the low aspect ratio device NSTX was found to be independent of beta in all cases examined [24]. If this transpires to be a generic result for spherical tokamaks it could significantly influence the choice of aspect ratio for future pilot plants and reactors. While these results are encouraging from a physics perspective the feasibility of such devices will depend on engineering considerations and as yet those have not been addressed in detail. The relatively low power will ease the engineering challenges. Relative to some of the large and powerful devices currently being considered as candidate reactors at R 0 6 m, P fus 1 GW, the P div /R 0 and n W values would be similar. P fus is reduced by about a factor of three and R 0 by about a factor of four so P div /R 0 increases by about 30%. The wall area is ∼ R 2 /A and so a reduction of P fus , R 0 and A by factors of three, four and two respectively would lead to n w increasing by about factor of 2.5. The absolute values of P div /R 0 and n w would be in the range of 40 MW/m and 3 MW m −2 respectively, which are challenging but in the region of those likely to have to be dealt with in much larger and more powerful devices. We can see this in figure 6 where at R ∼ 1.5 m P div /R 0 and n w approach possibly achievable values. Another critical area is the central column where high stresses have to be handled and, for steady state tokamaks, shielding has to be provided for the superconducting magnets. Possibly the use of high temperature superconductors, which require less space than low temperature superconductors and less shielding, would make an engineering solution feasible: we are investigating this possibility in our current work.
A step on such a path towards a low aspect ratio, relatively small reactor, could be to construct a pilot plant with Q fus ∼ 5. The parameters of a candidate device given by the system code are shown in table 1. A critical parameter is the H factor. We assume operation at an H factor relative to the beta-independent scalings of 1.5, which corresponds to H ∼ 1.8 relative to IPB98y2. Many of the smaller devices that have provided data Figure 6. The dependence of the key plasma parameters on R 0 corresponding to two scans in figure 5: Q fus = 30 and H (scaling)=1.5 for both the beta-dependent (blue) and beta-independent (red) scalings. For the beta-independent curves we calculate τ E (scaling) using equation (36) in [7]. V p versus R 0 is also shown.
for the ELMy H-mode database have demonstrated H factors relative to the beta-independent scaling at this level or higher (see figure 17 in [7]), albeit at much smaller absolute values of τ E . Further, the equivalent H factors found on the low aspect ratio devices, MAST and NSTX, indicate that such values may be possible at high field [25]. Both the beta-independent scaling and the H factor employed in this device therefore have a basis in experiment.
We include in the table estimates of the main dimensionless parameters that characterize the plasma; toroidal beta, β T , normalized collisionality, ν * , and normalized gyro-radius, ρ * . For the latter two we use the expressions for the global values of these parameters given in [26]. The calculation of ν * requires q 95 and we use equation (51) in [3]. Typical operating conditions for the low aspect ratio devices MAST and NSTX have values ∼ 10%/0.5/0.02 for these parameters respectively. The extrapolation is therefore mainly in terms of collisionality, which would be about an order of magnitude lower than currently achieved. However, values of ν * at this level (∼0.05) have already been obtained in existing devices such as JET [27] and JT60 [28], albeit at higher aspect ratio, with apparently no deleterious effects on confinement. It is noteworthy that the extrapolation is considerably less than for ST component test facilities such as ST-CTF, which have also been proposed [29]. These typically operate at higher temperatures and lower densities and would have a ν * ∼ 0.001.
More generally, the position of the pilot plasma relative to plasmas obtained in existing devices can be seen by evaluating β T , ν * , and ρ * for these devices and plotting versus aspect ratio and adding the points for the pilot. This is done in figure 7 where we have used the data in the international global H-mode confinement database [30] (specifically we use DB4V5 with the selection criteria described in table 4 of [31]. MAST data are not included because the volume averaged temperature is needed for the derivations of global ν * and ρ * and this is not included in the data base.) Representative points for ITER and the ST-CTF are included for comparison. From the figure we see that the pilot plant will occupy a position in dimensionless parameter space presently not occupied but close to existing devices. The figure illustrates the extrapolations in ν * , and ρ * that would be involved and confirms that they would not be large. On the other hand, the extrapolation in terms of the engineering parameters, especially of current and field, would be substantial and experiments at intermediate levels are therefore required. The on-going upgrades of NSTX and MAST should provide valuable results in this context. Tokamak Energy proposes to construct a small, high field (∼3 T) spherical tokamak that would also enable experiments at intermediate levels to be undertaken [32].

Summary and conclusions
In summary, a system code based on an established physics model has been developed and used to identify and quantify the main parameters that drive the performance of steady state tokamak pilot plants and reactors, and the main findings have been supported by analysis of the relevant tokamak physics equations. It has been found that for steady state tokamaks operating at fixed fractions of the density and beta limits, the fusion gain is dependent mainly on the absolute level of the fusion power and the energy confinement, and only weakly on the device size. The impact of the alternative, beta-independent scaling of the energy confinement has been investigated and found to be substantial, possibly leading to a reduction in the power needed for operation at high fusion gain by a factor of typically three to four, thereby potentially making possible relatively small pilot plants and reactors. Under the assumption of the beta-independent scaling, which has a basis in experiments on many devices, a relatively low power (P fus ∼ 185 MW), small size (R 0 ∼ 1.35 m) pilot plant can be envisaged from a physics perspective, and the main parameters of such a device are given. Many relevant details, for example details of power and particle control and current profile control, are not included in the system code and the next step in this investigative work is to examine those details especially for the low aspect ratio devices. These details will affect the numeric values of the parameters given by the system code but probably not the main qualitative findings. An obvious conclusion of these investigations is that it is crucial to resolve the remaining ambiguities on the beta dependence of the scaling of the energy confinement because the beta-dependent and beta-independent scalings lead to quite different fusion performance, and the beta-independent scalings potentially lead to pilot plants and reactors of lower power and smaller size.