2.1 Geometry

The first step is to construct a geometry for the PM array. For this section, we treat each PM as an ideal dipole point source. An ideal dipole has six degrees of freedom: three for position, two for spherical orientation and one for magnitude. Our optimisation fixes magnet positions for simplicity. If the position were a free variable, the optimiser would need to constrain neighbouring magnets from intersecting as they move. By fixing the positions, we can initialise an array without intersections and optimise only the magnetisation. To define the initial point cloud of dipoles without intersections, it is necessary to know the 3D build and orientation of the finite-volume magnets.

Some authors choose to align the magnets with a cylindrical grid (Hammond et al. Reference Hammond, Zhu, Korrigan, Gates, Lown, Mercurio, Qian and Zarnstorff2022). Others use a cartesian coordinate system (Lu et al. Reference Lu, Xu, Chen, Zhang, Chen, Ye, Guo and Wan2021a). For both cases, the coordinates do not align with a toroidal plasma, and multiple orientations are required to facilitate shaping. The MUSE magnets are chosen to align with a toroidal grid. All dipoles are aligned normal to a circular cross-section torus concentric with the VV. The magnetisation vector points either straight in or straight out, and this information is stored as a signed magnitude $\rho$ such that $\boldsymbol {m} = \rho \boldsymbol {m}_0$.Footnote ¹ Since the orientation and position are fixed by geometry, the optimiser manipulates only a scalar magnitude for each magnet.

It is desirable to pack the magnetisation volume as densely as possible. There are more magnets on the outboard side than on the inboard side due to toroidal curvature. Accounting for this, we divide the magnet volume into $N_\theta$ poloidal rings such that each ring has variable $N_\zeta = N_\zeta (\theta )$ magnet positions. At each surface coordinate $(\zeta,\theta )$ we initialise a magnet tower of $N_s$ radial slices. A ‘PM tower’ refers to the stack of magnets at a single ($\theta,\zeta$) position. In the context of optimisation, every tower consists of $N_s$ slices (see figure 4). In the context of assembly, however, the slices are merged into a single physical piece for each tower (see § 5). Note that $N_\zeta (\theta )$ is largest for the outer ring ($\theta = 0$) and smallest for the smallest ring ($\theta = {\rm \pi}$). Figure 2 shows the non-smooth function of the poloidal angle that arises from packing the magnets as densely as possible. This discrete function takes into account the circumference of a poloidal ring $C(\theta ) = 2 {\rm \pi}(R + a \sin \theta )$, the width of a magnet $w$, and the minimum structural gap which appears between magnets $d_m$. We must also take into account an assembly gap between support structures $d_s$. Since MUSE has two field periods, the stellarator is divided into four identical half periods. Let us define the boundaries of these ‘physical half periods’ to be our coordinate system's $x$ and $y$ axes. However, for the practical purpose of force balance, we will find it useful to let the PM structures straddle the half-period boundary symmetrically. As a result, the ‘PM holder quadrants’ are shifted $45^{\circ }$ out of phase, as shown in figure 3. This is also a useful feature for quality checks during assembly because the physical parts exhibit odd symmetry (see figure 19). We emphasise that the information in a two-field-period stellarator is entirely contained in a $90^{\circ }$ half period bounded between symmetry planes. We intentionally bisect the half-period into two pairs of distinct PM holders, each centred on one symmetry plane. As a consequence, the appropriate length for toroidal discretisation is actually an eighth of the circumference, $\varDelta = \tfrac {1}{2} C(\theta ) / (2 N_p) = C/8$. Then the expression for the number of toroidal positions in a poloidal ring is a floor function

(2.1)\begin{equation} N_\zeta (\theta) = \left \lfloor \frac{ \varDelta - d_s} {w + d_m} \right \rfloor. \end{equation}

This guarantees a gap of at least $d_m$ between every magnet. Due to toroidal curvature, the minimum minor radius $a$ differs for the inner and outer poloidal angles. On the inside the densest ring is determined by the base of the tower; on the outside it is determined by the top of the tower. Figure 3 highlights this feature, and it shows the relative size between the PM grid and the plasma shape. In the final design, $d_m = 1$ mm is the minimum gap between all magnet edges and $d_s = 2 d_m$. Because the space $d_s$ is removed from both sides, the physical gap between PM holders is at least 4 mm.

Figure 2.

Illustrations of the distribution of magnet positions and actual magnets. (a) The number of toroidal magnet positions as a function of poloidal angle $N_\zeta = N_\zeta (\theta )$. (b) The full 2D grid of magnet positions as a function of toroidal and poloidal angles $(\zeta,\theta )$, including the assembly gap and port holes. (c) The distribution of final MUSE magnets. Red and blue dots mark towers that have only positive or only negatively oriented magnets. The purple dots mark special locations where a single tower has both orientations. These are collected near the vertical port holes at toroidal angle $\zeta = 0$. An orthogonal coordinate grid is indicated in gold and grey contours. It is used to uniquely index each magnet position. Every 5 lines are labelled and bolded to facilitate counting. Dashed orange lines show the $\zeta = {\rm \pi}/2$ boundary that parts the $X$ and $Y$ quadrants and the $\theta = \pm {\rm \pi}/2$ boundary that approximately parts the inboard and outboard subsegments. Here, the poloidal axis is shifted to $(-{\rm \pi}, {\rm \pi})$ such that the densest region (the outboard midplane) is centred in the middle. Readers may wish to contrast this representation with alternatives in figures 8 and 10.

Figure 3.

(a) Top-down view of MUSE shows the PMs rendered as 3D blocks surrounding the VV. This plot highlights four identical physics half periods, surrounded by two pairs of identical PM holder quadrants, which each mirror a quarter period. (b) A constant toroidal angle cut, visualising the distance between plasma surfaces and the PM geometry. The translucent half-circles indicate the location of the magnet with tightest toroidal curvature, which is used to determine toroidal grid spacing in the magnet geometry.

