Propagation Modeling of Point Source Excited Magnetoinductive Waves Based on a New Plane Wave Expansion Approach

1Beijing Engineering Research Center of Optoelectronic Information and Instrument, Beijing Key Laboratory of Optoelectronic Measurement Technology, Beijing Information Science & Technology University, Beijing 100192, China 2Department of Biomedical Engineering, Oregon Health and Science University, Beaverton, OR 97006, USA 3School of Mechanical Engineering, University of Science & Technology, Beijing 100083, China


Introduction
Wireless underground sensor networks (WUSNs) become an active research in monitoring underground environments such as mineral contents, oil leakage, and formation stress [1][2][3][4].Since the propagation medium of underground environments is not air and contains soil, rocks, and water, wireless communication in WUSN faces severe challenges caused by the high path loss.
The classical wireless signal propagation for electromagnetic (EM) waves is not effective in working in underground environments [5][6][7].Since the waves with low frequency EM suffer from lower fading in loss media, this can be helpful to set the WUSN working at the VLF band similar to most of the telecommunication systems [8][9][10].However, large size antennae and high drive voltage are required to transmit EM waves in VLF band.This approach is usually impracticable in underground environments.The emergence of MI propagation technology provides a practical solution to this problem.
The MI technology has been extensively studied in the metamaterials field [11][12][13].In [14], Shamonina et al. applied the concept of MI waves for the loop current waves in one, two, and three dimensions.Syms et al. studied the positive and negative refraction of MI waves in two dimensions [15].In addition, they applied the concept of MI waveguide with the quantitative analysis to several kinds of waveguide devices [16].The MI waveguide was introduced by Sun and Akyildiz in [17] to the WUSNs field as a new proposed technology in the physical layer of sensor networks.The studies in [18][19][20][21] show that the MI waveguides have the capabilities of transmitting powers.Therefore, the MI waveguides are useful to recharge the WUSN nodes.In conclusion, the current studies of the MI structures mostly focused on the plane wave propagation.However, in a real WUSN, signals are always transmitted by a single loop, and the wave excited by the point source cannot be considered as a plane wave.
Since the frequency of MI networks only depends on the resonance frequency of the split ring resonators, this can be exploited to adjust the load capacitor and the number of windings to set the MI structure working at the VLF band, keeping the transmit antenna at a low size level.At the same time, the resonant capacitor will provide auxiliary power potential to drive the loop coils, which reduce the required driving voltage to a very low level.
Since the propagation model and wave patterns provide fundamental basis for link budget, topology controlling, and the design of communication protocol, the aim of this paper is to show the characteristics of the near field propagation of the MI waves that are excited by a point source.The frequent usage of square lattice MI structure is investigated.The simulation of time domain based on the fourth-order Runge-Kutta (RK) method is carried out to verify the mathematical derivation.
The rest of this paper is organized as follows.Section 2 presents the system model, including the dispersion equations of MI structure.The wave propagation analysis based on the plane wave expansion method is illustrated in Section 3. The simulation results and the verification using a time domain simulation are given in Section 4. The discussions are presented in Section 5, and Section 6 concludes this paper.

System Model
In classical WUSNs, nodes are distributed in monitoring areas.A square lattice with the same lattice constant "" can be considered a simple topological structure and suitable for WUSNs.As shown in Figure 1, every node in WUSNs is equipped with a capacitive loaded loop antenna.
In order to keep the isotropy, only the coplanar configuration of the network is discussed in this paper which is shown in Figure 2.
For the array of MI loops, only the mutual effects between nearest neighbours need to be considered [22].Thus, in a square lattice WUSN, each node can be considered to have relationships only with four neighbours.
In the coplanar loops, as shown in Figure 3, the flux is excited by a loop itself which has opposite direction compared with that excited by other loops.Thus, the mutual inductance  is negative and denoted by −. Figure 4 shows a simple equivalent circuit for one loop in the two-dimensional (2D) wave as shown in Figure 1.In Figure 4, each resonator is modelled by the connection of an inductance (), a capacitance (), and four mutual inductances of −, which are coupled from its nearest neighbors.In addition, in each loop, Kirchhoff 's voltage law is applied considering the relations with four neighbour loops' currents ( , ).
In the lossless case, we obtain where  , () is the voltage source of the node at th row and th column,  , () is the antenna current,  and  are the ( is the mutual inductance,  is the self-inductance, and  is the capacitor load.) inductance and capacitance of the antenna, respectively, and − is the mutual inductance.
In the frequency domain, (1) can be rewritten as where  =  + 1/ is the loop impedance and the coupling term is defined as −  = −.When the voltage source is vanished, we assume a travelling-wave solution for the simplest solution as where ⃗  =  ⃗  + V ⃗  is the wave vector and ⃗  =  ⃗  +  ⃗  is the position vector.Substituting (3) into (2), we obtain the constraint equation Equation ( 4) can be rewritten as where  0 = 1/ √  is the resonance frequency and  = 2/ is the coupling coefficient.Since (( 2 −  2 0 )/ 2 )/ < 2 and  is a very small value in WUSNs, the bandwidth will be very narrow.Figure 5 shows some curves of wave vectors.Considering the monochrome signal at the resonance frequency  0 , (5) can be rewritten as cos  + cos V = 0.
The resonant frequency wave vector curves are four straight lines that are joined end to end forming a square frame.In Figure 5, the curve of the ⃗  at lower frequencies is located inside the square frame of the resonance frequency curve, with higher frequency curves outside of the enclosure.

