Fully-staggered-array bulk Re-Ba-Cu-O short-period undulator: large-scale 3D electromagnetic modelling and design optimization using A-V and H-formulation methods

The development of a new hard x-ray beamline I-TOMCAT equipped with a 1-meter-long short-period bulk high-temperature superconductor undulator (BHTSU) has been scheduled for the upgrade of the Swiss Light Source (SLS 2.0) at the Paul Scherrer Institute (PSI). The very hard x-ray source generated by the BHTSU will increase the brilliance at the beamline by over one order of magnitude in comparison to other state-of-the-art undulator technologies and allow experiments to be carried out with photon energies in excess of 60 keV. One of the key challenges for designing a 1-meter-long (100 periods) BHTSU is the large-scale simulation of the magnetization currents inside 200 staggered-array bulk superconductors. A feasible approach to simplify the electromagnetic model is to retain five periods from both ends of the 1-meter-long BHTSU, reducing the number of degrees of freedom (DOFs) to the scale of millions. In this paper, the theory of the recently-proposed 2D A-V formulation-based backward computation method is extended to calculate the critical state magnetization currents in the ten-period staggered-array BHTSU in 3D. The simulation results of the magnetization currents and the associated undulator field along the electron beam axis are compared with the well-known 3D H-formulation and the highly efficient 3D H-{\phi} formulation method, all methods showing excellent agreement with each other as well as with experimental results. The mixed H-{\phi} formulation avoids computing the eddy currents in the air subdomain and is significantly faster than the full H-formulation method, but is slower in comparison to the A-V formulation-based backward computation. Finally, the fastest and the most efficient A-V formulation in ANSYS 2020R1 Academic is adopted to optimize the integrals of the undulator field along the electron beam axis by optimizing the sizes of the end bulks.


