Homotopy transitions and 3D magnetic solitons

Chiral magnets are a distinct class of materials, where the competition between the Heisenberg exchange interaction and the chiral Dzyaloshinskii-Moriya interaction^1,2 (DMI) gives rise to a wide variety of magnetic solitons—localized magnetic textures with particle-like properties. The most representative example of magnetic solitons in chiral magnets are magnetic skyrmions, which have been intensively studied both theoretically^3–5 and experimentally^6–12 in the last decades. In thick films and bulky samples, the magnetization vector field of chiral skyrmions resembles a filamentary structure composed of vortex-like tubes.

More recent studies have revealed a variety of 3D solitons in chiral magnets beyond skyrmions. The most prominent examples are chiral bobbers,^13,14 skyrmion bags,^15–18 heliknotons,¹⁹ and skyrmion braids.²⁰ Because of their distinct topology, the above-mentioned solitons possess different static, dynamic, and transport properties. In particular, the chiral bobbers are distinct by the presence of magnetic singularity—a Bloch point, while other magnetic textures are characterized by smooth magnetic vector fields. Because of this feature, contrary to other solitons, chiral bobbers are not topological solitons. Skyrmions and heliknotons, on the other hand, belong to different topological groups characterized by distinct topological invariants.

In this work, we present another type of magnetic solitons, which belong to a topological group of skyrmions and are stabilized in bulk crystals of chiral magnets. Because of the unique properties of these solutions, it is reasonable to classify them as a distinct class of magnetic solitons. Here, we will refer to them as hybrid skyrmion tubes. The cross sections of well-studied skyrmion tubes usually represent nearly identical configurations with small additional modulations due to the braiding effect²⁰ or the presence of a noncollinear background²¹ and/or free surfaces.²² We show that, besides such nearly homogeneous tubes, there are also solutions where the cross sections represent a continuous transformation between the skyrmions of different configurations but the same topological index—skyrmions of the same homotopy group. In a representative example, we illustrate the stability of such exotic spin textures and discuss their unique dynamic properties. In addition, we discuss an important consequence of applying such a homotopical concept to truly 3D localized magnetic spin textures. In particular, we present the solution representing a topologically trivial but statically stable and truly three-dimensional soliton—3D chiral droplet.

The model Hamiltonian for a chiral magnet with bulk-type DMI has the following form:

where n = M/M_s is the normalized magnetization |n| = 1 and $A$ and $D$ are micromagnetic constants of the Heisenberg exchange interaction and the DMI, respectively. V_m is the volume of the magnetic sample. The potential energy term U(n) includes an uniaxial (easy-axis or easy-plane) magnetic anisotropy, $Ku$⁠, and the interaction with the external magnetic field, B_ext∥e_z, is given by

Following the standard procedure,^3–5 the functional (1) can be written in the following dimensionless form that is more suited for analysis:

where $E=E/2A$ is the reduced energy and $u=Ku/MsBD$ and h = B_ext/B_D are the dimensionless values of the anisotropy constant and the strength of the external magnetic field, respectively. The integration in (3) is carried out over the reduced volume $V=Vm/LD3$⁠. The characteristic parameters $LD=4πA/D$ and $BD=D2/2MsA$ are the period of spin spiral and the saturation field for the isotropic case, $Ku=0$⁠.

The energy functional (3) remains valid for 2D (or quasi-2D) systems where the magnetization does not change along the z-axis and integration is carried out over dV = ldxdy, where l is the film thickness. The 2D model of the chiral magnet, besides ordinary axially symmetric π-skyrmions, provides a variety of magnetic solitons, e.g. skyrmion bags^15,16 and skyrmions with chiral kinks.¹⁷ The homotopy classification of localized solutions in 2D is provided by the following invariant:

The solutions with identical integer index Q are called the solutions of one homotopy class. The latter implies that one can continuously transform the vector fields of such solutions into each other. For instance, let us assume there are two localized magnetic textures n₁(x, y) and n₂(x, y) of identical charge Q. Then, one can introduce a vector field, n(x, y; s), which depends on an additional parameter s ∈ [0, 1]and possesses the property that for n(x, y; 0) = n₁(x, y) while for n(x, y; 1) = n₂(x, y). When n(x, y; s) is differentiable with respect to x, y, and s at any point, such transformation can be called a homotopy path. Since for any n₁(x, y) and n₂(x, y), there is an infinite number of homotopy paths, there is no unique method to construct such paths. To make the search more definite, we consider only those homotopy paths that also satisfy the minimum energy path (MEP) criterion. In particular, we use the geodesic nudged elastic band (GNEB) method^23,24 implemented in Spirit code.²⁵ The details of the MEP calculation are provided in  Appendix A.

