Analytical Model of 3-D Helical Solenoids for Efficient Computation of Dynamic EM Fields, Complex Inductance, and Radiation Resistance

In this article, we use the <italic>dynamic</italic> Green's function to produce a frequency-dependent magnetic vector potential <inline-formula><tex-math notation="LaTeX">$\vec{A}(\omega)$</tex-math></inline-formula> and derive expressions for the efficient (accurate and fast) computation of cylindrical components of the magnetic flux density vector <inline-formula><tex-math notation="LaTeX">$\vec{B}(\omega)$</tex-math></inline-formula> as a function of the solenoid's geometric and material parameters. <inline-formula><tex-math notation="LaTeX">$\vec{A}(\omega)$</tex-math></inline-formula> may be used to efficiently compute the frequency-dependent flux linkage <inline-formula><tex-math notation="LaTeX">$\Phi (\omega)$</tex-math></inline-formula>, the complex inductance <inline-formula><tex-math notation="LaTeX">$L(\omega)$</tex-math></inline-formula>, and the radiation patterns of the solenoid anywhere in space including both near-field and far-field regions, excluding the (source) regions of conducting wire. In addition, we propose the complex calibration coefficient <inline-formula><tex-math notation="LaTeX">$\chi (\omega)$</tex-math></inline-formula> to account for the finite-radius conductor. Several numerical examples are provided to validate the proposed helical model against the superposition of circular loops. The proposed model is demonstrated for a wide range of applications across the spectrum from 60 Hz to 170 GHz, representing low-frequency power systems to high-frequency mm-wave communication systems. A plan is being developed for experimental validation of the model.


I. INTRODUCTION
T HE ever increasing diversity and complexity of electrical systems that utilize a broad range of frequencies in the electromagnetic (EM) spectrum demands computationally efficient (accurate and fast) models and methods that facilitate system analysis, design, and optimization while ensuring electromagnetic compatibility (EMC).
The solenoid has a long history of usage in numerous and diverse electrical applications across the EM spectrum; e.g., 1) in the low-frequency and magneto quasi-static (MQS) regime, solenoids operating from tens of Hz to kHz may find application in power systems [1], [2], biomedical sensors [3], and data/power transfer for power electronics [4]; 2) solenoids operating in the tens of MHz to tens of GHz may find application in communications and remote sensing as spiral antennas [2], [5], [6], and in high-speed electronics as inductors and current injection loops [7], [8]; The author is with the Department of Electrical and Computer Engineering, University of Idaho, Moscow, ID 83844 USA (e-mail: azadehgol@uidaho.edu).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TEMC.2023.3298886.
Digital Object Identifier 10.1109/TEMC.2023.3298886 3) in the millimeter (mm) and sub-mm wave regime, solenoids operating near 95 GHz may find application in sheet-beam amplifiers and Gyrotrons, which may serve in electron cyclotron resonance heating, and magnetic resonance spectroscopy, etc. [9], [10].Therefore, an analytical model of the canonical solenoid that enables the efficient computation of dynamic EM fields/waves across the spectrum, in 3-D space, is indispensable.
Previous analytical models may lack in one or more ways [1]; for example: 1) some models may approximate the quasi-static mode using the static Green's function [1]; 2) some models may be valid only for particular configurations, such as infinitely long solenoids [11], tightly-wound wire winding with an infinitesimal wire pitch [11], [12], circularly symmetric loops [13], or a specific core material [14]; 3) other models be valid only in certain regions, such as only inside [14] or only outside the solenoid's volume, only along the solenoid's axes, and only in the far-field or only in the near-field regions, etc.On the other hand, computational models based on spatial discretization (e.g., finite element [15] or finite difference methods [2]) can offer a powerful tool in solving the EM fields of structures comprised of arbitrary geometry and material, however, they may face certain challenges in efficiently solving the fields of solenoids, which are often disparately multiscaled.
For example, Pollack et al. [16] represented the solenoidal magnetic fields via a series solution to Laplace's equation in helical coordinates using the separation of variables technique.As explained in [16,Sec.IV], the method requires determination of unknown separation constants of the series-expanded scalar potential through least-square fitting of the magnetic field B on a subset of 3-D spatial grid points inside the solenoid, where B is calculated through numerical evaluation of the Biot-Savart integral in discretized 3-D space.As described in [16, Sec.V], the computational cost (in memory and compute-time) of solving the disparately scaled solenoid's fields through 3-D spatial discretization is exorbitant even when parallelized across many processors.
Budnik and Machczynski [17] solved the solenoidal magnetic fields via computation of the Biot-Savart integral in the rectangular coordinate system, where the helical contour is discretized as straight line segments in 3-D space, this technique looks to be similar to the field-generation step described in [16, Sec.IV].Tominaka et al. [11] performed multipole expansion of a single helical current conductor to express the magnetic field of a helical dipole with an infinite length.Similar to [16], this model too relies on the separation-of-variables technique, where the unknown separation constants seem to be determined by application of periodic boundary conditions for a coil of infinite length, while their formulation presumes an infinitesimally thin helical current filament.
As will be shown in our work, although a model based on superposition of symmetric circular loops (M CL ) may reasonably approximate some of the field components of a helical structure with small winding pitch, the M CL is unable to quantify field variations due to helical asymmetry.In addition, one of the three orthogonal field components is nonexistent in the M CL , thus, it may be inadequate for applications requiring the aforementioned attributes.Moreover, the helical model (M H ) is computationally much more efficient compared to the M CL , as the computational cost of M H is nearly of order O(log(N T )), whereas that of M CL is nearly of order O(N T ).This is a significant advantage compared to many 3-D finite methods (e.g., [2] and [18]) that are approximately of the order O(N 3 E ), where N E is the number of discretized cells along each spatial axis and dictated by the solenoid's geometry, material, and excitation wavelength.Solenoids with multiple winding layers may be easily modeled through superposition of parametric integrals of each layer.As will be shown in this article, our proposed model seamlessly illustrates the solenoid's transitions across the reactive (energy storage) and the radiating modes, while identifying the resonant modes.
In contrast to [11] and [16], our proposed method: 1) Avoids the need to solve any separation-constants; 2) Avoids solving the Biot-Savart integral in 3-D discretized space, which is computationally expensive as discussed in [16] and [17]; 3) Enables efficient modeling of disparately-scaled cylindrical solenoids with a core of finite height and finite radius, wound by circular wires comprised of conductors and insulators of finite radii; 4) Provides direct access to the position of wires, thus enabling straightforward spatial control of the current; 5) Proposes a complex frequency-dependent calibration coefficient to conveniently account for the wire's finite radius.