Introduction
In 2004, Tanaka et al. first proposed the concept of the bulk high-temperature superconductor undulator (BHTSU) in which a dipole field was utilized to magnetize in-situ a series of rectangular high-temperature superconductors (HTS) [1]. In 2008, Kinjo et al. proposed the concept of the staggered-array BHTSU by using a superconducting solenoid to magnetize a series of staggered-array superconducting half-moon-shaped disks [2], and in 2013, they demonstrated for the first time an undulator field B 0 of 0.85 T in a 10-mm period, 4-mm gap BHTSU prototype subjected to an external field change ΔB of 4 T [3]. In 2019, Calvi et al. demonstrated an undulator field B 0 of ~0.85 T in a 10-mm period, 6-mm gap BHTSU prototype subjected to an external field change ΔB of 7 T [4]. Theoretically, the value of B 0 could be doubled for a 4-mm gap. Based on numerical calculations, an undulator field B 0 above 2 T can be obtained at 10-mm period, 4-mm gap with an external field change ΔB of 10 T, assuming no quenches occur during the field-cooled (FC) magnetization process. In late 2020, the Paul Scherrer Institute (PSI) scheduled the development of a 1-m long BHTSU with a 10-mm period and 4-mm gap for installation in the new I-TOMCAT microscopy tomography beamline planned for the upgraded Swiss Light Source (SLS 2.0) [5]. The x-ray source generated by the BHTSU will increase the brilliance of the beamline by well over one order of magnitude in comparison to state-of-the-art cryogenic permanent magnet undulator (CPMU) technology [6]. The new I-TOMCAT beamline would allow experiments to be carried out for very hard x-rays with photon energies in excess of 60 keV, thereby extending the photon range by almost a factor of two with respect to comparable instruments in medium-energy synchrotrons [7].
One of the key challenges for designing a 1-m long (100 periods) staggered-array BHTSU is the large-scale 3D modelling of the magnetization currents inside the 200 bulk superconductors comprising it. Considering the central superconducting bulks can generate a uniform sinusoidal undulator field, a simplified way to optimally design the 1-m long BHTSU is to retain five periods from both ends. The problem then becomes the optimization of the field integrals along the electron beam. In the case of using the finite-element method (FEM), a fine meshing of the superconducting bulks and the surrounding air subdomain is necessary to obtain smooth and accurate computational results for the undulator field. This, unfortunately, results in millions of elements, and thus degrees of freedom (DOFs). State-of-the-art FEM models implementing the E-J power law to simulate the superconductor"s nonlinear resistivity can solve such complex 3D HTS problems, but can take a significant amount of computation time [8] [9][10] [11]. Other available numerical methods for critical state modelling are reviewed in several references [9] [12][13] [14], but in the following, the current state-of-the-art in numerical methods for modelling the critical state in bulk superconductors with respect to large-scale 3D HTS modelling is summarized.
In 1994, Bossavit first proposed the variational formulation method for the generalized Bean model. The electric field E was treated as a subdifferential of a critical energy density; the value of E was set to either zero if the current density |J| was lower than the critical current density J c or infinity if |J| was larger than J c . This method was further developed by Prigozhin [25]. This approach employed a number of iteration steps to find the maximum penetration elements, therefore, minimizing the magnetic energy in the superconductor. In 1997, Nagashima et al.
adopted a sand-pile model, which was proposed by Tamegai et al. in 1993, for the numerical calculation of the magnetic field distribution in Y-Ba-Cu-O bulks [26] [27]. This analytical approach assumes the magnetization current flows along the periphery, namely in square current loops in rectangular conductors and circular current loops in cylindrical conductors [28][29] [30]. In Witzeling proposed the circuit method to compute the screening currents inside a superconducting cylinder, assuming the superconductors as an array of parallel wires and creating a relation between different current loops based on the law of induction [34]. In fact, this was the very first numerical method which solved the critical state currents in superconducting bulks.
Other circuit model based methods were developed by Morandi in 2014 and Nugteren et al. in 2015 [35] [36]. In 2001, Coombs et al. proposed the field-screened method for fast computation of the critical state in type-II superconductors [37]. This approach forced the current density in the element with the maximum vector potential to the critical current density J c after every iteration step until the external field is screened from the interior. Further studies based on this method were explored by Ruiz-Alonso et al. [38] [39]. In 2012, Vesstgarden proposed the fast fourier transform (FFT) based approximation method for modelling the electrodynamics in superconducting thin films [40]. This method has then been utilized mainly for simulating the thermal instabilities and flux avalanches in superconducting thin films [41] [42]. In 2015, Shen et al. showed the superiority of the FFT-method over fourier transform and Biot-Savart approaches in calculating the magnetic field from an array of fully magnetized bulks [43], and in 2018, Prigozhin et al. further developed the FFT method and solved 3D magnetization problems related to bulk superconductors [44]. In 2005, Gu et al. proposed the resistivity-adaptive algorithm (RAA) method to calculate the AC loss in HTS tapes [45]. The key concept of RAA is to perform a number of iteration steps to find a resistivity matrix to fulfill the Bean model [49]. The penetration current was updated after each iteration step with a field and angle dependent J c (B, θ). In 2020, the same authors extended the iteration method to compute the magnetization of bulk superconductors with both the critical state model and the E-J power law, for both ascending and descending magnetization stages [50]. In 2020, Zhang et al. proposed a fast and efficient backward computation method to calculate the critical state currents in a 2D periodical bulk HTS undulator [51]. The key concept of the backward calculation is to relax inwards the induced large surface currents gradually, obeying Maxwell"s equations and the critical state model. The approximate critical state solution is simulated by setting a large enough n-value in the E-J power law. In 2007, Campbell proposed a new method based on the force-displacement curve of the flux lines to determine the critical state in superconductors [53]. The key concept was to define a flux flow resistivity such that the relevant power law is 1/n rather than n to speed up the computation and obtain stable solutions. Other similar approaches were later employed by Gömöry [51]. So far, the numerical analysis of magnetization currents in the BHTSU has been carried out with the 3D H-formulation [80] [81], the MEM method [82], the RAA-combined T-Ω formulation [45] [83], and the recently proposed backward computation method [51]. In particular, the backward computation method has shown excellent performance in simulating the magnetization currents in a periodical BHTSU model with 1.8 million DOFs within 1.4 hours [51]. This paper extends the A-V formulation-based backward computation method to compute the critical state currents in a ten-period BHTSU in 3D. Both the magnetization currents and the associated undulator field are validated against the well-known H-formulation and the mixed H-φ formulation method, as well as with experimental results. Finally, we use the fastest and most efficient A-V formulation-based backward computation method to optimize the bulk sizes to minimize the integrals of the undulator field along the beam-axis.