In Fig. 1(a), we provide a representative example of the MEP between two skyrmions with Q = 1. We show only the spin texture corresponding to local minima and saddle points. Noticeably, the textures corresponding to the central minimum state (d) and saddle points (c), (e) contain a chiral kink¹⁷ while another minimum state is a skyrmion bag free of kinks (b). The MEP presented in Fig. 1(a) satisfies the above criteria for a homotopy path. The parameter s can be associated with the reaction coordinate, which has a meaning of the relative distance between the images (snapshots of the vector field) in the multidimensional parameter space.

To construct an initial configuration for the 3D skyrmion tube, one can take the stable 2D skyrmion configuration and place it at each xy-plane of the 3D simulated domain. A statically stable configuration of such a homogeneous skyrmion tube corresponding to the 2D skyrmion in Fig. 1(b) is provided in Fig. 2(a). To construct the initial state for a nonhomogeneous 3D skyrmion tube, we use a mapping from the homotopy path to the third spatial axis, s → z. In other words, to create a 3D magnetic texture, we sequentially lay down the spin texture of the images from MEP on top of each other along with the z-axis. The statically stable spin configuration obtained by the energy minimization of that initial state is provided in Fig. 2(b) and its cross sections are shown in (c)–(h). The intermediate region resembling a knot on the isosurface of the skyrmion string in Fig. 2(b) is well localized and, thus, can be thought of as a soliton that is hosted by a skyrmion string. This nonhomogeneous 3D skyrmion tube is a representative example of the solutions we refer to as hybrid skyrmion tubes. For conciseness, the localized intermediate region is called a knot in the following. The choice of such terminology is justified by pure visual analogy and has nothing to do with the topological knots. Whether one can continuously unwind that knot without the appearance of Bloch points represents an intriguing question that, however, goes beyond the scope of the present work.

It is worth mentioning that the inhomogeneities along the skyrmion tubes are always present near the free edges of the sample. This effect, known as the chiral surface twist,²² was discussed in detail in Ref. 18 for skyrmion bundles. In contrast to such edge localized inhomogeneities, the knots at hybrid skyrmion tubes have an additional degree of freedom to move along the skyrmion tube.

In a strong magnetic field, h ≥ 1 − 2u, the hybrid skyrmion tubes represent a metastable state embedded in the ferromagnetic vacuum. For h < 1 − 2u, the vacuum for such states is the cone phase. The range of fields and anisotropies where hybrid skyrmion tubes remain stable depends on the particular configuration—the type of 2D skyrmions in the tube cross sections. For instance, the hybrid skyrmion tube depicted in Fig. 2(b) at u = 0 remains stable at least in the range 0 < h < 0.35.

In general, the hybrid skyrmion tubes can be composed of a few intermediate states representing stable 2D skyrmions. Here, we consider the solution composed of two intermediate states only. Such configurations are easier to handle because they satisfy periodic boundary conditions. It is worth noting that the stabilization of hybrid skyrmion tubes does not require periodic boundary conditions along the z-axis. The latter makes the experimental observation of hybrid skyrmion tubes in thick films of chiral magnets quite promising. In the case of free boundaries along the z-axis, the tube has additional surface twist modulations. For sufficiently thick films (⁠$∼6LD=420$ nm for FeGe), these modulations do not decrease the tube stability. On the other hand, in the presence of the free boundaries, the hybrid skyrmion tube depicted in Fig. 2(b) can continuously unwind into the host skyrmion tube shown in Fig. 2(a). This process is illustrated by Movie 1 in the supplementary material, where each frame was obtained by a manual shift of the entire spin texture along the z-axis and following incomplete relaxation of the system. In the case of complete relaxation, the system converges to its initial state, which indicates the presence of an energy barrier for the knot to escape through the free boundary.

To demonstrate the unique dynamic properties of hybrid skyrmion tubes, we perform micromagnetic simulations based on the Landau–Lifshitz–Gilbert (LLG) equation,²⁶ 