To choose a specific geometry one should consider the total magnetisation thickness and the size of individual magnets. Giving the optimiser a larger region to place magnets expands the solution space, at the cost of larger assembly forces and more magnetic material. Larger footprints for each magnet tower yield fewer total parts, at the cost of fewer degrees of freedom and more non-magnetic gaps in the PM holder. The PM geometry is radially discretised into identical layers by splitting each magnet tower into unit slices. For a fixed $(\zeta,\theta )$ grid, one can propose multiple radial discretisations based on the minimum radius $a_0$, the maximum radius $a_N$ and the number of slices $N_s$. From this the slice thickness is determined, $\Delta a = (a_N - a_0) / N_s$. It is desirable to bring the tower base $a_0$ as close to the plasma vessel as possible, but a limit is set by the accuracy of the VV's radial extent, the need to leave an assembly gap and the desire to leave space for radiofrequency (RF) antenna straps. The top of the tower $a_N$ is constrained by the inner bore of the TF coils. We considered several geometries using 12–20 layers of $\Delta a \approx 1.6$ mm (to be precise, we specified in customary units 1/16 inch = 1.5875 mm). We considered cases where the fixed thickness magnetisation layer was shifted closer and farther from the plasma within a region that spans 92–124 mm radially. Ultimately, we selected $a_0 = 100$ mm and $N_s = 14$. In the numeric optimiser, the effect of radial displacement can be achieved by setting a lower magnet to 0 and an upper magnet to 1. Figure 4 illustrates this and other possibilities for a PM tower. Since the minimum distance between PM and plasma surface is 15 mm, and the dipole thickness is 1.6 mm, the plasma gap is approximately 10 times the magnet slice thickness, which justifies the ideal dipole approximation.

Figure 4.

Five types of magnet towers. It can be completely filled (a) or partially filled. This gives the effect of magnitude variation. Partially filled magnets can be flush against the PM holder material (b), or they can be recessed with an air gap after assembly (c). The difference between (b) and (c) gives the effect of position variation within one PM tower. At rare locations (about 30 out of 10 000), the optimiser decided that it is best to leave a gap between magnet slices (d) and (e). Undoing the gap was detrimental to our plasma solution. The gap can separate magnets of the same sign (d) or opposite sign (e). This is shown with purple shading and diagonal dashed lines. At construction, these gaps were filled by annular non-magnetic plastic spacers.

In MUSE, the PMs are constrained to fit outside the VV ($r_{{\rm outer}} = 85$ mm) and inside the circular-bore coils ($r_{{\rm inner}} = 127$ mm). This leaves a radial thickness of 42 mm. We must also take into account the thickness of the PM holder material (see § 3), the imperfect shape of our glass torus and necessary assembly gaps. As a result, our final magnet thickness is approximately 22 mm. There is space left on each side, part of which we used to insert a thin helicon wave antenna between the VV and PM holder. An area for future improvement could be concatenating neighbouring footprints in regions of strong magnetisation (see figure 2). In addition, it would be desirable to define the poloidal rings to include an assembly gap for poloidal subsegments. We did not do this and had to draw an ad hoc jagged partition in the helical space between magnetisation bands instead (further discussed in § 3).

2.2 Continuous optimisation

With the geometry fixed, we run the first pass continuous magnet optimisation for a target plasma equilibrium. We use the FAMUS code (Zhu et al. Reference Zhu, Hammond, Brown, Gates, Zarnstorff, Corrigan, Sibilia and Feibush2020a) which takes the following inputs: our magnet grid, a target plasma boundary and a set of TF coils. FAMUS is the PM optimisation branch of the FOCUS code for flexible coil optimisation (Zhu et al. Reference Zhu, Hudson, Song and Wan2018). A quasi-Newton solver minimises the objective function

(2.2)\begin{equation} f_B = \frac{1}{2} \int_S \left[ (\boldsymbol{B}_{\rm PM} + \boldsymbol{B}_{\rm TF}) \boldsymbol{\cdot} \boldsymbol{n} \right]^2 \,{\rm d}A, \end{equation}

where $\boldsymbol {n}$ is a unit vector normal to the plasma surface. This is a surface integral of the residual normal field from PM and TF coils on the plasma boundary. Since MUSE is a vacuum equilibrium, the contribution from plasma currents is not considered.

The objective is to find arrangements of dipole magnitudes that best produce the target equilibria, subject to magnet constraints. Because our geometry has just one orientation, we define the dipole moment to be everywhere parallel to the normal vector of a circular-cross section torus. The signed magnetisation is represented by a normalised scalar $\rho \in [-1,1]$. The optimiser allows $\rho$ to take decimal values in the continuous pass, but in the final configuration only integer values $(-1, 0, 1)$ are selected.

The target metric (2.2) is summed over the plasma surface. If there are $N$ dipoles and $M$ plasma surface elements, this makes an $N \times M$ matrix. Because the orientation is fixed, we can greatly improve the speed of iterated calculations by precomputing an induction matrix

(2.3)\begin{equation} G_{ij} = \frac{\mu_0}{4{\rm \pi}} \left[ 3 \left( \frac{\boldsymbol{r}_{ij} \boldsymbol{\cdot} \boldsymbol{n}_i } {|\boldsymbol{r}_{ij}|^5} \right) \boldsymbol{r}_{ij} - \frac{ \boldsymbol{n}_i} {|\boldsymbol{r}_{ij}|^3} \right] \boldsymbol{\cdot} \boldsymbol{n}_j. \end{equation}

In this expression $i$ indexes plasma surface elements, while $j$ indexes magnetic dipoles. Thus, $\boldsymbol {n}_i$ is the normal unit vector on each of $N$ plasma surface elements, $\boldsymbol {n}_j$ is the normal unit vector from each of $M$ dipole moments and $\boldsymbol {r}_{ij} = \boldsymbol {r}_i - \boldsymbol {r}_j$ is the pairwise position vector from magnetic dipole to plasma surface elements. This is particularly useful for MUSE because both sets of unit vectors $(\boldsymbol {n}_i, \boldsymbol {n}_j)$ are fixed. The distances and dot products can be precomputed once for the entire optimisation. Only the set of scalars $\{ \rho _j \}$ are updated such that the local field $(B\cdot n)_i = G_{ij} \rho _j$. This is analogous to electromagnetic induction in the sense that the matrix is determined entirely by geometry, but there is no actual induction. Compared with the similar quantity defined in Zhu et al. (Reference Zhu, Hammond, Brown, Gates, Zarnstorff, Corrigan, Sibilia and Feibush2020a), this induction matrix is dotted into fixed dipole orientations. This simplifies computations in the ‘auto-zot’ algorithm described in Section 2.3.

