A New Approach to the Modeling of Anisotropic Media with the Transmission Line Matrix Method

: A reformulation of the Transmission Line Matrix (TLM) method is presented to model non-dispersive anisotropic media. Two TLM-based solutions to solve this problem can already be found in the literature, each one with an interesting feature. One can be considered a more conceptual approach, close to the TLM fundamentals, which identiﬁes each TLM in Maxwell’s equations with a speciﬁc line. But this simplicity is achieved at the expense of an increase in the memory storage requirements of a general situation. The second existing solution is a more powerful and general formulation that avoids this increase in memory storage. However, it is based on signal processing techniques and considerably deviates from the original TLM method, which may complicate its dissemination in the scientiﬁc community. The reformulation presented in this work exploits the beneﬁts of both methods. On the one hand, it maintains the direct and conceptual approach of the original TLM, which may help to better understand it, allowing for its future use and improvement by other authors. On the other hand, the proposal includes an optimized treatment of the signals stored at the stub lines in order to limit the requirement of memory storage to only one accumulative term per ﬁeld component, as in the original TLM versions used for isotropic media. The good behavior of the proposed algorithm when applied to anisotropic media is shown by its application to different situations involving diagonal and off-diagonal tensor properties.


Introduction
The Transmission Line Matrix (TLM) method is a low-frequency, time-domain numerical method originally proposed by Johns and coworkers in the early 1970s [1] to solve electromagnetic wave propagation problems from a different point of view. The method was not a classical numerical one directly involved in solving the equations governing a certain phenomenon, as is the case with the Finite-Difference Time-Domain (FDTD) numerical method. In contrast, the TLM method is mainly based on finding an analogous transmission line circuit to describe the original problem by means of the study of the propagation of voltage and current pulses through it. The circuit (the TLM mesh) is formed by interconnecting elementary circuits (the TLM nodes), which are specifically designed to model a particular propagation or diffusion problem. Figure 1 shows the basic node used for the modeling of three-dimensional (3D) electromagnetic problems: the 3D Symmetrical Condensed Node (3DSCN). The node is said to be condensed because all the field components are defined at the node center and at the same time, which is a remarkable feature of the TLM method. The 3DSCN basically consists of twelve lines which are connected field components are defined at the node center and at the same time, which is a remarkable feature of the TLM method. The 3DSCN basically consists of twelve lines which are connected to similar lines of adjacent nodes. These twelve ports are termed the link lines, and they are mainly responsible for propagation. It is worth noting that connections at the node center are not actual connections; they are to be understood in a formal sense instead, in order to reproduce the coupling existing in the original differential equations governing the phenomenon, which is Maxwell's equations in this work [2]. On occasion, additional capacitive, inductive, or electric and magnetic loss lines are also formally connected at the node center to allow for a specific control of properties such as the electric permittivity, magnetic permeability, or electric and magnetic conductivity [3,4]. Although the electrical origin of the method has caused most of its applications to be of the electromagnetic type [5][6][7][8], TLM has also been successfully applied to particle or thermal diffusion problems [9] as well as to acoustic situations [10]. The fundamentals of TLM are similar, whatever the application. First, the medium to be modelled is discretized into elementary volumes, which are replaced with an analogous mesh of interconnected TLM nodes. Time is also considered at discrete time samples, The time-marching scheme is quite simple; for time iteration or time calculation, n, both voltage sets are connected through the scattering matrix, S  , characteristic of a particular node, in the following formula: Finally, these reflected pulses become the incident pulses for the next time step at the adjacent nodes for the link lines, and at the same node for the extra capacitive and inductive stubs. The process is repeated for the next time iteration. The singularity of the problem under study forces the particular geometry of the node and the scattering matrix elements, but the process is quite similar for different types of phenomena. The fundamentals of TLM are similar, whatever the application. First, the medium to be modelled is discretized into elementary volumes, which are replaced with an analogous mesh of interconnected TLM nodes. Time is also considered at discrete time samples, t n = n∆t. A set of incident pulses, n → V i , reaches the node center at each time step and produces a set of reflected pulses, n → V r , for that particular time step and node in the mesh. The time-marching scheme is quite simple; for time iteration or time calculation, n, both voltage sets are connected through the scattering matrix,S, characteristic of a particular node, in the following formula: Finally, these reflected pulses become the incident pulses for the next time step at the adjacent nodes for the link lines, and at the same node for the extra capacitive and inductive stubs. The process is repeated for the next time iteration. The singularity of the problem under study forces the particular geometry of the node and the scattering matrix elements, but the process is quite similar for different types of phenomena.
The time-marching scheme described above makes of TLM an alternative low-frequency numerical method with second-order accuracy, which provides appropriate results as long as the space sampling is below one-tenth of the minimum wavelength [11]. This is a common feature with other low-frequency numerical methods such as FDTD, a method with which certain equivalences have been found during the initial stages of TLM [12]. A detailed description of the TLM fundamentals and its most relevant applications, electromagnetic or not, may be found in [7].
An important feature of TLM is that it can be considered more of a conceptual approach than a purely mathematical numerical method, as is the case with the FDTD method, for instance. As might be expected, this property of the TLM method has advantages and drawbacks. From among the advantages, the most outstanding one is that the physical content of the different lines in the TLM node is clearly identified. This means that including a particular element in the node implies that all the aspects regarding this element are automatically implemented. As an example, exchanging the role or inductances and capacitances not only allows for the modeling of properties below the vacuum values or even negative values [13], but also reproduces the dispersive behavior to be expected if the corresponding metamaterial is manufactured in this way. However, some drawbacks must be mentioned. First, the algorithm devised for one type of problem is not directly extrapolated to another problem. A certain familiarization with transmission line concepts is required to take full advantage of the method and the ability to design appropriated TLM nodes for specific situations. Second and more importantly, obtaining the scattering matrix for a given node is usually a difficult task. Despite these difficulties and under this general scheme, the TLM method has proven to be a versatile and powerful numerical technique for solving different challenging electromagnetic situations. Some examples are the modeling of thin conducting plate situations, time-varying media, far-field calculations, dispersive and nonlinear materials, left-handed materials, just to name a few [6,14,15]. As regards the antenna analysis and design-a challenging problem for a low-frequency numerical method such as the FDTD or TLM methods-the addition of a specific sub-circuit representing the effects of the thin-wire antenna to the standard TLM node, shown in Figure 1, was first proposed in [16] to model antennas whose radius is significantly smaller than the node size. This was an imaginative and efficient solution which generated a number of subsequent works with new thin-wire models and even more complex situations [17][18][19][20]. Other applications of the TLM method to different antenna analysis and design problems can also be found in [21].
The conceptual TLM point of view favors the versatility of the method, but also that numerical procedures followed by the different research groups working on TLM are often quite independent, apparently different, and difficult to relate to. This may cause that some relevant advances from one group are not always completely exploited by other researchers, although the basic aspects of the TLM method are well-known within the TLM community. This is possibly the case with the modeling of anisotropic media for which the electromagnetic properties are now defined by tensor quantities. Although these media were already modeled in the early TLM works for the simple case in which principal directions are considered and thus only diagonal elements are present [7], the general case that includes off-diagonal elements introduces a coupling between the different field components, which makes the task considerably more difficult to deal with. Two qualitatively different TLM solutions have been reported to solve this difficulty in a different manner.
The first proposal is from Paul et al., who, in a brilliant set of three consecutive works [22][23][24], reformulated the TLM method to deal with general materials through a series of normalizations and schemes in which signal processing techniques, and specially the bilinear transform, played a fundamental role. This revised and elegant TLM scheme was then specified to model anisotropic and nonlinear media. The presented reformulation showed the important benefit that no matter how many effects had to be added to a series or a shunt sub-node, all these effects were comprised in one single accumulator quantity per field component. But this compact and elegant formulation separates the algorithm from the original version of the TLM and may complicate its application by other authors.
The second existing solution resulted later from a new reformulation of the TLM method [25]. This approach maintained the classical and conceptual origin of the TLM method in which an additional extra stub is identified for each component of the different tensor electromagnetic properties. However, despite this scheme being considered closer to the original TLM formulation and conceptually easier to understand for a wider TLM community, this independent treatment of tensor components necessitates the storage of the voltages at these extra stubs, which noticeably increases the memory and timecomputing requirements of the algorithm in general situations in which several off-diagonal components of the tensor quantities must be considered.
The purpose of this paper is to reformulate the TLM treatment of anisotropic materials in a new way that combines the advantages of the two existing solutions: the conceptually clear scheme of the original TLM described in [25], and the efficient reduction of memory storage and time-computing resources used in [22]. As a result, a relatively powerful but still simple algorithm is produced, which may help in making the treatment of anisotropic media with TLM an easier and more widely used methodology.
The paper is organized as follows: Section 2 describes a new TLM algorithm for the modeling of non-dispersive anisotropic media. Section 3 includes two numerical applications. The comparison of the obtained results with independent reference solutions shows the good behavior of the proposed model. Finally, conclusions are presented in Section 4.

