The time domain numerical method of three-dimensional conductors including radiation with lumped parameter circuit

We formulate a numerical method on the transmission and radiation theory of three-dimensional conductors starting from the Maxwell equations in the time domain. We include the delay effect in the integral equations for the scalar and vector potentials rigorously, which is vital to obtain numerically stable solutions for transmission and radiation phenomena in conductors. We provide a formalism to connect the conductors to any passive lumped-parameter circuits. We show one example of numerical calculations, demonstrating that the new formalism provides stable solutions to the transmission and radiation phenomena.

The time domain numerical method of three-dimensional conductors including radiation with lumped parameter circuit

Souma Jinno * , Shuji Kitora, Hiroshi Toki & Masayuki Abe
We formulate a numerical method on the transmission and radiation theory of three-dimensional conductors starting from the Maxwell equations in the time domain. We include the delay effect in the integral equations for the scalar and vector potentials rigorously, which is vital to obtain numerically stable solutions for transmission and radiation phenomena in conductors. We provide a formalism to connect the conductors to any passive lumped-parameter circuits. We show one example of numerical calculations, demonstrating that the new formalism provides stable solutions to the transmission and radiation phenomena.
In the circuit design with a desire to reduce wiring widths and gaps, the shape and size of the circuit conductor significantly influence the electromagnetic performance as the signal frequency increases [1][2][3] . Non-uniform current distribution caused by the proximity effect and propagation in a three-dimensional conductor creates the common-mode current 4,5 .
Furthermore, the time variation of potential becomes a source of ground bounce and simultaneous switching noise [6][7][8] . These electromagnetic effects due to circuit geometry create various noise phenomena that can lead to miss-design of circuit performance and the generation of unexpected electromagnetic noise. For required circuit performance, the circuit design in various fields requires a robust numerical method for rigorous full-wave analysis of three-dimensional circuit structures.
The Finite-Difference Time-Domain (FDTD) method, which directly solves the Maxwell equations for electromagnetic (EM) fields in time and space, is one of the most frequently used full-wave analysis tools [9][10][11] . Although the FDTD method with the EM fields is suitable for analyzing the electromagnetic field caused by radiation from a circuit, it cannot precisely treat the signals flowing in the conductors connected with the lumped-parameter circuit, where potential and current are used for the time development. We have to develop a robust algorithm to treat lumped-parameter circuits connected to conductors with any shape.
There are several papers to conduct full-wave analysis, which start from the Maxwell equations in one-dimensional transmission line systems using the integral-differential equations in terms of potential and current used in the circuit theory [12][13][14][15][16] . There are some differences among these studies, but essentially all works emphasize that the integral-differential equations include the Heaviside equation. However, these studies did not discuss the difficulties in solving the coupled differential-integral equations numerically in the time domain with the delay term.
One attractive formalism was proposed in the name of the PEEC (Partial Element Equivalent Circuit) method, which can conduct a time-domain full-wave analysis in a three-dimensional circuit system by considering the delay time 17 . However, even in the PEEC formalism, the divergence in the numerical calculation with delay time has been a long-standing problem [18][19][20][21] . One possible cause of the numerical problem is the approximate treatment of the delay time in the Galerkin formalism with the pulse function. The delay time is approximated with a location-independent constant within finite volumes by using the center-to-center approximation 22 .
Recently, our group has successfully performed numerical calculations in the time domain, including delay times in a one-dimensional two-wire circuit system 23,24 . However, for the analysis of various applications, it is necessary to connect the external circuits, i.e., lumped-parameter circuits represented by ordinary differential equations, and extend to the three-dimensional conductor system. We note that the three-dimensional conductor can be used for a multi-layer plane conductor system by separating a three-dimensional conductor into multiple planes. We are able also to describe a system with many line conductors with multiple layer conductors.