The space of PM solutions is mathematically ill-posed. Very different input magnet arrays can give the same plasma field output. Figure 5 shows an example of this degeneracy starting from three initial conditions: all 0, all 1 and random. Each subplot shows a half period of magnet geometry in three dimensions. We see that three very different solutions all give ‘optimised’ plasma metrics while some are clearly suboptimal in their use of magnets. This motivates adding target metrics for the total volume (an L-1 norm) of PMs

(2.4)\begin{equation} f_V = V \sum |\rho| \end{equation}

and the elimination of intermediary values, such that only ‘binary’ dipoles $\rho = \{-1,0,1\}$ instead of the continuous distribution $(-1,1)$ remain

(2.5)\begin{equation} f_D = \sum |\rho| (1 - |\rho|). \end{equation}

We use these to find a Pareto frontier between metrics ($f_B, f_V, f_D$), by comparing differently weighted optimisations. Although the latter two metrics are not differentiable at $\rho = 0$, we did not find convergence issues arising from solutions bouncing at the cusp. Instead, the magnets converged to intermediary solutions between whole numbers to satisfy $f_B$. Figure 5 also shows the port locations where magnets are forbidden. There are three ports that straddle the $X$-symmetry plane vertically and three ports that straddle the $Y$-symmetric plane horizontally. A photograph of the physical ports is given in figure 19.

Figure 5.

Mathematical ill-posedness of PM optimisation. (a,d), (b,e), (c, f) An initial seed and a solution converged to low $f_B$. Very different solutions all give the similarly low error, therefore the problem is under-constrained. In plot (b), we have highlighted five regions that are excluded from the optimisation for diagnostic ports. Both the initial and final states fix the magnetisation in these regions to zero. When reflected through stellarator symmetry, these five regions become 3 vertical ports and 3 horizontal ports per period, for 12 total port locations (cf. full torus photo in figure 19). Note that zero is blue for the plots (a,b) where the colour bar ranges from $[0,1]$. In contrast, zero is green for the plots (c–f) where the colour bar ranges from $[-1,1]$.

To find solutions with precisely $-1, 0, 1$ we devised additional discrete refinement methods. Our optimisation used the Fortran code FAMUS. This algorithm has been adapted into the Python code FICUS and the optimisation framework SIMSOPT (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021). Making a continuous first pass is quick computationally. This stage can be used to scope the trade-off between different magnet array geometries, target stellarator equilibria, magnetic field accuracy and PM volume. The difficulty is taking discrete steps to set $f_D = 0$ while maintaining a tolerably low $f_B$.

2.3 Discrete refinement

Discrete refinement uses multiple iterations to force the normalised dipole moment to be exactly $-1$, 0 or 1 while preserving the target field optimisation. Achieving $f_D = 0$ is necessary for constructing a set of physical magnets because $|\rho | = 1$ makes all magnet slices identical. We found, however, that it could not be achieved by a continuous optimiser. This is illustrated by the Poincaré plots in figure 6. Even starting from a fairly low value of $f_D \sim 10^{-4}$, naïve jumps to 0 yield intolerable degradation of $f_B$.

Figure 6.

(a–c), (d–f), (g–i) A triplet of plots resulting from a specific set of relative weights ($w_B,w_V,w_D$). The legend reports the resulting magnetic field error $f_B$ $({\rm T}^2\,{\rm m}^2)$, magnetised volume $f_V$ (${\rm cm}^3$) and the binary dipole metric $f_D$ (dimensionless). While the metrics retain physical units, FAMUS applies relative weights after normalising the metric by their value at the first optimisation step. The histograms show that the distribution of $|\rho |$ tends to collect at 0,1 as $f_D \to 0$. The 3D plots show the spatial distribution of the dipole set relative to a target plasma surface, with dipoles $|\rho | < 10^{-3}$ suppressed. The colouring is green for 0, blue for $-1$ and red for +1 (see the legend in figure 5). The Poincaré plots trace field lines at the $\zeta = 0$ toroidal cross section. Lowering $f_D$ trades off against maintaining good flux surfaces, which is quantified by low $f_B$. Let us note that the figure shown is not the final MUSE equilibrium, but an alternative NCSX-based target plasma candidate.

We address this problem by defining two discrete jump operators: (1) a COMPRESS step that sums the magnitudes of all dipoles of each tower and redistributes them into a stack of integer magnetisations plus a single non-integer residual; (2) a ROUND step that forces all dipoles below a threshold $|\rho | < \rho _*$ to 0 while setting all others to either 1 or $-1$. Both make a discontinuous jump in solution space. Different choices of rounding threshold $\rho _*$, each forcing $f_D$ to be exactly 0, lead to different neighbourhoods of local minima. The continuous optimiser is rerun to find the next generation of optimised solutions. A detail of the continuous optimiser is that dipoles get ‘frozen’ when $\rho = 0$ because $\partial _\rho f_B \propto \rho$. Vanishing magnitude implies a vanishing gradient and no change in optimisation. To reset the continuous optimiser, a small residue of $10^{-4}$ was added to all zero cells after compression, with exceptions for locations where the optimiser is intentionally disabled (e.g. port locations and unused geometry layers). In such locations, $\rho =0$ is set as a feature to permanently disable magnetisation.

Iterating between cycles of compressing and rounding followed by continuous optimisation improves the minimum $f_B$ when $f_D = 0$ over generations. This is reminiscent of annealing metal over repeated heat cycles. For each generation, we sampled multiple values of rounding threshold $\rho _*$ prior to compressing. Figure 7 shows how the best result improves across generations. The final configuration achieves $f_B = 1.29 \times 10^{-8}$ $({\rm T}\,{\rm m})^2$, $f_V = 3.78$ litres total magnetised volume and $f_D = 0$. To reach this, one more type of discrete refinement was used, which we dubbed ‘zotting’.

Figure 7.

Generations of $f_B$ are shown to progressively decrease while maintaining low $f_V$, subject to the constraint that $f_D=0$ after each discrete refinement round.

