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

Abstract

The signal fading in wireless underground sensor networks (WUSNs), which is caused by lossy media such as soil and sand, can be reduced by applying technology of magnetoinductive (MI) propagation. This technology can effectively establish a communication at very low frequency (VLF). In contrast to the previous studies in the literature, which mostly focus on the propagation of plane waves, we propose a new approach based on the plane wave expansion (PWE) to model the near field MI waves. The proposed approach is based on excitation of a point source, which is a common case in a practical WUSN. The frequent usage of square lattice MI structure is investigated. To verify the mathematical derivation, the simulation of time domain based on the fourth-order Runge-Kutta (RK) method is carried out. Simulation results show that the new model can provide a precise prediction to the MI wave's propagation, with the computation load being one-tenth of that of the time domain simulation. The characteristics of the propagation of the MI waves are presented and discussed. Finally, the reflection on the edge of the MI structure is reduced by analysing the terminal matching conditions and calculating a method for matching impedances.

Feng Liu 1 and Zhijun Zhang 2 and Lianqing Zhu 1 and Yajun Liu 3

Academic Editor:Hsuan-Ling Kao

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

Received 12 June 2014; Revised 7 October 2014; Accepted 11 October 2014; 4 March 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Wireless underground sensor networks (WUSNs) become an active research in monitoring underground environments such as mineral contents, oil leakage, and formation stress [1-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-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-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-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-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.

2. System Model

In classical WUSNs, nodes are distributed in monitoring areas. A square lattice with the same lattice constant " [figure omitted; refer to PDF] " 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.

Figure 1: Schematic representation of a capacitively loaded wire loop.

[figure omitted; refer to PDF]

In order to keep the isotropy, only the coplanar configuration of the network is discussed in this paper which is shown in Figure 2.

Figure 2: Geometry of the loops ( [figure omitted; refer to PDF] is the column number and [figure omitted; refer to PDF] is the row number).

[figure omitted; refer to PDF]

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 omitted; refer to PDF] . 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 ( [figure omitted; refer to PDF] ), a capacitance ( [figure omitted; refer to PDF] ), and four mutual inductances of [figure omitted; refer to PDF] , 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 ( [figure omitted; refer to PDF] ).

Figure 3: Mutual inductance between two coplanar loops. ( [figure omitted; refer to PDF] is the self-induced flux and [figure omitted; refer to PDF] is the mutual induced flux.)

[figure omitted; refer to PDF]

Figure 4: Equivalent circuit model of the loop in the lattice array. ( [figure omitted; refer to PDF] is the mutual inductance, [figure omitted; refer to PDF] is the self-inductance, and [figure omitted; refer to PDF] is the capacitor load.)

[figure omitted; refer to PDF]

