Advanced Pipeline for Designing Multi-Locus TMS Coils With Current Density Constraints

Objective: This work aims for a method to design manufacturable windings for transcranial magnetic stimulation (TMS) coils with fine control over the induced electric field (E-field) distributions. Such TMS coils are required for multi-locus TMS (mTMS). Methods: We introduce a new mTMS coil design workflow with increased flexibility in target E-field definition and faster computations compared to our previous method. We also incorporate custom current density and E-field fidelity constraints to ensure that the target E-fields are accurately reproduced with feasible winding densities in the resulting coil designs. We validated the method by designing, manufacturing, and characterizing a 2-coil mTMS transducer for focal rat brain stimulation. Results: Applying the constraints reduced the computed maximum surface current densities from 15.4 and 6.6 kA/mm to the target value 4.7 kA/mm, yielding winding paths suitable for a 1.5-mm-diameter wire with 7-kA maximum currents while still replicating the target E-fields with the predefined 2.8% maximum error in the FOV. The optimization time was reduced by two thirds compared to our previous method. Conclusion: The developed method allowed us to design a manufacturable, focal 2-coil mTMS transducer for rat TMS impossible to attain with our previous design workflow. Significance: The presented workflow enables considerably faster design and manufacturing of previously unattainable mTMS transducers with increased control over the induced E-field distribution and winding density, opening new possibilities for brain research and clinical TMS.

Unfeasible winding patterns are a critical problem when designing mTMS coils for electronic control of the stimulus location and orientation [3]. The electronic control in our mTMS method is achieved by adjusting the amplitude of simultaneously triggered pulses in a set of overlapping coils placed on the scalp. Each coil induces a specific E-field distribution on the cortical surface, enabling a fast and accurate control of the stimulation focus within a designated area without physically repositioning the transducer. Our method requires that each coil produces a predetermined E-field distribution that is often more complex than in conventional TMS. Consequently, the coil windings in an mTMS transducer can be complicated, as it is challenging to find energy-efficient yet manufacturable winding patterns that reproduce the desired E-field distributions. Further, as the stimulating E-field in mTMS is produced as a linear combination of the E-fields of the different coils, just one of the coils failing to accurately produce its desired E-field distribution would result in inferior performance unless redundancy is built into the system.
A powerful method for designing coils is to optimize for a minimum-energy surface current density that generates a desired E-field [8], [12], [16]. In a recent study, an efficient TMS coil for rat TMS was constructed by searching for the coil windings that require the least amount of energy for inducing the desired Efield, without constraining other properties of the E-field than the location, intensity, and direction of its maximum [17]. In a recent study on TMS coils for mice, manufacturability was considered in the optimization by including the maximum current density in the objective function [13]. Minimization of the maximum and total current density also has been previously applied in shim and gradient coil design for magnetic resonance imaging (MRI) [18], [19], [20], [21], [22], [23]. A similar technique has also been employed for stellarator coil design in fusion As |K| increases, more wires are required to match the total current. If |K| is too high, there will not be enough space for the required number of wires, making the winding pattern physically impossible.
reactors by parametrizing the engineering feasibility of the coil and including it in the objective function in optimization [24]. These methods, however, do not allow direct control of winding density, which may result in suboptimal designs in applications with tight requirements for both winding space and energy efficiency, such as mTMS.
In this paper, we present an mTMS design workflow that provides fine control over the winding density and the induced E-field distributions. To guide the optimization towards feasible windings without compromising energy efficiency, we implement a constraint that restricts the maximum surface current density in accordance with critical spatial and physical limits. We verify this approach by comparing mTMS transducer optimization with and without the surface current density constraint. We also manufactured a transducer to demonstrate the feasibility of our design method.