A ‘zot’ is a discrete change where the magnitude of a single magnet slice $\rho _j$ is set by hand. The surface integral $f_B$ is recomputed as a proxy to evaluate the effect on the plasma surface. We assume that lower $f_B$ corresponds to more faithful reproductions of the target equilibrium. Zotting was used to ‘annihilate’ neighbouring pairs that have opposite signs. In many cases, we found that zotting pairs decreased the integral $f_B$. This surprising observation can be explained by the fact that ROUND is a majority vote procedure. We chose the threshold that makes most magnets optimal as a group, but there remain local cases where specific magnets are suboptimal as individuals. A ‘wyrm’ is a single tower that houses magnets of both signs (named in anticipation that oppositely signed magnets in close proximity will launch out, like a Norse dragon). Several wyrms were removed completely through annihilation. Where this was not possible, we were able to guarantee at least three magnet thicknesses between oppositely signed regions for all towers.Footnote ² Zotting was also used for enlarging the unmagnetised region around ports, cutting a path for the assembly joint between poloidal subsegments of the PM holder (see § 3), and annihilating pairs across neighbouring towers. The final MUSE configuration is called ‘zot-80’ because it is the 80th hand-picked permutation. A ‘gap’ is a wyrm without sign flip (see cases (d,e) in figure 4). There are 24 wyrms and 7 gaps per half period in the zot-80 configuration.

After the design of MUSE, an algorithm called ‘Auto-Zot’ was explored to systematically test every possible magnet after a compress and round. This is a once-through non-accumulating search, where each of the $M$ magnets is perturbed independently relative to a base case. The ensemble of beneficial changes is combined and accepted together. The set of independently beneficial changes does not equal the superposition of those changes. An additional annealing loop is added to relax into a self-consistent solution. This approach is partially inspired by the MASTER algorithm from Lu et al. (Reference Lu, Xu, Chen, Chen, Zhang, Ye and Wan2021b).

Figures 8 and 9 shows the results of optimisation at three constant toroidal angle cross sections. The final configuration uses 12 674 ideal dipoles per half period. This corresponds to 3.245 l of magnetic material in the full torus. The solution uses about 12.7 % of the total volume available to the optimiser. To simplify the design, we concatenate the 50 696 total dipole slices into 9736 magnet towers by merging contiguous slices. While the dipole slices are all one fundamental unit, the magnet towers are distributed into 14 discrete sizes that are integer multiples of the fundamental unit. By volume, size 14 is the most populated bin, containing about a third of the total magnetic volume. By count, size 2 is the most populated bin, containing about a quarter of the magnet pieces (see figure 17).

Figure 8.

Optimised magnet distribution unfolded as a contour plot. For each toroidal and poloidal position, the slices in a magnet tower are summed to a signed integer between $(-14,14)$. Orange vertical lines indicate the partition between PM holder quadrants. The vertical ports leave a visible ring of magnets at toroidal angles $\zeta = 0,{\rm \pi}$. Unlike figures 2 and 10, this poloidal axis sets $\theta ={\rm \pi}$ in the middle to emphasise the region of thickest magnetisation, which appears near the inboard mid-plane. The $\zeta = {\rm \pi}/4$ cross section from figure 9 falls on the orange dashed line, whereas the $\zeta = 0,{\rm \pi} /2$ symmetry planes fall in between.

Figure 9.

(a,d), (b,e), (c, f) Cross sections at constant toroidal angle, $\zeta$. In (a–c), red and blue dots show the signs of slices in a PM tower. Two ‘wyrms’ can be seen at the $\zeta = 0$ plane. The quiver plot depicts the poloidal projection of the vector field $\boldsymbol {B}_{\rm PM}$. The magnets are up–down symmetric at symmetry planes $\zeta = 0, {\rm \pi}/2$. Panels (d–f) show Poincaré plots with good nested flux surfaces. There is a 2/10 magnetic island near the edge, but it lays outside the VV as illustrated by the surface outlined in red. In both (a–c) and (d–f), the dashed black line indicates the inner boundary of the VV at minor radius $a = 7.6$ cm.

4.1 Finite mu

Magnetic permeability adds a nonlinear interaction between PMs. Ideal PMs are linearly independent sources. The magnetisation of real magnets, however, is modified by the field produced by their neighbours. While neodymium PMs are an almost ideal source of magnetisation, they do have a small relative permeability $\mu _{\|} =1.05$ in the direction along the magnetisation. In addition, they have been measured to have a perpendicular permeability of $\mu _\perp = 1.15$ (Katter Reference Katter2005). For optimisation, we assume $\mu = 1$ such that the field from all magnet sources is linearly independent. This section examines the effects of $\mu > 1$ for both isotropic and anisotropic permeability cases. We include only PM–PM interactions. The TF–PM interactions are assumed to be negligible because, at the PM surface positions, the TF field (0.1–0.2 T) is small compared with the PM–PM fringe fields ($\geq$1 T). Moreover, the perpendicular interaction between $B_{\rm TF}$ and $\mu _\perp$ generates a field aligned dipole moment with characteristic scale $\sim 2$ cm. The magnitude of this field is small compared with the magnitude of the main TF field in the plasma.

MagTense is an open-source magnetics framework for Python and Matlab that computes finite $\mu$ effects (Bjørk et al. Reference Bjørk, Poulsen, Nielsen and Insinga2021). MagTense discretises a volume using tiles and computes the analytic demagnetisation field for finite elements. It assumes that every tile has uniform magnetisation. Since commercial magnets have nominally uniform bulk properties, it is valid to treat finite elements as uniform for our study. The magnetic permeability can be varied in parallel and perpendicular directions. We start with an optimised set of point-dipole slices from FAMUS, expand them into finite-volume magnetisation towers in FICUS, and use MagTense to compute the magnetic field on a spatial grid. The spatial grid is stored in cylindrical coordinates as an MGRID file, which is fed into VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983) to model the free-boundary MHD equilibrium. We report the effects of finite $\mu _{\|}$ and $\mu _\perp$ from PM–PM interactions on the rotational transform profile and the neoclassical effective ripple profile which is computed from the STELLOPT code NEO Nemov et al. (Reference Nemov, Kasilov, Kernbichler and Heyn1999).

Let us start with the isotropic case $\mu _\perp = \mu _{\|}$. Increasing the relative permeability is found to decrease rotational transform. This can be interpreted as decreasing the poloidal field at a fixed TF. Figure 12 shows that this effect can be compensated by correcting the TF coil current to decrease the TF. A $\delta I / I < 2.5\,\%$ change in current (10.4 A change over a 433.0 A base current) returns the original target equilibrium. A similar effect was found for the effective ripple. Fluxes from neoclassical transport are proportional to the metric $\epsilon _{{\rm eff}}^{3/2}$ in the asymptotically collisionless $1/\nu$ regime (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999). This proxy decreases in response to magnetic permeability. While having a lower ripple is desirable for particle confinement, we interpret this effect as a decrease in overall shaping, which is undesirable. This is supported by the fact that the iota profile also decreases. The same TF coil adjustment that corrects the iota profile also corrects the epsilon effective profile to that of our base case. Modulating TF current to tune rotational transform was also used, for example, on W7X to correct for elastic deformation under electromagnetic coil loads (Lazerson et al. Reference Lazerson2019)