A TLM Node for an Anisotropic Non-Dispersive Medium
In this section, a new reformulation of the TLM algorithm for dealing with nondispersive anisotropic media is considered. The description is carried out by following a scheme which resembles the procedures in the original TLM method. First, the topology of the node proposed is presented, identifying each term in Maxwell's equations with individual ports at the node, which may ease the choice of the circuit parameters. Second, the description of the Thevenin equivalent circuits, the modifications required by the presence of coupling elements caused by anisotropy and a novel unified treatment of the signal at the stubs are proposed and discussed to obtain the field and reflected pulses at each time step.

The TLM Node and the Parameters
Let us consider an anisotropic and non-dispersive material. The medium properties are defined by tensor quantities which couple the different electromagnetic field components. In the absence of coupling between the electric and the magnetic properties, the constitutive relations of such a material are given by the following formula: where → E stands for the electric field, → H is the magnetic field, → D is the electric displacement field, → B is the magnetic flux density, andĨ is the identity tensor. A more general case would include a magneto-electric tensor to describe the magnetic field effects on the electric displacement field and the electro-magnetic coupling tensor to describe the effect of the electric field on the magnetic flux density. For brevity's sake, these two effects have not been considered in this work, although they could be included without great difficulty, if required.
Regarding the constitutive parameters,ε andχ e denote the electric permittivity and electric susceptibility tensors, respectively, whileμ andχ m are the magnetic permeability and magnetic susceptibility tensors, respectively. Concerning the electric and magnetic current densities, , are the free electric current density and the free magnetic voltage density, respectively. It is worth noting that multiplications in Equations (2)- (5) should be replaced with time-domain convolution in the case of considering dispersive media, whose treatment is left for a subsequent work.
The phenomenon is governed by Maxwell's equations, Ampere's law, and Faraday's laws, given by the following: To help in the understanding of the procedure described in the following paragraphs, let us explicitly show the x component of Ampere's law and the x component of Faraday's law for the anisotropic material under study: For the usual case of an isotropic medium, these two equations define E x and H x and only require the first three terms on the right side. These are the uncoupled terms. However, for the anisotropic medium considered in this work, the last two pairs of extra coupling terms on the right side of (7) are to be considered in order to deal with the effect of E y and E z on E x or with an analogous coupling between the magnetic field components in (8).
The basic form of the TLM node for modeling this anisotropic medium is the 3DSCN shown in Figure 1, where only the main or link lines are represented. Extra lines, not represented, must also be included to model the deviation of the anisotropic medium from the vacuum behavior. In this work, the numbering of the ports corresponds to the original proposal in [26], with minor modifications to describe coupling effects. In this manner, link lines are ports 1-12, extra capacitive lines are lines 13-15, inductive lines are 16-18, electric loss lines are lines 19-21, and magnetic loss lines are numbered as lines 22-24. Appropriated subscripts for the extra lines will be added to this numbering scheme in order to denote the tensor quantity element with which a particular line is associated. The node length is ∆x, ∆y, and ∆z along the three Cartesian directions. The characteristic impedance of the link lines is, as usual, chosen to be the vacuum impedance, Z 0 = Y −1 0 = µ 0 /ε 0 , while the time step value will be specified later.
As regards the additional lines to adjust the properties of the medium, let us think of the node in Figure 1 as a vector circuit, since it must contain information on the full set of Maxwell's equations in their vector form expressed in (6). As it happens with Maxwell's vector equations, the 3DSCN may be split into coupled sub-circuits or scalar circuits: three shunt nodes to model the three components of Ampere's law, and three series nodes to describe the three components of Faraday's law. Figure 2a shows, for time t n , the node to describe the x component of Ampere's law, scalar Equation (7), while Figure 2b shows the sub-node corresponding to the x component of Faraday's law, scalar Equation (8). These nodes will be referred to hereinafter as the shunt node for E x and the series node for H x , respectively. shunt nodes to model the three components of Ampere's law, and three series nodes to describe the three components of Faraday's law. Figure 2a shows, for time tn, the node to describe the x component of Ampere's law, scalar Equation (7), while Figure 2b shows the sub-node corresponding to the x component of Faraday's law, scalar Equation (8). These nodes will be referred to hereinafter as the shunt node for Ex and the series node for Hx, respectively. The electromagnetic field components associated with this pair of scalar circuits are defined through the following analogy of the fields, with the common voltage at the shunt circuit and the common current at the series circuit: , .
It is worth noting that the minus sign used for calculating nHx, which is not strictly necessary but is maintained for historical reasons and used to make both analogies qualitatively similar, requires that the current at the series node is defined with the current opposite to the usual case, the right-hand rule, as shown in Figure 2b. Let us describe each one of these sub-nodes.
Concerning the shunt node for Ex, four main or link lines, ports 1, 2, 9, and 12, are those mainly responsible for the propagation of adjacent nodes through its connection to link lines. Ports 2 and 9 define Ex and Hy propagating along the z direction, while ports 1 and 12 define Ex and Hz propagating along the y direction. They model the curl term on the left side of (7), and depending on the choice of Δt, partially or totally describe the vacuum electric permittivity term associated with 0 ε through the capacitance Y0Δt/2 introduced by each link line to the circuit. In order to complete the model of the rest of terms on the right side of (7), an extra capacitive or open stub for each susceptibility term, e xj χ , is included, with j = x, y, z. These lines are denoted as ports 13xj, they have a characteristic admittance of Y0Yxj, and they introduce a capacitance Y0YxjΔt/2 to the circuit. An extra The red boxes at terminals AB denote a formal connection indicating that voltage V AB at common terminals AB is n V x for the link lines and for lines 13xx and 19xx, n V y for lines 13xy and 19xy, and n V z for lines 13xz and 19xz. (b) The series node for H x at time t n . Link lines, represented in dark blue, are connected to adjacent nodes. In pale blue: lines 16xj are short-circuited stubs, while lines 22xj are infinitely long stubs. The red boxes denote a formal connection indicating that the current is n i x for the link lines and lines 16xx and 22xx, n i y for lines 16xy and 22xy, and n i z for lines 16xz and 22xz.
The electromagnetic field components associated with this pair of scalar circuits are defined through the following analogy of the fields, with the common voltage at the shunt circuit and the common current at the series circuit: It is worth noting that the minus sign used for calculating n H x , which is not strictly necessary but is maintained for historical reasons and used to make both analogies qualitatively similar, requires that the current at the series node is defined with the current opposite to the usual case, the right-hand rule, as shown in Figure 2b. Let us describe each one of these sub-nodes.
Concerning the shunt node for E x , four main or link lines, ports 1, 2, 9, and 12, are those mainly responsible for the propagation of adjacent nodes through its connection to link lines. Ports 2 and 9 define E x and H y propagating along the z direction, while ports 1 and 12 define E x and H z propagating along the y direction. They model the curl term on the left side of (7), and depending on the choice of ∆t, partially or totally describe the vacuum electric permittivity term associated with ε 0 through the capacitance Y 0 ∆t/2 introduced by each link line to the circuit. In order to complete the model of the rest of terms on the right side of (7), an extra capacitive or open stub for each susceptibility term, χ e xj , is included, with j = x, y, z. These lines are denoted as ports 13xj, they have a characteristic admittance of Y 0 Y xj , and they introduce a capacitance Y 0 Y xj ∆t/2 to the circuit. An extra electric loss stub is also needed for each conductivity term, σ e xj . These lines are denoted as ports 19xj and have a characteristic admittance of Y 0 G xj , which is directly the admittance introduced by the line to the circuit since no reflected voltage is expected on this line. An independent current source, n I e x , takes into account the free electric current density term appearing in (7). Similar nodes for E y and E z can be defined with the appropriated link and source lines, together with extra capacitive ports, 14yj and 15zj, and electric loss lines, 20yj and 21zj, respectively.
As regards the series node for H x sketched in Figure 2b at time t n , it consists of the series connection of different lines. As regards the link lines, ports 4, 5, 7, and 8, they describe the curl terms on the left side of (8), and depending on ∆t, partially or totally define the vacuum permeability value, µ 0 , through the inductance Z 0 ∆t/2 introduced by each link line to the circuit. An inductive or short-circuited stub, 16xj, describes the magnetic susceptibility term due to χ m xj and a loss stub, 22xj, models the magnetic loss term due to σ m xj appearing in (8). Line 16xj has a characteristic impedance of Z 0 Z xj and introduces a total inductance of Z 0 Z xj ∆t/2 to the circuit, while line 22xj has a characteristic impedance of Z 0 R xj , which represents a resistance with this same value added in a series connection to the circuit for H x . In these lines, the index j takes the values x, y, and z. An independent voltage source, n V m x , takes into account the free magnetic voltage density term appearing in (8). The polarity of the link lines is fixed and shown in Figure 1; therefore, lines 5 and 8 contribute positively to i x , while link lines 4 and 7 do the opposite. The polarity of the rest of the lines is chosen to positively contribute to i x , as shown in Figure 2b.
In order to define the value of the characteristic admittances of the capacitive and loss lines at the circuit for E x , let us proceed as follows. As regards the capacitive lines 13xj, the circuit must define the three capacitive terms of the medium associated with the permittivity elements ε xj appearing in (7), with j = x, y, and z. The capacitances corresponding to index j is that of a parallel plate capacitor formed by a medium with permittivity ε xj between two parallel plates, with an area of S x = ∆y∆z and separated by a distance of ∆j. Taking into account that the link lines only contribute to the diagonal element, j = x, and that a capacitive line adds a capacitance equal to ∆t/2 times its characteristic admittance, the following equality must be met: Similarly, the loss line 19xj in the circuit of Figure 2a adds an admittance Y 0 G xj to the circuit that must define the admittance of a cubic volume with an area of S x and a width of ∆j, filled with a medium of conductivity of σ e xj ; thus, The more general case corresponding to the shunt node for E i , with i = x, y and z, can simply be derived by replacing S x in the previous equations, with S i = ∆a∆b, where ∆a and ∆b are the node dimensions along the two directions normal to the i direction. By doing so, the relative admittance on the capacitive and electric loss lines at the shunt nodes can be obtained through: An analogous calculation can be carried out for the inductive and magnetic loss lines at the series nodes. In this manner, the inductive lines corresponding to H i are defined by the following relative characteristic impedances: while the relative characteristic impedance of the magnetic loss lines is given by the following: Electronics 2021, 10, 2071