The Plane Wave Expansion of MI Waves Excited by Point
Source.We can consider the  − ⃗  ⃗  in (3) as a space signal whose Fourier transformation is known as wavenumber spectrum.As it is known, certain wave vectors are allowed in a 2D MI waveguide.In addition, we can describe the phenomena in other words that the MI waveguide is a space system and produces different response to the MI wave amplitude with different propagation constants; that is, MI waveguide is a system of space signal.It can be verified that the MI waveguide space system is a linear space invariant system (LSI) and is mathematically similar to a linear time invariant system (LTI).
Since the space parameters of the loop are the coordinates of the circle centre, a single loop exciting source can be considered as a point source and expressed with a pulse function.The wavenumber spectrum of point source is 1.As shown in Section 2, the MI waveguide is an ideal filter of wave vector with unit pass-band gain and the stop-band gain of 0, and there is a uniform system gain for every possible .As proven in the current studies, the possible form of wave vectors is four separate continuous curves in the loop array of square lattice.Thus, the space distribution of the loop current excited by a single loop source (point source) represents a single frequency signal and can be written using the inverse Fourier transformation as where () is the curve that is composed of the dots with coordinates (, V) as shown in (5). in ( 7) can be obtained considering the node at the origin of the space coordinates, which can be expressed by  where (0, 0, ) is the amplitude of the loop current at position of  = 0 and  = 0.In addition, (0, 0, ) could be a complex number.
Considering the current source is a periodic signal with a period of  0 , the loop current can be expressed by The current amplitude at any coordinates can be obtained from ( 6), (7), and (8) as When  0 → ∞, we can get the current at position (, ) excited by any nonperiodic point source: Because the propagation of vector curve, (), is symmetrical in both the  and  coordinates, ( 9) and ( 10) are rewritten as where V = arccos(( 2 −  2 0 )/ 2  − cos())/ as shown in (5).

Mirror source
Mirror source

Current source
Figure 7: The multimirror source.

The Reflection and Matching of Regular Boundaries.
In most applications, two kinds of regular edges are commonly used as shown in Figure 6.
In Figure 6, the loops in the shadow are the actual waveguide.For the convenience of analysis, the virtual loops are drawn with dotted lines and attached to the edge with zeros currents as the boundary conditions.
In an infinite 2D MI waveguide, an identical boundary condition is a mirror current source at the position symmetrical to the actual source for the straight line through the centers of attached loops.The MI waves with the reflection of the edge in the actual areas of the 2D MI waveguide could be considered as the superposition of the waves.In addition, these waves are excited by the two point sources.
As shown in Figure 7, in a finite 2D MI waveguide with more than one edges, the theory of geometrical optics can help us to find all of the mirror sources.In addition, this is useful to obtain the wave propagation that is resulted by the superposition of all PWE from the sources.
In some applications, the reflection is very harmful to the signal propagation.Thus, the edges should be designed to absorb any incident powers and the network should be terminally matched.
According to the situation in Figure 6(a), we first assume that there exists an infinite 2D MI waveguide.Furthermore, we remove half of the loops to form a line edge.As mentioned previously, the second loop (2 as shown in Figure 6(a)) has effect only on the first loop (1 as shown in Figure 6(a)).To keep everything normal in the first loop and remove the second loop, we should connect the impedance in series form as the second loop is still there.Thus, (2) can be rewritten in (13) with substituting   by the matching impedance  matching , where  , is not zero and  matching is obtained as Substituting ( 12) into ( 14) and considering the Fourier transform formula, we obtain where (, , ) is a real function of angular frequency and loop position.Thus, the matching device would be an inductor or a capacitor and its value should be calculated for every position and frequency.Note that it is impossible to realize perfect terminal impedance matching in a 2D MI waveguide.However, an approximate matching may be acceptable for a narrow band telecommunication with the consideration of the matching impedance of the resonance frequency.