where γ is the gyromagnetic ratio, α is the Gilbert damping, and $Heff=−1MsδEδn$ is the effective field. The last term in (5) is the Zhang–Li torque²⁷ arising from the presence of electric current,

where the vector $I=jμBp(1+ξ2)−1(eMs)−1$ is proportional to the current density j, ξ is the degree of non-adiabaticity, p is the polarization of the spin current, μ_B is the Bohr magneton, and e is the electron charge.

When a 3D soliton is embedded in the FM state, the current direction can be chosen arbitrarily. On the contrary, when the soliton is embedded in nonhomogeneous vacuum, e.g., cone phase, it is essential to choose the current direction perpendicular to the q-vector of the cone to prevent the excitation of the background. To simplify the following discussions, we consider the case of a hybrid skyrmion tube in ferromagnetic vacuum.

Movie 2 in the supplementary material illustrates the dynamics of the hybrid skyrmion tube when the current is along the z-axis. The simulations were performed with Mumax²⁸ code. The example of Mumax script and the initial state files are provided in the supplementary material. We used periodic boundary conditions in all spatial directions to simulate a bulk crystal. Under conditions when the electric current is parallel to the skyrmion tube, the knot on the isosurfaces moves along the tube in the direction opposite to the current. In this case, the position of the host skyrmion tube in the xy-plane remains fixed.

In Movie 3 in the supplementary material, we show the dynamics of the skyrmion tube when the current is perpendicular to the skyrmion tube, I∥e_x. In this case, both the host skyrmion tube and the knot move. As a result, all three components of the knot velocity are nonzero.

To quantify the knot velocity, we follow the approach used in our previous work on 2D skyrmion dynamics.²⁹ Extending this approach to 3D textures, the position of the knot $R=Rx,Ry,Rz$ can be defined as follows:

where $Nij(rk)=∫nz(r)dridrj$ is the magnon density averaged in the r_k-plane, where indices i, j, k ∈ {x, y, z}. The sign ± accounts for the direction of soliton motion: along the basis vector e_k (+) or in the opposite direction (−). The integer l_k is the number of times the soliton crossed the corresponding domain boundary. By tracing the position R at each moment in time, we can estimate the instant velocity of the knot, v = dR/dt, and compare it with the results of a semi-analytical approach based on the method of collective coordinates suggested by Thiele.³⁰ Assuming rigid motion of magnetic textures with velocity v, i.e., n(r, t) = n(r − vt), from (5) and (6), one can derive the Thiele equation, i.e.,

where the gyro-vector, G, and the dissipation tensor, $Γ̂$⁠, have components G_i and Γ_ij defined as follows:

where ɛ_ijk is the Levi–Cività symbol. The advantage of the Thiele approach is that the solution for the soliton velocity can be written explicitly. For instance, for I = Ie_x, the solution of (8) for uniform skyrmion tubes, Γ_xz = Γ_yz = Γ_zz = 0 takes a form identical to that of the 2D case,

For the hybridized skyrmion tube, all the entries of the dissipation tensor $Γ̂$ have nonzero values and the solution of (8) takes the following form:

It is worth noting two remarkable features of the solution (12). First, there is no continuous transition between (12) and (11) even when Γ_xz, Γ_yz, and Γ_zz tend to zero. Second, in the most general case of 3D solitons, the components of the gyro-vector G can be thought of as 2D topological charges in corresponding x, y, and z cross sections. For hybrid skyrmion tubes, however, similar to the 2D case, only the third component of the gyro-vector is nonzero and proportional to 2D topological charge in the xy-plane, G = 4πQL_ze_z.

To calculate the velocity components in (12), one has to find the components of the gyro-vector (9) and the dissipation tensor (10) for a particular magnetization distribution. For this, we use the magnetization vector field in static equilibrium obtained by the numerical energy minimization of the functional (1). It is worth emphasizing that the soliton velocity obtained from the solutions of the Thiele equation given by Eq. (8) corresponds to the velocities in a steady-state only. Thus, the velocities estimated from LLG simulations can be compared with the analytical solutions only when the transient state is over, and all the components of v have reached saturation.

