Time-domain Formulation of a Multi-layer Plane Circuit Coupled with Lumped-parameter Circuits using Maxwell Equations

We calculate electromagnetic phenomena in the multi-layer plane circuit starting from the Maxwell equations. We present a numerical method of potential and current density in two-dimensional conductors, where their time developments are treated as phenomena of wave propagation. We treat the plane conductors by dividing them into small finite-volume elements, similar to the case of the partial element equivalent circuit method, and the transport equations are then solved by the finite-difference time-domain method. Furthermore, we develop a calculation method for the boundary in a multi-layer plane by applying the method we have used in multi-transmission lines. We formulate the boundary conditions of a multi-layer plane coupled with lumped-parameter circuits and introduce an algorithm to reduce calculation costs that are largely associated with the two-dimensional extension from the multi-transmission-line case. We perform calculations of the wave propagation of potential, current density, and charge density in the time domain for a simple plane circuit. These calculations are presented as supplementary materials of the present paper.


time-domain formulation of a Multi-layer plane circuit coupled with Lumped-parameter circuits using Maxwell equations Souma Jinno * , Shuji Kitora, Hiroshi toki & Masayuki Abe
We calculate electromagnetic phenomena in the multi-layer plane circuit starting from the Maxwell equations. We present a numerical method of potential and current density in two-dimensional conductors, where their time developments are treated as phenomena of wave propagation. We treat the plane conductors by dividing them into small finite-volume elements, similar to the case of the partial element equivalent circuit method, and the transport equations are then solved by the finitedifference time-domain method. Furthermore, we develop a calculation method for the boundary in a multi-layer plane by applying the method we have used in multi-transmission lines. We formulate the boundary conditions of a multi-layer plane coupled with lumped-parameter circuits and introduce an algorithm to reduce calculation costs that are largely associated with the two-dimensional extension from the multi-transmission-line case. We perform calculations of the wave propagation of potential, current density, and charge density in the time domain for a simple plane circuit. these calculations are presented as supplementary materials of the present paper.
Electromagnetic noise is troublesome because it is difficult to understand the fundamental phenomena associated with noise. In recent years, circuit design that reduces the influence of electromagnetic noise has become increasingly important for reliable performance with the trend toward high-frequency and low-voltage electronics 1 . This trend is likely to result in electromagnetic noise, causing the malfunction of electric devices, and can be a major bottleneck in future circuit design 2 . There are many symptomatic treatments, such as filters, to suppress the electromagnetic noise 3,4 . However, these treatments require the experience of skilled engineers because the causes of electromagnetic noise phenomena are not clearly understood and cannot be described quantitatively using equivalent circuit modeling and impedance analysis in the frequency domain 5,6 . A fundamental understanding of the mechanisms of electromagnetic noise is needed for noiseless circuit design.
In previous studies, we elucidated the mechanism of electromagnetic noise phenomena in line circuit systems using the multi-conductor transmission line (MTL) theory derived from the Maxwell equations 7 . We calculated common-mode noise for a configuration of three-transmission-line theory and developed a noiseless circuit design 8,9 . We have also proposed a numerical calculation algorithm for hybrid systems, where the distributed-and lumped-parameter circuits are connected 10,11 .
For further analysis of electromagnetic noise in an actual circuit, such as a printed circuit board (PCB), we have extended the MTL theory to two-and three-dimensional conductors and formulated transport equations of potential and current. In the numerical calculation, we realized calculation of a multi-layer plane (MLP) circuit with a discontinuous ground using the finite-difference time-domain (FDTD) method 12 . Although we applied our algorithm to an MLP circuit to connect lumped-parameter circuits at boundaries, there are still major problems to be solved. First, we have to consider the directions of current according to the locations of the boundaries in the MLP. Second, the number of connection points of the lumped-parameter circuit greatly increases as compared to the case of the MTL theory. For these reasons, applying the previous calculation algorithm to arbitrary shapes of MLPs is difficult due to complicated boundary conditions and the associated heavy computer burden.
In the present paper, we would like to extend the MTL theory to the case of two-and three-dimensional conductors. To this end, we derive the transport equations in the form of finite-volume elements of potential and current density from the Maxwell equations. Using these equations, we are able to visualize the transitional phenomena as a function of time, which is important in order to understand the electromagnetic noise phenomena. In addition, we would like to resolve problems of complex calculations at boundaries of MLPs. For this purpose, we formulate the boundary equations to consider any direction of current at the boundaries. We use Kirchhoff 's current law (KCL) to consider the connective relation between lumped-and distributed-parameter circuits. We also use Kirchhoff 's voltage law (KVL) and branch constitutive equation (BCE) to consider the lumped-parameter circuit conditions.
In the present paper, a calculation algorithm for a two-dimensional conductor is presented in order to clarify the mechanism of electromagnetic noise occurring in the MLP circuit. We assume that two-or three-dimensional MLP circuit systems are connected with lumped-parameter circuits at the boundaries (an example is shown in Fig. 1). The proposed method provides information of transient phenomena of the potential and current density in the planar circuit system. In Section 2, two-and three-dimensional telegraphic equations expressed in terms of potential and current density are derived. Then, coupled differential equations are expressed in the form of discretized finite-volume elements for numerical calculation. In Section 3, we introduce the numerical method at the boundary. In Section 4, we describe how to reduce boundary calculation costs by considering the fact that most of the boundary of the MLP is open, by which the calculation cost of the boundary conditions can be significantly reduced.