Simultaneous discretization in space and time of delay-integral and partial-differential equations
In this paper, the potential U, the charge density q, the vector potential A , and the current density j are used as variables to describe electromagnetic phenomena in a three-dimensional conductor. We use the integral equations for the scalar and vector potentials named as retarded potentials obtained from the Maxwell equations in the Lorenz gauge 12-16 . where the distance between the two spacial positions is given as For stable numerical solutions, it is essential to treat the distance R(x − x ′ ) in the denominator and the delay time as rigorously as possible. Here, ε is the permittivity and µ the permeability defined by materials surrounding conductors. The use of the Lorenz gauge over other gauges as the Coulomb gauge may have an advantage for the description of potential and current in conductors 22 .
Besides, we use the continuity equation as the charge conservation for the charge and current density in conductor. We use the Ohm's law, E = ρj , represented by partial differential equations using the potentials and current density.
Here, ρ is the resistivity of the conductor. These four equations for the four quantities, U, A, q, j , are the fundamental equations to describe the transmission and radiation phenomena 16 . They are highly coupled equations, and we ought to develop a reliable algorithm for stable numerical solutions with suitable initial and boundary conditions. Here, the phenomenological Ohm's law represents the effect of the electromagnetic field on the particle motion with charged particles in conductors. A microscopic derivation of the Ohm's law from the many body quantum theory has been performed in the linear response theory 27 .
We start formulation of the numerical method for the stable solution of the above four coupled equations. First, we discretize the integral equations with the delay terms shown in Eqs. (1) and (2). Both the charge and the current exist in a three-dimensional conductor with finite volume. We use the following pulse functions to expand the charge density: These pulse functions f ′ s introduce small volumes in the conductor with the size of x y z.
∂q(x, y, z, t) ∂t + ∇ · j(x, y, z, t) = 0 , ∂A(x, y, z, t) ∂t = ρj(x, y, z, t) . www.nature.com/scientificreports/ The following pulse functions are used for the current density j α with α = x, y, z: Here, the current density's pulse functions are shifted from that of the charge density by a half-integer multiplying the discrete size, whose direction depends on the current density component. The subscript in the upper righthand of j, k, l indicates a half-integer deviation as j ± = j ± 1 2 δ(α, x) , k ± = k ± 1 2 δ(α, y) , l ± = l ± 1 2 δ(α, z) . In the above pulse function only the minus subscript has been used, while the plus subscript will be used later. The charge and current densities can be expressed as follows by using the above pulse functions: The unknown charge and current densities in the discretized space (j, k, l) and time (m) are denoted by q m (j,k,l) and j m+1/2 Using the collocation method, we take collocation points in space as the center of the rectangle, expressed by the pulse function. On the other hand, we noted that a collocation point of time should be the newest time defined by the pulse function. The following points define the collocation points for the charge density and the scalar potential: On the other hand, the current density and the vector potential are defined at a position shifted from the charge density and scalar potential by a half-integer multiplying the discrete size.
The variables defined by the above collocation points are represented as follows: The small volume V (j ′ ,k ′ ,l ′ ) is obtained by the pulse function f ′ s in the x, y and z directions. We further point out that it is important to introduce t (j,k,l) (x ′ , y ′ , z ′ ) , which is the delay time to propagate between the points (x j , y j , z k ) and (x ′ , y ′ , z ′ ) , which are to be integrated out in the spatial integration over where v is the velocity of the transmission of the signal, v = 1/ √ εµ.
To emphasis the delay, we use an integer n for the delay time from the discretized time m. Therefore, we rewrite m ′ = m − n and introduce the following pulse function: Here, the maximum value of n is written as N d , which expresses time to propagate for a signal to the farthest end of the conductor region. N x , N y and N z are numbers of divisions in x, y and z-direction. Finally, we can express Eq. (27) in terms of the delay time n.
A finite volume of potential defined by the pulse functions f j (x), f k (y) and f l (z) in Eqs. (5), (6) and (7) is expressed by solid line cube. The center black point denotes the potential at the collocation point. Red points represent current densities in x-direction, blue y-direction, and green z-direction. Collocation points of current density in x-, y-and z-direction are deviated by a half-integer multiple of the step size from the potential to use the leap-frog method. The deviations are expressed as www.nature.com/scientificreports/ Here, P n is called the "delay local potential coefficient, " which changes the integral range depending on the delay time n.
The delay local potential coefficient is used to calculate the effect of charge densities in the finite volume V (j ′ ,k ′ ,l ′ ) on the collocation point (x j , y k , z l ) , and the range of the effect is shaved in as shown in Fig. 3. In the conventional method, one uses the center-to-center approximation, where the x ′ , y ′ , z ′ coordinates in the delay time are approximated to be at the collocation point of the finite volume V (j ′ ,k ′ ,l ′ ) 22 . We are then able to take out the charge density with the delay time from the integral over x ′ , y ′ , z ′ , and the integrand becomes . Although the calculation of the local potential coefficient becomes simple, this center-to-center approximation leads to coupled delayed partial differential equations.
The vector potential can be also expressed in the same way using the equations (2), (14), (25), and (26): is the "delay local inductance", and the integral range depends on the delay time n.
(31) (10) and (11) is expressed by dashed line cube. Collocation points of potential are deviated by a half-integer multiple of the step size from the current density according to the direction. The deviations are expressed as j ± = j ± δ(α, x) , k ± = k ± δ(α, y) , l ± = l ± δ(α, z). www.nature.com/scientificreports/ Next, we discretize the continuity equation (3) using the FDTD method, which is the central difference method for the differential. The discretized space and time of the charge and current deviate from each other by a half integer in both the space and time directions.
Here, we define the collocation points shifted by an integer as j ± = j ± δ(α, x) , k ± = k ± δ(α, y) , l ± = l ± δ(α, z) . Using Eqs. (31), (33), (35) and (36), we can derive the unknown quantities, U, q, A α and j α at an advanced time using all the corresponding past quantities. Furthermore, we derive the delay difference equations in terms of the potential and current used in the circuit theory. From Eqs. (31) and (35), after eliminating the charge density, the updated equation for the potential is expressed as follows: Transferring the unknown potential at an advanced time to the left-hand side, the updating equation for the potential can be expressed as follows using the current densities with the continuity equation (35): Similarly, by eliminating the vector potential from Eqs. (33), (36) and transposing the new current density to the left-hand side, the updated equation for the current density is expressed as follows: Figure 3. The integral range of the delay local potential coefficient P n (0,0,0)(j,k,l) at z = 0 . The black dot represents the potential's collocation point, and the solid square surrounding the black dot represents the finite volume of the potential for the case of x = y and tv = 1 2 x . The colored circle centered at the collocation point (x 1 , y 1 , z 1 ) has a radius of (n + 1)�tv and represents the integral range in (x ′ , y ′ , z ′ ) of the n delay effect on the colocation point (x 1 , y 1 , z 1 ). www.nature.com/scientificreports/ Here, in the case �α/2 < v�t , the non-delay local mutual inductance L 0 becomes finite for neighboring cells, and simultaneous equations are necessary to be solved. On the other cases, when the neighboring non-delay local inductance is zero for �α/2 ≥ v�t , and the current density at each location can be solved directly: Using Eqs. (38) and (40), we are able to obtain the most advanced potential U and current density j α using known past potentials and current densities.