In the most general case of 3D solitons, when all three components of the velocity are nonzero, one can introduce two deflection angles: $β1=arctanvy/vx$ and $β2=π/2−arctanvx2+vy2/vz$⁠. The angle β₁ is also known as the skyrmion Hall angle, because of which we refer to β₂ as the second skyrmion Hall angle. Table I demonstrates good agreement between the values of β₁ and β₂ estimated from LLG simulations and those calculated with the Thiele approach for parameters j = −5 × 10⁸e_x A/m², α = 0.01, and ξ = 0.05. Noticeably, for the particular tube depicted in Fig. 2(b), the angle β₂ is larger than β₁. As follows from the Thiele equation (8), when I = Ie_z, in agreement with the micromagnetic simulations, only the z-component of the knot velocity is nonzero, v_z = −Iξ/α.

The homotopy concept used for constructing hybrid skyrmion tubes can be extended further for solitons localized in all three dimensions. For this, one can consider the homotopy transformation of any 2D skyrmion with Q = 0 into the FM phase. Here, we examine such transformation on the example of chiral droplet shown in Fig. 3. The MEP in Fig. 3 is quite similar to that shown in Fig. 1. Two stable states representing the isolated soliton and the saturated FM state are separated by the saddle point. In this transformation, the collapse of the droplet goes by its shrinking, which has been previously discussed in Refs. 31 and 32.

By mapping the obtained images of the corresponding homotopy path onto the spatial z-axis, we constructed the initial state of the 3D soliton. However, it turned out that such an initial guess does not lead to a stable configuration in the FM phase but instead becomes stable in the cone phase vacuum only. The isosurfaces and cross sections of the stable configuration of the 3D chiral droplet obtained by numerical energy minimization are shown in Fig. 4. As seen in (a), the red color in the surface is missing, which means that by projecting the vector field of the 3D chiral droplet onto an S₂ sphere, one cannot cover it entirely. The latter means that both the topological invariant (4) and Hopf invariant are zero in this case. The cross sections in (b) and (c) show well-localization of the soliton. We estimated the characteristic size of the 3D chiral droplet equal to $∼1LD$ in all spatial directions.

The stability range for the 3D chiral droplet in terms of the magnetic field, h, and anisotropy, u, is shown in Fig. 5. Note that the phase transition lines between the phases are reproduced from Ref. 33. As follows from the diagram, the 3D chiral droplet stability region is in the range where the skyrmion lattice is the lowest energy state, while the cone phase is a meta-stable state. The blurred edges at h ≲ 0.28 and u ≲ 0.15 denote that we did not succeed in identifying the stability range below these values reliably. At these values of magnetic field and anisotropy, the helicoid competes in energy with the cone phase. A more precise estimation of the lower bound for the 3D chiral droplet stability requires much larger sizes of the simulated domain above the limit of our computational resources.

In isotropic chiral magnets with weak magnetocrystalline anisotropy, the uniaxial anisotropy required for 3D chiral droplet stability can be induced by demagnetizing fields when the sample has a prolate shape of an ellipsoid or a wire, in the limiting case. The synthesis of such nanowires of B20-type crystals was reported earlier for MnSi³⁴ and Fe_1−xCo_xSi.³⁵ 

The 3D chiral droplet motion can be induced by applying different external stimuli. First, let us consider the motion of the 3D chiral droplet induced by the external magnetic field slightly tilted and rotating about the z-axis, h = 0.34(sin θ_h cos 2πνt, sin θ_h sin 2πνt, cos θ_h). We performed the micromagnetic simulations assuming the tilt angle θ_h = 0.1 and the small frequency of rotation ν = 5 MHz, which is close to a quasi-stationary process. Movie 4 in the supplementary material demonstrates such motion of the 3D chiral droplet consisting of two types of motion—translation along the z-axis and rotation about the z-axis. Remarkably, similar dynamics has been reported for heliknoton¹⁹ and this seems to be a common property of all 3D localized solitons hosted by a conical phase. The motion of the 3D chiral droplet free of rotation can be induced by the cone phase rotation under a static external field. To demonstrate this, we performed the simulations where we manually rotated a small number of spins located far from the droplet. As illustrated in the supplementary material, Movie 5, the translation motion of the droplet along the z-axis, in this case, occurs without its rotation.

In the case of Zhang–Li torque, the gyro-vector of the 3D droplet equals zero and the Thiele equation (8) provides a trivial solution, v = ξI/α. On the other hand, due to the presence of the cone phase, one has to choose I ⊥ q to suppress the excitation of the cone phase itself. An analysis of the cone phase dynamics induced by I∥q is provided in  Appendix B.