of 22
A maximum allowable time step may be obtained by imposing non-negative values for each diagonal element Y ii and Z ii in (12) and (14). This maximum value is usually chosen as the time step for the TLM algorithm. With this choice, the link lines define the vacuum properties if the node length is identical for all directions, while the capacitive and the inductive added stubs describe the susceptibility, which distinguishes the medium properties from the vacuum ones.
The analogy is completed by including the source elements through an electric source current at the shunt node and a magnetic source voltage at the series node, which are obtained simply by multiplying each source density quantity by the corresponding transversal area, S i :

The Field Definition and the Scattering Process
The propagation of voltage pulses through a transmission line is a dynamic process in which, at a specific moment, incident pulses disappear and become reflected ones. The situation becomes difficult when several transmission lines are to be considered simultaneously. The Thevenin Theorem is a useful tool to ease this task. The situation is depicted at the top of Figure 3. At the incident stage, an incident voltage pulse, n V i , propagates towards terminals AB, reaches these terminals at a certain time, t n , and produces a reflected voltage, n V r . This process is described at t n by the Thevenin equivalent circuit shown at the bottom of Figure 3 [7]. Despite the fact that the reflected pulse does not appear in the Thevenin circuit, it may be obtained in terms of the incident pulse and the total voltage or current on the line as follows: Electronics 2021, 10, x FOR PEER REVIEW 9 of 23 Let us now consider Figure 4, which shows the Thevenin equivalent circuit corresponding to the parallel node for Ex in Figure 2a, when all the incident pulses reach the node center at the n-th time step. The elements associated with nVx are shown against a white background. As regards the coupling terms, two lines represented in the pale blue region stand for the terms εxy and e xy σ , lines 13xy and lines 19xy. The connection with the lines associated to nVx must be understood in a formal sense, since they are actually current sources which depend on the voltage nVy. This fact is described in Figure 4 by specifying a different voltage, nVy, at this pair of lines. A similar situation is found for the lines associated with εxz and e xz σ , ports 13xz and lines 19xz, depicted in the pale red region and with a different voltage, nVz, at their terminals. Let us now consider Figure 4, which shows the Thevenin equivalent circuit corresponding to the parallel node for E x in Figure 2a, when all the incident pulses reach the node center at the n-th time step. The elements associated with n V x are shown against a white background. As regards the coupling terms, two lines represented in the pale blue region stand for the terms ε xy and σ e xy , lines 13xy and lines 19xy. The connection with the lines associated to n V x must be understood in a formal sense, since they are actually current sources which depend on the voltage n V y . This fact is described in Figure 4 by specifying a different voltage, n V y , at this pair of lines. A similar situation is found for the lines associated with ε xz and σ e xz , ports 13xz and lines 19xz, depicted in the pale red region and with a different voltage, n V z , at their terminals. region stand for the terms εxy and e xy σ , lines 13xy and lines 19xy. The connection with the lines associated to nVx must be understood in a formal sense, since they are actually current sources which depend on the voltage nVy. This fact is described in Figure 4 by specifying a different voltage, nVy, at this pair of lines. A similar situation is found for the lines associated with εxz and e xz σ , ports 13xz and lines 19xz, depicted in the pale red region and with a different voltage, nVz, at their terminals. The current flowing through line 1 is given by the following equation: As regards the rest of the lines, analogous expressions can be obtained by simply using the appropriate admittance and eliminating the incident voltage pulse in the case of loss lines. According to Kirchhoff's current law, the sum of currents along all the lines connected to the ground terminal must be zero. By doing so and after some simple calculations, the following is produced: where, Figure 4. Thevenin equivalent circuit of the shunt node for E x at time t n . Green boxes denote a formal connection, allowing for the voltage at the upper terminals to be n V x for the link lines, the source, and lines 13xx and 19xx, while the voltage is n V y for lines 13xy and 19xy (in the pale blue region) and n V z for lines 13xz and 19xz (in the pale red region).
The current flowing through line 1 is given by the following equation: As regards the rest of the lines, analogous expressions can be obtained by simply using the appropriate admittance and eliminating the incident voltage pulse in the case of loss lines. According to Kirchhoff's current law, the sum of currents along all the lines connected to the ground terminal must be zero. By doing so and after some simple calculations, the following is produced: where, The terms are organized according to their conceptual content. Thus, the first term, n L e x , refers to the contribution to n V x of the incident voltage pulses at the link lines with E x polarization. The n S e x term is due to the signal stored at the stubs associated or coupled with E x at previous time steps. This term corresponds with the accumulator term in [22], although in this work, it requires neither a node with equal length for all Cartesian directions nor a time step fixed by ∆t = ∆l/(2c), with c being the vacuum speed of light. Finally, the term n V e x models the excitation term contributing to n V x . At this point it is worth noting that the calculation of n E x using (19) requires storing pulses at three stubs per field, while the formulation presented in [22] only uses one accumulator term per component. The situation would be even worse, requiring three more stubs, if magneto-electric coupling had to be considered.
A new method of solving this important drawback is proposed in this work in order to optimize the scheme and reduce the memory requirements to one term per field component. The solution is based on the fact that stub voltages remain at the same node, which allows for the calculation of this term in a unique and global form. To do so, let us remember that for any capacitive stub, the incident voltage at the n-th time step is equivalent to the reflected voltage at the previous time step, at the same node. Since the line of relative admittance Y xk connects to voltage n V k , time step n results in the following equation: Electronics 2021, 10, 2071 10 of 22 Substituting (21) into n S e x and using (19) for time step n − 1, the term due to the stubs can be expressed as: which reduces the effect of the three stubs to a single accumulator parameter that only depends on the previous values of this parameter and the voltage.
In order to describe the information contained in the three shunt nodes, or equivalently, the three voltage components, the following compact form involving the column vector and tensor quantities can be obtained as follows: where each link line contributes to n → L e according to the electric field component it defines, that is to say: and the rest of vector and matrix terms are evident. The situation defining the magnetic field components is analogous, but the series circuits must now be used. The Thevenin equivalent circuit for the series node to define the x-oriented current at the n-th time step is shown in Figure 5. The lines directly associated with n i x are represented against a white background. As regards the off-diagonal elements of the permeability and magnetic conductivity tensors, an inductive stub and a magnetic loss line are shown against a blue background, lines 15xy and 22xy, to describe the coupling terms due to µ xy and σ m xy . Analogous lines 15xz and 22xz model the effect of µ xz and σ m xz , represented against a red background. The connection to the lines associated with n i x are again formal, since the current flowing through these lines is not n i x , but n i y and n i z . The situation defining the magnetic field components is analogous, but the series circuits must now be used. The Thevenin equivalent circuit for the series node to define the x-oriented current at the n-th time step is shown in Figure 5. The lines directly associated with nix are represented against a white background. As regards the off-diagonal elements of the permeability and magnetic conductivity tensors, an inductive stub and a magnetic loss line are shown against a blue background, lines 15xy and 22xy, to describe the coupling terms due to μxy and m xy σ . Analogous lines 15xz and 22xz model the effect of μxz and m xz σ , represented against a red background. The connection to the lines associated with nix are again formal, since the current flowing through these lines is not nix, but niy and niz. Figure 5. Thevenin equivalent circuit of the series node for Hx at time tn. Green boxes denote a formal connection, allowing for the current to be nix at the link lines, the source, and lines 16xx and 22xx¸while the current is niy for lines16xy and 22xy (in the pale blue region) and niz for lines16xz and 22xz (in the pale red region). Now, imposing the Kirchhoff's voltage law on the circuit in Figure 5 and rearranging the terms, the following is produced: where, Figure 5. Thevenin equivalent circuit of the series node for H x at time t n . Green boxes denote a formal connection, allowing for the current to be n i x at the link lines, the source, and lines 16xx and 22xx, while the current is n i y for lines16xy and 22xy (in the pale blue region) and n i z for lines16xz and 22xz (in the pale red region). Now, imposing the Kirchhoff's voltage law on the circuit in Figure 5 and rearranging the terms, the following is produced: where, As with the shunt circuits, the first term describes the contribution of the link lines, the second term is due signal stored at the inductive stubs, the third term stands for the source, and the right side of (26) describes the coupling between the different magnetic field components. The sign with which each link line contributes to n L m x depends on the voltage polarity of each line: according to Figure 5, lines 5 and 8 generate a current oriented along n i x , while lines 4 and 7 have the opposite effect.
As with the shunt circuit for E x , the need to store incident voltages at the stubs in (27) implies an undesired increase in memory storage and time-computing resources. An optimized solution for this inconvenience is also proposed for this circuit. To do so, let us proceed as with the shunt circuits. Since an inductive stub is a short-circuited stub, the incident field is related with the reflected voltage for the same line at the previous time step. Concretely and using the scheme of Figure 3, the following formula is produced: Substituting (28) into n S m x and using (26) for time step n − 1, the term due to the stubs can be reduced to the following: which means that an optimized scheme with only one storage term per field component is derived and makes that memory storage and time computing requirements for modeling anisotropic media are almost identical to those of the original TLM method when applied to isotropic media. Considering the three series circuits for the magnetic field, a compact form using the column vector and tensor quantities may be derived with the following: where the vector and matrix terms in (30) are evident, except for n → L m , in which the contribution of a certain link line depends on the current or magnetic field component it defines and its voltage polarity, which defines the appropriated sign for this current. Namely, for this case: Once the voltage and current components at the node are obtained through (24) and (30), the reflected pulse at each link line must be obtained. This can be performed by considering two opposite lines and the common voltage and current they define at the n-th time step, n V c and n i c . The situation is shown in Figure 6. The box at the center represents the effect of other lines coupled with this pair by means of the shunt and series circuits involved with these two opposite lines.
Once the voltage and current components at the node are obtained through (24) and (30), the reflected pulse at each link line must be obtained. This can be performed by considering two opposite lines and the common voltage and current they define at the n-th time step, nVc and nic. The situation is shown in Figure 6. The box at the center represents the effect of other lines coupled with this pair by means of the shunt and series circuits involved with these two opposite lines. Figure 6. Circuit involving two opposite link lines associated with a common voltage, nVc, and a common current, nic. The box at the center denotes a formal coupling with other lines also present in the shunt node associated with nVc, and the series node associated with nic.
The voltage on the line contributing positively to nic is denoted as nVpos, and the voltage on the opposite line is denoted as nVneg. It can be seen that the reflected pulses can be obtained using the following simple equations [27]: As an example, this expression applied to link lines 1 and 12 yields the following: Figure 6. Circuit involving two opposite link lines associated with a common voltage, n V c , and a common current, n i c . The box at the center denotes a formal coupling with other lines also present in the shunt node associated with n V c , and the series node associated with n i c .
The voltage on the line contributing positively to n i c is denoted as n V pos , and the voltage on the opposite line is denoted as n V neg . It can be seen that the reflected pulses can be obtained using the following simple equations [27]: As an example, this expression applied to link lines 1 and 12 yields the following: The rest of the link lines can be treated in a similar way. It is worth noting that the optimization scheme proposed in this work makes the voltage pulses reflected at the stubs unnecessary, since their global contribution per field component is considered through the electric and magnetic accumulators, n S e x and n S m x , respectively, which are updated, as indicated by Equations (23) and (30).
Summarizing, the procedure is as follows:

At the variable initiation stage:
The node size is chosen according to the desired maximum valid frequency of the results; The time step, together with the characteristic admittance or impedance of each line in the node, are calculated using (12), (13), (14), and (15); Matrices A e and A m , together with their corresponding inverse matrices, are obtained from the last equation in (20) and (27); Incident voltages on the link lines are set to zero; Initial voltage and current vectors, together with the terms due signals previously stored at the stubs, are set to zero: 0 ; The number of time step calculations, n max , is calculated according to the desired frequency resolution, ∆f =(∆t n max ) −1 .

2.
The time-marching process is started. For the calculation at time step n: The column vector terms due to the incident pulses at link lines, n i , respectively. The output fields are obtained and stored from these vectors using Equation (9) or its equivalent; The reflected pulses at the link lines are obtained using (33); Reflected pulses are transformed into incident pulses for the next time step calculation; Finally, index n is updated to n + 1, and the process is repeated until n reaches the final time step, nmax.
As mentioned before, the procedure is a reformulation intended to combine the advantages of two existing solutions: (i) The conceptual clarity provided by the individual treatment of each term in Maxwell's equations in a way that is similar to that present in the original TLM [25], and (ii) the optimization of memory storage and time calculation provided by [22]. The first aspect may help to provide a better understanding of the method, its dissemination, and the subsequent improvement for researchers more familiar with the initial TLM works, while the second aspect avoids the important drawback concerning memory and time-calculation requirements.
As a final comment, the procedure presented here does not include magneto-electric coupling, thus resulting in two independent calculations-one for voltage calculation and the other for current calculation-which involve matrix and vector quantities with three components. The inclusion of a more general situation in which magneto-electric coupling is present would be qualitatively similar, but all the six field components would be simultaneously required for each calculation. As a consequence, for instance, the field summations in (19) and (26) should be evaluated for the 6-field components, using only one 6 × 6 element matrix, A, instead of two independent 3 × 3 matrices, A e and A m .