A. Surface Current Density Constraint
A TMS coil with a desired E-field distribution can be designed with a two-stage process [8]. First, we find the optimal current density K (A/m) on a surface that covers the maximum possible extent of the TMS coil. We then discretize this surface current density into a coil winding. The optimal K must produce the desired induced E-field spatial distribution [3], [8], [12] whilst minimizing, e.g., the energy, resistive power loss, or maximum current density [16]. Other design constraints, such as maximum coil current and minimum wire spacing [6], can be included at this stage. The TMS coil can be virtually represented by a flat or curved surface triangle mesh with each flat triangular element containing a constant surface current density [12]. To ensure a physically valid current distribution, the optimized surface current distribution is expressed as a scalar stream function Ψ [9], [12], [13]; a coil winding path template is obtained by discretizing Ψ into its isolines. This method for TMS coil design is described in detail in [12].
To obtain manufacturable TMS coil winding patterns, we need to account for winding density during the computational design. Discretizing K into a winding [6], [10], [14] yields a coil that reasonably approximates K if the wires are distributed such that the total surface current in a region is matched by the total current in the wires within the region. The number of winding turns required in an area depends on the required surface current density and the current in the wire, whereas the maximum possible winding density depends on the wire diameter. Thus, with a wire of a given diameter and maximum current, it is impossible to accurately reproduce surface current patterns in regions where the total required surface current exceeds the current-carrying capacity of the wires, as illustrated in Fig. 1. To prevent optimization results where |K| is too high to be realizable, a constraint on K is applied: |K| ≤ I max /d wire everywhere, where I max is the maximum current that can be driven into the coil and d wire is the wire diameter. This approach to enforcing coil feasibility has been previously explored in the context of MRI coils with minimum power dissipation [25].
Let us represent K with a piecewise linear scalar stream function Ψ in a triangle mesh consisting of flat elements: where r is the position andn the unit normal vector of the mesh [8]. To express the surface current density constraint in terms of Ψ, we write where the threshold value k = I max /d wire .
As Ψ is piecewise linear, in mesh triangle j we have where Ψ j,i is the stream function value at node i of triangle j, e j,i the vector between the other two nodes, and G j the area of the triangle ( Fig. 2(a)). Thus, to satisfy (2), we implement constraints to ensure |K| ≤ k in every triangle. To do so, we apply the approach used in Ref. [12] for constraining E-field magnitude on the cortex: we constrain |K| in each triangle with 16 linear constraints, which together closely approximate the nonlinear constraint |K| ≤ k ( Fig. 2(b)). The density constraint can thus be expressed as Ax ≤ b, where x is a vector containing the values of Ψ at the N nodes mesh nodes, b is a vector of 16 times N tri (number of constrained triangles in mesh) elements, each entry being equal to k, and A is a 16 N tri × N nodes matrix formed so that Ax gives K along the constraint directions for each triangle. Each row of A thus corresponds to a single linear constraint, such as the one highlighted in blue in Fig. 2(b), by preventing the magnitude of K (as defined by Ψ) in a particular triangle, in one of the 16 directions, from exceeding the threshold value k.
The TMS coil optimization problem is quadratic and can be solved with the interior-point method [8]. To accurately obtain the K that minimizes the cost function, a high-resolution mesh is required. To apply the surface current density constraint on each triangle of such a mesh would require a very high number of linear constraints, of the order of 60,000. Each constraint adds a dimension to the logarithmic 'barrier function' present in the interior-point method. Based on our tests, common solvers cannot reliably handle this problem due to various numerical reasons, such as numerical instability, even when most constraints end up irrelevant as |K| k in most triangles. Instead of immediately constraining every triangle, we iteratively add triangles to A and b. After an initial round of optimization without the density constraint, successive rounds are performed where constraints are added to all triangles where |K| > k in the previous run. This process is continued until the condition |K| ≤ k is met everywhere. Typically, most triangles are left unconstrained.