Figure 12.

(a,c) Effect of finite $\mu > 1$ on rotational transform $\iota$ profiles and the neoclassical transport metric $\epsilon _{{\rm eff}}^{3/2}$ profiles, using surfaces equally spaced in toroidal flux as the radial coordinate. (b,d) Both sets of profiles can be reduced to the $\mu =1$ case, provided a small decrement in TF coil current. The inset bounds the magnitude of this change $\delta I / I < 2.5\,\%$.

Anisotropic permeability effects are shown in figure 13, which varies $(\mu _{\|}, \mu _\perp )$ independently and plots the resulting rotational transform profiles. The effect of $\mu \ne 1$ is dominated by $\mu _{\|}$, and $\mu _\perp$ has a comparatively small effect. The dominance of parallel permeability $\mu _{\|}$ can be explained by figure 9. Since the radially oriented magnets exist to produce poloidal field, most of the field points parallel to the PM. There is not much field perpendicular to the PM. Thus, the effect of $\mu _\perp$ is subdominant. A different PM stellarator, one which allows magnets to be orientated orthogonal to each other, may reach another conclusion. These results demonstrate that $\mu > 1$ magnetic permeability effects do not pose a challenge for the MUSE PM stellarator. There is a shift in the $\iota$ profile, which can produce islands, but the effect can be compensated by adjusting the TF by $2.5\,\%$. This justifies using the $\mu = 1$ model for optimisation.

Figure 13.

Expansion and inclusion of the $(\mu _{\|}, \mu _\perp )$ combinations from figure 12. (a,b) We show 96 profiles, where $\mu _\perp$ ranges from 1.00 to 1.15, and $\mu _{\|}$ ranges from 1.00 to 1.05. Observe that the curves varying $\mu _\perp$ at fixed $\mu _{\|}$ are clustered together. (c,d) We compare Poincaré plots between the base case (blue) and the most perturbed case (red). Without correction, the axis is shifted, and some flux surfaces are removed from inside the VV volume. (e, f) The most perturbed case with correction (green) matches the base case.

4.2 Sensitivity

The fields in MUSE are designed assuming that every PM has a precise position, magnetisation magnitude and orientation. In reality, the magnetic material itself will have some finite tolerance magnetisation strength and orientation. There will also be imprecisions from assembly that affects position and orientation. In this section, we consider both stochastic and correlated sensitivities, and we set bounds on assembly tolerance.

The effect of uncorrelated perturbations is studied by generating $(N=600)$ magnet sets with stochastic perturbations in orientation, magnitude and position. Then we compute the surface error field metric $f_B$ for each. Figure 14 shows the distributions for four sets of perturbations. Each perturbation is sampled from a uniform distribution. The $\delta M$ varies magnitude by up to 3 %. The $\delta \boldsymbol {r}$ set perturbs the centre-of-mass position by up to 1 mm. The $\delta \boldsymbol {n}$ set varies orientation by up to $5^{\circ }$. Each of the above independently varies all dipoles in the system by a random perturbation drawn from a uniform distribution. A fourth set combines the above by perturbing every magnet by $\delta M$, $\delta \boldsymbol {r}$ and $\delta \boldsymbol {n}$. The threshold for good Poincaré plots and STELLOPT metrics was just above $\tilde f_B = 1.5$, a 50 % perturbation. Even the outliers of the combined perturbation lie below this threshold. The histograms show that these upper bounds on magnetisation, orientation and position perturbation magnitudes yield acceptable stellarator metrics. Therefore, we used this set of tolerances for magnet purchases.

Figure 14.

Histograms showing the relative change for the error field metric $\tilde f_B = f_{B,{\rm perturbed}} / f_{B,{\rm nominal}}$. Here $N=600$ sets of perturbed magnets were generated for each stochastic sample.

Laser metrology was performed on the magnet holders (see figure 11). The first round of PM holders found significant deviations up to 2.5 mm. These localised around the edge of the PM holder, and the deviations tended to shrink the printed PM holder relative to the CAD. The distortion is a cumulative error that peaks at the PM holder's edges. It appears to accumulate over the PM holder body and is significantly less over most of the volume. A series of correlated perturbation studies were done to determine what magnitude of distortion is tolerable. Stellarator flux surfaces were found to be relatively insensitive to distortions up to 10 mm. However, we found difficulty inserting test magnets into square cavities where the 3D print distortion exceeds 1.3 mm. Therefore, the PM holder tolerance is not set by the physics of magnetic sensitivity. Instead, the tolerance is set by a mechanical limit on magnet spacing. It was hypothesised that the shrinking distortion results from uneven cooling, as the first print was fan-cooled over a few hours. To address this, the subsequent magnet holders were inertially cooled over multiple days. Cross braces were also added within the torus to distribute thermal stress more evenly. The final 3D prints had bulk deviations less than 1.2 mm, while deviations less than 2.5 mm were tolerated in isolated pockets. For both physical and mechanical limits, we conclude that the final PM holders are satisfactory.

4.3 Magnetostatic forces

To evaluate forces in a PM stellarator, one must use the near-field interaction between neighbouring magnets to determine internal stress and the net force on the support structure. The dipole approximation is insufficient because it diverges as $1/r^3$ near the source. An exact analytic formula that is valid everywhere in space for the magnetic field from a rectangular prism magnet is presented in Appendix A. Given the magnetic field, forces on the PM can be evaluated using the method of magnetostatic charges. Magnetisation is a magnetic surface charge density, $\sigma _m = M$. One can discretise a charge element $\Delta q_m = \sigma _m \Delta A$, where $\Delta A$ is an area element such that the net force is $F = \sum B \Delta q_m$.

This force is analogous to the electrostatic Coulomb force $F = qE$. The magnetic monopole $q_m$ does not physically exist, but it is a valid computational tool in magnetostatics. The magnetic surface charge density, being equal to material magnetisation $M$, has units of A/m. Then the magnetic charge $q_m$ has units ${\rm A}\,{\rm m}$, which is consistent with the magnetic dipole moment $m = q_m d$ having units ${\rm A}\,{\rm m}^2$. Since all our magnet towers share a common footprint (40 mm$^2$), they all carry the same total magnetic charge, $q_m \approx 750$ ${\rm A}\,{\rm m}$ or, equivalently, ${\rm N}\,{\rm T}^{-1}$. Since the total charge scales like surface area $R^2$ (for photographic enlargement), we find that pressure and stress density is scale-independent for a PM stellarator.