A. Background and Basic Definitions
In [1], we presented an analytical model for estimating the 3-D MQS fields of cylindrical solenoids with helical wire winding, using the static Green's function and relying on the Schur vector product of the electric current source vector and its tangent to obtain a frequency-independent magnetic vector potential A. Herein, we use the dynamic Green's function to obtain a frequency-dependent A(ω).Maxwell's differential equations in time-harmonic form (with e jωt dependence) express the dynamics of EM fields/waves in media, as follows [19]: where the angular frequency ω = 2πf (rad/s), f (Hz) is the cyclic frequency, the imaginary number j = √ −1, E (V/m) is the electric field intensity vector, B (Wb/m 2 ) is the magnetic flux density vector, H (A/m) is the magnetic field intensity, D (C/m 2 ) is the electric flux density, J c (A/m 2 ) is the electric conduction current density, J i is the electric impressed (source) current density, and M i (V/m 2 ) is the impressed (source) magnetic current density.
Considering the cylindrical symmetry of the structure in Fig. 1 [1], [2], the cylindrical coordinate system is a reasonable choice for expressing both the EM fields/waves and the geometry.In the Cartesian (rectangular) coordinate system, the vector ν = {ν x x, ν y ŷ, ν z ẑ}.In the cylindrical coordinate system, the vector ν = {ν ρ ρ, ν φ φ, ν z ẑ}.The unit-vector ν = ν ν , where ν = ν = √ ν • ν * is the magnitude of the vector ν according to the two-norm for real vectors [20] and for complex vectors [21], where x * is the complex conjugate of x, and p • q designates the dot product of p and q [19].The rectangular and cylindrical coordinate systems are related via the transformation rules in [19].