B. Equality-at-FOV Constraint
The designed mTMS coils should accurately produce the desired (target) E-fields. To ensure that the E-field induced by a coil matches the target E-field, we constrain the maximum difference between the E-fields. As the E-field far from the targeted cortical region is typically not important, we define a stimulation field of view (FOV), a region within which the E-field of the optimized coil should closely match the target E-field. We refer to the associated constraint set as the equality-at-FOV constraint. The stimulation FOV can be manually determined or automatically defined from the desired stimulation E-fields.
The equality at FOV is defined by a pair of linear constraints for each E-field degree of freedom per cortex point within the FOV. Two pairs of constrains are enough in a spherical model where only the E-field components tangential to the surface can be nonzero ( Fig. 3(a)). At each point, one constraint pair limits the magnitude of the induced E-field in the direction of the target E-field (E target ) in that point: |E || − E target | ≤ αE max , where α is a constraint parameter, E || the magnitude of the induced E-field component parallel to E target , and E max the maximum magnitude of the target E-field within the FOV. The other constraint pair restricts the magnitude of the component perpendicular to E target (E ⊥ ): |E ⊥ | ≤ αE max . Thus, the equality-at-FOV constraint limits deviations from the target E-field to √ 2αE max at maximum ( Fig. 3(b)).
The constraint can be expressed in matrix form Ax = b, where x is an N nodes -element vector of stream function values, A is a 4N FOV × N nodes matrix where N FOV is the number of points in the FOV, constructed such that Ax yields induced E-field component magnitudes E || , −E || , E ⊥ and −E ⊥ in each FOV point, and b is a 4N FOV -element vector with corresponding magnitudes E target + αE max , −E target + αE max , αE max and αE max , respectively, for each FOV point.

C. Optimization Workflow
Multi-locus TMS transducers are designed to generate a range of focal stimulation E-fields via linear combinations of a small number of coil-induced E-fields. In a previously presented optimization workflow [3], transducers were designed by first determining a set of target stimulation E-fields, optimizing current patterns to reproduce them, and then decomposing those current patterns to get suitable transducer coils. In this study, we present a refined workflow (Fig. 4) that provides finer control over the target stimulation E-fields while ensuring that the coil winding densities and required currents remain within desired limits.
In our previous workflow [3], we computed the set of target stimulation E-fields by using a commercial coil model in different locations and orientations near the cortex. Here, we replace the model with a surface current density pattern that induces an E-field distribution with a focal region (the area where ) of custom length and width [26], giving us more control over the stimulation E-field shape. To find a suitable current pattern, we optimize for the minimum-energy surface Fig. 4. The optimization workflow with density-and equality-at-FOV constraints for finding manufacturable mTMS coils. Only a single, downsampled coil mesh is shown for visual clarity; in reality, each coil has its own high-resolution mesh at a different distance from the head. current distribution [8] that induces an E-field distribution on the cortical surface with specific magnitude and orientation in five control points subject to the equality-at-FOV constraint. The E-field directly below the center of the coil mesh must be of the magnitude E max , pointing in any arbitrary direction, and the E-field in four points around it must have magnitude E max / √ 2 and the same direction as in the peak point. Two of the points are at distance a (width) along the spherical surface from the maximum in directions perpendicular to its E-field direction, and the other two at distance b (length) in directions parallel to the E-field, as shown in Fig. 5. The set of desired stimulation E-fields is then found by moving the optimized surface current pattern into different positions and orientations and determining the resulting E-field distributions.
Next, we determine the target E-fields of the transducer coils from the stimulation E-fields. We first express the stimulation E-fields as a set of vectors, for each E-field, by defining two components per cortex point corresponding to the two tangential directions illustrated in Fig. 3(a). We collect the E-fields into a matrix and apply singular value decomposition (SVD) to obtain singular vectors that can approximate the stimulation E-fields via linear combination. Using the E-fields described by these singular components as targets, we optimize for a set of minimum-energy surface current distributions that generate these transducer E-fields, applying the surface current density and equality-at-FOV constraints. The resulting current distributions are discretized into sets of stream function isolines, providing winding templates for the transducer coils.
As constraining the maximum current density restricts the maximum E-field output of a coil, high target E-field magnitudes The five points defining the focal region and the resulting E-field, induced by a minimum-energy current pattern that reproduces specific E-field vectors in each of the five points. The direction is the same in each point, but the points indicated by the red and purple circles have magnitude 1/ √ 2 times that of the maximum point (black circle). The parameters a and b customize the shape of the E-field, allowing us to select an energy-efficient E-field with a shape that suits our needs as the target stimulation E-field.
can result in the optimization problem having no feasible solution. To alleviate this problem, we estimate the required maximum output for each coil by determining the linear combinations of transducer E-fields that generate each target stimulation Efield. By defining the target magnitudes as the maximum E-field among those linear combinations, we avoid over-constraining the problem by demanding coils to generate stronger E-fields than operating the transducer requires.