Numerical Result and Time Domain Simulation
Since it is difficult to obtain the analytical solution of the PWE in close forms, the numerical evaluation is carried out.Table 1 shows the conditions of numerical analysis.We consider node array with 121 by 121.The source is located in the array centre.The time domain simulation is carried out to verify the numerical results.As shown in (1) and O p e r a t i n g f r e q . 1 0 K H z using a loop at the th column and th row, the differential equations can be set up with the form as For the loops at the edge and corner, some of the items are removed.Matching impedance is also inserted according to (15).
Assuming that  , ()/ =  , (), ( 16) is rewritten as The second-order differential equation can be transformed into two equations with one order.Thus, the whole number of differential equations in this test is 29282.With the fourth-order RK method, we can obtain the time domain numerical solution.To limit the bandwidth of the signal, a rising cosine window is used to modulate the carrier as where   = 10 and  0 = 20000.Both PWE and simulation results reveal that the narrow MI wave packet propagates in four main directions.Figure 8 shows a snapshot of the propagation at the time of 0.9 seconds.
We redraw the current amplitude of the loops along the diagonal line in Figure 9.As can be seen in Figures 8 and  9, the PWE results match with the results of time domain simulations, except that there are more fluctuations in the PWE results.In addition, we see that due to the symmetry of the loop array, the propagation of the MI wave is also  symmetrical at 61th row and 61th column, which divides the loop array into four identical parts.Figure 8 illustrates that the propagation speeds and the wave patterns of the wave packets (obtained with two different methods) are almost the same along the diagonal line.The wave packet patterns also show that there is serious dispersion with the signal transmission.
Figure 8 shows that the MI wave packet propagates along several different paths.This is clearly seen in Figure 10, where the amplitudes of the wave front at 0.9 s are given for different azimuth angles form 0 to 90 degrees.In Figure 10, the PWE results contain more chaotic ups and downs than the time domain simulation results.However, the outlines of the pattern are coincident.

Discussion
The PWE method is based on the plane wave superposition.Thus, at any time, the space is already full of steady plane waves.As can be seen in Figures 8(a) and 8(b), before arriving of the narrow wave packets, there are plenty of wave distributions in the spare area.Moreover, in the real and negative times, the future wave packets introduce fluctuations at presence of time and area (without physical meaning).This explains why the PWE results contain more ripples than the simulations of time domain.
The simulation of time domain is based on solving the differential equations and the Runge-Kutta method.These are commonly used as tools for this purpose.The loop array of 121 by 121 includes solving 29282 differential equations and introducing a huge amount of calculations.This occupies a great deal of computing resources.Therefore, in this paper, the running time takes about 70 hours to accomplish the MI wave simulation using a computer with a 3.4 GHz 4-Core CPU and 32 GB memory.However, it takes only lower running time than half an hour to accomplish identical simulation with the PWE method.Therefore, even though the PWE method is only an approximation of the real propagation, it is a very efficient method to analyse the MI waves.
To explain the formation of the pattern of the MI wave propagation, we review the periodic structure about the bandwidth and the energy flow.Using (5), we can obtain the bandwidth of different propagation directions.The normalized bandwidth of the simulated loop array is given in Figure 11.
As can be seen in Figure 11, the diagonal direction has the maximum bandwidth.
By considering the center node (, ) and its four neighbors in Figure 2, the neighbors' power which is transmitted and received by the center node is expressed as Substituting ( 3) into (19) and assuming that  =  , , we can obtain the dispersion equation as where  left to center is the power flow in  direction and is marked with   and the  down to center is marked with   ; thus (20) can be rewritten as The result in ( 21) is also deduced through the group speed calculation.Figure 12 shows the normalized energy propagating speed of the simulated loop array.
As illustrated in Figures 11 and 12, the main part of the signal power is bounded into the diagonal directions and this is called the main propagation direction (MPD).In other directions, the propagation is more complicated.To illustrate this issue, an example in Figure 13 is presented.
In Figure 13, a narrow wave packet, called W1, departs from the source into the direction of D1 and the energy flowing direction will be D3 towards the MPD.The secondary wave packet, called W2, which is excited by W1, will eventually turn to the direction of D2.This is due to the energy flow and the new propagation path which is parallel to the MPD.This path is called auxiliary propagating direction (APD).As shown in Figure 8, the wave packet propagates along the MPD directions; thus new secondary waves will be excited and resulted in several trajectories of discrete parallel wave propagation.

Conclusion
The numerical results and time domain simulations are presented.The results and simulations show that the proposed PWE method is efficient and accurate to model the MI wave propagation in the near field of a point source.It can speed up the analysis dozens of times.From the theoretical derivations and simulations, we can conclude that the MI wave propagation in a 2D resonance array is not similar to that of an EM wave in some continuous mediums.The MI waves, which are excited by a point source, mainly propagate along the diagonal directions (MPDs).Secondary waves form multipropagation trajectories (APDs) are parallel to the MPDs.Since the main power is transmitted along the MPDs, the 2D array should employ a switching and relay mechanism when setup of the communication between two nodes is not in the same MPD.For future works, the effects of other topologies, the communication protocols, and beam forming technologies will be analysed in the MI 2D waveguides.

Figure 2 :
Figure 2: Geometry of the loops ( is the column number and  is the row number).

Figure 3 :
Figure 3: Mutual inductance between two coplanar loops.(  is the self-induced flux and   is the mutual induced flux.)

Figure 4 :
Figure 4: Equivalent circuit model of the loop in the lattice array.( is the mutual inductance,  is the self-inductance, and  is the capacitor load.)

Figure 6 :
Figure 6: The line edge of a 2D MI waveguide.

Figure 8 :
Figure 8: Snapshot of the propagation of MI narrow wave packet excited by a point source at the 61th row and the 61th column (the grey level represents the amplitude of the current).

Figure 9 :
Figure 9: The current amplitude along the diagonal line.