Theory of the A-V formulation-based backward computation and its extension to 3D modelling
Large eddy currents can be induced on the surface of a FC-magnetized superconducting bulk when assuming the magnetic flux pinning force is infinitely large. In such a case, the superconducting bulk acts as a permanent magnet, having a uniformly magnetized internal field equal to the initially applied value. However, in reality this situation can never be realized since the flux pinning force in a type-II superconductor is always limited to a finite value. Assuming one superconducting bulk is FC-magnetized slowly under isothermal conditions, the eddy currents will gradually penetrate inwards inside the bulk following a quasi-static critical state model. In the end, the bulk superconductor is magnetized as much as possible with the minimum electro-magnetic entropy production [14], but without generating higher magnetic field than the initially applied [37]. This is indeed the key theory of the backward computation method recently proposed for critical state modelling of type-II superconductors [51]. During the backward iterations, the large surface currents induced by rapid FC magnetization relax inwards step-by-step obeying Maxwell"s equations and the critical state model with isothermal assumptions, as shown in the schematic in Figure 1. The relaxation stops when no more superconductor elements can be penetrated. Any additional penetration trying to further minimize the electro-magnetic entropy production will result in the field-screened phenomenon [37]. J refers to the magnetization current density, i refers to the element number, and P refers to the penetration sign.
The algorithm of the backward computation method can be implemented in any programmable FEM software with different forms of Maxwell"s equations. Here we employ the popular FEM software ANSYS using its default A-V formulation method, as described by Eqs. (1)-(4).

∇ × A = B
(1) superconductor elements is ~10.20 T, larger than the initial 10 T, which is impossible as stated by the field-screened method [37].
To conclude, we can draw conclusion that no more superconductor elements can be further penetrated and the simulation result obtained by the A-V formulation-based backward computation method agrees with the MEMEP method.
Studies on the critical state in a superconducting cubic bulk were recently carried out by Pardo et al. using a MEMEP 3D variational method and compared with the 3D H-formulation method [14][84] [85]. An obvious off-plane (x-y) bending effect for the magnetization currents is observed where far from the mid-plane, since the high induction field at the diagonal is compensated by a J z component which has opposite signs at the diagonal. However, the assumption of square current loops (i.e., the sand-pile model) can still provide a good approximation of the magnetic moment inside a cuboid bulk HTS [86]. In fact, the current bending effects disappear and the current follows in square loops when the applied magnetic field is above the full penetration field. In this section, we will extend the A-V formulation-based backward computation method to calculate the critical state magnetization currents in a 3D cuboid bulk HTS based on the sand-pile model assumption. .
Taking the i th circuit in Figure 3 as an example, the resistivity is set to infinity in both the yand z-directions and zero in x-direction for elements whose center coordinates are "y > x"; the resistivity is set to infinity in both the xand z-directions and zero in the y-direction for elements whose center coordinates are "y < x"; the resistivity is set to infinity in the z-direction and zero in both the xand y-directions for the corner element whose center coordinate is "y = x". The critical current density J c,i in the i th circuit is updated after every backward calculation and its value equals the averaged J c (B), as expressed by Eq. (8). The averaged current densities of J x,i ̅̅̅̅ (y>x) and J y,i ̅̅̅̅ (y<x) are calculated for the top elements (M i ) and the right elements (N i ), respectively, according to Eqs. (9)- (10).
As shown in the schematic of the numerical algorithm in Figure 4, when J x,i ̅̅̅̅ (y>x) is lower than −J c,i and J y,i ̅̅̅̅ (y<x) is greater than J c,i , the penetration sign P i of the i th circuit is set to 1 and the absolute value of the magnetization current density in the i th circuit is forced to J c,i ; when the value of the penetration sign P i is 1, the magnetization current density in the i th circuit is updated with the latest J c,i . The current direction in the corner element is assumed to be 45 degrees against the xor y-direction.
The relaxation stops when no more current loops can be penetrated. Assuming the HTS bulk is FC magnetized at 10 K, its critical current density J c (B) is expressed by the following equation where the values of J c1 , J c2 , B L , B max and y are 1.0x10 10 A/m 2 , 8.8x10 9 A/m 2 , 0.8 T, 4.2 T and 0.8, respectively [51].