D. Validation
To validate the workflow, we designed a 2-coil mTMS transducer that can translate the peak E-field over a 1-cm line segment (E-field translation line) on a spherical rat cortex model of radius 13.7 mm. Due to the significantly smaller size of the rat brain compared to the TMS coil, our previous workflow [3] would result in a coil that is impossible to manufacture, making this a suitable application to validate the method developed in this study. We defined 11 stimulation E-fields evenly spaced by 1 mm to guarantee a smooth shift along the 1-cm translation line segment. The focal region width a was set at 4.8 mm, slightly below half the length of the translation line, so that the focal regions at the extremes of the line segment did not overlap. Length b was set at 7.8 mm for the E-field proportions to resemble those of a typical figure-of-eight coil [27]. In comparison, the commercial coil MC-B35 Butterfly (MagVenture A/S, Denmark) is slightly less focal, with a = 5.4 mm and b = 9.8 mm.
An SVD was performed on the set of 11 E-fields, of which the first two singular vectors, explaining 99% of the total variance, were used for the coil optimization. To approximate the target E-fields with E max = 100 V/m, the first and second components required 100 and 65.5 V/m, respectively, determined by leastsquares fitting the transducer E-fields to the stimulation E-fields within the FOV. Importantly, the optimizer would fail to find a feasible solution if the second coil were required to output 100 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. rather than 65.5 V/m. The linear rise time of the current was set at 60 μs. The coils were modeled as two rectangular meshes of 10 cm × 20 cm, the longer side along the E-field translation line. The meshes were placed 4.6 and 6.6 mm away from the cortical surface to account for coil-layer stacking. Despite the second SVD component E-field requiring the lower peak magnitude, reproducing it required the higher-intensity current pattern, thus the mesh closer to the cortex was used with that component as the target E-field. Notably, this order is the inverse of the similar transducer presented in [3].
The maximum wire current I max was set at 7 kA and wire diameter d wire at 1.5 mm, resulting in surface current density threshold 4.67 kA/mm for the density constraint. For the equality-at-FOV constraint, the tolerance parameter α was set to 0.02 corresponding to a maximum error of 2.8%. We applied an additional exact equality constraint at E-field maxima to ensure that the peak strength of the induced E-field matches the target. The FOV was defined as the set of points where at least one of the 11 stimulation E-fields had a magnitude equal to or above 0.8E max . Using the full focal region (E > √ 2E max ) as the FOV prevented the optimizer from finding feasible solutions when the density constraint was applied. Therefore, we reduced the FOV size rather than increasing α to prioritize E-field fidelity near the E-field translation line. The chosen FOV was sufficient for the optimized coils to accurately reproduce the primary features of the desired transducer E-fields.
Two optimization rounds were performed for both coils. The first round used only the equality-at-FOV and the second one both the equality-at-FOV and the density constraint. To discretize the surface current distributions into manufacturable winding paths that best approximate them, the numbers of isolines were chosen as the largest numbers that yielded current paths at least d wire apart from each other. All computations were performed with custom-made scripts written in MATLAB 2021a (MathWorks Inc, USA), and optimizations used MATLAB's built-in implementation of the interior point method algorithm. The workstation used for the computations had an Intel Xeon W-2133 CPU with 32 GB RAM and ran on Windows 10.
Flat coil formers with 1.5-mm-wide grooves following the isolines of the optimized current patterns were 3D-modeled with SolidWorks 2021 (Dassault Systèmes SE, France) and manufactured with a Formlabs Form 3 (Formlabs, USA) 3D printer using Formlabs Tough 1500 Resin. In the 3D model design, the isolines were joined into continuous winding paths by manually adding grooves connecting them. The coils were wound with a single layer of 1.3-mm-diameter litz wire (Rupalit Safety with 3-layer Mylar insulation, Rudolf Pack GmbH, Germany) and potted with epoxy glue (Bison 2-Component Universal Epoxy). This wire was chosen due to the unavailability of suitable 1.5mm-diameter wire. The individual E-fields of the coils, as well as an example mTMS-shifted E-field induced by pulsing both coils together, were measured on a 70-mm-radius hemisphere with our TMS characterizer [27]. The transducer was placed on a platform 15 mm away from the E-field sensor of the characterizer; the center of the bottom coil winding was approximately 16.5 mm away, and that of the top coil approximately 18.5 mm away. The E-fields of the manufactured coils were characterized by delivering trapezoidal pulses [28] with our custom power electronics at 200 V and measuring the average E-fields induced by the rising current. Pulse rise time was 30 μs, hold period 10 μs and fall time 15.7 μs (top coil) and 14.8 μs (bottom coil). The example mTMS-shifted E-field was generated with voltages of 173 and 207 V for the top and bottom coils, respectively, a combination that keeps peak E-field magnitude equal to that of the 200-V pulse in the top coil.
The E-fields of the optimized current patterns were calculated at the measurement points to compare the computed surface current density patterns and the manufactured coils. The 200-V pulses ramped the coil current in 30 μs from 0 to 1.25 kA in the top coil and to 1.05 kA in the bottom coil, whereas in the optimization, the pulse duration was 60 μs and the peak current 7 kA. As the induced E-field is proportional to the rate of change of the coil current, the calculated E-field of the top coil was scaled down by 0.36 and that of the bottom coil by 0.3, to make the measurements and calculations comparable. We measured each E-field at a set of 1000 points centered around the middle of the coil. We also measured the E-field profiles along the E-field translation line (121 points; 7.3-cm line segment).