Derivation of Coupled Difference Equations of Potential and Current Density in MLP circuits
In previous studies, we derived equations of potential and current in a multi-transmission-line (MTL) system from the Maxwell equations 7 . We have reported discretized equations in an MLP circuit system and compared the results of simulations and experiments 12 . In the present paper, we discuss the derivation of discretized equations in detail in order to further understand the calculations. We start with the retarded scalar U t r ( , ) and vector A t r ( , ) potentials obtained by applying the Lorenz gauge from the Maxwell equations.  where q t r ( , ) and j t r ( , ) are the charge and current densities in the conductors. We assume that, in addition to U t r ( , ) and A t r ( , ), both Ohm's law and charge conservation hold in the planar conductors.
Here, ρ denotes the resistance in the plane conductors. These four equations satisfy the necessary and sufficient conditions for solving the scalar potential, vector potential, charge density, and current density in conductors. Here, we realize a numerical instability in the presence of the delay time in Eqs.. (1) and (2). This stability problem is a large challenge and a long-standing problem for circuit theories such as the partial element equivalent circuit (PEEC) method 13 . At present, we avoid this numerical difficulty by neglecting the delay time in the coupled differential equations 14 . For numerical Figure 1. Example of a two-dimensional multi-conductor plane circuit system described in the present paper. We regard MLP circuits as distributed-parameter circuits. Lumped-parameter circuits, such as voltage source, resistance, capacitance, or inductance circuits, are connected by distributed-parameter circuits at the boundaries, which provide the boundary conditions. The potential U and current density j depend on the position r and time t.
calculation, we use finite-volume elements by taking an average in the small finite volume. In the present paper, we assume a multi-layer planar conductor circuit, as shown in Fig. 2, and we use a cubic finite-volume element to divide the conductor planes into finite-volume elements. We refer to one of the divided elements as cell i k l ( , ), where i is a conductor plane, k is the x coordinate, l is the y coordinate, and the mesh dimensions are Δx and Δy.
Here, we consider the distribution in two directions (x and y) because the thickness of the plane conductor, denoted Δz, is sufficiently thin compared with Δx and Δy. The volume of cell i k l ( , ) is represented by . We can express all variables using quantities of cell elements as: We apply the averaging procedure to cell i k l ( , ) for Eqs. (1) and (2): , , ( , ) ( , ) ( , ) Equations (1) and (2) can be expressed as the sum of cell elements. Here, P and L represent the local potential coefficient and the local inductance coefficient of the conductors: The local coefficients P and L have two sets of subscripts, which represent the relation between cells i k l ( , ) and i k l ( , ) ′ ′ ′ . Next, we derive spatial difference equations by taking the average of Eqs. (3) and (4) in cell i k l ( , ). Here, Eq. (3) includes the spatial differentiation, which is the inclination in the x, y, and z directions in cell i k l ( , ). For instance, spatial differentiation of any function f x y z ( , , ) in the x direction can be converted into the cell spatial difference using the following relation: i k l ( , ) 1 2 1 2 Figure 2. We introduce small finite-volume elements in the multi-plane conductors. The x y ( , ) position of the i-th plane is discretized as i k l ( , ), where k and l correspond to numbers in the x and y directions, respectively. The volume of each small element is V x y z = Δ Δ Δ , where Δz is the thickness of the plane and is considered to be very thin.
We express spatial differentiation with respect to x as a central difference shifted from x by half a cell, and take the average in the x direction. The average spatial differentiation in cell i k l ( , ) is converted into the cell spatial difference between cells, which is positively or negatively shifted by half a cell in the x direction from cell i k l ( , ). Furthermore, the divergence of the space is included in Eq. (4) and can be expressed using the inner product: Here, we do not consider the thin z direction of a rectangular cell for simplicity. The average divergence in cell i k l ( , ) is converted into the sum of the spatial difference in each direction. Finally, Ohm's law (3) and the continuous Eq. (4) can be given in terms of the cell spatial differences: , , , , Here, Eqs. (13) and (14) express the vector components of Eq. (3). The difference of the half integer of space between the potential and the current densities comes from Eqs. (11) and (12). We used the scalar potential and the current density as variables in circuit theory and eliminated the vector potential and charge density cell elements from Eqs. (7), (8), (13), (14) and (15). Finally, the transport equations are expressed in the form of recurrence formulas by discretizing the time direction: x y 1 2 x y 1 2  (18), respectively. Moreover, if the current densities can be calculated, we can also quantify the vector potential and charge density using Eqs. (8) and (15), respectively. The suffix m at the upper right represents the discretized time. The difference of the half integer of time between the potential and current densities comes from the FDTD method. Here, N represents the number of plane, and N x and N y represent the numbers of divisions in the x and y directions, respectively. The calculation of boundaries expressed by gray cells in Fig. 3 is explained in the next section.