Figure 15 shows the PM–PM forces in our system. As in the finite permeability study, this plot neglects the TF–PM forces because $B_{\rm TF}$ is an order of magnitude smaller than $B_{\rm PM}$ at the surface of the PMs. Let us recall the fact that, while the numerical magnet optimisation decomposed towers into dipole slices, physical construction used single magnets for all towers. We examine the force on each magnet individually by summing the $B$ field from all other magnets and applying the surface charge method. The resulting force peaks at approximately 20 N in the inboard region where magnetisation is thickest. The PM holders straddle the symmetry plane to minimise net force because individual force vectors cancel pairwise. This stores energy inside the PM holder as internal stress instead of external forces or torques. The magnetic pressure inside individual magnet towers (0.5 MPa from 1 T field intensity) is negligible compared with the tensile strength of the magnets ($\sim$80 MPa). The shear stress from TF field is even less because the TF field is of order 0.1–0.2 T. Therefore, our analysis is only concerned with the force and stresses on the PM holder, and not stresses on the magnets themselves. The remaining radial force is of order 1000 N (the weight of approximately 100 kg), which necessitates special assembly equipment. Section 5.3 describes how a hand-turned crankshaft is deployed to aid construction.

Figure 15.

Net force (colour bar) on each magnet tower ($\max \sim 20$ N). The legend sums over towers to compute the net total force on each PM holder quadrant (${\sim }1000$ N). There are approximately 2500 towers in each quadrant. Due to symmetry, only the radial component of the net force on each quadrant remains. The peak forces localise to the quadrant boundaries because internal forces cancel (they appear as stress instead) on the interior of magnetisation.

5.1 Magnet stencil

PM assembly is streamlined when magnets of a specific size and orientation are inserted as a batch. The stencil is designed to identify holes that contain the same magnet size and orientation while covering all others. Magnet size is distinct from cavity depth because some magnets are recessed into their cavities by the optimisation algorithm (see figure 4). We first prototyped 3D-printed plastic covers that conformed to the toroidal shape and later pivoted to 2D paper stencils that can be cut from ordinary card stock. The latter requires an algorithm to map the torus and its magnet holes onto flat 2D sheets, as illustrated in figure 16.

Figure 16.

The stencils are made by projecting a torus $(\zeta,\theta )$ onto a cylindrical segment $(u,v)$, then unrolling the cylinder into a plane $(x,y)$. (a) 3D perspective. The depicted angle ($\Delta \zeta = 40^{\circ }$) is exaggerated for illustration (the actual stencils used $15^{\circ }$). (b) 2D projection of the trapezoidal cross section. The cylinder is defined such that its intersection with the const $\zeta$ planes are circular. This sets $a_{{\rm cyl}} = a_{{\rm torus}} \cos \phi$ and shifts the centre of the cylinder to $R_{{\rm torus}} \cos \phi$, where $\phi = \Delta \zeta / 2$. As a result, outboard poloidal angles radially project ‘in’ from the torus to the cylinder, whereas inboard poloidal angles project ‘out’. (c) Curvature of the cylinder unrolled into planar coordinates. All four corners are right angles. The vertical extent covers only $180^{\circ }$ of the poloidal extent. The height of the unfolded map as a function of cylindrical–poloidal angle $v$ is an elliptic integral. (d) Photo showing how the stencil selectively reveals holes of a common type on the PM holder.

A torus cannot be unwrapped into flat paper because it has non-zero Riemann curvature, so we make a series of approximations. The first step is to partition each toroidal quadrant into straight cylinder wedges. Looking from above, these cylindrical wedges project into trapezoids in the $x$–$y$ plane. This removes curvature from the toroidal direction, and it introduces a shift in the poloidal angle coordinate. Next, we unroll each cylinder into a flat sheet, which maps straight lines into curved boundaries. It is helpful to define a chain of one-to-one maps

(5.1)\begin{equation} \left(\begin{matrix} \zeta \\ \theta \end{matrix} \right) \to \left(\begin{matrix} u \\ v \end{matrix} \right) \to \left(\begin{matrix} x \\ y \end{matrix} \right) \end{equation}

that carries points from toroidal coordinates $(\zeta,\theta )$, to coordinates on a cylinder $(u,v)$, to flat stencil coordinates $(x,y)$. Each step removes a degree of curvature. The cylinder is parametrised by a sector size $\Delta \zeta$. In our case, each stencil covers a sixth of a quadrant, so $\Delta \zeta = {\rm \pi}/12$. The $(\zeta,\theta )$ and $(u,v)$ coordinate systems share the same toroidal angle everywhere, but their poloidal coordinates coincide only at the ends of the cylinder. For this reason, it is convenient to centre the toroidal coordinate $u=0$ such that it bisects the cylindrical axis of each stencil segment. We define the poloidal coordinate $v=0$ on the outboard midplane where the major radius is largest. The first map is given by

(5.2)\begin{gather} u = \zeta, \end{gather}

(5.3)\begin{gather}v = \theta - \sin(\theta) \left[ \cos (\zeta) - \cos \left( \frac{\Delta \zeta}{2} \right) \right] \frac{ R + a \cos(\theta) }{a}, \end{gather}

where $R$ and $a$ are the major and minor radii of the torus. The second term in (5.3) corrects for the shift in poloidal angle between a toroidal wedge and its cylindrical approximation. It vanishes at either end where $\zeta = \pm \Delta \zeta /2$ because the mapping $v=\theta$ is exact on the boundary. The shift in angle is largest at $\theta = \pm {\rm \pi}/2$, and it is zero on the $x$–$y$ plane where $\theta = 0,{\rm \pi}$. The amplitude of this correction vanishes in the limit $\Delta \zeta \to 0$. Next, the cylinder is mapped to a flat 2D stencil by

(5.4)\begin{gather} x = [R + a \cos(v)] \cos \left( \frac{\Delta \zeta}{2} \right) \tan(u) , \end{gather}