Numerical Results
The TLM algorithm presented in the previous sections was applied to different situations to test its capability in successfully modeling non-dispersive anisotropic media.
The first problem considered the reflection of a plane wave, which, when propagating through vacuum in the x direction, reaches a carbon-fiber composite material. This is an anisotropic conducting material used for aircraft construction. The situation is sketched in Figure 7. The three slabs had identical widths, d = 3.75 mm, were normally oriented to the x-axis, and were supposed to be of infinite extension in the y and z directions, which turned it into a one-dimensional (1D) situation. The three slabs presented identical and isotropic values for the dielectric permittivity, 43ε 0, and magnetic permeability, µ 0 . The magnetic conductivity tensor was zero. As regards the electrical conductivity, it increased along the carbon fiber directions and was almost negligible for normal directions. The fibers were contained in the x plane in all cases, but they were oriented differently for each slab: along the y-axis, at 45 • from the yand the z-axis, and along the z-axis, for the first, second, and third slabs, respectively. The only non-zero elements of the electric conductivity tensor at each layer were σ e1 yy = 12 S/m, σ e2 yy = σ e2 yz = σ e2 zy = σ e2 yy = 8 S/m, and σ e3 zz = 12 S/m. The slab was surrounded by infinite vacuum regions at both ends. A z-polarized incident wave was excited and the reflected fields, E z and E y , were calculated to obtain the corresponding reflection coefficients for the layered material, Γ zz and Γ zy . Nodes of equal length of 93.75 µm for the three Cartesian directions were used, which means that each slab width was 40 nodes. This node size ensured proper modeling up to frequencies slightly below 50 GHz [11]. The maximum allowable time step was fixed by the vacuum regions, ∆t = ∆l/(2c) = 156.36 fs. fibers were contained in the x plane in all cases, but they were oriented differently for each slab: along the y-axis, at 45° from the y-and the z-axis, and along the z-axis, for the first, second, and third slabs, respectively. The only non-zero elements of the electric conductivity tensor at each layer were 1 12 S/m, e yy σ = The slab was surrounded by infinite vacuum regions at both ends. A z-polarized incident wave was excited and the reflected fields, Ez and Ey, were calculated to obtain the corresponding reflection coefficients for the layered material, zz Γ and . zy Γ Nodes of equal length of 93.75 µm for the three Cartesian directions were used, which means that each slab width was 40 nodes. This node size ensured proper modeling up to frequencies slightly below 50 GHz [11]. The maximum allowable time step was fixed by the vacuum regions, 5 / ( 1 6.36 2 ) t l c Δ = Δ = fs. The plane wave situation allowed for imposing symmetry conditions at the y and z limits, thus only one node was used for those two directions. As regards the x direction, the medium was truncated at distances dl and dr, large enough so as to separate and better identify the incident from the reflected fields. In our case, dl = 2000 and dr =380. The global mesh size was of 2500 × 1 × 1 nodes. Absorbing boundary conditions were imposed at The plane wave situation allowed for imposing symmetry conditions at the y and z limits, thus only one node was used for those two directions. As regards the x direction, the medium was truncated at distances d l and d r , large enough so as to separate and better identify the incident from the reflected fields. In our case, d l = 2000 and d r = 380. The global mesh size was of 2500 × 1 × 1 nodes. Absorbing boundary conditions were imposed at both limits in the x direction by simply connecting the link lines reaching the boundaries, x = 0 and x = 2500, to a discrete load with the vacuum impedance value. This is a remarkable simple condition available for the TLM method, which shows perfect absorption for normal incidence.
A Gaussian shaped electric current along the z direction was used to excite a zpolarized electric field with spectral information up to around 10 GHz: with g = 15.34 × 10 9 s −1 , t m = 0.3127 ns, and the value of I 0 was chosen to obtain a unitary maximum amplitude for the generated electric field. A total of 2 15 = 32768 time steps were calculated, which produced a frequency resolution of 0.195 GHz. Figure 8 is a plot of the time evolution of the incident of the z-polarized electric field and the two relevant components of the reflected electric field (calculated 10 nodes away from the left numerical limit). The y-component of the reflected field had a significantly lower amplitude than the z-component; hence, it was scaled by a factor of 10. both limits in the x direction by simply connecting the link lines reaching the boundaries, x = 0 and x = 2500, to a discrete load with the vacuum impedance value. This is a remarkable simple condition available for the TLM method, which shows perfect absorption for normal incidence.
A Gaussian shaped electric current along the z direction was used to excite a z-polarized electric field with spectral information up to around 10 GHz: with g = 15.34 × 10 9 s −1 , tm = 0.3127 ns, and the value of I0 was chosen to obtain a unitary maximum amplitude for the generated electric field. A total of 2 15 = 32768 time steps were calculated, which produced a frequency resolution of 0.195 GHz. Figure 8 is a plot of the time evolution of the incident of the z-polarized electric field and the two relevant components of the reflected electric field (calculated 10 nodes away from the left numerical limit). The y-component of the reflected field had a significantly lower amplitude than the z-component; hence, it was scaled by a factor of 10.  Figure 9 shows the magnitude of the reflection coefficients versus frequency. An independent FDTD solution presented in [28] was included for comparison. A good qualitative agreement was clearly seen between both independent results, were almost coinci- Figure 8. Incident and reflected electric fields at a three-layered carbon-fiber material. The y-oriented reflected field was magnified by a factor of 10. Figure 9 shows the magnitude of the reflection coefficients versus frequency. An independent FDTD solution presented in [28] was included for comparison. A good qualitative agreement was clearly seen between both independent results, were almost coincidental, which validated the presented algorithm. Although the proposed TLM results and the FDTD ones are almost indistinguishable in Figure 9, a quantitative comparison was carried out by evaluating the relative difference between both reflection coefficient magnitudes, defined as the difference between the magnitude of both reflection coefficients normalized to the FDTD magnitude used as reference. This relative error versus frequency is plotted in Figure 10, which shows differences below 1%, thus corroborating the validity of the TLM results. The increase in the difference at higher frequencies is explained by the different dispersion behavior of the TLM and FDTD methods [11]. The second example is not an actual medium with anisotropic behavior at certain properties, but a hypothetical one in which a more difficult case is considered since all the Although the proposed TLM results and the FDTD ones are almost indistinguishable in Figure 9, a quantitative comparison was carried out by evaluating the relative difference between both reflection coefficient magnitudes, defined as the difference between the magnitude of both reflection coefficients normalized to the FDTD magnitude used as reference. This relative error versus frequency is plotted in Figure 10, which shows differences below 1%, thus corroborating the validity of the TLM results. The increase in the difference at higher frequencies is explained by the different dispersion behavior of the TLM and FDTD methods [11]. Although the proposed TLM results and the FDTD ones are almost indistinguishable in Figure 9, a quantitative comparison was carried out by evaluating the relative difference between both reflection coefficient magnitudes, defined as the difference between the magnitude of both reflection coefficients normalized to the FDTD magnitude used as reference. This relative error versus frequency is plotted in Figure 10, which shows differences below 1%, thus corroborating the validity of the TLM results. The increase in the difference at higher frequencies is explained by the different dispersion behavior of the TLM and FDTD methods [11]. The second example is not an actual medium with anisotropic behavior at certain properties, but a hypothetical one in which a more difficult case is considered since all the properties are simultaneously anisotropic. Two situations were considered, which only differed in the choice of co-ordinate axes. The first one considered Cartesian co-ordinate The second example is not an actual medium with anisotropic behavior at certain properties, but a hypothetical one in which a more difficult case is considered since all the properties are simultaneously anisotropic. Two situations were considered, which only differed in the choice of co-ordinate axes. The first one considered Cartesian coordinate axes that coincided with the principal directions, while the second one used rotated co-ordinate axes, which would require using non-diagonal tensor properties. Since the situation was only a mathematical artifact, both results had to coincide and would thus serve as a test for the TLM method solution proposed in this work.
For the co-ordinate axes along the principal directions, the medium had vacuum properties, except for the following parameters: ε xx = 4ε 0 , σ e xx = 2.65 · 10 −5 S/m, µ zz = 4µ 0 , and σ m zz = 3.763 Ω/m. The values of both conductivities were chosen to meet the condition σ e xx /ε 0 = σ m zz /µ 0 . This relation between the electric and the magnetic conductivity ensured propagation without frequency dispersion, which allowed for the prediction of a common attenuation factor for all frequencies in the case of plane wave propagation [2,29]. In this example, the specific values of both conductivity terms had been chosen for an attenuation factor of exp(−1) after travelling 100 nodes. It is worth noting that this attenuation factor is only approximated since the medium considered here has permittivity and permeability values which differ from the vacuum ones, the medium was not isotropic, and the wave was not a plane wave but an elliptical one.
A z-oriented magnetic field was excited inside this medium. With the aid of symmetry conditions for the z direction, the medium could be considered as two-dimensional (2D) contained in the xy plane. As in the first example, discrete loads with the vacuum impedance were connected to the link lines reaching the boundaries at the x and y limits. This absorbing condition is not so efficient for arbitrary incidence but this fact is not relevant for the example considered, since the mesh is large enough to avoid pulses reflected at the boundaries.
A mesh of 600 × 600 × 1 TLM nodes was used, with a uniform size of 1 m in the three Cartesian directions. The maximum allowable time step for this node size was again used, which was ∆t = 1.6678 ns. Under these conditions, the E x , E y , and H z fields were expected to propagate with an elliptical front wave with a different phase velocity for the x and y directions. Namely: which means that a wave spends four time steps to travel a node distance along the x direction, while it takes eight time steps for the y direction. A relevant property of this medium is that the wavelength for propagation along the y direction at a given frequency is one half of the wavelength for the x direction at that frequency. Since the numerical results are expected to be valid for frequencies whose wavelength is sampled at least ten times, the node size and v y limit the validity of the results to frequencies approximately below f max =7.5 MHz. Bearing this in mind, a Gaussianshaped source similar to that in (35) was used to excite the H z field at the mesh center, with g = 14.99·10 6 s −1 and t m = 219.2 ns, in order to not excite frequencies above f max . Figure 11 is a plot of H z at time step 1000. The elliptical shape is evident, with the vertical axis length one-half of the horizontal axis length, as expected from the velocity values in (36). Disregarding the delay term t m at around 130 time steps, the field propagated a distance of around 250 nodes and 125 nodes along the x and y directions, respectively, as expected. In addition, the signal attenuated in the x direction due to the effect of σ m zz , but attenuation was higher for the y direction since the effect of σ m zz was combined to that of σ e xx . Electronics 2021, 10, x FOR PEER REVIEW 18 of 23   Figure 11: the source point, A, and two points, B and C, separated by 104 nodes along the x and y directions, with coordinates (300, 300, 1), (404, 300, 1), and (300, 404, 1), respectively. The field at the source point was scaled to one-fifth of its value to better appreciate the attenuated field at points B and C. The figure clearly shows the expected higher attenuation for the field propagating along the y direction. The time delay between the field at points A with respect to the field at points B and C is 419 and 828 time steps, respectively, which is in good agreement with the expected delays of 416 and 832 time steps, according to the distance travelled and the phase velocities shown in (36). Independent TLM results obtained with the method proposed in [22] are also represented to validate the method proposed in this work. A qualitative good agreement is clearly observed in the Figure.    Figure 11: the source point, A, and two points, B and C, separated by 104 nodes along the x and y directions, with coordinates (300, 300, 1), (404, 300, 1), and (300, 404, 1), respectively. The field at the source point was scaled to one-fifth of its value to better appreciate the attenuated field at points B and C. The figure clearly shows the expected higher attenuation for the field propagating along the y direction. The time delay between the field at points A with respect to the field at points B and C is 419 and 828 time steps, respectively, which is in good agreement with the expected delays of 416 and 832 time steps, according to the distance travelled and the phase velocities shown in (36). Independent TLM results obtained with the method proposed in [22] are also represented to validate the method proposed in this work. A qualitative good agreement is clearly observed in the Figure Figure 11: the source point, A, and two points, B and C, separated by 104 nodes along the x and y directions, with coordinates (300, 300, 1), (404, 300, 1), and (300, 404, 1), respectively. The field at the source point was scaled to one-fifth of its value to better appreciate the attenuated field at points B and C. The figure clearly shows the expected higher attenuation for the field propagating along the y direction. The time delay between the field at points A with respect to the field at points B and C is 419 and 828 time steps, respectively, which is in good agreement with the expected delays of 416 and 832 time steps, according to the distance travelled and the phase velocities shown in (36). Independent TLM results obtained with the method proposed in [22] are also represented to validate the method proposed in this work. A qualitative good agreement is clearly observed in the Figure.  The quantitative comparison of the TLM results proposed in this work with the TLM reference results is shown in Figure 13, where the difference between both solutions normalized to the reference solution is plotted. The difference is considerably lower than in Figure 10 since the proposed and the reference solutions are now based on the TLM method. The good agreement validates the obtained results when the principal axes coincide with the co-ordinate axes.
Electronics 2021, 10, x FOR PEER REVIEW 19 of 23 The quantitative comparison of the TLM results proposed in this work with the TLM reference results is shown in Figure 13, where the difference between both solutions normalized to the reference solution is plotted. The difference is considerably lower than in Figure 10 since the proposed and the reference solutions are now based on the TLM method. The good agreement validates the obtained results when the principal axes coincide with the co-ordinate axes. The capability of the algorithm to deal with non-principal axes was tested by rotating the medium around an axis, which passes through point A and is parallel to the z direction. Although this rotation considerably complicated the description because the off-diagonal elements of the properties had to be considered, the results should be clearly comparable and only be affected by this same rotation. A counter-clockwise rotation of 120° was chosen to test possible undesired effects of negative stubs. The medium parameters with respect to these axes were: For this case, and to test its versality in dealing with nodes of different sizes, the length of each node was chosen to be 1 m for the x and z directions, and 0.5 m for the y direction. The equivalent TLM mesh was then of 600 × 1200 × 1 nodes and the time step was one-half of the value used in the previous example. The maximum frequency of the valid results was slightly more difficult to predict due to the axes rotation but was still around or slightly below the previous value of 7.5 MHz. As before, symmetry conditions were applied in the z direction and the connection of link lines to load impedances at the x and y limits. After the rotation, points A, B, and C were approximately described by coordinates (300, 600, 1), (248, 780,1), and (210, 496, 1), respectively. Figure 14 shows Hz at time step 2000, equivalent to time step 1000 represented in Figure 11, for a counter-clockwise rotation of 120°. As expected, the results are qualitatively coherent with those sketched in Figure 11. The capability of the algorithm to deal with non-principal axes was tested by rotating the medium around an axis, which passes through point A and is parallel to the z direction. Although this rotation considerably complicated the description because the off-diagonal elements of the properties had to be considered, the results should be clearly comparable and only be affected by this same rotation. A counter-clockwise rotation of 120 • was chosen to test possible undesired effects of negative stubs. The medium parameters with respect to these axes were: ε xx = 1.75ε 0 , ε yy = 3.25ε 0 , ε xy = −1.299ε 0 , σ e xx = 6.662 · 10 −6 S/m, σ e yy = 19.875 · 10 −6 S/m, σ e xy = −11.475 · 10 −6 S/m, µ zz = 4µ 0 , σ m zz = 3.763 Ω/m, while the remaining parameters were symmetrical or corresponded to the vacuum values.
For this case, and to test its versality in dealing with nodes of different sizes, the length of each node was chosen to be 1 m for the x and z directions, and 0.5 m for the y direction. The equivalent TLM mesh was then of 600 × 1200 × 1 nodes and the time step was onehalf of the value used in the previous example. The maximum frequency of the valid results was slightly more difficult to predict due to the axes rotation but was still around or slightly below the previous value of 7.5 MHz. As before, symmetry conditions were applied in the z direction and the connection of link lines to load impedances at the x and y limits. After the rotation, points A, B, and C were approximately described by co-ordinates (300, 600, 1), (248, 780,1), and (210, 496, 1), respectively. Figure 14 shows H z at time step 2000, equivalent to time step 1000 represented in Figure 11, for a counter-clockwise rotation of 120 • . As expected, the results are qualitatively coherent with those sketched in Figure 11. Electronics 2021, 10, x FOR PEER REVIEW 20 of 23 Figure 14. Magnetic field Hz at time step 2000 for a medium with anisotropic electric permittivity and conductivity, and magnetic permittivity and conductivity. The principal axes were rotated by 120° with respect to the Cartesian ones. Node size is 1 m for the x and z directions and 0.5 m for the y direction. The time step used is one-half of that used for the example in Figure 11; therefore, both Figures represent the same instant.
The comparison of the magnetic field at points B and C for the case in Figure 11 (using uniform mesh and principal axes) with results in the case considered now (with rotated axis and non-uniform mesh) is shown in Figure 15. Both figures are almost indistinguishable from each other in the time domain.  Figures 11 and 14-obtained using a uniform mesh with principal axes oriented along the Cartesian axes, compared with results obtained using a nonuniform mesh with principal axes rotated by 120° with respect to the Cartesian ones.
The relative difference at points B and C versus frequency between the results in Figure 15 is shown in Figure 16. Each absolute difference was normalized to the magnitude of the solution obtained for the uniform mesh and for the Cartesian axes at its respective point. This plot is actually a good example with which to test the TLM limits fixed by the phase velocities and the node size; in our case, it was a maximum valid frequency of around 7.5 MHz. As expected, results show good agreement with the reference for low frequency values, but present an increasing degradation as the limit frequency is reached. Figure 14. Magnetic field H z at time step 2000 for a medium with anisotropic electric permittivity and conductivity, and magnetic permittivity and conductivity. The principal axes were rotated by 120 • with respect to the Cartesian ones. Node size is 1 m for the x and z directions and 0.5 m for the y direction. The time step used is one-half of that used for the example in Figure 11; therefore, both Figures represent the same instant.
The comparison of the magnetic field at points B and C for the case in Figure 11 (using uniform mesh and principal axes) with results in the case considered now (with rotated axis and non-uniform mesh) is shown in Figure 15 The comparison of the magnetic field at points B and C for the case in Figure 11 (using uniform mesh and principal axes) with results in the case considered now (with rotated axis and non-uniform mesh) is shown in Figure 15. Both figures are almost indistinguishable from each other in the time domain.  Figures 11 and 14-obtained using a uniform mesh with principal axes oriented along the Cartesian axes, compared with results obtained using a nonuniform mesh with principal axes rotated by 120° with respect to the Cartesian ones.
The relative difference at points B and C versus frequency between the results in Figure 15 is shown in Figure 16. Each absolute difference was normalized to the magnitude of the solution obtained for the uniform mesh and for the Cartesian axes at its respective point. This plot is actually a good example with which to test the TLM limits fixed by the phase velocities and the node size; in our case, it was a maximum valid frequency of around 7.5 MHz. As expected, results show good agreement with the reference for low frequency values, but present an increasing degradation as the limit frequency is reached.  Figures 11 and 14-obtained using a uniform mesh with principal axes oriented along the Cartesian axes, compared with results obtained using a non-uniform mesh with principal axes rotated by 120 • with respect to the Cartesian ones.
The relative difference at points B and C versus frequency between the results in Figure 15 is shown in Figure 16. Each absolute difference was normalized to the magnitude of the solution obtained for the uniform mesh and for the Cartesian axes at its respective point. This plot is actually a good example with which to test the TLM limits fixed by the phase velocities and the node size; in our case, it was a maximum valid frequency of around 7.5 MHz. As expected, results show good agreement with the reference for low frequency values, but present an increasing degradation as the limit frequency is reached. Electronics 2021, 10, x FOR PEER REVIEW 21 of 23 Figure 16. Relative difference in Hz between the TLM results proposed in this work using a uniform mesh with principal axes oriented along the Cartesian axes, compared with results obtained using a non-uniform mesh with principal axes rotated by 120° with respect to the Cartesian ones. The difference was normalized to the solution obtained for principal axes oriented along the Cartesian axes.
The following is a summary of the results of this section: • A first example of an actual medium exhibiting anisotropic behavior of the electrical conductivity shows that the proposed algorithm provides results in good agreement with independent FDTD results; • A second example considers a more challenging case with simultaneous anisotropic contributions of all the electromagnetic properties under study: o First, a mesh with nodes of equal length for all the Cartesian co-ordinate direction was tested when these co-ordinate axes were oriented along the principal axes of the electromagnetic properties, therefore making all the tensor properties diagonal. The comparison with an independent TLM solution validated the proposed method; o The final test considered a more complicated situation in which the node had a different length for each direction and the principal axes were rotated with respect to the co-ordinate axes, thus requiring the use of non-diagonal tensors; o This last situation clearly showed the expected limits of the method: good performance for low frequencies, with poorer results as the frequencies were high enough, implying the need to take approximately less than ten samples per wavelength.