current source I [19]
where μ = μ 0 μ r core for 0 ≤ ρ ≤ r c ∧ 0 ≤ z ≤ h c and μ = μ 0 for otherwise, and μ r core is the relative permeability of the core material inside the solenoid.The primed symbols designate the source coordinates, the unprimed symbols designate the observation field coordinates, d is the magnitude of the differential length vector over the contour C of the line source, and the Green's function where R (m) is the magnitude of the distance between source and observation positions, and the phase constant β = ω √ μ (rad/m).The magnetic flux density vector B (Wb/m) may be computed by taking the curl of A, as follows [19]: We may express the magnetic flux linkage Φ(z) as a function of position z along the solenoid's height [19], [13] where ds is the differential area vector normal to the surface S, and d is the differential length vector around the contour C that encloses the surface S. We may obtain the total flux Φ solenoid by integrating Φ(z) over the solenoid's height along z.We may define the inductance L (H) as the ratio of Φ solenoid over the source current I, as follows [19]: In the rest of Section II, we develop the analytical expressions for computing in 3-D space the frequency-dependent magnetic flux density in (5), magnetic flux linkage in (6), and complex inductance for the cylindrical solenoid with helical wire winding.

B. Parametric Helix
Referring to the Fig. 1, we define the solenoid's geometry by the core radius r c (m), core diameter d c = 2r c , wire conductor radius r w (m), wire conductor diameter d w = 2r w , wire insulation radius r i (m), wire insulation diameter d i = 2r i , and the winding layer radius r L = r c + r i .We designate the number of winding turns-per-layer along the core's height (along ẑ) by N T ∈ N * (natural numbers without zero: +1, +2, +3, . . .), the core's height h c = N T × d i (m), and we designate the number of winding layers in the direction of the core's radius (along ρ) by N L ∈ N * .Furthermore, we assert the constraint: The helical wire represents the contour of the electrical current flowing around the cylindrical core of the solenoid.We express the path of the contour in the rectangular coordinate system using (8), which is in parametric form [22] C where the parameter τ spans the range 0 ≤ τ ≤ 2π × N T .Our choice for the z-component in (8) means that the outer boundaries of the insulation region of adjacent wires touch each other as they wind around the helical solenoid.Using the coordinate transform relations in [19], the vector C R (τ ) in rectangular coordinates may be transformed to the cylindrical coordinates, to yield We may obtain the tangent vector T [22] by taking the derivative of ( 9) with respect to (w.r.t.) τ , in the following: and magnitude of the tangent vector T is We define the unit-vector û in the direction of the tangent vector, as

C. Distance Vector, R: From Source to Observation
In the cylindrical coordinate system, we define the field observation position P = {ρ ρ, φ φ, z ẑ}, and the source position P = {ρ ρ, φ φ, z ẑ}.Furthermore, we define two position vectors ν and ν w.r.t. the coordinate system's origin, where ν and ν , respectively, describe the position of P and P .The distance vector R = ν − ν , and its magnitude R = R may be expressed as Transforming the position of the point P located on the helical contour in (8) from rectangular to cylindrical coordinates [19], we obtain the position vector R, as follows: We evaluate R on the contour of the helix by substituting the coordinates of ( 14) for P into (13) to yield and substitute ( 15) into (4), to obtain G( ν, ν ).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

D. Differential Arc Length Vector, da
We use (9) to define the differential arc length vector da [22] of the helical contour, as then, by relying on (10) we rewrite (16) in terms of the parameter τ and the tangent vector T , in the following: and using (17) together with (11), we find the differential arc length vector's magnitude da , as

E. Differential Length Vector, d
We may obtain the tangent vector v at the position P by taking the derivative of ( 14) w.r.t.τ , to yield and we obtain the unit-vector v in direction of the tangent where We define the differential vector dR for the helical contour path of (14) in the cylindrical coordinate system, as We define the differential length vector d [23] for the helical contour of ( 14), as and insert the differential elements of ( 21) into (22) to yield where Λ = {ρ, ρ φ, ẑ}, and p • q designates the Schur vector product of the vectors p and q.Using ( 14), we substitute ρ → r L in (23), and obtain the magnitude of d in (24), which is equivalent to da in ( 18)