One has to note that contrary to the Thiele equation, which predicts that the 3D chiral droplet should move without deflection, v∥−I, in numerical simulations, we still observe a slight deflection in the soliton motion. Since, in the case of the hybrid skyrmion tube, we got pretty good agreement between numerical simulations and the Thiele approach, we attribute this effect to the numerical artifact caused by the excitation of the cone phase, which, in turn, prevents reaching the regime of a steady motion. Movie 6 in the supplementary material shows the dynamics of the droplet caused by the in-plane electric current I∥e_x. To prevent the excitation of the cone phase, in this simulation, we pinned the spins on the bottom plane of the simulated domain.

In this work, we have discussed new types of 3D magnetic solitons stabilized in chiral magnets—hybrid skyrmion tube and 3D chiral droplet. We show that the magnetic textures of these solutions can be explained in terms of homotopy transitions between various 2D skyrmions.

We have studied the static and dynamic properties of new solitons. The parameters of the external magnetic field and magnetic anisotropy at which the presented solutions remain stabilized are different and characterized by the various states of the vacuum. Hybrid skyrmion tubes can be stable in the FM state and cone phase, while the 3D chiral droplets are stable only in the conical phase surroundings. It is shown that the dynamics of the hybrid skyrmion tube can be induced by the electric current modeled by the Zhang–Li torque in the LLG equation. The electric current applied along the tube causes the motion of the knot on the skyrmion tube in the direction opposite to that of the current. The electric current applied perpendicular to the skyrmion tube leads to more complicated dynamics, which can be described in terms of two skyrmion Hall angles. The results of LLG simulations agree well with the analysis based on the Thiele equation.

For the case of 3D chiral droplets, we show that their motion besides the electric current can be induced by a rotating external magnetic field. Since the cone phase is easily excited when the electric current is not perpendicular to the cone wave vector, an analysis based on the Thiele approach is difficult in this case. The analytic solution for the dynamics of the cone phase is provided. Movies provided in the supplementary material illustrate numerical simulations for the motion of 3D solitons.

See the supplementary material for additional movies illustrating the dynamics of 3D solitons discussed in the main text. Supplementary Movie 1 illustrates the escape of the knot on the hybrid skyrmion tube through the free boundary of the plate. Supplementary Movies 2 and 3 show the dynamics of the hybrid skyrmion tube induced by current I∥e_z and I∥e_x, respectively. Supplementary Movie 4 illustrates the dynamics of a 3D chiral droplet induced by a small precession of the external magnetic field about the z-axis with low frequency. Supplementary Movie 5 shows the dynamics of a 3D chiral droplet under the cone spiral rotation induced by a weak electric current I∥e_z applied in the bounded volume far from the droplet. Supplementary Movie 6 shows the dynamics of a 3D chiral droplet induced by current I∥e_x when the spins in the plane z = 0 are pinned.

The magnetization in the movies is presented by the isosurfaces n_z = 0 and the standard color code used in Mumax for visualization of the unit vector fields. The exceptions are Movies 4–6 where we use red–blue color code for the n_y-component of magnetization: n = (0, −1, 0), red and n = (0, 1, 0), blue.

We thank Filipp Rybakov, Bernd Schroers, and Bruno Barton-Singer for fruitful discussions. We are also grateful to Moritz Sallermann for technical support in minimum energy path calculations with Spirit. We acknowledge financial support from the Deutsche Forschungsgemeinschaft through SPP 2137 “Skyrmionics” Grant No. KI 2078/1-1 and the European Research Council under the European Union’s Horizon 2020 Research and Innovation Programme (Grant No. 856538 - project “3D MAGiC”).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The MEP in Fig. 1(a) represents a homotopy path between two states (b) and (d). To get the MEP without singularities, one has to make a reasonable assumption for the initial path. The straightforward and often used approach based on the interpolation between two configurations by simple rotation of the spins is not appropriate for this purpose. It is easy to check that the topological index Q is often not conserved along the MEP with this approach. To avoid such a behavior of the GNEB solver, we have constructed the ansatz for the initial path using the energy minimization and drag-and-drop function implemented in Spirit. In particular, as seen in Figs. 1(b)–1(e), the intermediate states represent the elongation and twisting of the magnetic texture near the chiral kink. Starting with the state depicted in Fig. 1(b) and using the drag-and-drop option, we enforce the chiral kink to elongate and bend to form a state similar to that in (c) and (e). We intentionally stop the relaxation process before the system converges to one of two local minima as in Fig. 1(b) or Fig. 1(d). Then, we collect such unrelaxed snapshots of the spin texture with the different levels of elongation of the part containing the chiral kink and use them as an initial guess for MEP.