III. RESULTS
The target E-fields and the stimulation FOV are depicted in Fig. 6. Optimizing for the two transducer E-fields with the density and equality-at-FOV constraints took 20 minutes. With the previous workflow, it would have taken 45 minutes longer due to optimizing current patterns for the 11 stimulation fields.
The optimized current patterns and the resulting windings are shown in Fig. 7. Using only the equality-at-FOV constraint, the maximum surface current densities were 6.59 and 15.4 kA/mm for the top coil (target E-field Fig. 6(b) and bottom coil (target E-field Fig. 6(c), respectively, both exceeding the density constraint threshold 4.67 kA/mm. With optimal winding patterns of wire width 1.5 mm, the coil currents required to reproduce these patterns would have been 9.88 and 23.2 kA, respectively. With the density constraint applied, both current density patterns were more spread out with maximum magnitudes equal to the threshold value, reducing the coil current to 7 kA for both patterns.
In both optimization runs, the target E-fields were replicated with a maximum difference in FOV equal to √ 2 · 0.02 · 100 V/m Without density constraint, the windings become extremely tightly packed in regions of high current density. As the isolines form a set of closed paths, the final design must combine the isolated loops into a single winding. This is generally done by manually adding connecting paths in a 3D modeling software.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. ≈ 2.8 V/m. The shapes of the E-fields of all coils can be seen in Fig. 8. With the top coil, the density constraint had minimal effect on the E-field shape. With the bottom coil, the optimization without density constraint resulted in E-field magnitudes remaining higher in a wide area surrounding the FOV compared to the optimization with density constraint. Without density constraint, the highest feasible numbers of isolines were 10 and 3 for the top and bottom coils, respectively, whereas with the constraint, 16 and 12 isolines could be used. Using 16 and 12 isolines, without density constraint the minimum distances between adjacent lines were 0.9 and 0.38 mm for the top and bottom coils, respectively, whereas with density constraint the distances were 1.55 and 1.50 mm.
Two coils were manufactured based on the isoline patterns shown in Fig. 7(f), with the coil plates ( Fig. 9) assembled into an mTMS transducer. The measured E-fields of both coils closely match the calculated E-fields (Fig. 10), with only minor discrepancy: In the top coil, the field is attenuated faster along the negative x direction; in the bottom coil, the focal point at positive x is slightly stronger than its negative counterpart. This discrepancy propagates to the E-field combining stimulation with both coils (Fig. 10(c)). The focalities of the measured and simulated E-fields were similar in all three measurements. Fig. 9. The two 3D-printed coils based on the isolines shown in Fig. 7(f). A single layer of wire was wound into the grooves, with the complicated design of the bottom coil requiring a detour through the top coil. The winding paths were found by determining grooves for each isoline and then manually inserting connecting paths joining them into a single winding. The connections were placed such that current in each loop flows in the direction dictated by the surface current pattern.

IV. DISCUSSION
In this study, we present an mTMS coil design workflow with four major complementary advancements that enable the design of manufacturable mTMS coils with a high degree of control over E-field shapes and coil-winding density. Our method introduces 1) an improved surface current density constraint that accounts for desired wire width and maximum current, enabling the design of previously unattainable, manufacturable coils. 2) A method for defining custom target E-fields that enforces specific dimensions for their focal regions, as well as a constraint set that ensures the induced E-fields match the optimization targets more accurately than attainable with the method in [8]. 3) Coil-specific E-field strength estimation that ensures feasible combined use of the individual coils in mTMS. 4) A refined application of SVD that significantly reduces the computational time of coil design compared to a previously published pipeline [3].
We verified our optimization algorithm by designing and manufacturing a compact 2-coil mTMS transducer that can electronically shift the stimulus location along a 1-cm line segment on a spherical rat cortex. We showed that the method with the density constraint can result in a manufacturable design where earlier methods would otherwise result in coils that require either impossibly tight windings or high currents that cannot be achieved by conventional TMS power electronics. Using density and equality-at-FOV constraints together resulted in current patterns that reproduced the desired target E-fields within the FOV as accurately as the optimization with only the equality-at-FOV constraint. 3D-printed coils based on the optimized current patterns successfully replicated their E-fields at the measurement distance, with minor mismatch in the E-field distributions for both coils. As both measured E-fields were stronger on the right side, this mismatch can in part be due to millimeter-scale alignment errors between the coil and the measurement probe.
As the coil characterizer in our measurement setup had been designed to model stimulation of the human cortex [27], the measurement points were roughly 1-2 cm further away from the coils than the focal region on the spherical cortex model used in the optimization. The measurement points were also sampled from a spherical surface with a radius considerably larger than that of the rat cortex model. Thus, while the E-fields resemble their calculated counterparts at measurement distance, the measurements do not directly demonstrate how the transducer performs at the distance it is intended to be used at. However, it can be expected that the E-fields similarly match elsewhere in space as well, though high-spatial-frequency features that have already been attenuated at measurement distance could differ at closer ranges. Demonstrating that the E-fields are accurately reproduced at the rat cortex distance would require constructing a new sensor setup tailored for rat TMS characterization.