F. Source Current Vector, I
A nonzero excitation frequency may produce a nonuniform current distribution in the finite-conductivity and finite-radius wire, due to the skin-effect [24], the proximity effect [25], etc.However, to limit the complexity of formulation in this work, we will assume a uniform current distribution across the wire cross-section, then use (20) to yield where I 0 is a constant designating the input current's peak amplitude.The magnitude of components of I c are plotted in [1, Fig. 2] for I 0 = 1.0.We now define the current in the helical contour in terms of the parameter τ , to obtain III. REDEFINING (3), USING THE SCHUR PRODUCT Note the integral in (3) assumes that I and d are functions of the source coordinate variables {ρ , φ , z }.However, I in (26) has been defined as a parametric vector function of the parameter τ , thus, considering the equivalence between d in (24) and da in (18), we propose to modify the integral in (3) by using I from (26), da from (17) or d from (23), and û from ( 12) to rewrite the vector potential as [1] A Note that the form of the integral in (27) echoes that of the standard line integral of F along the contour C (i.e., C F • T da [22]).However, in (27) instead of the vector dot product F • T we have used the vector Schur product F • T [26], also known as the element-wise product or the Hadamard product.The physical interpretation of the Schur (Hadamard) vector product was discussed in [1, Sec.III-A].

IV. SOLUTION OF A
We obtain the following expression of the magnetic vector potential for the helical current source, by substituting the Green's function from (4), the tangent vector from (10), and Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the parametric current from (26), into (27) where R is given in (15).The number of turns-per-layer may be specified by the difference between the upper and lower limits of the definite integral, where τ 1 is the starting position and τ 2 is the ending position on the helix, e.g., τ 1 = 0 and τ 2 = 2πN T .

VI. SOLUTION OF THE MAGNETIC FLUX LINKAGE, Φ
We may obtain Φ through substitution of (28) for A in (6) and evaluation of a line integral (see [1,Sec. VI]).We may compute Φ(z) at a specified position z, by using ( 28) and (23) while evaluating A at ρ = ρ 1 = r L − r i , which corresponds to the wire insulation's boundary toward the core's center, as where 15) was used in setting R b = R| ρ=ρ 1 , and (23) was used to make the substitution d | ρ =ρ 1 = {0 ρ, ρ 1 φ, d i 2π ẑ}dϕ.We set φ = 0 based on the periodic nature of cosine.We used the differential variable dϕ for the outer integral, whereas dτ was used for the inner integral.We may obtain the total flux of the solenoid Φ solenoid by integrating Φ(z) w.r.t.z over the helix's height h c , given in the following equation.Modeling multiple winding layers may proceed as described in [1, Sec.VIII]

VII. REDUCTION TO THE CIRCULAR LOOP
If d i = 0, I z = 0 and the helix collapses to a circular loop with a single turn-per-layer centered at z .Consequently, 0 ≤ τ ≤ 2π, A z = 0, I ρ = 0, and (29) reduces to 1 ρ dτ (32) where Furthermore, the magnetic flux linkage for the circular loop may be expressed as where R b was defined in (30) and now d i = 0. Note that the helical expressions ( 29)-( 31) support a natural formulation for the efficient modeling of multiple turnsper-layer, whereas the fields resulting from the circular-loops expressions (32), (33) must be summed via superposition to model a solenoid with N T > 1.Furthermore, in addition to the ρ and ẑ field-components, the M H provides expressions for the φ field-component while including contributions of I z , whereas the M CL is unable to quantify the aforementioned attributes.