In the lossless case, we obtain [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the voltage source of the node at [figure omitted; refer to PDF] th row and [figure omitted; refer to PDF] th column, [figure omitted; refer to PDF] is the antenna current, [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the inductance and capacitance of the antenna, respectively, and [figure omitted; refer to PDF] is the mutual inductance.

In the frequency domain, (1) can be rewritten as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the loop impedance and the coupling term is defined as [figure omitted; refer to PDF] . When the voltage source is vanished, we assume a travelling-wave solution for the simplest solution as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the wave vector and [figure omitted; refer to PDF] is the position vector. Substituting (3) into (2), we obtain the constraint equation [figure omitted; refer to PDF]

Equation (4) can be rewritten as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the resonance frequency and [figure omitted; refer to PDF] is the coupling coefficient. Since [figure omitted; refer to PDF] and [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] , (5) can be rewritten as [figure omitted; refer to PDF]

Figure 5: The propagation vector curves.

[figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] at lower frequencies is located inside the square frame of the resonance frequency curve, with higher frequency curves outside of the enclosure.

3. Propagation of MI Waves

3.1. The Plane Wave Expansion of MI Waves Excited by Point Source

We can consider the [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the curve that is composed of the dots with coordinates ( [figure omitted; refer to PDF] ) as shown in (5). [figure omitted; refer to PDF] in (7) can be obtained considering the node at the origin of the space coordinates, which can be expressed by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the amplitude of the loop current at position of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . In addition, [figure omitted; refer to PDF] could be a complex number.

Considering the current source is a periodic signal with a period of [figure omitted; refer to PDF] , the loop current can be expressed by [figure omitted; refer to PDF]

The current amplitude at any coordinates can be obtained from (6), (7), and (8) as [figure omitted; refer to PDF]

When [figure omitted; refer to PDF] , we can get the current at position ( [figure omitted; refer to PDF] ) excited by any nonperiodic point source: [figure omitted; refer to PDF]

Because the propagation of vector curve, [figure omitted; refer to PDF] , is symmetrical in both the [figure omitted; refer to PDF] and [figure omitted; refer to PDF] coordinates, (9) and (10) are rewritten as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] as shown in (5).

3.2. The Reflection and Matching of Regular Boundaries

In most applications, two kinds of regular edges are commonly used as shown in Figure 6.

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

(a) Side edge

[figure omitted; refer to PDF]

(b) Hypotenuses

[figure omitted; refer to PDF]

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.

Figure 7: The multimirror source.

[figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] by the matching impedance [figure omitted; refer to PDF] , [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is not zero and [figure omitted; refer to PDF] is obtained as [figure omitted; refer to PDF]

Substituting (12) into (14) and considering the Fourier transform formula, we obtain [figure omitted; refer to PDF] where [figure omitted; refer to PDF] 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.

4. 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.

Table 1: Numerical analysis conditions.

	Item	Value
1	Shape	Circle
2	Loop radius	0.1 m
3	Turns	20
4	Loop space	1 m
5	Operating freq.	10 KHz

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 using a loop at the [figure omitted; refer to PDF] th column and [figure omitted; refer to PDF] th row, the differential equations can be set up with the form as [figure omitted; refer to PDF]

For the loops at the edge and corner, some of the items are removed. Matching impedance is also inserted according to (15).

Assuming that [figure omitted; refer to PDF] , (16) is rewritten as [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] .

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.

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).

(a) The PWE result

[figure omitted; refer to PDF]

(b) The time domain simulation

[figure omitted; refer to PDF]

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 9: The current amplitude along the diagonal line.

[figure omitted; refer to PDF]

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.

Figure 10: The wave front amplitude.

(a) The PWE result

[figure omitted; refer to PDF]

(b) The time domain simulation

[figure omitted; refer to PDF]

5. 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.

Figure 11: The bandwidth of different propagation directions.

[figure omitted; refer to PDF]

As can be seen in Figure 11, the diagonal direction has the maximum bandwidth.

By considering the center node [figure omitted; refer to PDF] and its four neighbors in Figure 2, the neighbors' power which is transmitted and received by the center node is expressed as [figure omitted; refer to PDF]

Substituting (3) into (19) and assuming that [figure omitted; refer to PDF] , we can obtain the dispersion equation as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the power flow in [figure omitted; refer to PDF] direction and is marked with [figure omitted; refer to PDF] and the [figure omitted; refer to PDF] is marked with [figure omitted; refer to PDF] ; thus (20) can be rewritten as [figure omitted; refer to PDF]

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.

Figure 12: The energy flow speed.

[figure omitted; refer to PDF]

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.

Figure 13: The established process of APD.

[figure omitted; refer to PDF]

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.

6. 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.

Acknowledgments

The authors wish to thank National Natural Science Foundation of China (fund no. 61273071), Beijing Nature and Science Foundation (fund no. 8142028), and National High Technology Research and Development Program (fund no. SS2013AA060204) for their financial support of the work in this paper.

Conflict of Interests

The authors have declared that no conflict of interests exist.

References

[1] M. Li, Y. Liu, "Underground structure monitoring with wireless sensor networks," in Proceedings of the 6th International Symposium on Information Processing in Sensor Networks (IPSN '07), pp. 69-78, Cambridge, Mass, USA, April 2007.

[2] M. C. Vuran, I. F. Akyildiz, "Cross-layer packet size optimization for wireless terrestrial, underwater, and underground sensor networks," in Proceedings of the 27th IEEE Communications Society Conference on Computer Communications (INFOCOM '08), pp. 780-788, Phoenix, Ariz, USA, April 2008.

[3] A. R. Silva, M. C. Vuran, "Communication with aboveground devices in wireless underground sensor networks: an empirical study," in Proceedings of the IEEE International Conference on Communications (ICC '10), pp. 1-6, Cape Town, South Africa, May 2010.

[4] Z. Sun, I. F. Akyildiz, "Connectivity in wireless underground sensor networks," in Proceedings of the 7th Annual IEEE Communications Society Conference on Sensor, Mesh and Ad Hoc Communications and Networks (SECON '10), pp. 1-9, Boston, Mass, USA, June 2010.

[5] X. Dong, M. C. Vuran, "A channel model for wireless underground sensor networks using lateral waves," in Proceedings of the 54th Annual IEEE Global Telecommunications Conference: "Energizing Global Communications" (GLOBECOM '11), Houston, Tex, USA, December 2011.

[6] E. P. Stuntebeck, D. Pompili, T. Melodia, "Wireless underground sensor networks using commodity terrestrial motes," in Proceedings of the 2nd IEEE Workshop on Wireless Mesh Networks (WiMESH '06), pp. 112-114, Reston, Va, USA, September 2006.

[7] L. Li, M. C. Vuran, I. F. Akyildiz, "Characteristics of underground channel for wireless underground sensor networks," in Proceedings of the IFIP Mediterranean Ad Hoc Networking Workshop (Med-HocNet '07), Corfu, Greece, 2007.

[8] A. Munoz, A. Mediano, "Output design considerations in wireless portable Through-The-Earth communications system using current injection," in Proceedings of the 42nd European Microwave Conference (EuMC), Amsterdam, The Netherlands, 2012.

[9] M. R. Yenchek, G. T. Homce, N. W. Damiano, J. R. Srednicki, "NIOSH-sponsored research in through-the-earth communications for mines: a status report," IEEE Transactions on Industry Applications , vol. 48, no. 5, pp. 1700-1707, 2012.

[10] L. Yan, J. A. Waynert, C. Sunderman, "Measurements and modeling of through-the-earth communications for coal mines," IEEE Transactions on Industry Applications , vol. 49, no. 5, pp. 1979-1983, 2013.

[11] M. Eleftheriou, N. Lazarides, G. P. Tsironis, Y. S. Kivshar, "Surface magnetoinductive breathers in two-dimensional magnetic metamaterials," Physical Review E: Statistical, Nonlinear, and Soft Matter Physics , vol. 80, no. 1, 2009.

[12] O. Sydoruk, A. Radkovskaya, O. Zhuromskyy, "Tailoring the near-field guiding properties of magnetic metamaterials with two resonant elements per unit cell," Physical Review B , vol. 73, no. 22, 2006.

[13] P. Gay-Balmaz, O. J. F. Martin, "Electromagnetic resonances in individual and coupled split-ring resonators," Journal of Applied Physics , vol. 92, no. 5, pp. 2929-2936, 2002.

[14] E. Shamonina, V. A. Kalinin, K. H. Ringhofer, L. Solymar, "Magnetoinductive waves in one, two, three dimensions," Journal of Applied Physics , vol. 92, no. 10, pp. 6252-6261, 2002.

[15] R. R. A. Syms, E. Shamonina, L. Solymar, "Positive and negative refraction of magnetoinductive waves in two dimensions," European Physical Journal B , vol. 46, no. 2, pp. 301-308, 2005.

[16] R. R. A. Syms, E. Shamonina, L. Solymar, "Magneto-inductive waveguide devices," IEE Proceedings: Microwaves, Antennas and Propagation , vol. 153, no. 2, pp. 111-121, 2006.

[17] Z. Sun, I. F. Akyildiz, "Magnetic induction communications for wireless underground sensor networks," IEEE Transactions on Antennas and Propagation , vol. 58, no. 7, pp. 2426-2435, 2010.

[18] F. Zhang, S. A. Hackworth, W. Fu, M. Sun, "The relay effect on wireless power transfer using witricity," in Proceedings of the 14th Biennial IEEE Conference on Electromagnetic Field Computation (CEFC '10), Chicago, Ill, USA, 2010.

[19] O. Jonah, S. V. Georgakopoulos, "Wireless power transfer in concrete via strongly coupled magnetic resonance," IEEE Transactions on Antennas and Propagation , vol. 61, no. 3, pp. 1378-1384, 2013.

[20] H. Lim, K. Ishida, M. Takamiya, T. Sakurai, "Positioning-free magnetically resonant wireless power transmission board with staggered repeater coil array (SRCA)," in Proceedings of the IEEE MTT-S International Microwave Workshop Series on Innovative Wireless Power Transmission: Technologies, Systems, and Applications (IMWS '12), pp. 93-96, 2012.

[21] Y. Narusue, Y. Kawahara, T. Asami, "Impedance matching method for any-hop straight wireless power transmission using magnetic resonance," Proceedings of the IEEE Radio and Wireless Symposium (RWS '13) , pp. 193-195, 2013.

[22] M. Eleftheriou, N. Lazarides, G. P. Tsironis, "Magnetoinductive breathers in metamaterials," Physical Review E , vol. 77, no. 3, 2008.