Boundary conditions of a three-dimensional conductor with a lumped-parameter circuit
In this section, we discuss the treatment of potentials and current densities at the boundaries between a conductor and a lumped-parameter circuit. We consider a case where the current density flows to the surface A of a conductor in the vertical direction through a lumped parameter circuit as shown in Fig. 4. There are other surfaces B, C, D, E, F in Fig. 4, where a lumped parameter circuit may be connected. We derive a boundary condition to connect the lumped-parameter circuit to the delay difference equations we have derived in "Simultaneous discretization in space and time of delay-integral and partial-differential equations" section. We use a nodal potential and a branch current as unknown variables in the lumped-parameter circuit, considering the nodal potential to be equivalent to the potential of the conductor and the branch current equivalent to the conductor's current at the boundary. In a previous study, based on the sparse-tableau method, we have derived an update equation in the time domain for the nodal potential U l and branch current I l of the lumpedparameter circuit 26 . We can summarize the branch constitutive equation (BCE), Kirchhoff 's voltage law (KVL), and Kirchhoff 's current law (KCL) in a matrix and vector form 26 .
Here, U l and I l represents the nodal potential and the branch current, which are unknowns in the lumpedparameter circuit. A l is an incidence matrix representing the node-to-branch connection relationship and the direction. Their components are 0, 1 and −1 , where the row element corresponds to the node, and the column is the branch. For example, when the node p and the branch q are not connected, the matrix element at the pth row and qth column is 0. When the node and branch are connected, and the direction is outward from the node, the matrix element is 1, and when the direction is inward from the node, the element is -1. Z l is the time-domain impedance matrix, which is a diagonal matrix. The matrix element in the qth row and qth column are R when the branch q has a resistor R, �t/(2C) for a capacitor C, and 2L/(�t) for an inductor L. Besides, V s represents the voltage source vector, and the qth element has the value of the voltage source in the branch q. I s represents the current source vector, and qth element has the value of the current source in the branch q.
The matrices ǫ and δ are diagonal matrices for the sign dependence of the electric elements in the branch, and the diagonal β th element is represented as follows: Next, we write the equation that the delay difference equation (38) satisfies at the boundary. There is no concept of space in the lumped-parameter circuit, and we identify the potential of the lumped-parameter circuit to the potential of the three-dimensional conductor at the boundary. In the example shown in Fig. 4, U m l is set equal to U m (j,k,l) . On the other hand, the current density is defined at a half integer time in the delay differential equation Here, γ α(j,k,l) is a coefficient for the sign dependence of the direction of the current density, which depends on the location where the lumped-parameter circuit is connected as shown in Fig. 4 and its value is shown in Table 1.
We note that the current density defined at the integer time is zero at the boundary where no lumped parameter circuit is connected, since there is no current flowing in and out through lumped-parameter circuit. By summarizing the potential and current density at the collocation point where the lumped-parameter circuit is connected, we write Eq. (44) in the matrix form using vectors U d and j d .
Here, A d is an incidence matrix representing the connection relationship between the potential U d and the current j d at the boundary. Since the direction of the current is determined by γ α , the matrix elements of A d are 1 or 0. The rows correspond to the potential and the columns to the current density. For example, when the potential p and the current density q are connected to the same lumped-parameter circuit, the element in the pth row qth column will be 1. On the other hand, when it is not connected, it will be 0. Z 0 d is a non-delay local impedance matrix, and the element in the pth row qth column is 1 2 γ α(q) �t �α P 0 pq where p and q are collocation points of potential. Besides, Q U d is a delay term vector, which is a contribution to the boundary potential and represents the third and fourth items on the right side of Eq. (44). Figure 4. The connection between the potentials and current densities in lumped-parameter and threedimensional circuits at the boundary is shown. The node potential of the lumped-parameter circuit U m l and the potential of the three-dimensional circuit at the boundary U m (j,k,l) are equivalent, and the branch current of the lumped-parameter circuit I m l and the current of the three-dimensional circuit at the boundary j m x(j − ,k,l) are equivalent. The direction of the current is perpendicular to the surface of the boundary. There are six different surfaces from A to F as the boundary. In the figure shown here, as an example, the lumped-parameter circuit is connected at the surface A. www.nature.com/scientificreports/ Finally, we use Kirchhoff 's current law to describe the connection between the lumped-parameter circuit and the three-dimensional conductor.
Here, γ α S α is a diagonal matrix whose components are multiplication of γ α and cross-sections of finite volume in the α(= x, y, z) direction. The above equations (41), (45) and (46) are the boundary conditions. To solve these equations, we write a matrix expression of the same form as the expression (41).
First, the unknown variables U d , U l , I l and j d are summarized below.
The incidence matrixes for the lumped-parameter and three-dimensional circuits are summarized below.
The impedance matrices of the lumped-parameter and three-dimensional circuits are summarized below.
The lumped-parameter circuit's voltage source and the delay term in the delay difference equation are summarized as follows.
The boundary conditions for connecting the lumped-parameter and three-dimensional circuits can be expressed as follows from the above expressions.
Here, ε and δ are matrices for changing the sign depending on the components of the lumped-parameter and three-dimensional circuits.
Furthermore, γ S is the matrix for converting the units of current density in the three-dimensional circuit to current in Kirchhoff 's current law.