(5.5)\begin{gather}y = a \cdot E\left[ v , \sin \left( \frac{\Delta \zeta }{2} \right) \right], \end{gather}

where $E$ is an incomplete elliptic integral of the first kind

(5.6)\begin{equation} E(\alpha, m) = \int_0^\alpha \sqrt{ 1 - m^2 \sin^2 x } \, {{\rm d} x} \end{equation}

that measures the height $y(v)$ for lines of constant $u$ (see figure 16c).

With this mapping, we project the boundary of the PM holders and the four edges of each square magnet cavity onto a 2D CAD image, which can be cut from light card stock. We used a computer-controlled paper cutter (the Cricut Maker 3), to automate the process. There are 29 magnet types (14 of each sign, plus 1 for holes that have both), and 4 unique magnet holder segments (2 for $X/Y$ times 2 for outboard/inboard). Therefore, there exist at most 116 sets of stencils. In practice, 81 stencil sets were used because not every option was populated. The total stencil area is further reduced by using stellarator symmetry to half the number of printed stencils.

5.2 PM assembly

MUSE uses 14 types of magnets that are composed from a single fundamental unit, as shown in figure 17. The polarity of these magnets can not be distinguished by eye. We define ‘north’ to be the direction that points parallel to a compass aligned with geomagnetic north. Stellarators and tokamaks have a net dipole moment. It is defined by the direction of net plasma current or the equivalent rotational transform. MUSE's net dipole moment points up out of the table. To match conventions, we arrange the PM such that this net dipole moment agrees with the ‘north’ of the individual PMs. Figure 9 shows the north field at the $\zeta = 0$ symmetry plane. Since stellarators are helical, north points down at the $\zeta = {\rm \pi}/2$ symmetry plane. Nonetheless, the net dipole moment north points up on average. With magnetic polarity fixed, we use reference magnets to orient the otherwise indistinguishable PM towers. A piece of magnetic material is used to hold and organise magnets during assembly.

Figure 17.

The distribution of magnet sizes in MUSE is peaked at size 2 and 14. The inset shows a photograph of MUSE magnet types, which are all composed from a single fundamental unit.

To pick up a magnet during assembly, we used a hollow cylinder guide magnet that slides along a non-magnetic wooden stick, as shown in figure 18. The guide magnet selects a polarity that remains fixed when the assembly tool twists and turns in hand. We used superglue (ethyl cyanoacrylate adhesive) to fix the magnets in the PM holder. Glue was preapplied in the square cavity. It takes 5–15 s for the glue to be dry enough to hold a magnet in the presence of background magnetic forces. The bond continues to strengthen over 24 h. Acetone was used to undo glue bonds where necessary. We note that this solvent does not affect the 3D-printed nylon material. Some magnets did not sit squarely into their cavity after insertion. We used a hammer and punch to level the magnet in place. The punches are made from non-magnetic brass. A 3/16-inch diameter is most useful for leveling and 1/16-inch diameter is useful for removing incorrectly placed magnets from behind using the drain hole. We used calipers to measure cavity depth to detect when a magnet was placed incorrectly. We also used calipers to verify when a magnet was positioned correctly in a recessed cavity (see case (c) in figure 4). Non-magnetic titanium pliers are useful for removing incorrect magnets after they have been partially dislodged by applying a thin punch through the drain hole. Dental tools are useful for removing debris such as dried glue from the corners of the square holes. They are also useful for removing residual powder from the powder-based 3D print. For the strongest regions of magnetisation, the cylindrical guide magnet was not sufficient to bring a target magnet into its cavity. The magnets tend to pivot and deflect towards already-glued neighbours. In these cases, we used the square hole of a titanium non-magnetic socket wrench to rigidly guide a magnet into its cavity.

Figure 18.

A cylindrical guide magnet is used to insert magnets. The guide magnet gives a distinct polarity to the otherwise indistinguishable magnet in hand. It can be detached while normal force is applied to secure the target magnet while glue dries. Also shown are the escape holes, which are used to eject misplaced magnets from the PM holder cavities.

As for the assembly procedure, we found that working in pairs to assemble each poloidal subsegment improves efficiency and decreases errors. Three common errors were identified: inserting a magnet into the wrong hole, inserting a magnet with the wrong polarity and having the glue dry while a partially inserted magnet is tilted. Prior to inserting a set of magnets, we used the stencil to mark all the corresponding holes and then removed the stencils during actual assembly. Since the magnets themselves have unmarked polarity, the guide cylinder is helpful for identifying orientation throughout each stencil set. It is strategic to start with the smallest magnets because these thin slices tend to flip in the presence of larger magnets. Next, we do all the ‘wyrms’ because they contain thin magnets that are susceptible to flipping. Then we go through the remaining magnets in order from smallest to largest. After inserting all the magnets in one stencil set, checks are conducted to verify that the set is indeed complete and that all magnets are oriented correctly and leveled.

Inevitably, some magnets were inserted incorrectly and needed to be removed. This occurred roughly 100 times or on 1 % of the total magnets. A magnet can be driven out by punching a thin rod through the drain hole on the back of the PM holder. When a magnet is removed, the hole is cleaned using acetone and dental tools before reinsertion. On rare occasions, it was necessary to apply acetone carefully to weaken the glue. Also on rare occasions, a rotary tool was required to clean residual glue from a cavity after the magnet was removed. This occurred fewer than 10 times. After all the PM holders were completed, we used an axial Hall probe to check for incorrect magnitudes and orientations. Tilted, displaced and significantly demagnetised magnets were detected by cross-checking the field strength at corresponding locations on identical pairs of PM holders (see figure 19). Approximately 10 magnets were found with the wrong polarity. The resulting error field has potentially large effects, so these magnets were removed and corrected. This check discovered some demagnetised size-1 (thinnest) magnets. Using a Hall probe we found that approximately 30 % of these smallest magnets were demagnetised by a magnitude of approximately 20 %, in spite of a formal quality assurance document from the vendor stating the spread in magnetisation is less than ${<}3\,\%$. A measurement of extra unused size-1 magnets found the perturbation occurring in similar proportions. This suggests that the ‘demagnetisation’ is actually a quality control issue from the vendor. We assessed the effect of this unexpected error using uncorrelated perturbation calculations similar to those described in § 4.2. Fortunately, it was found that for perturbations applied to the size-1 magnets, even a 50 % change in magnitude was inconsequential. A study of $N=600$ random draws showed that the nominal $f_B = 1.29 \times 10^{-8}$ was degraded to $1.45 \times 10^{-8}$ in the worst case, while the threshold for acceptance is $< 2.2 \times 10^{-8}$. Moreover, we found that for the 20 % magnitude perturbations on the size-1 magnets, the error field metric was even improved on average. The surprising observation could be explained by the fact discretising the magnet thickness makes size 1 the smallest non-zero unit. It would appear that the stellarator would gladly take ‘size 0.8’, even if randomly distributed. This fortuitous outcome retires the risk from these partially demagnetised magnets.