3D H-formulation-based methods
Two 3D H-formulation-based methods are also used to simulate the BHTSU problem: the well-known, traditional Ampere"s (Eq. (2)) and Faraday"s laws. µ is assumed to be the permeability of free space, µ 0 . The E-J power law [87] is used to simulate the nonlinear resistivity of the superconductor, where E is proportional to J n and n = 100 is assumed to reasonably approximate the critical state. The E-J power law takes into account the assumption made for J c , which may be assumed constant or J c (B) as described above (Eq. (11)). FC magnetization is simulated as described in [51], by setting an appropriate magnetic field boundary condition on the outer boundaries of the air subdomain such that, for 0 ≤ t ≤ 100 s, µ 0 H z (t) = 10 -t/t ramp , where t ramp = 100 s. Thus, the initial condition is µ 0 H z (t = 0 s) = 10 T and the magnetic field is ramped linearly down at a rate of 0.1 T/s to µ 0 H z (t = 100 s) = 0 T. In the 3D H-φ formulation [78], the bulk superconducting subdomains are still modelled using the H-formulation, but in the nonconducting, air subdomain, the model solves the magnetic scalar potential φ. This significantly reduces the number of DOFs, since it is a scalar rather than a 3D vector in the H-formulation case. In addition, no artificial resistivity is needed in the air subdomain (usually set to 1 Ω·m or similar in the H-formulation). In the φ formulation, Ampere"s law states that when neglecting displacement currents, and the magnetic scalar potential is defined as In addition, the divergence free condition of B results in the governing equation in the air subdomain One must also take care to ensure appropriate coupling at the interface between the two formulations: the tangential components of the magnetic field are equated on the H-formulation side and the perpendicular components are equated on the φ side as described in [88]. FC magnetization in this case is simulated in the same way as the H-formulation for the applied field, and the initial condition is the same for the bulk superconducting subdomains, but in the air subdomain, an initial condition of φ = -10· z/µ 0 (from Eq. (13)) is set.