Effect of the delay time in numerical results
In this section, we discuss the treatment of delay time for numerical calculations in the time domain. Compared to the rigorous method presented in this paper without approximation of delay time, we use the center-to-center approximation, which has been adopted in the PEEC method 17 .
Here, x i , y k , z l and x j ′ , y k ′ , z l ′ are the collocation points. Exploiting this approximation, we can discretize the charge and current densities in time, including the delay time, which does not influence the integral over space (52) ǫ ββ = −1 for capacitor and three dimensional circuit , 1 otherwise .
(55) t (j,k,l)(j ′ ,k ′ ,l ′ ) = R(x j ′ − x j , y k ′ − y k , z l ′ − z l ) v . Table 1. Coefficients γ x , γ y , and γ z are depending on the lumped-parameter circuit's connection between the nodes of potential (j, k, l) and current density (j ± , k ± , l ± ) at the boundaries A through F shown in Fig. 4. www.nature.com/scientificreports/ in the local impedance. Since we can pull out the charge and current densities from the space integral, we obtain a simplified delay local impedance: Here, in the center-to-center approximation, the discretized delay time n is uniquely defined by the pulse function in Eq. (29) regardless of the space integral range, which is equivalent to the approximation using floor function 28 .
For the demonstration of the two treatments of the delay time, we use a single-plane circuit consisting of a plane conductor as shown in Fig. 5. The step size in space is x = y = 0.2 mm and the step size in time is 0.5 × �x/v . Figure 6 shows the numerical results of delay local self-and mutual-impedances in this numerical condition. Here, the calculation methods of delay local impedance are described in the Supplementary Information. When we rigorously treat the delay time, some delays occur even in the local self-impedance since the delay time influences the integral over positions in the finite volume. We observe that the delay local impedances smoothly vary as positions of the cells and the delay time are increased. In contrast, when we use the center-to-center approximation, all the coordinates within the finite volume simultaneously contribute to the local impedance. Furthermore, when we define the delay time in terms of a pulse function, the delay local mutual-impedance of different positions is defined at the same delay time. Figure 7(a) compares the numerical results using two different methods in handling the delay time, and (b) compares the numerical results with and without the delay time in the rigorous method. As shown in Fig. 7(a), we see that the numerical results are stable when the delay time is treated rigorously. On the other hand, the calculations using the center-to-center approximation provides unstable results. Therefore, the rigorous treatment of the space-dependent delay time is vital for stable numerical results when we take into account the delay time. We further show numerical results by changing spacial mesh sizes x, y and α for the rigorous and approximate cases in Sec. B of Supplementary Information. We find stable solutions in various parameter ranges for the rigorous case, but no stable solutions for the case of the center-to-center approximation in the parameter range explicitly calculated. Figure 7(b) compares the numerical results with and without the delay time. The current in the single-plane circuit with the delay time is attenuated, caused by the external radiation. On the other hand, when the delay time is not taken into account, the current is not attenuated, and moves back and forth between the two ends of a single-plane antenna. From these results, it is possible to treat the radiation and transmission phenomena from the circuit by taking explicitly the delay time into account. We have made calculations of two plane circuits using the rigorous method together with the center-to-center approximation for the delay term. We have obtained stable solutions in the rigorous method.