Figure 19.

(a) The 12 completed PM holder subsegments. Each column is a triplet that makes one quadrant. $Y$is on the left and $X$ is on the right. (b) The 16 planar, circular TF coils are positioned inside the water-jet cut support structure. (c) The glass VV is joined by 3D-printed low-thickness VV couplers. Glass ports were hot welded to the torus. The PM quadrants accommodate up to three ports each. Seven are installed at present on the north, west and south quadrants.

Opportunities to improve PM holder designs for future experiments include the following.

1. (i) Increasing the thickness of the wall behind the magnets (currently 1 mm). During insertion, we punctured a few cavities when tapping short magnets into deep holes. This was repaired by gluing the punctured fragment back onto the PM holder, and it did not leave any permanent deformation.

2. (ii) For parting the poloidal subsegments it is better to make cuts along the minor radius direction, rather than cuts that are horizontal or vertical along the ($R, Z$) direction because the former mates poloidal subsegments at $90^{\circ }$ faces. We unintentionally used radial cuts for the $Y$ quadrants but not for the $X$ quadrants. The slanting face redirects the compression force of the bolt joining two subsegments and introduces about 1 mm of sliding motion, even when the bolts are tightened completely. We plan to resolve this by using Teflon tape to increase friction, but the issue could be avoided altogether by using radial cuts to ensure a $90^{\circ }$ mating surface.

3. (iii) The $X$ quadrants have jagged toroidal cuts, whereas the $Y$ quadrants have a comparatively simple cut at constant poloidal angle. The latter is more favourable. It is an optimisation coincidence that we were able to find such a cut for $Y$ but not for $X$. It should be possible to predesign poloidal gaps in the dipole geometry since similar gaps are already included for the toroidal cut that separates $X$ from $Y$.

5.3 Torus assembly

MUSE torus assembly consists of three systems: the PM holders, a glass VV and the TF coils. A photograph of each system is given in figure 19.

The MUSE VV is made from four borosilicate glass quadrants. Each $90^{\circ }$ turn (7.5 cm minor radius, 30 cm major radius) was supplied by Schott-Kimax as off-the-shelf parts designed for use in chemical process plants. The advantage is low cost compared with custom metal fabrication. The transparent glass also allows direct visualisation of the stellarator plasma. Colleagues at the University of Wisconsin–Madison used a similar product to build a VV for the HELICOTOR experiment (O. Schmitz, private communication). Our glass quadrants did not immediately fit into a torus due to the bends being too square, the openings not forming exactly $90^{\circ }$, and the openings tilting out of plane such that two flanges would not meet at a flat surface. We bought 10 pieces and could not find a combination of 4 that fit into a torus. We hired professional glass blowers to align the flanges, smooth the bends and add ports. This process successfully produced a VV that has achieved $7 \times 10^{-8}$ Torr. The completed VV has seven ports: three vertical ports to match the $X$ pipe, three horizontal ports to match the $Y$ pipe and a single port pipe that can fit in either $X$ or $Y$ positions.

The TF coil system uses legacy ‘L-2’ coils (Bonanos Reference Bonanos1964), which were constructed for Francis Chen's linear experiments (Chen Reference Chen1995). They are 33-turn water-cooled copper coils wound in two layers, 16 turns in each layer with 1 turn crossing over. MUSE uses 16 coils. The resistance of each coil is about 23 $m\varOmega$. At nominal operation, the TF system produces a 1.5 kG magnetic field on axis, while drawing 70 kW of electrical power (433.0 A, 160 V). We purchased a pair of MagnaPower DC power supplies (500 A, 100 V).

To assemble the torus, we first installed three helicon antennae on two of the VV quadrants for ion heating. Then the three poloidal subsegments were fit around the VV quadrant. Four TF coils were mounted on each quadrant, and the four quadrants were brought together radially using an assembly jig. At the final stage of toroidal mating, conducting grounding bands were inserted inside the glass VV flanges, and VV clamps were installed outside the VV flanges. The VV clamps are custom made to have a low radial profile so that they can fit over the flange and below the PM holders. The torus support structure is water-jet cut from plates of G-10 fiberglass laminate. The TF coils butt against an hourglass-shaped central column. Figure 20 shows the assembly structure, and figure 21 shows the completed device.

Figure 20.

MUSE is assembled in four quadrants of VV, PM holders and TF coils. Only the VV and PM are shown in this photo. The assembly structure brings the four quadrants together using a hand-turned threaded shaft. This is necessary to overcome the mutual repulsion between PM holders.

Figure 21.

The MUSE PM stellarator is mounted on an optical table. Simple planar coils surround the circular-cross-section PM holder and VV. The torus support structure is made from water-jet cut G-10 plates. Lexan panels are installed between TF coils to add torque resistance. The TF power cables and water cooling manifold are above.

5.4 Cost and timeline

The physics design of MUSE's stellarator equilibrium and PM optimisation was completed in March 2021. The engineering Final Design Review was held in December 2021. Construction was completed in January 2023 with first plasma achieved in February 2023. The entire stellarator was built in approximately 1 year, including procurements. To the best of the authors’ knowledge, this is the fastest build for an optimised stellarator of any size.

In terms of material expense, the PM cost $16 5000 (2022 USD), which included a 30 % contingency for extra magnets. The VV cost approximately $5000 for raw material, but $10 000 for modifications by expert glass blowers because the original parts were out of round. This endeavor was published on the front cover of Fusion, a coincidentally named journal of the American Glassblowers Society (Paris & Souza Reference Paris and Souza2022) The TF coils were freely available from PPPL, but the warehouse inventory valued them at $1000 each. The 3D-printed PM holders cost about $2500 per quadrant. Therefore, the total material cost of the MUSE PM stellarator is just under $60 000. This sum includes only the material cost of stellarator components. It does not include the costs of employee salary, administrative overhead or auxiliary equipment such as power supplies.