Large-scale 3D modelling of the BHTSU critical state currents
As described in [2][4], a sinusoidal undulator field B y along the electron beam axis (z-axis) can be generated in staggered-array superconducting bulks after FC magnetization with a superconducting solenoid. In the models, we approximate this by applying a uniform background magnetizing field. The past prototypes and FEM simulations were all based on staggered-array half-moon-shaped bulk superconductors acting as a compact insertion device [3][4][80] [81]. In this section, we study the undulator field generated by staggered-array cuboid bulks, to which the pre-stress can be applied more homogeneously because of their regular shape. Figure 5  travel along the z-axis they are forced to oscillate in the xz-plane and generate hard x-rays under the action of the Lorentz force associated with the sinusoidal field B y . As the magnetization currents in each cuboid bulk are symmetric based on the assumption of square loops, the ANSYS A-V formulation model of the ten-period BHTSU can be simplified as shown in Figure   5(d), with applied flux parallel boundaries at x = 0 and y = ±9 mm. During FC magnetization from 10 T to zero, the temperature inside the ten-period BHTSU is kept constant at 10 K.
The A-V formulation-based backward computation method is implemented in ANSYS 2020R1 Academic. To obtain accurate simulation results, the entire FEM model shown in Figure 5 The H-and H-φ models were built in COMSOL 5.6: the H-formulation is implemented using the Magnetic Field Formulation (mfh) interface and the φ-formulation is implemented using the Magnetic Fields, No Current (mfnc) interface. In addition to the general settings outlined in Section 2.2, we make use of the geometric symmetry of the problem to model 1/8 th of the problem, as described in [89]. Thus, ¼ of each of the bulks is modelled around the ab-plane and ½ of the total undulator is modelled along its length, mirrored across the xy-plane along the z-axis at z = 0. In the H-formulation parts, this is achieved using the "Magnetic Insulation" node (n x E = 0) for the yz-plane at x = 0 and xz-planes at y = ±9 mm, and using the "Perfect Magnetic Conductor" node (n x H = 0) for the xy-plane at z = 0, respectively. In the φ-formulation, the "Magnetic Insulation" node is also used, but in this case, n· B = 0. For the half-symmetry along the undulator length, the "Zero Magnetic Scalar Potential" node is utilized such that φ = 0 for the xy-plane at z = 0. A mapped mesh of element size 0.25 mm x 0.25 mm x 0.25 mm is used in the region 0 ≤ x ≤ 10 mm, -9 mm ≤ y ≤ 9 mm, 0 < z < 150 mm and the rest of the air subdomain is meshed as "Free Tetrahedral" (using COMSOL"s predefined "finer" mesh setting). The model has a total of 0.13 million elements in the HTS bulks and 1.22 million elements in the air subdomain. All elements are linear (first-order) edge (curl) elements. There are then approximately 3.2 million and 1.3 million DOFs for the H-and H-φ formulations, respectively.

Comparison with experimental results
A five-period Gd-Ba-Cu-O bulk undulator with a period length of 10 mm and the magnetic gap of 6 mm was tested at the University of Cambridge [4]. As a first attempt to verify the BHTSU concept, FC magnetization experiments were conducted with a conservative background solenoid field of 6 T, rather than the simulated 10 T. Figure 9 shows

Optimal design of the BHTSU
In synchrotron radiation light sources or x-ray free electron lasers, the electron beams are expected to follow their original direction of motion after experiencing the field of an undulator, without neither an offset nor an angle with respect to the orbit defined by the accelerator magnets. It is therefore essential to minimize the first and second integrals of the undulator field B y along the z-axis according to Eqs. (15)- (16). Figure 9. Comparison between the measured undulator field B y (the central period) after FC magnetization from 6 T at 10 K [4], the simulation results in 2D [51] and the simulation results in 3D obtained by using the A-V formulation-based backward computation method.
The measurement B y is taken from the central period of the Gd-Ba-Cu-O bulk undulator prototype.
The influence of the bulk sizes on the undulator field B y was investigated by running a series of 3D BHTSU models. As compared in Figure 10, the amplitude of B y increases as the length of the bulks" sides increases from 10 mm to 14 mm. Further increasing the bulk length does not appreciably increase the undulator field B y . Therefore, the size of the main superconducting cuboid bulks is set to "14 mm x 14 mm x 4 mm". The integral values (IB y and IIB y ) of the undulator field B y are minimized by optimizing the sizes of the end bulks in the ten-period BHTSU and the associated magnetic gaps. First attempts are made by merely optimizing the sizes of the two end bulks on the top, but a smooth transition of the undulator field B y cannot be obtained since two field peaks at the ends always exist. Further studies indicate that slightly reducing the sizes of the two end bulks at the bottom can eliminate the two field peaks at the ends. Finally, the size of the two end bulks at the bottom is set to "13 mm x 13  magnetization current density J s in the staggered-array HTS bulks and the induced undulator field B y along the z-axis. It should be pointed out that the magnetization currents in the HTS bulks are not fully symmetric with respect to the xy-plane because of some small numerical error induced by the iterative analyses. The numerical errors can be corrected by retaining the undulator field B y in the left-half model and extending the point-symmetric B y values based on the fact that a symmetric undulator geometry can always generate a point-symmetric undulator field B y , having a zero value of the first field integral IB y (−∞, +∞) along the z-axis, as shown in Figure 12(e). Thus the main task is to minimize the second field integral IIB y (−∞, +∞) as much as possible to reduce the offset of the electron beams passing through the HTS undulator. Figure 12(f) plots the integral values of IIB y (z) along the z-axis for the given bulk sizes shown in Figures 12(a)-(b). The optimized IIB y (−∞, +∞) is ~1.2 x 10 -6 Tm 2 , satisfying the requirements for an insertion device commissioned in the synchrotron radiation light source.
It is worth mentioning that the optimization of the undulator field B y is based on each of the HTS bulks having the same J c characteristics (as given by Eq. (11)). The material properties often differ between different HTS bulks, as demonstrated in [4], even within the same batch. In such a case, to improve the accuracy of the models even further, we need to first evaluate the J c of each HTS bulk by referring to the measured undulator field B y and then tune the amplitude of B y and the associated phase errors through moving the HTS bulks up and down [92].

Conclusion
In this work, the theory of the recently-proposed 2D A-V formulation-based backward computation method is extended to calculate the critical state magnetization currents in a ten-period staggered-array BHTSU in 3D. The numerical algorithm is fastest and most efficient A-V formulation is then adopted to obtain the optimal design of the ten-period BHTSU. After more than twenty iterations of the sizes of the end HTS bulks, the second field integral of the ten-period BHTSU is minimized to 1.2 x 10 -6 Tm 2 . The numerical algorithm of the 3D backward computation is fast, reasonably straightforward and adaptable to other programmable FEM software packages, allowing it to be readily extended to simulate other large-scale HTS magnetization problems in the future.