Summary
In this paper, a new numerical calculation method was developed for full-wave analysis in the time domain, considering transmission and radiation in the three-dimensional conductors. The characteristic feature of present numerical method is the simultaneous discretization of time and space using the collocation method, which enables us to treat the delay time rigorously.
We compared numerical results for two cases where we treat the delay time rigorously and approximately using the center-to-center approximation. The results showed that rigorous treatment of the delay time presented in this paper gave a stable numerical result, while the approximate treatment, which uses center-to-center approximation, gave unstable results. For stable numerical calculation of delay time, it was essential to consider the space dependence in the delay term rigorously. We also compared numerical results with and without delay time. The results showed that, in the single-plane circuit, the current significantly decreased in the propagation process due to radiation from the single-plane circuit.
(56) Z n (j,k,l)(j ′ ,k ′ ,l ′ ) = 1 4π dx ′ dy ′ dz ′ Figure 5. A single-plane circuit consists of a conductor plane of 1 mm width and 20 mm length. In numerical calculations, the step size in space is set to x = y = 0.2 mm and the step size in time is set to �tv = �x/2 , where v is the velocity of the signal. A current source is connected to the left end of the single-plane circuit, where pulse current is used as input to the boundary. The boundary was calculated numerically using the method described in "Boundary conditions of a three-dimensional conductor with a lumped-parameter circuit" section.  Figure 6. Delay local impedances between two finite volumes after discretizing the single-plane circuit as functions of discretized delay time n. The solid line represents the delay local self-impedance Z n (1,1,1)(1,1,1) . The dashed and dotted lines are the delay local mutual impedances Z n (1,1,1)(1,2,1) and Z n (1,1,1)(2,2,1) between the neighboring cells of the origin (x 1 , y 1 , z 1 ) , respectively. The red lines denote the result of the rigorous method presented in this paper, the green lines the center-to-center approximation method. In the Supplementary Information, we explain the detailed calculation method of delay local impedance in the plane circuit.