A. Equality-at-FOV and Current Density Constraints
By applying the density constraint, we get the tightest feasible windings that produce the desired E-fields without the need for cumbersome fine-tuning of optimization parameters. Sparser windings resulting from reducing the maximum current density increases the inductance and resistance of the coil [18], [19]; in [18], optimizing for minimax current density yielded a surface current pattern with inductance and resistance over twice those of the patterns found by minimizing those quantities. By limiting the maximum current density by the smallest degree necessary to ensure manufacturable windings, we avoid trading energy efficiency for winding sparsity where addition of winding space is not of significant benefit. This can be especially beneficial in applications demanding high-power pulses and complex coil designs, such as mTMS.
Though it would likely be possible to attain the tightest feasible winding with the minimization approach as used, e.g., in the rat coil studies by Cobos Sánchez et al. [13], [16], doing so would require precise balancing of weight parameters in the loss function, and any change in the optimization setup would likely require rebalancing the weights. An algorithm for balancing the weights to attain a desired winding density has been proposed by Harris et al. [29]. However, this algorithm does not account for the required current magnitude in the final coil design, only the winding distance. The current patterns obtained with our approach meet the design requirements without a need for cumbersome fine-tuning, making it an efficient method when there is a known maximum current density that should not be exceeded.
The density constraint operates under the assumption that wires will be only d wire away from each other in regions of peak surface current density. If the wire distance significantly exceeds that value, the coil may require currents higher than I max to operate. Although determining windings from evenly spaced stream function isolines does not guarantee the distance being exactly d wire , we can consistently get results close to it by selecting the largest feasible number of isolines. In this study, the minimum winding distances for both coils were within common machining tolerance of d wire , meaning the assumption was satisfied to the degree that can be practically achieved in manufacturing.
Although feasible maximum currents are successfully enforced by our method, it has some notable limitations as well. Isolines with sharp turns and small loops can be impossible to reproduce accurately in physical windings, and the density constraint does not prevent such features. Additionally, applying the equality-at-FOV constraint can result in an unsolvable optimization problem when the tolerance threshold is strict and FOV is wide. This problem could be alleviated in the future by varying the error threshold depending on the distance from key regions. This way, we could have more control over the peripheral E-field without sacrificing fidelity in critical areas.