Method of Boundary calculation of MLp with a Lumped-parameter circuit
In this section, we explain the method of connecting the two-dimensional distributed-parameter circuit and the lumped-parameter circuit at any boundary position. We use the node potential and branch current in the lumped-parameter circuit and the cell potential and cell current density in the distributed-parameter circuit as unknown variables. In order to treat these circuits, we use KVL and the BCE in the lumped-parameter circuit and Eqs. (16), (17) and (18), which are cell constitutive equations (CCE), in the distributed-parameter circuit. Assuming that the node potential is equivalent to the cell potential and branch current is equivalent to the cell current, we can express current conservation by KCL, which is the connection relationship between the lumpedand distributed-parameter circuits, as shown in Fig. 4.
We formulated the time-domain incident-potential equation (IPE), which is a coupled equation of a lumped-parameter circuit and an MTL, in a previous study 11 . As for the lumped-parameter circuit, we can use time-domain IPE by using the incidence matrix of directed graph A and time-domain impedance matrix Z:   . Connection between the lumped-and distributed-parameter circuits at boundaries. In order to use KCL, we assume that the cell potential is equivalent to the node potential and that the cell current is equivalent to the branch current. In the Supplementary information, we present the time variations of potentials and charge densities as numerical calculation results (see Supplementary Information, Movies 1 and 2). The upper plane has a narrower width of W 1 , and the bottom plane has a wider width of W 2 . The length of both planes is L, and the distance between them is H. The independent voltage source V S , which has internal resistance R S , is connected at the left side (a-a'), and approximate impedance matching is realized using an appropriate value of R L at the right side (b-b'). (2019) 9:17891 | https://doi.org/10.1038/s41598-019-53288-x www.nature.com/scientificreports www.nature.com/scientificreports/ connections between nodes and branches, and its entries are 0, 1, or −1. If the component of the α-th row and the β-th column is 1, then the α-th node is connected to the β-th branch, and the current direction is outward from the node. If the component of the α-th row and the β-th column is −1, then the current direction is inward to the node. If the component of the α-th row and the β-th column is 0, then the nodes and branches are not connected. In addition, Z l is a time-domain impedance matrix of size × p p, which has diagonal elements depending on β-th lumped-parameter circuit elements. From a previous study 10,11 , the time-domain impedance elements are derived as follows. The time-domain impedance elements for the resistor, capacitor, inductor, and independent voltage source are R, t C /(2 ) Δ , L t 2 /( ) Δ , and 0, respectively. Vector V s is for an independent voltage source, and vector I s is for an independent current source. Then,  and δ are × p p diagonal matrices, which are used to change signs depending on the β-th lumped-parameter circuit elements: for other elements (20) Next, we derive the time-domain IPE in the same form as Eq. (19) for the distributed-parameter circuit. We consider boundary equations of lumped-parameter circuits and the MLP circuit at boundaries A through D shown in Fig. 5. Since the size of the current density cell is halved at these boundaries and results in a complicated numerical calculation of the impedance matrix, we replace the boundary current density with the time average of the current density flowing in a potential cell i k l ( , ) at the boundary, which is expressed as follows:   www.nature.com/scientificreports www.nature.com/scientificreports/ We use the recurrence formula of the potential in Eq. (16) as a CCE for the boundary equation. Unknown variables in the distributed-parameter circuit are the boundary potential and the current density flowing in the potential cell and are extracted from the summation in Eq. (16), which is brought to the left-hand side: 1 2 x y We introduced factors xi k l  Table 1 to consider the shift in the current density by half in Eq. We can express the CCE as a matrix equation, as used in Eq. (19): Here, matrix A d is a unit matrix because each cell is not connected to other cells at boundaries. Vector U d and j d are the cell potential and the cell current density, respectively. Moreover, Z d is the time-domain cell impedance matrix expressed as . We write ∼ U for all terms in the second, third, and fourth summation terms in the right-hand side of Eq. (22) and treat this term as a source term in the CCE.
Next, we can use KCL to derive the connection relation between distributed-and lumped-parameter circuits using the incidence matrix: ). Finally, we can express the time domain IPE for the coupled lumped-and distributed-parameter circuit by combining Eqs. (19), (23) and (24). We express the variables of the boundary, which are U l , U d , I l , and j d , as follows: The incidence matrices of lumped-and distributed-parameter circuits are joined horizontally: The time-domain impedance matrices in both lumped-and distributed-parameter circuits are arranged diagonally: Here, the matrices  and δ are diagonal matrices with the signs chosen according to the β-th element of the lumped-parameter circuit and the distributed-parameter circuit, as in Eqs. (30)  The matrix γ S is also a diagonal matrix to fit the unit of current in KCL and consider the direction of current in the β-th distributed-parameter circuit.