VIII. IMPEDANCE Z, RESISTANCE R σ , AND E
We may define the impedance, Z, across the coil's terminals as Z = R σ + jωL, where {Z} = R σ represents the losses due to dissipation in a finite-conductivity electrical conductor.R σ for the solenoid may be approximated as [2] where R pul (Ω/m) is the resistance per-unit-length (pul) of the wire and may be found through (12) on [24, p. 185].If E is needed, then (29) may be used in (2) to yield and the voltage along contour C may be obtained via Maxwell-Faraday's relation; i.e., V C = C E. d = −jωΦ [19].
The formulation in Section II is based on a single-phase constant β [in (4)], which represents EM propagation in homogeneous and lossless media covering both inside and outside the solenoid, and it readily provides the EM fields for a solenoid whose core has the same material properties as its surroundings, for example, an air-core solenoid placed in air (with μ r core = 1, r core = 1, and σ core = 0) leads to the propagation constant γ core = γ air = jβ = jω √ μ 0 0 .
However, a magnetic core (which is typically a good conductor (σ 2 core /(ω core ) 2 1) has relative permeability μ r core > 1.0 and electrical conductivity σ core = 0, thus, γ core = jω √ μ core ˆ core = γ air , where the core's complex permittivity ˆ core = 0 r core + (σ core /(jω)).This implies EM propagation across heterogeneous and lossy media, and prompts us to propose the following five-step algorithm for finding the EM fields inside and outside a magnetic core solenoid placed in the air.
Step 1: In ( 29), set μ = μ core , set β = ω √ μ core ˆ core , and express the magnetic flux density for inside and on the surface of the core, as Step 2: Use B in to compute the equivalent surface current density, J s on the surface of the cylinder, leading to where the surface current density on the hollow cylindrical tube is , and on the bottom disk it is Step 3: Apply the surface equivalence principle [19] to replace the solenoid structure with its equivalent surface current density J s , then use J s together with the free-space Green's function to compute the magnetic vector potential ( A out ) outside of the core, as where G( r, r ) = e −jβR /R, β = ω √ μ 0 0 , and R was defined in (13).
Step 4: Compute the magnetic flux density vector ( B out = ∇ × A out ) outside of the core.Note that the curl operator is w.r.t the observation (unprimed) coordinates, whereas the surface integral is w.r.t. the source (primed) coordinates, thus we may take the curl inside the integral to obtain where ι = ((e −jβR (1+jβR))/R 1/3 ){ρ κ 1 , φ κ 2 , ẑ ((κ 3 +κ 4 )J s φ + κ 5 J s ρ )/(ρ(1 + jβR))}, and the coefficients Step 5: Finally, the magnetic flux density ( B out ) outside the core may be found by evaluating the surface integral in (39) on the cylindrical surface of the core, per (40).For a given solenoid, the current densities need to be computed only once, then may be reused to accelerate the evaluation of the surface integrals at different observation points, as needed (40)

X. RESULTS AND DISCUSSION
Unless stated otherwise, for the plots in this section, we assume μ r = 1, I 0 = 1, r c = 2 mm, AWG 51 wire with r w = 0.0112 mm, r i = 2 × r w = 0.0224 mm, N T = 100, and N L = 5.Note that the solenoid's height is defined as h c = N T × d i , thus, h c increases linearly with increasing N T , however, h c is unaffected by increasing N L that also increases the number of turns.We place the above described solenoid along z between z = 0 and z = h c , and along ρ centered at ρ = 0. Unless stated otherwise, the field plots in this section are evaluated at ρ = r c /2, φ = 0, and z = h c /4.

A. Calibration Coefficient
To account for the finite radius of conductor, we propose a calibration strategy based on normalization of A by the computed total current (comprised of impressed, conduction, and displacement currents) around a single wire of finite radius r w .Note that ( 28) and ( 29) are contour integrals over the infinitesimally thin helical line expressed in (9), and they do not consider the finite radius of the wire's conductor depicted in Fig. 1(b), where S 1 is the region expected to be occupied by the conductor of the finite-radius wire.Prior to calibration, the infinitesimally thin wire's conductor may produce frequency-dependent and complex displacement current through the surface S 1 , however, as we assumed that the wire's conductor is PEC (for simplicity), S 1 must be void of any displacement current.Also, the conduction current in PEC is zero, thus, the total computed current in S 1 must be solely comprised of the impressed (source) current.Referring to Fig. 1(b), the constant calibration coefficient χ may be computed using Ampere's law, as where I computed is obtained through contour integration (in spherical coordinates) of the magnetic field H θ around the contour C 1 enclosing the wire's conductor surface S 1 , H θ is obtained through transformation of B in (29) from cylindrical coordinates to spherical coordinates [19], A cal is the calibrated vector potential, and I φ and A were defined, respectively, in ( 26) and (27).Fig. 2 shows χ for the abovementioned example, where Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the deviations are both frequency-dependent and complex, as expected.An alternative calibration strategy based on experimental measurement of the terminal inductance L was described in [1, Sec.XII-E].

B. Anatomy of H
The proposed model may be used to efficiently solve the EM fields/waves across the frequency spectrum from dc to high frequencies.Since the static model was previously discussed in [1], herein we focus on the high-frequency regime.As an example, we choose the source frequency at f s = 70 GHz and plot H in Fig. 3 and compare the magnetic field intensity of the helical model (M H ) in ( 29) versus the aggregated circular-loops model (M CL ) in (32), where we assume tight wire winding with d i = 2r w .We observe that M H and M CL are in excellent agreement for H ρ (ρ), H ρ (z), H z (ρ), and H z (z).For H ρ (φ) and H z (φ), M H is able to decipher field variations due to slight helical asymmetry embedded in the example solenoid, whereas M CL is φ-invariant due to the symmetric circular loops.Next, we increase the pitch between adjacent wires to d i = 20r w and show the results in Fig. 4. We observe that the agreement between M H and M CL is not as good compared to the tightly-wound case (due to limitations of M CL ) and the field's φ-variance is now commensurate with variance along the other two axes.In Fig. 5, M H predicts H φ , whereas M CL cannot do so, note that the φ-variance increases with increasing pitch d i .In Fig. 6, we assume a magnetic core comprised of pure iron with μ r = 5000 and σ = 10 7 S/m, set f s = 60 Hz, and use M H according to the five-step algorithm outlined in (36)-(40) to predict B, as expected, in all cases the flux density is relatively larger in the core compared to outside.Note the discontinuities in B at transitions between core and air, due to presence of J s on the core's surface.

C. Computational Complexity
Fig. 7 demonstrates the computation time for the M H expressed by (29) and for the M CL expressed by (32), where CPU time is reported for computing each field component at a single point in 3 and the fit for B φ falls in between.In contrast, the computation time for M CL is approximately of order O(N T ), where fit ρ 0.1103 + 0.0051N T for B ρ , and fit z ≈ 0.0326 + 0.0067N T for B z .For both the M H and M CL , the computation time is expected to increase linearly with N L , i.e., O(N L ), based on (43) in [1].
The timing estimates are based on numerical evaluation of the integrals [27] that are computed serially on one core of a modest desktop computer with a single CPU [28] and 16 GB RAM.If field values are desired at a large number of spatial points, then the corresponding field integrals may be conveniently distributed across multiple processors in parallel for faster computation.

D. Complex Inductance L, and Radiation Resistance R r
For this section and the following section, let us assume N T = 500, and N L = 1.We use the flux Φ in (6) to compute the inductance L in (7) at terminals of the coil along its height between z = 0 and z = h c , as shown in Fig. 8(a).As can be observed, the inductance is complex and given by the following equation: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where L ∈ C, with real part {L} = L r and imaginary part {L} = L i .There is excellent agreement between the helical and aggregated circular models.
If the inductance is purely real (L ∈ R) then the imaginary part {Z} = ωL represents the stored magnetic energy.However, as shown in Fig. 8(a) and (b), the inductance is complex (L ∈ C) and frequency-dependent.Thus, we may redefine the impedance where R r = −ωL i represents the radiation resistance [29].Chou et al. [30] developed a theory on the relation between radiation resistance in frequency-domain partial element equivalent circuit models and complex inductance.As shown in Fig. 8(b), L r and L i span over positive, zero, and negative values.A positive L r implies that the coil acts as a magnetic energy storage device (inductor), whereas a negative L r implies the coil acts as an electric energy storage device (capacitor).A negative L i leads to R r > 0 and implies energy is dissipated (lost) as radiation, while a positive L i leads to R r < 0 and implies energy is generated by the coil.Note that in Fig. 8(b), L i is negative for most frequencies but there are a few frequencies at which L i ≥ 0, this implies a nonpassive model [31] of the coil system and it occurs likely because we assumed lossless (nonphysical) material properties through { , μ} ∈ R leading to β ∈ R in the Green's function of (4).If R σ + R r ≤ 0, then using the Z model of (43) in time-domain simulations may lead to instability, this may be addressed by using lossy material ({ , μ} ∈ C) in ( 4).

E. Radiation Patterns and Total Radiation
Using (29), the radiation pattern in the Fraunhofer region (ρ > 2h 2 c /λ [29]) is plotted for several interesting frequencies, as shown in Fig. 9.For all frequencies, the current source I 0 = 1 in (25).The resonant frequency, f res , of a structure is defined as the frequency at which the total stored (magnetic and electric) energy is zero, leading to {Z(f res )} = 0, which happens when L r = 0 in (43).In our example coil, these resonant frequencies occur at f res {52, 92, 125, 168} GHz.
At the MQS frequency of f = 60 Hz where the inductance is real (L i = 0) as shown in Fig. 8(b), the radiation level is below approximately 10 −4 (A/m) 2 across all angles [see Fig. 9(a)], resulting in negligible total radiation [see Fig. 9(b)].As frequency increases to 1 GHz, L i becomes increasingly negative, corresponding to an increase in the radiation level, which remains below roughly 1 (A/m) 2 across all angles.
At f = f res 52 GHz where |L i | is larger than at all other resonant frequencies, the maximum radiation across all angles reaches nearly 4.2 × 10 5 (A/m) 2 , and the total radiation reaches its maximum value of nearly 3.5 × 10 6 (A/m) 2 .Such results may contribute to assessment of the solenoid's EMC, e.g., radiated emissions [32].At f = f res 92 GHz, L i = L r = 0 and the coil neither stores nor radiates energy, thus, all input energy is dissipated in the conductor represented by R σ .As shown in Fig. 8(a), note that L varies sharply for N T < 100 and asymptotically converges to a constant for N T 200.At N T ≈ 45, L r = 0 where the solenoid experiences resonance at f 70 GHz.Note that for N T 45, L r > 0 implies stored magnetic energy dominates, whereas for N T 45, L r < 0 implies stored electric energy dominates in the solenoid.

XI. CONCLUSION
We extended the model in [1] to enable the efficient computation of the dynamic EM fields/waves of 3-D finite-height cylindrical solenoids with helical winding of wires with finite radius.Parametric integral expressions, based on the Schur vector product of I in (26) and T in (10), were derived for all three cylindrical field components of the helical solenoid.The dynamic EM fields/waves of an example solenoid were examined in detail using two different models: (a) the proposed helical model M H , and (b) the model based on superposition of circular loops M CL .Results obtained from the two models were shown to be in very good agreement for the tightly-wound solenoid.Several items may be worth noting: 1) M H could detect field variations due to helical asymmetries, whereas M CL was invariant w.r.t.angle φ; 2) M H provided detailed information about the angular component of the dynamic EM fields/waves; whereas the angular field-component was missing in M CL ; 3) as shown in Fig. 7, the computational complexity of M H was considerably better than M CL , with approximately O(log(N T )) for the former and O(N T ) for the latter.We outlined a five-step algorithm through (36)-(40) to predict B inside and outside a solenoid with magnetic core, placed in air.A complex scaling coefficient χ(ω) in (41) based on the computed total current was proposed for calibration of the helical model to account for the finite radius of the wire.
The complex inductance of (42) was discussed and its imaginary component was correlated to the radiation resistance R r in (43), the radiation patterns, and the total radiation, these results may enable the assessment of the solenoid's EMC, e.g., radiated emissions.
A plan is being developed for experimental validation of the analytical model and results through laboratory measurement of fabricated solenoid hardware.Current limitations in the proposed model, which motivate some future topics of research include the following: 1) modeling the frequency-dependent and nonuniform current distribution across the wire's cross-section; 2) regularization of the singular integrals for inclusion of the source (wire) regions; 3) modeling nonlinear core material; 4) developing an inverse model that takes the measured H from the solenoid as input and computes the current components in the solenoid's wire winding as output, e.g., for analysis and design-optimization of current probes.

Fig. 1 .
Fig. 1.(a) Cylindrical solenoid with core (gray cylinder), helical wire winding (blue contour), core height h c , and core radius r c .(b) Cross-section view of the wires [1].C 1 designates the closed contour, which encloses the conductor surface S 1 ; these are defined in spherical coordinates { R, θ, φ}.

Fig. 4 .
Fig. 4. Field plots for wide pitch with d i = 20r w .(a)-(c) Comparing H ρ for helix versus circular loops.(d)-(f) Comparing H z for helix versus circular loops.

Fig. 7 .
Fig.7demonstrates the computation time for the M H expressed by(29) and for the M CL expressed by(32), where CPU time is reported for computing each field component at a single point in3-D space versus N T , for N L = 1.The computation time (in seconds) for 100 ≤ N T ≤ 10 k is approximately of order O(log(N T )) for M H , where fit ρ −2.4359 − 0.0156 √ N T + 0.5018 log(N T ) for B ρ , fit z ≈ −4.0493 − 0.0211 √ N T + 0.7948 log(N T ) for B z ,

Fig. 9 .
Fig. 9. (a) Polar plot of the radiation pattern in Fraunhofer region at several frequencies, where the radial axis is in Log 10 scale.(b) Total radiation summed overall θ versus frequency.