B. Improved Coil Optimization Workflow
The introduced workflow advances mTMS transducer design methodology in several ways. First, the custom target E-field definition allows direct control over the desired E-field shape, removing reliance on existing commercial coil models [26]. As the required coil energy increases with stimulation focality [8], [26], choosing a suitable target E-field shape is important. If necessary, the E-field distribution could be controlled even more precisely by adding more control points or other constraints such as restricting the maximum magnitude outside the focal region. The presented target E-field definition approach can also be used for designing transducers capable of previously unattainable stimulation, such as TMS where the spatial shape of the stimulation E-field changes between pulses. Such transducers could be used, for example, in studying how stimulation focality affects TMS efficacy, a poorly understood topic with significant importance to coil design [26].
Estimating the required E-field strength for each coil separately as opposed to requiring each coil to match the desired stimulation strength allows us to decide better the order of the coils in the transducer. Generally, the coils with shapes leading to the least energy-efficient E-field production are placed closest to the scalp, but if those coils need comparatively weak output strengths, it can be optimal to place them further away. Furthermore, knowing how strong a field each coil needs to be able to produce can be of vital importance when applying the surface current density constraint. Without estimation of the coil-specific output strength, the optimizer found no solution when applying the density constraint in our validation design problem, whereas using our method resulted in a feasible transducer.
Another advancement of our workflow is the considerably reduced computation time thanks to applying SVD directly to the stimulation E-fields rather than at a set of optimized coils. This approach enables rapid design iteration even with modest computational resources. Our workflow also removes the need for optimizing current patterns to reproduce the stimulation E-fields, which saved approximately two-thirds of the total optimization time in our validation test. The time saving becomes more relevant for mTMS transducers with more than two coils. For such transducers, the ratio of the number of target stimulation E-fields to the number of transducer E-fields is much higher-when optimizing a 5-coil transducer, Nieminen et al. utilized 8964 target E-fields to cover a 30-mm-diameter region on the cortex [4]; with the same computational setup as in this study, optimizing current patterns for such targets would have taken over 600 hours, whereas performing the equivalent task with the method proposed in the present study would take less than an hour. For transducers with more than five coils, the difference is even greater. Thus, removing the need to optimize current patterns for target E-fields allows designing complex mTMS transducers within a reasonable time without requiring the use of a high-performance computer cluster. The time savings can be especially large if this optimization step would need to be performed multiple times, for example, to find the optimal balance between the focality of the target E-field and the energy efficiency of the transducer.
Determining the transducer E-fields directly from the SVD of the stimulation E-fields also makes for a closer representation of the stimulation E-fields compared to those generated with our previous method. In some applications, imperfect representation of the desired stimulation E-fields may result in the stimulation range of the design falling short of the required one.

V. CONCLUSION
We demonstrated that constraining the maximum surface current density in mTMS coil optimization can result in manufacturable coil designs where a previously introduced approach would yield unfeasible results. We introduced a design workflow that incorporates a current density constraint in addition to improvements resulting in faster optimization and increased E-field control compared to previous methods. Coil design involving high current densities is a critical part of many technologies ranging from TMS to MRI to nuclear fusion; the constraint implementation demonstrated in this study provides a useful tool for finding optimal designs in various complex applications.

DECLARATION OF COMPETING INTEREST
J. O. Nieminen, L. M. Koponen, and R. J. Ilmoniemi are inventors on patents and patent applications on TMS technology, including technology relevant to mTMS discussed in this work.