The calculations were performed using the following parameters: L_x = L_y = 4L_D with L_D = 64a.

In terms of spherical angles (Θ, Φ), the magnetization reads as follows: n = (sin Θ cos Φ, sin Θ sin Φ, cos Θ)^T. Assuming that at t = 0, the magnetization profile is given by the cone phase Θ = Θ_c, Φ = 2πz + ϕ₀, we can find an analytic solution for the LLG equation with the current I = Ie_z turned on at t > 0. The uniformity of the magnetic texture in the xy-plane leads to the following system of equations for (Θ, Φ):

Here, τ = γB_Dt/4π² is the dimensionless time and i = 4π²I/γL_DB_D is the parameter of the electric current. Although the system (B1) represents two coupled nonlinear partial differential equations, its solution can be written for the case when Θ = Θ(τ).

We found that strong currents can lead to a transition from the cone phase to the FM state. The critical current at which this happens is given by

The formula (B2) is meaningful only at ξ ≠ α. In the case ξ = α, the cone phase cannot be excited by the current. Below the critical regime, i.e., $i∈(ic−,ic+)$⁠, the angle Θ monotonically changes from Θ_c—the equilibrium cone angle without the electric current to Θ_I—the equilibrium cone angle in dynamical steady state in the presence of the current:

For $i∉(ic−,ic+)$⁠, the angle Θ_I equals to 0 or π depending on the sign of the current, i. At $i∈(ic−,ic+)$⁠, the solution of (B1) for Θ can be written as follows:

where a₁ = h − i(ξ − α)/2πα and a₂ = 1 − 2u. The solution for Φ is given by

As follows from (B4), at τ → ∞, one has cos Θ = cos Θ_I = a₁/a₂. In this case, Eq. (B5) describes the rotation of the cone phase with constant velocity Iξ/α, which agrees with the Thiele prediction for topologically trivial configurations. Noticeably, the dynamics described by the LLG equation reaches a dynamically steady state exponentially fast. Note also that solutions (B4) and (B5) remain valid not only for the bulk systems but also for the films with free boundaries along the z-axis. The latter remain valid for any h and u where the surface modulations¹³ do not perturb the cone phase. This follows from the fact that the solutions obtained satisfy the boundary conditions: ∂_zΘ = 0, ∂_zΦ = 2π on the free surface, z = const, for any τ.

To verify the solutions obtained, we compare them with the results of LLG simulations performed for different values of the current but fixed values of h = 0.34, u = 0.26, α = 0.01, and ξ = 0.05. The critical current values (B2) for these parameters are $ic−≈−0.22$⁠, $ic+≈1.29$⁠. We have performed simulations for two current values within the critical range: i₁ = −0.15, i₂ = 0.15, and for one outside this range i₃ = 2. The results are provided in Fig. 6. For all the three cases, we see a good agreement between numerical and analytical solutions at least for the chosen simulation time. At longer times, the z-component of the magnetization tends to the limit value accordingly to (B3). For the cases shown in (a) and (b), these limit values are about 0.51 and 0.91, respectively. For the case shown in (c) corresponding to the current above the critical value, the magnetization tends to −1. As one can deduce, applying the current I∥q allows manipulating the cone phase angle and its dynamics in a controllable way.

The obtained analytical solutions, besides their pure academic interest, can be used also for testing the accuracy of numerical schemes for solving the LLG equation. In particular, the solutions of (B1) in a special case of zero damping, α = 0, can be written in a more compact form as follows:

The stability of LLG solvers typically requires the presence of nonzero damping. Therefore, the solutions (B6) can be useful for testing new LLG solvers, which are free of this limitation. Both cases α = 0 and α ≠ 0 can be generalized for the case of time-dependent currents I = I(t) easily. A detailed analysis of this case is beyond the scope of the present study and will be provided elsewhere.