Conclusions
A new reformulation of the TLM method for modeling non-dispersive anisotropic media is presented in this work. The algorithm is proposed in a way that resembles the original TLM formulation, based on the assignment of a specific transmission line to each term in Maxwell's equations. Although this scheme may be considered by some authors as conceptually simpler than other approximations and may therefore help to provide a better understanding, use, and possibility for further development by other authors, it presents the important disadvantage of requiring a priori the storage of an individual voltage pulse for each element of the electromagnetic tensor quantities. This is not a minor drawback in the case of anisotropic media, in which an important number of off-diagonal tensor properties are to be considered. To solve this difficulty, an optimized treatment of the signal stored on the capacitive and inductive lines is proposed, which reduces the memory storage requirements to one term per field component. As a result, the algorithm Figure 16. Relative difference in H z between the TLM results proposed in this work using a uniform mesh with principal axes oriented along the Cartesian axes, compared with results obtained using a non-uniform mesh with principal axes rotated by 120 • with respect to the Cartesian ones. The difference was normalized to the solution obtained for principal axes oriented along the Cartesian axes.
The following is a summary of the results of this section:

•
A first example of an actual medium exhibiting anisotropic behavior of the electrical conductivity shows that the proposed algorithm provides results in good agreement with independent FDTD results; • A second example considers a more challenging case with simultaneous anisotropic contributions of all the electromagnetic properties under study: First, a mesh with nodes of equal length for all the Cartesian co-ordinate direction was tested when these co-ordinate axes were oriented along the principal axes of the electromagnetic properties, therefore making all the tensor properties diagonal. The comparison with an independent TLM solution validated the proposed method; The final test considered a more complicated situation in which the node had a different length for each direction and the principal axes were rotated with respect to the co-ordinate axes, thus requiring the use of non-diagonal tensors; This last situation clearly showed the expected limits of the method: good performance for low frequencies, with poorer results as the frequencies were high enough, implying the need to take approximately less than ten samples per wavelength.

