Particle in Cell (PIC) simulation has been used to design a quasi-spherical, net power, p+B11 fueled, fusion reactor incorporating high temperature superconducting (HTS) magnets. By burning a fully thermalized plasma, our proposed MET6 reactor uses the principles of Magneto-Electrostatic Trap design of Yushmanov (1980) to improve the classic Polywell design, resulting in a predicted power-balance Q ≈ 1.5. Because the input power consumed by the reactor barely balances the waste bremsstrahlung radiation, future research will focus on reducing the bremsstrahlung output to reach practical net power levels.
Research Article
A Quasi-Spherical Fusion Reactor Burning Boron-11 Fuel
https://doi.org/10.21203/rs.3.rs-4993024/v1
This work is licensed under a CC BY 4.0 License
Version 1
posted
You are reading this latest preprint version
This paper reports on the computer design aimed to produce a net power fusion reactor constructed with six high-temperature-superconducting (HTS) electromagnets. Computer simulation results were analyzed using Yushmanov’s magneto-electrostatic-trap (MET) theory [1]. Our spherical design, which we named MET6, has a more efficient shape for plasma confinement than cylindrical designs [2,3]. Yushmanov showed that, considering only particle losses, power-balance (Q) scales linearly with the volume-to-surface ratio of a confined plasma core. Since a sphere has a bigger volume-to-surface ratio than a cylinder of the same volume, our spherical design might be expected to produce larger Q-values than recent cylindrical ones [2,3], but only if its proton-boron bremsstrahlung can be kept to manageable levels.
Yushmanov’s analysis was specialized to reactors burning deuterium plus tritium (DT) fuel. Because of widespread concerns [4] for the long-term availability of the rare isotope tritium, we have chosen to design a reactor burning proton plus boron-11 (pB11) fuel. Compared to tritium, terrestrial boron-11 is so plentiful that supplying it to future generations will not substantially increase the operational costs of the MET6 reactor. Natural boron is 80% comprised of the isotope boron-11. In our analysis we compute reactor power-balance to serve as a figure of merit in comparing various reactor designs. Power-balance is defined as the ratio of useable power output to the drive power input required to confine and heat the fuel plus power to compensate for heat lost to bremsstrahlung radiation. Because of the relatively high 5+ charge state of the boron in the fuel, our power output is nearly overpowered by background bremsstrahlung radiation, leaving very little useable output power from fusion.
Reactors burning DT fuel are little troubled by bremsstrahlung radiation. In our case, we approximate the input power requirement to be just that needed to balance the bremsstrahlung losses. This simplifies the calculation of power input which can now be computed from a textbook formula for radiated power. The resulting power-balance is predicted to be 1.5, just barely break-even. The design still holds merit as a contender for commercial power because the predicted power-balance has a large uncertainty due to limitations of the simulation. A power-balance of 2 would be adequate for plasma ignition. Methods of improving the power-balance in this range are suggested in the closing remarks.
MET6 evolved naturally from its predecessor, the hexahedral Polywell reactor, invented by Robert Bussard in 1989 [5]. Polywell was criticized by Rider [6] and also in the 2018 Jason Report summarizing the U.S. Department of Energy’s ARPAE funding program [7]. The Jason Report stated that Polywell had been “ruled unworkable” referencing Rider [6], claiming that Polywell’s ions’ non-thermal energy distribution would require too much input energy to maintain and the device therefore could not reach unity power-balance. At that time Rider’s general conclusions had been challenged by Rostoker et al. who showed that Rider had assumed operating conditions that contradicted all known reactor designs, making his conclusions arguably useless for practical reactor design.
The design of our MET6 sidesteps the controversy over Polywell’s non-thermal energy distribution. In contrast to Polywell, the plasma in our MET6 reactor is predicted to have a Maxwellian energy distribution, with only a slight departure predicted due to the extreme high energy tail of electrons up-scattering to energies exceeding the height of the trapping potential. This minor departure from Maxwellian distribution was ignored by Yushmanov [1], Dolan [2] and Sporer [3] and likewise in this paper. Pure Maxwellian energy distribution is an assumed feature of the MET6 design.
Seven year ago, one of us (JGR) published an OOPIC designed simulation for a net-power Polywell reactor burning aneutronic pB11 fuel [9]. That design suffered from the same flaw as all previous Polywell designs, in assuming to maintain a non-thermal particle energy distribution. In addition, the design incorporated resistive magnets instead of more compact superconductors, resulting in a reactor size bigger than ITER. The present design of MET6 simulates a plasma in thermal equilibrium and incorporates compact super-conducting magnets to make a more compact break-even design.
Hardware Design – MET6 is a fusion reactor design comprised of six cylindrical, HTS coil-magnets mounted on the faces of a regular hexahedron (i.e. cube). Fig. 1 (a-c) shows PIC [9,10] simulated snapshots of the plasma in the central plane of the cube. Individual macro-particles’ positions are plotted as black dots, electrons in (a), protons in (b), and boron-11 ions in (c). When the simulation program runs, it displays a new version of the particle positions at each time step, i.e. every 73 picoseconds of plasma time. Analysis of the particle movements guides our optimization of the design. Two important parameters of the design are the sizes of the magnets and vacuum tank. The fusion power scales with this size and is proportional to the volume of the vacuum tank, labeled ① in Fig. 1(a). The tank is required to be big enough to leave clearance outside the magnets to accommodate insulating, hollow legs supporting the magnet from the tank. The legs are not expressly simulated, but the simulation predicts clear spaces left to accommodate the legs without causing plasma losses. The invention of the hollow support legs and clear-space to accommodate them is the subject of a 2015 patent on a proposed Polywell reactor [11]. Although our MET6 is not a Polywell, the support legs are the same as described in the patent. The magnets in Fig. 1 appear to be floating in space but 8 spaces clear of plasma are predicted to accommodate support legs in the shadow of the magnets. Legs would reach from the tank wall to each of eight magnet boxes shown.
Fig. 1(a-c) shows the major internal components of the reactor labeled with circled integers as follows: ① labels the square cross section of the cubic vacuum tank, of diameter 4.79m and shown with tick marks. The metal tank is simulated as an electrical conductor held at zero volts. The tank serves as an equipotential for injecting electrons, protons, and boron fuel ions. ② labels the location of an electron emitter, envisioned as a hot filament biased to the same zero voltage as the tank. Starting at time t=0, 15A of zero-energy electrons were generated from the point source depicted as an orange dot located 1 cell width inside the righthand tank wall in panels (a-c). Electrons are accelerated inward by an electric field formed by applying a 2.4 megavolt (MV) DC voltage to 8 conducting magnet boxes serving as an electrical anode. The square boxes labeled ③ and ④ represent the righthand one of the four magnet coils mounted on the four side faces of the cube. Three other pairs of boxes represent the other three magnets surrounding the central plane of the cube. Each coil intersects the central plane in two places, outlined by square conducting boxes. Inside the boxes are simulated HTS wires carrying current appropriate to generate the magnetic field required to confine electrons as shown in Fig. 1(a).
Starting at time t=0, electrons flow continuously into the center of the tank through the open bore of the righthand magnet coil i.e. between squares labeled ③ and ④. The trapped electrons’ density increased with time until it reached the spatially uniform density shown in snapshot (a), made at a plasma time of t=0.38ms. This time was chosen as a stopping point for the simulation. It is the earliest time that the plasma density reached a stable, steady-state condition. As long as the electron drive current was held at 15A the plasma density was constant, independent of time.
This density is still many orders-of-magnitude too low to produce net power. In real-world start up, the electron drive current would be further increased to raise plasma density to net-power levels. In the simulation, the diagnostics at this early time are enough to predict the net-power conditions, making it unnecessary to explicitly simulate beyond t=0.38ms. A key finding is that the shape of the plasma core is constant in simulated time so that the net-power density could be determined by the β=1 condition that defines the same surface at net-power as we see in Fig. 1. With β = (4e-11)nT/B2 (from the Plasma Formulary [21] pg. 29), β is a function of electron density n, temperature T, and magnetic field B, all in cgs units and evaluated at the surface of the plasma core indicated by the double-headed arrow in (a). The β=1 equation was inverted to find the surface electron density n, which also must equal the surface fuel density due to the requirement for quasi-neutrality. The fuel density determines the reactor’s power output and, dividing it by the drive power, the power-balance, which is an important measure of reactor performance.
In a practical reactor, the magnetic field internal to the plasma sheath [throughout the region marked by the arrow in Fig.1 (a)] would become zero as the density rises toward β=1. In the simulation, at the time of Fig. 1, the plasma density was still too low to generate any significant diamagnetic field. To help speed the calculation, the OOPIC simulation was run using the “electrostatic” mode, in which the magnetic field is static in time and computed cell by cell from programmed expressions. To simulate diamagnetism, the t=0 field from the HTS wires was replaced cell-by-cell with a value B=0 for all cells inside a square region of cells spanning the center of the tank. Outside the region the field was computed from textbook formula for field from 8 infinitely-long wires. The width of the field free region is indicated by the double-headed arrow in Fig. 1(a). The diameter of the region was used as an adjustable parameter, adjusted to match the width of the recirculating electron beam to just fill the aperture at the position of the vertical arrow in Fig. 1(a). The wider was the field free region selected, the wider was the beam width at the aperture. The wider the beam at the aperture, the greater the current of electrons scraped off by the aperture.
Ions enter the tank from two point-sources, protons from a source located at position ⑤ in Fig. 1(b) and boron-11 ions from a source at position ⑥ in Fig. 1(c). Positions of the point-sources were fixed at the center of the lefthand tank-wall and the velocities/currents of the ions specified by adjustable parameters in the input file. The velocities were adjusted by trial and error for the ions to pass freely into the core of the reactor along the cusp line through the bore of the lefthand magnet. The currents were adjusted to stabilize the densities of the two species of ions to be equal to each other and each matching half the electrons’ charge density.
The design of ion sources for tokamaks is the subject of others of our papers [12,13]. Injection into tokamaks is more demanding than into MET6 due to the tokamak’s transverse magnetic field, which tends to deflect charged-particle beams. To avoid such deflection, tokamaks require neutral or neutralized beams whereas MET6’s cusp injection works well with simpler charged particle beams, i.e. those simulated here. As seen in Figs. 1(b-c) the ions incoming from the sources ⑤ and ⑥ are focused into narrow beams by the cusp field in the lefthand magnet.
On entering the tank, newborn ions travel through the open bore of the lefthand magnet. The initial energies of the ions, selectively imparted by the sources, were adjusted to match the accelerating potential they will find on entry. Fig. 2(d) shows the steady-state potential along the horizontal midline of the reactor, the path followed by ions on their first pass through center. The force on a newborn ion is proportional to the negative gradient of the potential shown. The ions’ initial velocities are slowed by the rising slope of the potential until they reach the lefthand peak of the potential, where the force is zero. The ions residual momentum carries them into the core, down the falling slope of the potential well. They accelerate to the tank-center where the force is again zero. To the right of the center, they decelerate until they reach the righthand peak, marked “U-ΔU = 2.2MV” in Fig. 2(d). At this point, they reverse direction and turn back to the left. To avoid loss of newborn ions to the righthand tank wall, their initial kinetic energies were adjusted to be slightly less than the potential energy e(U- ΔU), which subtracts from their kinetic energy at the top of the righthand peak.
To visualize the ions’ early history, their positions in velocity-position phase space are displayed in Fig. 2(e); here, gray dots represent protons and black dots represent boron ions. The upper arrow marks the incoming paths of protons and the middle arrow marks the incoming path of the boron ions, respectively. Due to their different charge-to-mass ratio, boron ions follow a different path than protons of the same energy. Due to their higher charge state, boron ions scatter more and become trapped more often than the protons on their first pass. Most protons are seen to reflect from the rising potential they experience in the righthand magnet and then exit the tank to the left after a single pass forward and back through center. The bottom arrow shows the exit path of protons; it is more densely populated than the ions’ exit path just above it because more protons survive trapping on their first pass than ions.
The stopping conditions of the simulation leading to the diagnostics at this t=0.38ms were determined by the condition that the particles’ densities equilibrated in time and the ions’ positive charge density approximately cancel the electrons’ negative charge density at the center of the reactor. The ion injection currents were adjusted to make the sum of the 2 ions’ central density values [graphs (e) and (f)] approximately equal to the electrons’ central density value [in (d)]. This condition made the central plasma approximately charge-neutral, with only a slight excess of negative charge. The visible discrepancy between the level of the proton charge density in (e) and boron charge density in (f) is accidental, due to random statistical variations in cell counts. The curves in (e) and (f) are subject to relatively large statistical noise because they represent charges in a horizontal row of cells only one cell wide in the vertical direction. Each cell contains only a few macroparticles, and consequentially the exact number is subject to relatively large statistical uncertainty. The apparent discrepancy between proton and boron densities does not affect the predicted power-balance of the reactor which is determined by the ratio of bremsstrahlung power to fusion power at the final electron temperature. Making the proton and boron charge densities equal was adopted as a first try for simplicity, later found to be non-optimum for power-balance, as described in the discussion of equations (1) and (2) in the next section of this paper.
Figs. 1(d-f) show the central densities of electrons, protons and boron-ions, plotted on linear scales. The vertical scale of the boron-ions in (f) is shown expanded by a factor of 5 compared to the protons’ scale, to compensate for the 5+ charge state of the fully stripped boron ions. The matching of the charge densities between (e) and (f) is not perfect, only our best effort to match central proton charge density to central boron charge density, both subject to statistical variations.
The electrons formed the 2D electrostatic potential color coded in Fig. 2(a). The (online published) colors code the potential voltage; magenta codes the applied voltage of 2.4MV, found uniformly inside the 8 magnet boxes. Green codes the tank-central value of potential, 1.2MV, half the applied voltage. The color red codes voltages near zero, outside the ring of magnet boxes and extending to the tank walls. Fig. 2(d) shows a 1D graph of the potential voltage made along the horizontal midline of the 2D plot (a). The central value is half the applied voltage, a common feature found in the MET reactor concepts [1]. Although the individual electrons in Fig. 1(a) are in constant motion, the potential they form is static in time from t=0.38ms on. The stopping condition of the simulation was chosen to make the ions’ central density slightly less than that of the electrons. This makes the central potential nearly the same with or without ions injected. At steady state the central charge density is partially neutralized by the trapped ions, but a slight negative charge remains to keep the attractive potential approximately the same as it was before injected ions began to accumulate.
The magnetic field from our six-magnet reactor is dominated by cusps. Fig. 2(b) shows a rendering of the magnetic field’s equal-magnitude surface in a six-magnet Polywell constructed at the University of Sydney [14]. The size of our reactor design is bigger than the Sydney ones, but the shape of the magnetic field is the same. The shape is dominated by 26 cusp lines, a typical three of which are marked by arrows in (b). Charged plasma particles travel freely along the cusp lines, undeflected by the parallel magnetic field vector. This feature allows controlled injection of charged particles from electron, proton, and boron sources mounted at either ends of the horizontal cusp line.
Unfortunately, the 26 cusp lines also provide paths for hot plasma particles to escape the core and become lost on the tank walls. The simulation software is, by design, two-dimensional which requires selecting a plane through the cube in which to simulate and intentionally capture the main features relevant to a 3D plasma. The plane selected is the cube’s central plane, outlined by a square of solid and dashed lines in Fig. 2(b). This plane contains face and edge cusps, typical ones of which are indicated by the two lower arrows in (b), but none of the cube-corner-vertices, a typical one of which is indicated by the upper, left-pointing arrow. The cube’s 8 corner-vertices are out of the central plane and therefore their losses are not simulated. By this omission we ignore the losses through the cube-corner-vertices but include the losses through the more lossy face and edge, i.e. in-plane, cusps. Luckily, the rate of loss of these omitted cusps, proportional to the area of the circle at the head of the left-pointing arrow, are smaller than those of the in-plane cusps, proportional to the areas of the cusps pointed to by the lower two arrows. By including the more lossy face and edge cusps in the simulation, while omitting the less lossy cube-corner-vertex cusps, we tend to overestimate the cusp losses compared to what they would be in a more realistic, i.e. 3D, simulation. Because the bremsstrahlung losses are found to overpower the cusp losses, the differences between losses rates through face, edge, and corner cusp losses do not add to the uncertainty of our estimated power-balance.
For 2D simulation, we omit the field from the top and bottom coil magnets. The remaining four in-plane magnets were simplified by replacing the curved conductors by straight wires positioned at the same intersection points as the actual curved coils. Our simplified magnetic field is simulated as if generated by electric currents in eight straight wires mounted perpendicular to the simulated plane and intersecting the plane at the centers of the 8 magnet boxes. The magnet boxes are simulated as metal conductors biased to 2.4MV. The 8 biased, conducting boxes provide the electric field which combines with the magnetic field from the wires to confine particles.
The polarity of the currents is arranged to generate magnetic field vectors pointing inward at the center-bores of each of the 4 in-plane coil magnets. The magnitude of the bore-field was adjusted to be 25 Teslas (T) at the center of an isolated magnet. This is in keeping with the maximum field available in commercial HTS magnet coils [15]. When four magnets are assembled into a square, opposing fields from adjacent coils partially cancel to reduce the bore field in each magnet to 18T. The spacing of the magnets was selected to produce the same field magnitude in the corners between the coils as in the bores but pointing outward from the tank’s center instead of inward. Static magnetic field at each PIC cell location inside the tank was computed from a textbook formula [16] for the in-plane component fields from 8 infinitely-long straight wires. Horizontal and vertical field components were computed and stored in two 50x50 arrays, one for x-field and one for y-field. The eight biased, conducting magnet boxes shown, plus eight HTS wires simulated (but not shown), combine to simulate the electric and magnetic fields from the four superconducting coils mounted on the side faces of the hexahedral reactor. Fields from the top and bottom coils do not contribute in the central plane because the opposing fields from these two coils cancel each other in the median plane due to the up-down symmetry of the cube.
Fig. 2(c) shows the x-component of the magnetic field from the 8 wires, plotted along the horizontal midline of the tank. As expected the magnitude of the field is maximum (18T) in the bores of the two magnets on the midline and zero at the center of the tank. The field shown is the applied field which would exist prior to the injection of plasma. At steady state the plasma modifies the central field by adding a diamagnetic current which suppresses the applied field in the central region, expanding the cavity volume that contains the fusion fuel.
As a check that the top and bottom coils could realistically be omitted from the 2D simulation, a separate simulation of applied magnetic field was made with the public domain 3D simulation packages WarpX [17] and MagPyLib [18]. The 3D simulation included six circular coil magnets on the six faces of the cube, including top and bottom coils the same size and field-strength as those in the MET6 design. In the simulation plane the magnetic field from the 3D simulation was found to be qualitatively the same as the field from the eight straight wires, approximately matching the midline field shown in Fig. 2(c). This gives us confidence that our simulated magnetic field from the eight wires captures the essential details of the magnetic field that would be found in a real cubic reactor. Although the WarpX simulation contained a realistic magnetic field, the results have not yet progressed to a level that could be used for further quantitative analysis. Our WarpX simulation so far lacks the essential anode bias which forms an electric field that recirculates most of the escaping electrons back into the core whenever they happen to exit the central region through a cusp. Without the bias potential the plasma would have an unrealistic shape and excessive losses.
Fig. 1(d) shows a graph of the electron’s particle density plotted along the midline of the 2D distribution in (a). The units of the vertical axis are electrons per square meter. Two peaks of concentrated density appear, centered in the bores of the left and righthand magnets. The righthand peak is marked by an arrow. The peaks show concentrations of trapped low-energy electrons which have down-scattered in energy to become trapped in local minima of the electric potential generated by the bulk of electrons. A dense concentration of these trapped electrons would be detrimental to the efficient operation of the reactor. By their negative charge concentration, they repel and rob energy from the newborn electrons entering from the source at position ②. To reduce the trapped electrons’ concentration a special electrode was simulated in the bore of the righthand magnet, as shown at the head of the arrow in (a). This electrode was simulated as a thin metal disk with a circular aperture drilled in its center. The electrode was invented to shape the electric field in the bore of the magnet and to extract a portion of the trapped electrons which hit it. The current of electrons hitting the electrode was adjusted to be a significant fraction of the injected electrons current, the minimum fraction that was found to halt the growth of the peaks.
The fraction of electrons hitting the extra electrode was adjusted by varying the width of the B=0 imposed field-free region; the larger the width the larger the fraction. Increasing the fraction of electrons hitting the aperture was found to reduce the height of the peaks in (d), but at the cost of reducing the power-balance. The larger the current selected, the smaller was the fraction of trapped electrons populating the peaks in Fig. 1(d). The extra electrode provides an adjustible mechanism of “pumping-out” trapped electrons, as described by Yushmanov in his 1981 paper:
“To keep the gaps 'empty' it is necessary to create in them an artificial pumping-out mechanism which has a number of specific characteristics including selectivity with respect to the trapped electrons and a carefully adjusted intensity within fixed narrow limits. No specific technical ideas for methods of creating such a pumping-out mechanism have as yet been put forward.” [19]
Our placement of the aperture in the righthand magnet satisfies the need to actively pump low-energy electrons, perhaps the first solution proposed yet. The steady and uniform density of trapped electrons and ions in Fig. 1 shows that this solution works well enough. Little effort was made to optimize the adjusted electron current loss on the aperture. Adjusting the current and power lost by the electrons hitting the electrode would impact the predicted power-balance Q, but only slightly compared to the power consumed by compensating for the bremsstrahlung losses. Such optimization of the pumping current is a detail of the design we chose to leave for future fine tuning of the design parameters.
In addition to providing for the pumping out of slow electrons, the extra electrode has a second important function. It shapes the potential along the path of the newborn ions to increase their probability of being trapped on their first passage through the core. The potential function shown in Fig. 2(d) exhibits two peaks. The range of heights of the peaks define a range of injection energies suitable to choose for the incoming ions. Their injected energies need to be in a range higher than the eU value at the lefthand peak and lower than the eU value at the righthand peak. Ions in this energy range pass freely in through the lefthand cusp and then are reflected by the potential at the righthand cusp. This gives each ion the maximum probability of being trapped in passing twice through the bulk of the trapped plasma. Those ions not trapped on their round-trip passage exit the core to the left to be lost on the lefthand tank wall. The bottom-most of the three arrows in Fig. 2(e) marks their exit path. The density of ions along the bottom-most arrow in Fig. 2(e) differs for protons and boron-ions. This is due to the different probabilities of capture of ions on their first pass forward and back through center. Because most of the newborn ions exit after a single pass, the efficiency of injection of ions from the sources on the lefthand wall will improve as the central plasma density grows with time. At this early time (t=0.38ms) the injected currents of protons and borons were adjusted to be 2.5A and 2.75A, respectively. These currents were set by trial and error to just balance the loss rate of the newborn ions which happen to survive an initial round-trip passage through the trapped plasma.
The distribution of the ions in Fig. 2(e) reflects the history of fueling the reactor over the first 0.38ms after time zero. As simulated time increases, the fueling process fills in the central region of Fig. 2(e) and the thermalization process degrades the ions’ initial velocities to smaller values than those shown in 2(e). After many milliseconds the dots would accumulate near the center of the phase-space (velocity vs. position) distribution, providing a way to estimate the final temperature of the plasma after it reaches thermal equilibrium. Rather than undertaking an impractically long OOPIC simulation run to estimate this temperature, we relied on the analysis of Yushmanov [1], which found that thermalization lowered the plasma temperature from the applied eU to an equilibrium temperature equal to “about 1/16 of e(U- ΔU),” where ΔU is the potential depression in the gaps due to electrons’ space-charge accumulating in the gaps. From the horizontal arrow in our Fig. 2(d), the temperature of the thermalized plasma can be estimated as T = e(2.2MV)/16 = 140keV. The quantitative effect of thermalization is to lower the plasma temperature by a factor of 16. This estimate of thermalization loss is a principal contribution of Yushmanov’s theory, which distinguishes our MET6 reactor design from the earlier Polywell designs. In Polywell, ions’ energies were mistakenly assumed to maintain their original eU temperature until they fused.
Estimating the Power-balance Q from the PIC Diagnostics –The standard definition of power-balance is the ratio of fusion power-output to drive-power input. The fusion power-output is calculated from the textbook [20] formula:
PpB = np nB <σv> EpB V (1)
where np is the protons’ particle density, nB is the B11 ions’ particle density, <σv> is the fusion reactivity averaged over the thermal velocities (v) of the T =140keV thermal plasma, EpB = 8.7Mev, the energy yield of one pB11 fusion event (from the Plasma Formulary [21] pg.44), and V is the volume of the fusing plasma. From the diameter of the cubic plasma core shown in Fig. 2(f), V = (1.53m)3 = 3.6m3 = 3.6e6cm3.
To simulate the plasma’s natural quasi-neutrality, the magnitude of charge density of the ions must approximately equal the magnitude of the charge density of the electrons. We initially chose to divide equally the positive charge density between protons and boron-ions. (This choice was subsequently adjusted to reduce the bremsstrahlung contribution to the power-balance as simulated.) To make equal contributions from protons and borons, the 5+ charge state of the fully stripped boron ions implies particle densities nB = np/5 = n/10, where n is the electrons’ density derived from the β=1 condition at the surface of the plasma core. This density was calculated by inverting the equation β=1 (from the plasma formulary [21] pg. 29) to yield the following expression for the electron density at the surface:
n = Bs2 / [(4e-11) Ts], (2)
where Bs is the magnitude of the magnetic field vector at the surface of the plasma core and Ts the electron-energy, both expressed in cgs units. The resulting units of n are also cgs units, cm-3.
In Fig. 1(a), the double-headed, gray arrow points to selected points on the surface of the plasma core. These points were selected where the surface of the plasma core is parallel to the y-axis of the simulation plane. In the online version of Fig. 1(a) electrons at the surface were seen to move parallel to the y-axis, indicating that the magnetic vector points parallel to the y-axis at these points. Thus, the By component of field is the only non-zero component of the field’s three-vector. As one of its standard diagnostic displays, the simulation software displays the cell-by-cell values of the 3 components of magnetic field, including By as a 2D array. Fig. 2(f) shows a graph of By as a function of x, extracted from the 2D array along the line marked by the double-headed arrow. The double-headed arrow from Fig. 1(a) is reproduced in (f) and marked “1.53m.” A second, left-pointing arrow in (f) marks the value of By at the tip of the double arrow, i.e. at the surface of the plasma core. This surface value is Bs = By = 8T, which equals 80 kilogauss (8e4) in cgs units. The electron energy at the surface is taken to be the most probable energy of the thermal distribution, given by the plasmas’ temperature, 140keV. Inserting these values into eq. (2) yields our estimate of electron density:
n = Bs2 / [4e-11 T] = (8e4)2 / [(4e-11) (1.4e5)] cm-3 = 1.3e15cm-3. (3)
The remaining factor needed to evaluate the fusion power from eq. (1) is the average reactivity, <σv>, at the plasma temperature of T = 140keV. Fig. 3 shows a graph [22] of this experimental quantity as a function of temperature for different fuels. The curve for p+B11 fuel is colored red in the online publication and pointed to by the downward diagonal arrow. Heavy (blue) arrows mark the temperature 140keV and the associated reactivity <σv> is 1.0e-16 cm-3 s-1. Substituting into eq. (1) yields the simulated estimate of the reactor’s nuclear fusion output:
PpB = np nB <σv> EpB V = (n/2) (n/10) <σv> EpB V (4)
= [(1.3e15)2/20)] (1.0e-16) (8.7e6) (3.6e6) eV/s = 2.6e26 eV/s = 42 megawatts,
where the symbols are defined previously at the introductions of eq. (1) and eq. (2). This power-output
estimate forms the numerator of the power-balance ratio. The denominator is the bremsstrahlung power output which dominates the requirement for drive power input. The bremsstrahlung power estimate is taken from Section 2.64 of the Glasstone and Lovberg textbook [23],
Pbrem = 5.35e-31 n (np+Z2nB) (TkeV )1/2 V, (5)
5.35e-31 n (n/2+25n/10) V = 5.35e-31 (1.3e15)2 (0.5+2.5) (140)1/2 (3.6e6) = 116 megawatts.
where n, np, nB, V are as in eq. (4) above, Z = 5 is the charge of the fully stripped boron ions, and TkeV = 140 is the electrons thermal energy in units of keV. This drive power exceeds the fusion power output in eq. (4), leading to an estimate for power-balance less than the desired minimum value of unity:
Q = 42MW / 116MW = 0.36 (6)
Improving the predicted power-balance can be accomplished by changing the fuel mixture to have a reduced percentage of boron relative to hydrogen. The bremsstrahlung power output contains the factor n(np+Z2nB), which is subject to the constraint np + 5nB = n, coming from the required balance between negative electron charge-density n and positive ion charge-density comprised of protons plus Z=5 boron ions. Reducing the boron density nB has a magnified effect on the bremsstrahlung power due to the Z2 factor (=25) multiplying it. Assume we reduce the boron density from n/10 to n/25. This would reduce the n(np+Z2nB) factor in the bremsstrahlung power from 3n2 to n2, a factor of 3. Reducing nB would also reduce the (np nB) factor in the numerator of power-balance from n2/20 to n2/25. The net effect on power-balance would be to increase it by the factor (20/25) / (1/3) = 60/25 = 2.4, which yields the adjusted power-balance as:
Q1= 2.4 x 0.36 = 0.9 (7)
still too small for practical net power.
Possible Recovery of Bremsstrahlung Emitted Energy – As pointed out in Section 2.73 of Glasstone and Lovberg [23], it is possible in principle to recover the energy emitted as radiation if the confined plasma is optically thick to the emitted photons. In this case x-rays of the characteristic bremsstrahlung wavelength convert to heat before they exit the core. X-rays of a characteristic frequency would be attenuated exponentially with a mean free path in centimeters estimated by Spitzer [24]:
mfp ≡ 7e-5 T1/2 ν3 / n2, (8)
where T is the plasma temperature in keV, ν is the x-ray frequency in s-1, and n is the plasma’s particle density in cm-3. From Fig. 2.8 of Glasstone and Lovberg [23], the most probable x-ray wave length is 0.1 angstroms or 1e-9cm. Converting this photon wave-length to photon frequency yields ν = c/0.1e-8 = (3e10cm/s)/(1e-9cm) = 3e19s-1. Substituting the simulated values of T, ν, and n into eq. (8), yields the mean free path of the bremsstrahlung photons estimate to be:
mfp = 7e-5 (140)1/2 (3e19)3 / (1.3e15)2 = 7e-5 (12) 27e57 / (1.7e30) = 1.3e25 cm.
This path length is astronomically long compared to the 5m size of the reactor. In other words, the confined plasma is optically thin to the x-rays, not optically thick, meaning that a negligible amount of the emitted bremsstrahlung energy can be captured in the plasma.
We next considered the possibility of converting the escaping bremsstrahlung energy into electricity, and then inputting this generated electrical power as heat to compensate for the cooling due radiation emission. For our plasma temperature of T=140keV, the emitted bremsstrahlung radiation is in the form of a broad x-ray spectrum (as shown in Fig. 2.8 of [23]) with median energy 100keV. X-rays of these energies can be converted into heat by capturing them in a boiling-water jacket surrounding the vacuum tank. The boiling water would then be used to drive a steam turbine to generate electricity via the Rankine cycle with an efficiency around 40% [25]. The detailed design of the water jacket and turbine system are beyond the scope of this paper, but as a reality check we can estimate how thick the water jacket would need to be to convert the emerging bremsstrahlung x-rays to heat. X-rays are attenuated exponentially in normal materials. The fraction of emitted radiation absorbed in a water jacket can be estimated by the formula:
(I0 – I) / I0 = exp(-μx), (9)
where I and I0 are the intensities of radiation before and after the water absorber, μ is the known radiation length of 100keV x-rays in water, and x is the thickness of the water along the path of the x-ray i.e. the radial direction from the reactor’s center. The radiation length of 100keV x-rays in water, taken from standard tables [26] is μ = (1.7e-1cm2/g) (1.0g/cm3) = 0.17cm-1. A modest water jacket just 24cm thick would produce μx = .017x24 = 4 and would therefore absorb a fraction of the emitted x-rays 1-e-4 = 0.98. Virtually all the emerging x-rays would be recycled by a water jacket 24cm thick. In order to reach the water to be recycled, x-rays must pass through the magnets and vacuum tank walls. Assuming these losses can be made negligible compared to the water, the water’s Rankine efficiency converts 40% of the x-ray energy to electricity, which is then is injected to heat the plasma. This would decrease the bremsstrahlung loss which is the denominator of the power-balance. The adjusted power-balance estimate then becomes:
Q2 = Q1 / (1-0.40) = 0.9 / (1-0.40) = 1.5. (10)
A reliable PIC simulation has been coupled to the classic MET analysis method of Yushmanov [1] to design a compact, quasi-spherical reactor, MET6. Table 1 shows the final simulation parameters of MET6 compared to Sporer’s cylindrical design [3].
Table 1 – Comparison of Sporer’s HTS reactor with our MET6 reactor
The two reactors have very different power balances Q, the main metric for performance. Although the MET6 design is predicted to have poorer power-balance, it burns a more plentiful fuel than Sporer’s design, giving it a massive advantage in sustainability. Plus it does not produce neutrons. MET6 burns hydrogen plus boron fuel, economically available from terrestrial sources. Sporer’s design, like most of the many tokamak designs (ITER, SPARC, EAST, etc.), would burn tritium, scarcely available and very expensive from a few aging Candu reactors [27].
As the plasma density rises during startup, electrons’ space charge accumulates in the cusps unless suppressed by some active means [19]. The accumulation of space charge would be fatal, leading to the loss of the potential-well due to the reduced energy of incoming electrons after plowing through the trapped charge. A major discovery of this work is that scraping the exterior surface of the bounding sheath can eliminate this loss of the potential-well. Suppression of the space charge trapped in the cusps required 40% of the incoming electrons’ current be lost by scraping the gap-confining aperture. The surface layer of the sheath was scraped off by the edges of the aperture. The fraction of electrons scraped was adjusted by adjusting the diameter of the aperture or alternately adjusting the diameter of the diamagnetic exclusion region. Adjusting the diameter of the diamagnetic regions gave finer control than adjusting the diameter of the aperture because the aperture was already set to be its minimum nonzero value of 2 PIC cells.
Varying the width of the diamagnetic region to adjust the level of scraping current utilizes a defocusing effect, which increases with the length of the beam’s transit through the field-free region. The passage of the recirculating electron beam through the wider or narrower B=0 region defocused the beam more or less depending on the width of the region. A defocused beam resulted in more loss of electrons at the point of the arrow in Fig. 1(a). The larger the diameter of the exclusion region the larger the fraction of current scraped. The amount of scraping was adjusted using the feedback from the simulation’s diagnostic shown in Fig. 1(d). The trapped electrons form two peaks located inside the lefthand and righthand magnets’ bores. [For ease of identification we note here that the righthand peak is marked by an arrow in Fig. 1(d)]. With no scraping the ratio of the peak height to the density at the center continued to grow with simulated time, gradually reducing the height of the confining peaks in the potential in Fig. 2(d) which in turn released the trapped ions. When the amount of scraping was adjusted to be above a certain threshold the ratio of peak height to the central density became constant in time. In other words, the scraping halted the build-up of trapped electrons, in effect pumping out low energy electrons from the cusp traps at the same rate they are trapped.
The amount of power required for scraping depends on the detailed structure of the plasma inside the thin sheath marking the boundary between field-free interior and plasma-free exterior to the central ball of plasma. The details of the spatial structure in the sheath is not well simulated by our choice of PIC cell size much larger than the electron gyro-radius. A PIC technique with better resolution is published by Park et al. [28], which reports a massive PIC simulation consuming more than a million CPU hours to run. The device simulated by Park et al. is a picket fence device, not a candidate for a net-power reactor because it lacks the essential blocking electric field required to return electrons that leak out the cusps. In our design the blocking electrode is the grounded tank. The grounded tank plus the positive magnet bias voltage create an electric field that repels electrons as they approach the tank from inside. Without such repulsion, the electrons that leak out the cusps would hit the tank and be lost, which would ruin the power-balance. Although Park et al. do not claim to present a workable reactor design, they do demonstrate that resolving the detailed structure of the simulated sheath would require choosing the PIC cell size to be at least as small as the electron gyro-radius in the sheath.
With our relatively large cell size, we do not resolve the detailed spatial structure in the sheath. But such good resolution is not required to validate the MET6 reactor’s design using our simulation. Our limited simulations suggest a surface layer thicker than the real-world layer would be. This may affect the level of the current required to pump out the trapped electrons from the cusps, but not the temperature of the plasma on which the power-balance depends. The sufficient level of the pumping current was set by trial and error adjustment of the density of trapped electrons in Fig. 1(d). Because bremsstrahlung dominates the power requirement, the level of power invested in electron pumping has negligible effect on the reactor’s predicted power-balance.
To see if our conclusions depend on PIC cell size, a separate simulation run was made with the same device parameters shown in Table 1 but with PIC cell size ½ the cell-size used so far. The higher-resolution produced the same confining potential (U-ΔU), the same electron temperature T, and the same electron density n, as the lower-resolution simulation. With these factors held constant, repeating the analysis would suggest the same power-balance as predicted with the lower resolution.
The improved-resolution simulation, using 100x100 cells to span the simulation plane, was not used for routine refinement of the MET6 design because it challenged the workflow timescales needed to test different parameters. Each regular run of the OOPIC software, using 50x50 cells, required about 50h of desk-top processor time to produce a result. The test run with 2x higher resolution took sixteen times more processor time per run. The higher resolution run required computing the fields in four times more cells to cover the simulated area, and also required the weights of the macroparticles be decreased by a factor of four to maintain the same statistical accuracy of simulated density in each cell. Individual cell counts are used to produce the 1D sections of 2D density displays, like in Figs. 1(d-f). The processor time necessarily increased as the inverse 4th power of cell size.
Particle in Cell (PIC) simulation has been used to design a quasi-spherical, net-power, aneutronic, fusion reactor incorporating high temperature superconducting (HTS) magnets. The reactor is predicted to have power-balance Q = 1.5, with a predicted power-balance uncertainty estimate to be less than 10. Thus, the MET6 design is optimistically considered as a viable candidate for possible commercial fusion power application.
Before bringing this design to the stage of a prototype, several engineering challenges must be faced in theory. The first is the design of superconducting ring magnets and their mounting hardware, which is well assisted by the construction of ITER in France, where the central solenoid will contain a stack of six ring magnets of approximately the same size and magnetic field as the magnets proposed here for MET6. The first few magnets for ITER’s central solenoid have been delivered and successfully tested [29]. The successful development process for these magnets can be replicated for the MET6 magnets. Both maximum field and supporting materials’ stresses are of common concern both in ITER and in MET6.
The power-balance of Q=1.5 may be too small to predict useable net power. The question of how big the simulated power-balance needs to be has different answers depending on the level of uncertainty in the calculation and whether the reactor operates as an ignition device or as a power amplifier. Moreover, any practical concept that can attain near unity power balance will also have the effect of stimulating new scrutiny, especially related to improved methods of analysis and design. As a power amplifier the conversion of gross power into net power is subject to efficiency factors for converting nuclear energy into electrical energy outside the reactor and then back to heat the plasma inside the reactor. The efficiency losses to and from electrical energy would require a power-balance on the order of Q=5-10 [19]. On the other hand, for an ignition device, the inefficiencies of electrical conversion are avoided. In this case, a power-balance as little as Q=2 might be useful to obtain ignition. Here, the MET6’s quasi-spherical shape and cusp-confinement are an advantage over devices of other shapes like tokamaks and field-reversed-confinement (linear) reactors [30]. The longer path length for the fusion products inside our spherical confinement volume makes for more efficient conversion of the emitted alpha’s power into heat without losses. Lossless conversion of fusion energy into plasma heating is the essential advantage of ignition devices over power amplifiers.
A power-balance of Q=2 could arguably be obtained by reducing the electron energy to reduce the bremsstrahlung radiation relative to the fusion energy. This is theoretically possible by forming a deeper potential well than that shown in Fig. 2(d). Our simulations show that the central potential voltage and the resulting confined electron energy T is smaller for larger injected electron current. The T1/2 dependence of emitted bremsstrahlung power on electron temperature in eq. (5) should allow reducing bremsstrahlung by lowering the trapped electrons’ energy.
The path to simulating improved power-balance in MET6 is the subject of ongoing research by our group. Currently the public domain simulation package WarpX [17] is being pursued for this quest and has 3D capability, improving on our OOPIC software which is fundamentally 2D.
The ultimate proof of the conceptual design will require continuing investigations leading to a scale-model, prototype reactor of the MET6 design.
Author Contribution
J.G.R. wrote the main manuscript text and prepared figures 1-3. All authors reviewed the manuscript.
Acknowledgment –
The authors thank Brendan Sporer for valuable discussions.
- E.E. Yushmanov, The power gain factor Q of an ideal magneto-electrostatic fusion reactor, Nucl. Fusion 20 3 (1980)
- T.J. Dolan, Magnetic electrostatic plasma confinement, Plasma Phys. Control. Fusion 3, 1539 (1994)
- B.J. Sporer, Analysis of Two Fusion Reactor Designs Based on Magnetic Electrostatic Plasma Confinement, J. Fusion Energ 41 5 (2022)
- R. Pearson, A.B. Antoniazzi, W.J. Nuttall, Tritium supply and use: a key issue for the development of nuclear fusion energy, Fusion Engineering and Design 136 (2018)
- R.W. Bussard, Method and apparatus for controlling charged particles, US patent 4826646 (1989)
- T.H. Rider, Fundamental limitations on plasma fusion systems not in thermodynamic equilibrium, Physics of Plasmas4, 1039 (1997)
- JASON The MITRE Corporation, Prospects for Low Cost Fusion Development, report JSR-18-01 (2018)
- N. Rostoker, A. Qerushi, and M. Binderbauer, Colliding Beam Fusion Reactors, Journal of Fusion Energy 22 2 (June 2003)
- Joel G. Rogers, Aneutronic Polywell Reactor, 19th US-Japan Workshop on Inertial Electrostatic Confinement Fusion, Kansai University, Osaka, Japan, downloaded from: https://www.researchgate.net/publication/321757410_Aneutronic_Polywell_Reactor, July 11 (2024)
- Panels of Fig. 1 were produced by a particle-in-cell (PIC) simulation program OOPIC Pro, licensed from Tech-X Corporation of Boulder, Colorado, USA in 2016. OOPIC Pro is no longer supported by Tech-X, but an equivalent public domain version of the software is available from the University of Michigan under the name XOOPIC [10]. We compared XOOPIC with OOPIC Pro and found the results identical within statistical uncertainties.
- Prof. J.P. Verboncoeur, XOOPIC, The Plasma Theory and Simulation Group, http://ptsg.egr.msu.edu/ (2023)
- Joel G. Rogers, Modular Apparatus for Confining a Plasma, U.S. Patent 9,082,517B2 (2015)
- F.J. Wessel, J.J. Song, H.U. Rahmana), G. Yur, N. Rostoker, and R.S. Whitea, Propagation of Neutralized Plasma Beams, Physics of Fluids B: Plasma Physics, 2 6 (1990)
- Frank J. Wessel, Andrew Egly, and Joel Rogers, Neutralized, Intense-Ion Beams for Fusion, Plasma Phys. Control. Fusion66 065028 (2024)
- M. Carr, D. Gummersall, S. Cornish, and J. Khachan, Low beta confinement in a Polywell modeled with conventional point cusp theories “Fig. 12,” Physics of Plasmas 18 11 (2011)
- Commonwealth Fusion Systems, 117 Hospital Rd., Devens, MA 01434. https://cfs.energy/technology/hts-magnets
- John David Jackson, Classical Electrodynamics, John Wiley and Sons, Inc., New York, (1962)
- This research used the open-source particle-in-cell code WarpX https://github.com/ECP-WarpX/WarpX, primarily funded by the US DOE Exascale Computing Project. Primary WarpX contributors are with LBNL, LLNL, CEA-LIDYL, SLAC, DESY, CERN, and TAE Technologies. We acknowledge all WarpX contributors.
- Michael Ortner and Lucas Gabriel Coliado Bandeira, Magpylib: A free Python package for magnetic field computation, https://www.sciencedirect.com/science/article/pii/S2352711020300170 (2020)
- E.E. Yushmanov, The influence of electron capture in gaps on the efficiency of the magneto-electrostatic trap, Nucl. Fusion 21 329 (1981)
- Samuel Glasstone and Ralph R. Lovberg, Controlled Thermonuclear Reactions, an Introduction to Theory and Experiment eq. (2.10), D. van Nostrand Company, Inc., Princeton, New Jersey (1960)
- J.D. Huba, NRL Plasma Formulary, Beam Physics Division, Plasma Physics Branch, Naval Research Laboratory, Washington, D.C. (2009)
- Christian Hill, Nuclear fusion cross sections and reactivities, downloaded 2024/4/18 from website “https://scipython.com/blog/nuclear-fusion-cross-sections/” (2024)
- Samuel Glasstone and Ralph R. Lovberg, Controlled Thermonuclear Reactions, an Introduction to Theory and Experiment, D. van Nostrand Company, Inc., Princeton, New Jersey (1960)
- Lyman Spitzer, Physics of Fully Ionized Gases, Interscience Publishers, Inc. pg. 80 (1956)
- https://en.wikipedia.org/wiki/Rankine_cycle
- https://physics.nist.gov/PhysRefData/XrayMassCoef/ComTab/water.html
- R. Pearson, A.B. Antoniazzi, W.J. Nuttall, Tritium supply and use: a key issue for the development of nuclear fusion energy, Fusion Engineering and Design 136 (2018)
- J. Park, G. Lapenta, D. Gonzalez-Herrero and N. Krall, Discovery of an Electron Gyroradius Scale Current Layer: Its Relevance to Magnetic Fusion Energy, Earth’s Magnetosphere and Sunspots, 6 74 (2019)
- Nicolai Martovetsky, Kevin Freudenberg, Graham Rossano, Kenneth Khumthong, Nikolai Norausky, Ed Ortiz, Zbigniew Piec, Kurt Schaubel, Jeffrey Sheeron, Gordon Kaskin, John Smith, Roberto Bonifetto, Roberto Zanino, Andrea Zappatore, Thierry Schild, Florent Gauthier, Cornelis Jong, Yuri Ilyin, Leonard Myatt, Kristen Cochran, Marco Breschi and Alan Langhorn, Testing of the ITER Central Solenoid Modules, https://www.osti.gov/servlets/purl/1923232 (2022)
- R.M. Magee, K. Ogawa, T. Tajima, I. Alfrey, H. Gota, P. McCarroll, S. Ohdachi, M. Isobe, S. Kamio, V. Klumper, H. Nuga, M. Shoji, S. Ziaei, M.W. Binderbauer, and M. Osakabe, First measurements of p11B fusion in a magnetically confined plasma, Nature Communications 14:955 (2023)
No competing interests reported.