calculation Algorithm of Boundary conditions to Reduce calculation cost
In this section, we explain the appropriate calculation algorithm for an MLP circuit. In the previous algorithm shown in Table 2, we calculated all boundaries using Eq. (29), although the size of the matrix is proportional to the number of cells in the MLP circuit. The calculation is burdensome in the case of a complicated circuit. In the new algorithm, we divide the boundary into two cases, considering that most of the boundaries in the MLP circuit are open. One is the case with a lumped-parameter circuit (W), and the other is the case without a lumped-parameter circuit (WO), as shown in Table 2. We can reduce the size of the matrix in Eq. (29) considerably and obtain a simplified boundary equation in the case of no lumped-parameter circuit.
In the case without a lumped-parameter circuit, the KCL at the boundary is derived from Eq. (29): The current density in the distributed-parameter circuit flowing perpendicular to the boundary is always zero when a lumped-parameter circuit is not connected. As for the condition of the distributed-parameter circuit, the following equation can be obtained by substituting Eq. (33) into Eq. (22): ( )

Discussion with numerical Results
We performed numerical calculation using the present method and the corresponding experiment. We then compared the numerical results with the experimental results together with the results calculated by the PEEC method 15,16 . We used the equations and space discretization methods used in the PEEC method. There are differences in the calculation method of boundary conditions and the updating time step between our formulation and the PEEC method. We calculated the boundaries by connecting the lumped-parameter circuit and the transport equations derived from the Maxwell equations, as discussed in the present paper, while the PEEC method converts all finite-volume cells into equivalent circuit models and solves both the boundary and no-boundary parts at the same time. As for time discretization, we used the leapfrog time discretization method, which is the concept of the FDTD method, for updating the potential and the current density 17 , while the PEEC method uses a lumped-parameter circuit simulation algorithm 18,19 .
We consider a two-layer plane circuit with two bends used in the simulation and experiment (see Supplementary Information, Figs 1 and 2). We compare our simulation results with sub-nanosecond time domain reflectometry (TDR) experiments and PEEC method simulation. We can see good agreement between the numerical results of the proposed method and the experimental results (see Supplementary Information,  Fig. 3). We perform simulation with both the proposed method and the PEEC method with the same numerical calculation parameters, such as the space and time steps. The proposed method can reproduce the reflected wave with higher accuracy because of the calculation method for boundary conditions and leapfrog time updating.
In addition, we present electromagnetic noise simulation in a simple circuit, which has a micro-strip line configuration, as shown in Fig. 4. The ground noise voltage, which causes external radiation and common mode noise, is generated in various directions according to the position of the signal line 20 (see Supplementary  Information, Movies 1 and 2).
At the moment, the delay time in the scalar and vector potentials in the charge and current densities has been dropped in the numerical calculations. We are currently studying the inclusion of the delay time effect in the present formulation 14 . In addition, multi-layer-plane conductors exist in a uniform medium, and the skin effect is not considered. In order to precisely consider the change of dielectric properties in the conductors and surrounding materials as well as the skin effect, it is necessary to extend our calculation to three dimensions. If three-dimensional calculation is realized, then the propagation speed will change inside and outside of the conductor plane, which means that we can consider the variation of dielectric properties in each plane, as well as the skin effect. We hope to consider these factors in the near future.
This calculation method enables us to visualize transient basic electromagnetic phenomena, such as propagation, reflection, and coupling with other conductors. We can also clarify the relationship between complex circuit shapes and electromagnetic phenomena observed in high-frequency filters, antennas, and meta-materials [21][22][23] . This calculation method will be important not only in elucidating electromagnetic noise phenomena but also in developing future advanced circuit design.

Summary
In the present paper, we introduced a method to quantify the transient phenomena of the potential and current density in a multi-layer plane (MLP) circuit. The theoretical equations were obtained from the Maxwell equations. In addition, we used Ohm's law and continuity equations in order to describe the transport phenomena. In the numerical method, we use the spatial difference method used in the PEEC method and the central difference approximation used in the FDTD method. We also developed a calculation method for the boundary in the MLP and proposed an appropriate calculation algorithm in order to greatly reduce the calculation cost. We applied the present theoretical method for a two-dimensional circuit system and demonstrated the power of the newly developed theory.