Conclusions
A new reformulation of the TLM method for modeling non-dispersive anisotropic media is presented in this work. The algorithm is proposed in a way that resembles the original TLM formulation, based on the assignment of a specific transmission line to each term in Maxwell's equations. Although this scheme may be considered by some authors as conceptually simpler than other approximations and may therefore help to provide a better understanding, use, and possibility for further development by other authors, it presents the important disadvantage of requiring a priori the storage of an individual voltage pulse for each element of the electromagnetic tensor quantities. This is not a minor drawback in the case of anisotropic media, in which an important number of off-diagonal tensor properties are to be considered. To solve this difficulty, an optimized treatment of the signal stored on the capacitive and inductive lines is proposed, which reduces the memory storage requirements to one term per field component. As a result, the algorithm proposed combines and exploits the advantages of two previous existing solutions: a direct and conceptual approach close to the original TLM method, together with an optimized treatment of the signals stored on the stub lines, which reduces the memory storage and time-computing requirements of those used in the original TLM scheme for modeling isotropic media.
The good behavior of the proposed method is tested through a pair of examples involving non-dispersive anisotropic media. Their comparison with independent existing solutions proves the good behavior of the proposed method at the usual frequency limit for which the wavelength has been sampled at least ten times. Results in the time domain are virtually indistinguishable from the reference works, although some differences may be observed in the frequency domain, specially at frequencies near the expected limit of the method.
The extension to deal with frequency-dependent properties bridging the difficulties inherent to signal processing models seems to be the next problem to be considered and is left for a subsequent work.