Stored energies and Q-factor of two-dimensionally periodic antenna arrays

TL;DR

This paper derives a shape-independent, current-based expression for the Q-factor of two-dimensionally periodic antenna arrays, enabling accurate bandwidth prediction and beam steering considerations.

Contribution

It introduces a novel, shape-independent Q-factor formulation for periodic arrays that directly relates to current density and can be rapidly computed.

Findings

Numerical results show good agreement between the proposed Q-factor and fractional bandwidth.

The method accurately predicts Q-factor for arrays with Q ≥ 5.

The approach includes the input-impedance Q-factor by Yaghjian and Best, with consistent results.

Abstract

The Q-factor for lossless three-dimensional structures with two-dimensional periodicity is here derived in terms of the electric current density. The derivation in itself is shape-independent and based on the periodic free-space Green's function. The expression for Q-factor takes into account the exact shape of a periodic element, and permits beam steering. The stored energies and the radiated power, both required to evaluate Q-factor, are coordinate independent and expressed in a similar manner to the periodic Electric Field Integral equation, and can thus be rapidly calculated. Numerical investigations, performed for several antenna arrays, indicate fine agreement, accurate enough to be predictive, between the proposed Q-factor and the tuned fractional bandwidth, when the arrays are not too wideband (i.e., when $Q\geq 5$). For completeness, the input-impedance Q-factor, proposed by Yaghjian and Best in 2005, is included and agrees well numerically with the derived Q-factor expression. The main advantage of the proposed representation is its explicit connection to the current density, which allows the Q-factor to give bandwidth estimates based on the shape and current of the array element.

Figures

Equations

$$\mathbf{J}(\mathbf{r}+{\mathbf{\zeta }}_{mn})=\mathbf{J}(\mathbf{r}){\mathrm{e}}^{\mathrm{j}{\mathbf{k}}_{\mathrm{t}00}\cdot {\mathbf{\zeta }}_{mn}}.$$

$$Q=\frac{2\omega \max \{W_{\mathrm{e}},W_{\mathrm{m}}\}}{P_{\mathrm{tot}}}.$$

$$W_{\mathrm{e}}=\frac{ϵ}{4}{\int }_U|\mathbf{E}{|}^2-|{\mathbf{E}}_{\mathrm{p}}{|}^2\mathrm{d}v,$$

$$W_{\mathrm{m}}=\frac{\mu }{4}{\int }_U|\mathbf{H}{|}^2-|{\mathbf{H}}_{\mathrm{p}}{|}^2\mathrm{d}v,$$

$$-\nabla \cdot (\mathbf{E}\times {\mathbf{H}}^{\ast })+\mathrm{j}\omega (ϵ|\mathbf{E}{|}^2-\mu |\mathbf{H}{|}^2)=\mathbf{E}\cdot {\mathbf{J}}^{\ast },$$

$$\begin{array}{rl}\frac{1}{2}\mathrm{Re}{\int }_{\delta U_d}& \left[\mathbf{E}(\mathbf{r})\times {\mathbf{H}}^{\ast }(\mathbf{r})\right]\cdot \hat{\mathbf{n}}(\mathbf{r})\mathrm{d}S\\ & =-\frac{1}{2}\mathrm{Re}{\int }_{\Omega }\mathbf{E}(\mathbf{r})\cdot {\mathbf{J}}^{\ast }(\mathbf{r})\mathrm{d}v,\end{array}$$

$$P_{\mathrm{r}}=\mathrm{Re}\left\{-\frac{1}{2}{\int }_{\Omega }\mathbf{E}\cdot {\mathbf{J}}^{\ast }\mathrm{d}v\right\}=\mathrm{Re}{P}_{\mathrm{c}},$$

$$({\nabla }^2+k^2)G({\mathbf{r}}_1,{\mathbf{r}}_2)=-\sum _{(m,n)\in {\mathbb{Z}}^2}\delta ({\mathbf{r}}_1-{\mathbf{r}}_2-{\mathbf{\zeta }}_{mn}){\mathrm{e}}^{-\mathrm{j}{\mathbf{k}}_{\mathrm{t}00}\cdot {\mathbf{\zeta }}_{mn}}.$$

$$G({\mathbf{r}}_1,{\mathbf{r}}_2)=\frac{1}{2\mathrm{j}S}\sum _{(m,n)\in {\mathbb{Z}}^2}\frac{1}{k_{\mathrm{z}mn}}{\mathrm{e}}^{-\mathrm{j}{\mathbf{k}}_{\mathrm{t}mn}\cdot ({\mathbf{\rho }}_1-{\mathbf{\rho }}_2)}{\mathrm{e}}^{-\mathrm{j}{k}_{\mathrm{z}mn}|z_1-z_2|},$$

$$\nabla \cdot \mathbf{A}+\mathrm{j}\omega ϵ\mu \varphi =0.$$

$$\mathbf{E}=-\nabla \varphi -\mathrm{j}\omega \mathbf{A}.$$

$$\varphi ({\mathbf{r}}_1)=\frac{1}{ϵ}{\int }_{\Omega }\varrho ({\mathbf{r}}_2)G({\mathbf{r}}_1,{\mathbf{r}}_2)\mathrm{d}{v}_2,$$

$$\mathbf{A}({\mathbf{r}}_1)=\mu {\int }_{\Omega }\mathbf{J}({\mathbf{r}}_2)G({\mathbf{r}}_1,{\mathbf{r}}_2)\mathrm{d}{v}_2.$$

$$\begin{array}{rl}|\mathbf{E}(\mathbf{r}){|}^2=& |\nabla \varphi {|}^2+{\omega }^2|\mathbf{A}{|}^2\\ & +2\mathrm{Re}\{\mathrm{j}\omega \nabla \cdot (\mathbf{A}{\varphi }^{\ast })\}-2k^2|\varphi {|}^2.\end{array}$$

$$\begin{array}{rl}|{\mathbf{E}}_{\mathrm{p}}(\mathbf{r}){|}^2=& |\nabla {\varphi }_{\mathrm{p}}{|}^2+{\omega }^2|{\mathbf{A}}_{\mathrm{p}}{|}^2\\ & +2\mathrm{Re}\{\mathrm{j}\omega \nabla \cdot ({\mathbf{A}}_{\mathrm{p}}{\varphi }_{\mathrm{p}}^{\ast })\}-2k^2|{\varphi }_{\mathrm{p}}{|}^2.\end{array}$$

$$\begin{array}{r}{\varphi }_{\mathrm{p}}({\mathbf{r}}_1)=\frac{1}{ϵ}{\int }_{\Omega }\varrho ({\mathbf{r}}_2)G_{\mathrm{p}}({\mathbf{r}}_1,{\mathbf{r}}_2)\mathrm{d}{v}_2,\\ {\mathbf{A}}_{\mathrm{p}}({\mathbf{r}}_1)=\mu {\int }_{\Omega }\mathbf{J}({\mathbf{r}}_2)G_{\mathrm{p}}({\mathbf{r}}_1,{\mathbf{r}}_2)\mathrm{d}{v}_2,\end{array}$$

$$G_{\mathrm{p}}({\mathbf{r}}_1,{\mathbf{r}}_2)=\frac{1}{2\mathrm{j}S}\sum _{(m,n)\in \mathcal{P}}\frac{1}{k_{\mathrm{z}mn}}{\mathrm{e}}^{-\mathrm{j}{\mathbf{k}}_{\mathrm{t}mn}\cdot ({\mathbf{\rho }}_1-{\mathbf{\rho }}_2)}{\mathrm{e}}^{-\mathrm{j}{k}_{\mathrm{z}mn}|z_1-z_2|}.$$

$$\begin{array}{rl}& W_{\mathrm{e}}=\frac{ϵ}{4}{\int }_U(|\nabla \varphi {|}^2-k^2|\varphi {|}^2)+{\omega }^2(|\mathbf{A}{|}^2-|{\mathbf{A}}_{\mathrm{p}}{|}^2)\\ & -k^2(|\varphi {|}^2-|{\varphi }_{\mathrm{p}}{|}^2)+2\mathrm{Re}\{\mathrm{j}\omega \nabla \cdot (\mathbf{A}{\varphi }^{\ast }-{\mathbf{A}}_{\mathrm{p}}{\varphi }_{\mathrm{p}}^{\ast })\}\mathrm{d}v.\end{array}$$

$$\begin{array}{rl}{\int }_{U}2\mathrm{Re}& \{\mathrm{j}\omega \nabla \cdot (\mathbf{A}{\varphi }^{\ast }-{\mathbf{A}}_{\mathrm{p}}{\varphi }_{\mathrm{p}}^{\ast })\}\mathrm{d}v\\ & =2\mathrm{Re}\{\mathrm{j}\omega {\int }_{\delta U}\hat{\mathbf{n}}\cdot (\mathbf{A}{\varphi }^{\ast }-{\mathbf{A}}_{\mathrm{p}}{\varphi }_{\mathrm{p}}^{\ast })\mathrm{d}S\}.\end{array}$$

$$\begin{array}{rl}{W}_{\mathrm{e}}=\frac{ϵ}{4}& {\int }_U(|\nabla \varphi {|}^2-k^2|\varphi {|}^2)\\ & +{\omega }^2(|\mathbf{A}{|}^2-|{\mathbf{A}}_{\mathrm{p}}{|}^2)-k^2(|\varphi {|}^2-|{\varphi }_{\mathrm{p}}{|}^2)\mathrm{d}v.\end{array}$$

$$\begin{array}{rl}\frac{\omega }{4}{\int }_U& \mu |\mathbf{H}{|}^2-ϵ|\mathbf{E}{|}^2\mathrm{d}v\\ & =-\frac{1}{4}\mathrm{Im}{\int }_{\Omega }\mathbf{E}(\mathbf{r})\cdot {\mathbf{J}}^{\ast }(\mathbf{r})\mathrm{d}v=\frac{1}{2}\mathrm{Im}{P}_{\mathrm{c}},\end{array}$$

$$W_{\mathrm{m}}=W_{\mathrm{e}}+\frac{1}{2\omega }\mathrm{Im}{P}_{\mathrm{c}}.$$

$$\begin{array}{r}|\nabla \varphi {|}^2-k^2|\varphi {|}^2=\mathrm{Re}\{\nabla \cdot (\varphi \nabla {\varphi }^{\ast })+\varphi {\varrho }^{\ast }/ϵ\}.\end{array}$$

$$\begin{array}{r}{\int }_U\mathrm{Re}\{\nabla \cdot (\varphi \nabla {\varphi }^{\ast })\}\mathrm{d}v={\int }_{\delta U}\mathrm{Re}\{\hat{\mathbf{n}}\cdot (\varphi \nabla {\varphi }^{\ast })\}\mathrm{d}S,\end{array}$$

$$\begin{array}{r}{\int }_U|\nabla \varphi {|}^2-k^2|\varphi {|}^2\mathrm{d}v=\frac{1}{ϵ}\mathrm{Re}\{{\int }_{\Omega }\varphi {\varrho }^{\ast }\mathrm{d}v\}.\end{array}$$

$$\mathrm{j}\omega \varrho =-\nabla \cdot \mathbf{J}$$

$$\varphi ({\mathbf{r}}_1)=\frac{\mathrm{j}}{\omega ϵ}{\int }_{\Omega }G({\mathbf{r}}_1,{\mathbf{r}}_2){\nabla }_2\cdot \mathbf{J}({\mathbf{r}}_2)\mathrm{d}{v}_2.$$

$$\begin{array}{rl}& W_{\mathrm{e},1}=\frac{ϵ}{4}{\int }_U|\nabla \varphi {|}^2-k^2|\varphi {|}^2\mathrm{d}v\\ & =\frac{\mu }{4k^2}\mathrm{Re}{\int }_{\Omega }{\int }_{\Omega }({\nabla }_1\cdot {\mathbf{J}}^{\ast }({\mathbf{r}}_1))G({\mathbf{r}}_1,{\mathbf{r}}_2){\nabla }_2\cdot \mathbf{J}({\mathbf{r}}_2)\mathrm{d}{v}_2\mathrm{d}{v}_1\\ & =⟨\mathbf{J},{\mathcal{W}}_{\mathrm{e},1}\mathbf{J}⟩,\end{array}$$

$$\begin{array}{rl}|\mathbf{A}(\mathbf{r}){|}^2=& {\mu }^2{\int }_{\Omega }{\int }_{\Omega }{G}^{\ast }(\mathbf{r},{\mathbf{r}}_1)G(\mathbf{r},{\mathbf{r}}_2){\mathbf{J}}^{\ast }({\mathbf{r}}_1)\cdot \mathbf{J}({\mathbf{r}}_2)\mathrm{d}{v}_2\mathrm{d}{v}_1.\end{array}$$

$$\begin{array}{rl}|{\mathbf{A}}_{\mathrm{p}}& (\mathbf{r}){|}^2\\ =& {\mu }^2{\int }_{\Omega }{\int }_{\Omega }{G}_{\mathrm{p}}^{\ast }(\mathbf{r},{\mathbf{r}}_1)G_{\mathrm{p}}(\mathbf{r},{\mathbf{r}}_2){\mathbf{J}}^{\ast }({\mathbf{r}}_1)\cdot \mathbf{J}({\mathbf{r}}_2)\mathrm{d}{v}_2\mathrm{d}{v}_1.\end{array}$$

Full text

I Introduction

The Q-factor of an oscillating system is an indirect measure of the width of the system’s resonance. In many cases, the Q-factor accurately predicts the bandwidth at which an external excitation can be efficiently applied to the system [1, 2]. This practically important property has encouraged a development of Q-factor representations for electric circuits [3], electrically small antennas [4, 5, 6, 7, 8, 9], and finite radiating structures [10, 11, 12]. For periodic structures, the existing Q-factor expressions are restricted to a few geometries [13, 14, 15]. This paper proposes the Q-factor formulation for periodic structures with arbitrary element shape.

The Q-factor (also called quality factor, $Q$, radiation Q, antenna Q [1, 8, 6]) is defined to be proportional to the ratio between cycle means of the stored and the loss energies [16]. In electromagnetics, the relation to energies links the Q-factor with distributed quantities of the system, such as radiation modes and electric current density. This relation paved way to many powerful Q-factor-based results in antenna theory.

For small antennas, the fundamental bounds, relating size, bandwidth and efficiency, are important examples of Q-factor results. An essential step in establishing such bounds is the development of a suitable representation for the Q-factor. Chu [5], Wheeler [4] and Harrington [7] obtained classical lower bounds on the Q-factor for canonical geometries by analyzing Q-factor representation, based on circuit models for spherical and cylindrical wave expansions of radiating fields. Collin and Rothschild [6] used the Q-factor expressed directly in reactive (stored) and radiated energies to establish bounds for spherical and cylindrical volumes. The current density representation of the Q-factor, proposed by Vandenbosch [10] (see also the earlier work of Geyi [17]) and later refined by Gustafsson and Jonsson [11, 12], is the key element to establish $Q$-bounds for desired geometries using optimization methods [18, 19, 20, 21, 22, 23, 24]. A practically convenient Q-factor formula, proposed by Yaghjian and Best for antennas with a defined port, is expressed in terms of input impedance and its derivative [1], and has been used for a genetic-algorithm-based optimization of an antenna structure [22].

For larger radiating arrays, which can be accurately described by a unit cell representation, much fewer results exist. Tomasic and Steyskal [13, 14] proposed a lower Q-factor bound for a 1D-array of infinitely long cylinders with either magnetic or electric current sources in free space and over a ground plane. Kwon and Pozar [15] proposed a Q-factor for a 2D-flat infinite array of strip dipoles in a plane located in free space, over a ground plane, and over a dielectric substrate with a ground plane. The key element of their derivation is a Fourier transform of the surface current density on the dipoles; the current density is there assumed to be in one direction only, and a single propagating mode is permitted. Kwon and Pozar also suggest, that their method can be applied to other flat geometries with one-dimensional currents, however, no extension of this idea has been published.

In this paper, we propose and derive the Q-factor expression for three-dimensional radiating arrays with two-dimensional periodicity. The proposed expression is a step towards developing fundamental bounds, obtainable by current density optimization, similarly to non-periodic case [18, 19, 20, 21, 22, 23, 24]. The expression takes into account the exact shape of the radiating elements, and is formulated in terms of electric current density in the array’s unit cell. There is no restriction on the number of propagating Floquet modes that are permitted, and beam steering is allowed. The current density is in general three-dimensional with no restrictions on directions of the currents. The expression involves integration over a finite volume of current density supported within a unit cell, and can thus be rapidly calculated numerically. The Q-factor is calculated as a function of frequency for several antenna arrays. To validate the proposed Q-factor expression, we compare our numerical results with the tuned bandwidth, determined from data computed by a commercial full-wave solver, and also with the input-impedance-$Q$ formula of Yaghjian and Best [1]. Note that our expression, in distinction to the input-impedance-$Q$, does not require a specified input port. Additionally our expression generalizes the result of Kwon and Pozar [15], which we demonstrate by restricting our result by their assumptions and geometry.

The paper is organized as follows. Section II defines the periodic array geometries, considered in this paper, and introduces the stored energies and the radiated power for periodic structures. The periodic Q-factor expression, which is the main theoretical result of this paper, is derived in Sections III–IV: in Section III, we represent the stored energies in terms of electric scalar and magnetic vector potentials. This is convenient for the subsequent derivation of the stored energies in terms of current density in Section IV. The complex and the radiated powers in terms of the current density are included for completeness in the end of Section IV. Section V contains comparison of the proposed Q-factor with the result of Kwon and Pozar [15] for an array of dipoles. In Section VI we approximate the currents by a finite basis, in order to obtain explicit matrices to facilitate numerical calculations of the Q-factor. Numerical examples for a few unit cell antenna-designs are presented in Section VII. We end the paper with conclusions in Section VIII.

II Periodic arrays and stored energies

In this section, we define the family of periodic arrays, that are considered in this paper. The Q-factor, electric and magnetic stored energies, and radiated power are defined in a periodic setting. The periodic Green’s function is stated together with its associated quantities.

We consider a three-dimensional array of lossless perfectly-electric-conducting (PEC) elements of arbitrary geometry. The array is two-dimensionally periodic on a rectangular grid. The array elements are of finite extent and sufficiently regular to support a solution to Maxwell’s equations. The unit cell is enclosed in a box $U_{d}=\{{\bm{r}}=(x,y,z)\in\mathbb{R}^{3}:x\in[0,a],y\in[0,b],z\in[-d/2,d/2],a>0,b>0,d>0\}$, with the box surface denoted as $\delta U_{d}$ and composed of the top surface $\delta U_{\rm top}$ at $z=d/2$, the bottom surface $\delta U_{\rm bottom}$ at $z=-d/2$, and the side walls $\delta U_{\rm sides}$. Fig. 1 shows an example of such a structure. Here the PEC regions are shown in gray, and a unit cell $U_{d}$ is shown as a blue column.

The array elements support an electric current density ${\bm{J}}$, which is periodic up to a phase shift:

[TABLE]

Here ${\bm{r}}\in\mathbb{R}^{3}$, ${\bm{\zeta}}_{mn}=am\hat{{\bm{x}}}+bn\hat{{\bm{y}}}$; $m,n\in\mathbb{Z}$, and ${\bm{k}}_{{\rm t}00}$ is a coordinate-independent transverse phasing wave-vector. The equivalent magnetic sources are neglected in this paper. The inclusion of magnetic sources into consideration can be done similarly to the finite case, see e.g. [12]. The goal is to find an expression for the Q-factor in terms of the unit cell electric current density ${\bm{J}}$. The time convention $\mathrm{e}^{\mathrm{j}\omega t}$ is used but suppressed throughout the paper, and $\omega$ and $t$ denote the angular frequency and time respectively.

The Q factor in terms of stored electric $W_{\rm e}$ and magnetic $W_{\rm m}$ energies per unit cell and the total power $P_{\rm tot}$ (the radiated power and ohmic heating) per unit cell is commonly defined [16, 25] as

[TABLE]

For the lossless case the total power is equal to the radiated power $P_{\rm r}$.

We define the stored energies per unit cell as a difference between the total electric (magnetic) energy and radiated electric (magnetic) energy in each unit cell

[TABLE]

[TABLE]

where $U$ is the box $U_{d}$ as $d\to\infty$, ${\bm{E}}_{\rm p}$ is the component of the ${\bm{E}}$-field, which includes the electric field of propagating Floquet modes only, ${\bm{H}}$ is the magnetic field, and ${\bm{H}}_{\rm p}$ is the corresponding propagation modes component of the ${\bm{H}}$-field. Here we note that the radiated energy is contained in the propagating Floquet modes, as they do not decay at $z\to\infty$. This definition resembles the far-field subtraction approach for single finite-sized radiating elements [11, 12], and for a single propagating mode it is similar to the stored energies definition in [15].

A convenient representation of the radiated power $P_{\rm r}$ can be obtained from the differential form of Poynting’s theorem for harmonic electromagnetic fields, which is [26]

[TABLE]

where ${\bm{E}}$, ${\bm{H}}$ is the total electric and magnetic fields respectively, and the superscript $*$ indicates the complex conjugate. Here $\epsilon$ is the free-space permittivity, and $\mu$ is the free-space permeability. The real part of the volume integral of (5) over the unit cell provides the relation

[TABLE]

where $\hat{{\bm{n}}}$ is the outward-normal unit vector of the boundary surface $\delta U_{d}$ of the unit cell, $\Omega$ is the region of the metal inclusions within $U_{d}$. The right-hand side of (6) can be identified as the radiated power, and thus

[TABLE]

where we denoted the complex quantity in brackets as the complex power $P_{\rm c}$.

To express stored energies explicitly in terms of the unit-cell current density, we utilize the free-space Green’s function for Helmholtz equation. Here we include the periodic up to a phase-shift requirement (1) on solutions from Floquet theorem by a phase correction of the point source in the unit cell:

[TABLE]

Here $\delta({\bm{r}})$ is the Dirac delta-function, $k=\omega/c_{0}$ is the free-space wavenumber; $c_{0}=1/\sqrt[]{\epsilon\mu}$ is the speed of light in free space. The solution of this equation is the 2D-periodic Green’s function. Throughout this paper, the spectral form of the 2D-periodic Green’s function is used [27]

[TABLE]

where ${\bm{k}}_{{\rm t}mn}={\bm{k}}_{{\rm t}00}+2\pi\frac{n}{a}\hat{{\bm{x}}}+2\pi\frac{m}{b}\hat{{\bm{y}}}$, $S=ab$, $k_{{\rm z}mn}=\sqrt[]{k^{2}-{\bm{k}}_{{\rm t}mn}\cdot{\bm{k}}_{{\rm t}mn}}$, and the coordinates are separated into transverse and longitudinal components ${\bm{r}}_{i}={\bm{\rho}}_{i}+\hat{{\bm{z}}}z_{i}$, $i=\{1,2\}$. The transverse phasing wave-vector is ${\bm{k}}_{{\rm t}00}=k\sin\theta_{0}\cos\phi_{0}\hat{{\bm{x}}}+k\sin\theta_{0}\sin\phi_{0}\hat{{\bm{y}}}$ with the polar $\theta_{0}$ and azimuthal $\phi_{0}$ angles, which define a direction of the fundamental-mode plane wave. Whenever $(k^{2}-{\bm{k}}_{{\rm t}mn}\cdot{\bm{k}}_{{\rm t}mn})$ is negative, we choose the negative branch of the square root to evaluate $k_{{\rm z}mn}$, so that $\mathrm{e}^{-\mathrm{j}k_{{\rm z}mn}|z_{1}-z_{2}|}$ is exponentially damped.

III Stored energies in terms of potentials

In this section we derive a potential representation of the stored energies. To obtain the current-density representation of the stored energies, we first express the electric and magnetic stored energies in terms of the electric scalar potential $\phi$ and the magnetic vector potential ${\bm{A}}$. A similar potential representation approach was used in [11, 12] to derive the refined expressions for the stored energies in terms of sources for single antennas.

The potentials are defined by the choice of the Lorenz gauge [26, 28, 29]

[TABLE]

The choice of gauge implies that $(\nabla^{2}+k^{2})\phi({\bm{r}}_{1})=-\varrho({\bm{r}}_{1})/\epsilon_{0}$ and $(\nabla^{2}+k^{2}){\bm{A}}({\bm{r}}_{1})=-\mu_{0}{\bm{J}}({\bm{r}}_{1})$. Here $\varrho$ is the electric charge density. We also recall the connection between the electric field and potentials

[TABLE]

The solution to the Helmholtz equation yields the potentials expressed in terms of the sources

[TABLE]

[TABLE]

Note that $\varrho,{\bm{J}}$ have support in $\Omega$, which is of finite size. Here, the 2D-periodic Green’s function ${G}({\bm{r}}_{1},{\bm{r}}_{2})$ is given in (9). We observe that periodicity up to a phase-shift carries over also to the charge density, due to the continuity equation, i.e. $\varrho({\bm{r}})=\varrho({\bm{r}}+{\bm{\zeta}}_{mn})\mathrm{e}^{\mathrm{j}{\bm{k}}_{{\rm t}00}\cdot{\bm{\zeta}}_{mn}}$.

The first term of the integrand of (3), the total electric energy density, is expressed using (10), (11) and the vector identity ${\bm{A}}\cdot\nabla\phi^{*}=\nabla\cdot({\bm{A}}\phi^{*})-\phi^{*}\nabla\cdot{\bm{A}}$:

[TABLE]

Similar treatment of the electric energy density of the propagating modes gives

[TABLE]

The potentials ${\bm{A}}_{\rm p},{\phi}_{\rm p}$ correspond to the propagating modes, and are defined as follows

[TABLE]

where the propagating-modes part of the periodic Green’s function is given as

[TABLE]

Here the summation is performed over the set of propagating modes $\mathcal{P}=\{(m,n):k^{2}-{\bm{k}}_{{\rm t}mn}\cdot{\bm{k}}_{{\rm t}mn}\geq 0\}$, i.e. the modes for which $k_{{\rm z}mn}$ is purely real.

The substitution of electric energy densities (14) and (15) into the stored energy definition (3) gives the representation of the electric stored energy in terms of potentials

[TABLE]

Here the terms are grouped for their further convenient evaluation. The integral of divergence term in this expression vanishes. To show that, we apply the divergence theorem [26, 28, 29]

[TABLE]

A direct substitution of the Green’s function (9) and its propagating-modes part (17) in potentials (12), (13), (16) asserts that $\{{\bm{A}}({\bm{r}}),\phi({\bm{r}}),{\bm{A}}_{\rm p}({\bm{r}}),{\phi}_{\rm p}({\bm{r}})\}\mathrm{e}^{\mathrm{j}{\bm{k}}_{\rm t00}\cdot{\bm{r}}}$ are periodic with lattice vectors $a\hat{{\bm{x}}}$ and $b\hat{{\bm{y}}}$. Evidently, $({\bm{A}}\phi^{*}-{\bm{A}}_{\rm p}{\phi}_{\rm p}^{*})$ is also periodic, and the integral over $\delta U_{\rm sides}$ in the right-hand side of (19) vanishes as the outward normal vectors $\hat{{\bm{n}}}$ are directed oppositely at the opposite side-walls. Integration over $\delta U_{\rm top}$ and $\delta U_{\rm bottom}$ vanishes as $d\to\infty$, since each term of $({\bm{A}}\phi^{*}-{\bm{A}}_{\rm p}{\phi}_{\rm p}^{*})$ contains at least one evanescent-mode factor, exponentially decaying with growth of $d$. Thus, the electric stored energy in terms of potentials is

[TABLE]

A similar approach as above can be used to find the potential representation of the magnetic stored energy (4). Here, however, we instead integrate the expression (5) of Poynting’s theorem over the unit cell volume and take the imaginary part to obtain the following relation for the imaginary part of the complex power $P_{\rm c}$

[TABLE]

where the imaginary part of the integral of $\nabla\cdot({\bm{E}}\times{\bm{H}}^{*})$ is vanishing. Indeed, application of the divergence theorem turns the volume integral into surface integral of ${\bm{E}}\times{\bm{H}}^{*}$ over the surface of $U$. Due to the periodicity of ${\bm{E}}\times{\bm{H}}^{*}$, integral over side walls $\delta U_{\rm sides}$ vanishes. The integral over unit cell’s top $\delta U_{\rm top}$ and bottom $\delta U_{\rm bottom}$ surfaces is purely real, which is shown in Appendix A. Combining (21) with the plane-wave behavior $\int_{U_{d}}(|{\bm{H}}_{\rm p}|^{2}-|{\bm{E}}_{\rm p}|^{2}/\eta^{2})\mathrm{d}v\rightarrow 0$ as $d\to\infty$ and the definitions of the stored energies (3) and (4), we obtain the connection between the electric and magnetic energies

[TABLE]

Here, $\eta=\sqrt[]{\mu/\epsilon}$ is the free space impedance for the propagating modes. The complex power, $P_{\rm c}$, is formulated in terms of sources in the end of Section IV.

IV Stored energies in terms of sources

In the potential formulation of the electric stored energy (20), we have combined the integrand terms in pairs. The pairs are chosen such that the integral of each pair is finite. In this section, we derive from these pairs the representation of the stored energy in terms of the electric current density ${\bm{J}}$.

We begin with the evaluation of the integrand’s first pair, $(|\nabla\phi|^{2}-k^{2}|\phi|^{2})$, which is a weak form of a Helmholtz equation. Using the vector identity $|\nabla\phi|^{2}=\nabla\cdot(\phi\nabla\phi^{*})-\phi\Delta\phi^{*}$ and the Helmholtz equation for the scalar potential, we obtain that

[TABLE]

The divergence term in this expression vanish after an integration over $U$, to see this note that the divergence theorem gives

[TABLE]

where $\delta U$ is the boundary of $U$. Periodicity of $(\phi\nabla\phi^{*})$ follows from periodicity of $\phi({\bm{r}})\mathrm{e}^{\mathrm{j}{\bm{k}}_{\rm t00}\cdot{\bm{r}}}$. This, in turn, means that the integrals over the opposite side walls of $U$ cancel each other. Additionally, the integral over the top and the bottom of $U$ vanishes, since $(\phi\nabla\phi^{*})$ is purely imaginary for large $|z|$. For the details of this derivation, see Appendix B. Thus we have

[TABLE]

A combination of the Green’s function representation of the scalar potential (12) and the continuity equation

[TABLE]

for current density yields the expression for the scalar potential in terms of currents

[TABLE]

The continuity equation (26) gives us the electric current representation of the charge density $\varrho$. Substitution of $\varrho$ and $\phi$ into (25) yields the source representation of the first pair as a quadratic form

[TABLE]

where the scalar product $\langle{\bm{J}},C({\bm{J}})\rangle=\int_{\Omega}{\bm{J}}^{*}\cdot C({\bm{J}})\mathrm{d}v_{1}$ and $\mathcal{W}_{\rm e,1}$ is the linear operator acting on ${\bm{J}}$. The above divergence form (28) is equivalent to the operator $\mathcal{W}_{\rm e,1}$ acting on the set of current densities ${\bm{J}}$ considered here.

To evaluate the second pair in the stored energy representation (20), we use (13) to express the vector potential squared magnitude

[TABLE]

We can similarly express the contribution due to the vector potential of the propagating modes. We utilize (16) to get

[TABLE]

Direct substitution of (29) and (30) into the second pair of (20) together with an interchange of integration order, see Appendix C, yields

[TABLE]

where the kernel is

[TABLE]

In order to simplify $g$, we start with utilizing the structure of the integrand and note that the subtraction of sums in the respective Green’s functions reduce the integrand to summation over $(m,n,p,q)\in\mathbb{Z}^{4}\setminus\mathcal{P}^{2}$ of terms of the form

[TABLE]

Furthermore, we utilize the Floquet modes orthogonality under integration in the $xy$-plane over the unit-cell (71) to obtain a sum over $\mathbb{Z}^{2}\setminus\mathcal{P}$. The subsequent integration in $z$ reduces the kernel to

[TABLE]

as shown in Appendix C. The rather remarkable reduction in complexity in going from (32) to (34) facilitates an effective numerical computation of $g$, in particular with the exponential decay in $m,n$ of the terms in the sum. From the conjugate symmetry $g({\bm{r}}_{1},{\bm{r}}_{2})=g^{*}({\bm{r}}_{2},{\bm{r}}_{1})$ it follows immediately that $W_{\rm em,1}$ is purely real.

When $z_{1}\neq z_{2}$, the convergence of the double sum in (34) is ensured by the exponential damping of the terms with growth of both $m$ and $n$. For the special case with flat thin antennas parallel with the periodicity plane, we have $z_{1}=z_{2}$, and the kernel reduces to

[TABLE]

where the terms decay at least with the cubic power of both $m$ and $n$. We observe that both (34) and (35) are independent of the choice of the coordinate system. This is markedly different from the finite support antenna case in which the far-field subtraction results in a weakly coordinate dependent stored electric energy [30].

We apply a similar procedure to the last pair in the electric stored energy representation (20), which involves the scalar potential and its propagating-mode part. Substitution of (27) yields that

[TABLE]

In the same manner we can find the expression for the scalar potential contribution due to propagating modes

[TABLE]

Combining the latter two expressions yields the third pair in (20) which reduces to the same kernel $g$ in (34) above. Thus

[TABLE]

From the conjugate symmetry of the kernel $g$ it follows that $W_{\rm em,2}$ is real. Finally, we combine (28), (31) and (38) to write the current representation of the electric stored energy

[TABLE]

To find the current density representation of the radiated power (7) and the magnetic stored energy via (22), we need to evaluate the complex power

[TABLE]

We proceed with evaluating the integrand by substituting the potential form of the electric field (11) and consequently the potentials in terms of currents (13) and (27)

[TABLE]

By applying the vector identity

[TABLE]

and the fact that

[TABLE]

we rewrite the complex power as

[TABLE]

Note that the kernel of the operator $\mathcal{Z}$ is the same as the electric field integral equation (EFIE) kernel in the method of moments (MoM) for EM-simulations [31]. We notice that the imaginary part of the first term of (44) is $-2\omega W_{\rm e,1}$. Equation (22) can thus be used to obtain the magnetic stored energy

[TABLE]

where

[TABLE]

Thus

[TABLE]

For the case of PEC structures, only the currents on the surfaces of the elements exist. The adjustment of the expressions for stored energies and complex power to surface currents is straightforwardly done by replacing the volume current density ${\bm{J}}$ with the surface current density ${\bm{J}}_{\rm S}$ and the volume elements $\mathrm{d}v$ with the surface elements $\mathrm{d}S$ in the integrals in (28) (31), (38), (44) and (46).

V Array of dipoles: comparison with Kwon and Pozar

In this section, we examine our proposed expressions for stored energies, radiated power and the Q-factor by considering a particular array element shape – narrow, infinitely thin rectangular dipoles situated in the array plane. This case was considered earlier by Kwon and Pozar [15] for a single propagating mode. There, they assume only $x$-directed currents on rectangular dipoles contained in a plane, parallel to the $xy$-plane:

[TABLE]

with $f({\bm{r}})=f(x,y)$ an electric current density profile. Below, we illustrate that our Q-factor expression, under these assumptions on the current density and number of propagating modes, coincides with the result of [15].

Using vector calculus identities, we can rephrase (28) as

[TABLE]

where $\overline{\overline{\nabla}}_{1}\nabla_{1}{G}({\bm{r}}_{1},{\bm{r}}_{2})$ should be interpreted as a Jacobian of $\nabla_{1}{G}({\bm{r}}_{1},{\bm{r}}_{2})$. The substitution of (48) for the surface current densities yields

[TABLE]

By inserting the Green’s function (9) with $|z_{1}-z_{2}|=0$, and taking the second order derivative, we obtain

[TABLE]

where $k_{{\rm x}m}={\bm{k}}_{{\rm t}mn}\cdot\hat{{\bm{x}}}$, and $F_{mn}=\int_{\Omega}f({\bm{r}}_{1})\mathrm{e}^{\mathrm{j}{\bm{k}}_{{\rm t}mn}\cdot{\bm{\rho}}}\mathrm{d}S$ are the Fourier transform coefficients. We recall that $k_{{\rm z}mn}\in\mathbb{R}$ for propagating modes, and hence the propagating-modes terms $(m,n)\in\mathcal{P}$ in the sum are purely imaginary and do not contribute to $W_{\rm e,1}$. For the evanescent modes, $k_{{\rm z}mn}=-\mathrm{j}|k_{{\rm z}mn}|$, which follows from the branch choice. This transforms (51) into

[TABLE]

where Kwon and Pozar assume [15] that only one mode propagates, i.e. $\mathcal{P}=\{(0,0)\}$. A treatment of (38) with similar steps yields

[TABLE]

By identifying the Fourier coefficients as done in (51), we easily obtain the other two contributions (28), (46) of the stored energies

[TABLE]

[TABLE]

Combining the stored-energies terms according to (39) and (47), we obtain

[TABLE]

[TABLE]

where $k_{{\rm y}n}={\bm{k}}_{{\rm t}mn}\cdot\hat{{\bm{y}}}$.

The radiated power is found similarly by combining (7) and (44) with the derivation steps used in this section

[TABLE]

The expressions (56)-(58) are exactly coinciding with the result of Kwon and Pozar [15]. Note that a difference by a $1/2$ factor, caused by a difference in definitions, is both in the stored energies and the radiated power, and is eliminated when the Q-factor is computed.

VI Matrix formulation of stored energies and radiated power

In this section we discuss how the source representations of the stored energies (28), (31), (38), (47) and the complex power (44), together with the EFIE-based MoM solver, can be used to compute the Q-factor for a 2D periodic array of an arbitrary shape. We start with introducing a set of basis functions $\{{\bm{f}}_{n}({\bm{r}})\}$, which approximates the surface current density ${\bm{J}}_{\rm S}$ on $\Omega$

[TABLE]

where ${\rm\bf I}=(I_{1},I_{2},...,I_{N})^{\rm T}$ is a vector of the current density coefficients, $N$ is the number of the basis functions. We proceed with substituting this basis expansion in the components of the electric stored energies (28), (31) and (38). This yields the matrix representation of the electric stored energy (39)

[TABLE]

where the kernel matrix elements are given by ${\rm\bf W}_{\rm e}^{(m,n)}=\langle{\bm{f}}_{m},(\mathcal{W}_{\rm e,1}+\mathcal{W}_{\rm em,1}-\mathcal{W}_{\rm em,2}){\bm{f}}_{n}\rangle$, and the superscript $(.)^{\rm H}$ denotes the Hermitian transpose. We continue in this fashion obtaining the matrix form of the complex power

[TABLE]

where the impedance matrix elements are given as ${\rm\bf Z}^{(m,n)}=\langle{\bm{f}}_{m},\mathcal{Z}{\bm{f}}_{n}\rangle$. In accordance with (22), we find the kernel matrix for the magnetic stored energy

[TABLE]

Combination of the above matrix representations with (2) and (7) gives the matrix forms for the electric and magnetic Q-factors, and the total Q-factor respectively

[TABLE]

To calculate the $Q$-factor for a given geometry and excitation, the current density coefficients ${\rm\bf I}$ can be computed e.g. by solving the EFIE by MoM, as done in the next section. When the goal is to find the optimal-bandwidth currents on a given spatial support, the Q-factor (63) can be minimized with respect to the current density coefficients ${\rm\bf I}$, see finite antenna examples in e.g.  [18, 19, 24].

VII Numerical results

The goal of this section is to illustrate the practical use of the proposed Q-factor expression. We numerically validate the Q-factor representation by comparison with the actual tuned-impedance bandwidth, and the input-impedance Q-factor by Yaghjian and Best [1]. Additionally, we compare our Q-factor expression with the Q-factor from [15] for a simple dipole array. For the case of narrow-band arrays, we illustrate that all Q-factors give similar results and that they predict the fractional bandwidth.

We start with recalling how the actual tuned bandwidth is obtained and converted into the equivalent-tuned-bandwidth Q-factor $Q_{\rm B}$; and also how the input-impedance Q-factor $Q_{\rm Z}$ is calculated.

The approach to find $Q_{\rm B}$ follows Yaghjian and Best [1]: given an input impedance $Z(\omega)$ of an unmatched antenna, we tune it at the angular frequency $\omega_{0}$ with a series element of impedance $\mathrm{j}X_{\rm s}(\omega)$ such that $\mathop{\mathrm{Im}}Z(\omega_{0})+X_{\rm s}(\omega_{0})=0$. The series element can either be an inductance $L$ with $X_{\rm s}(\omega)=\omega L$ or a capacitor $C$ with $X_{\rm s}(\omega)=-1/(\omega C)$ depending on a sign of $\mathop{\mathrm{Im}}Z$. The overall tuned input impedance is thus $Z_{0}(\omega)=Z(\omega)+\mathrm{j}X_{\rm s}(\omega)$, and the characteristic impedance $Z_{\rm ch}$ is set as a constant equal to $\mathop{\mathrm{Re}}Z(\omega_{0})$. The reflection coefficient is then calculated as $\Gamma=(Z_{0}-Z_{\rm ch})/(Z_{0}+Z_{\rm ch})$ and the fractional bandwidth $B$ is estimated at the threshold level $\Gamma_{0}$. Finally the fractional bandwidth is converted to an equivalent Q-factor [1] by

[TABLE]

The Yaghjian and Best formula provides an estimate of the Q-factor, matched at a given angular frequency, and involves the antenna’s input impedance $(R+\mathrm{j}X)$

[TABLE]

where $(.)^{\prime}$ denotes the angular frequency derivative.

Below we compare the here derived Q-factor with the $Q_{\rm B}$ and $Q_{\rm Z}$, computed from an input impedance given both by our in-house MoM code and CST MW Studio simulation in the frequency domain. The implementation of the proposed Q-factor expression is based on and is an extension of our in-house MoM code [32], where RWG basis functions [31] were used. In all the examples, both CST and the in-house-MoM-code models have voltage-gap excitation. We used Ewald’s method [27] to efficiently evaluate the Green’s function (9), for alternative acceleration methods see e.g. [33]. The convergence of our MoM results are validated by comparison of input impedances with solutions obtained from CST MW Studio 2017.

VII-A Dipole array

We first consider an infinite planar array of strip dipoles. This example has been treated in Section V, where we analytically imposed the assumptions of [15] by Kwon and Pozar to compare our expressions with their result. Here, however, we do not apply the assumptions of [15] when calculating Q-factor with our expression (63).

The dependence of the Q-factor on the phasing angle (beam steering) is presented in Fig. 2a for an array of center-fed dipoles, length $l=0.45\lambda$, width $w=0.02\lambda$, and periodicity $p=0.5\lambda$, see the inset in Fig. 2b. The dipoles are aligned along the $x$-axis; the scan angle $\theta_{0}$ is given in the E-plane ($\phi_{0}=0^{\circ}$) and the H-plane ($\phi_{0}=90^{\circ}$), where $\phi_{0}$ is counted from the positive $x$-semiaxis. The here derived Q-factor in expression (63) is given by solid lines, and the Q-factors by our implementation of Kwon and Pozar method [15] are given by dashed lines. Good agreement between the two methods is observed, and the small discrepancies are the outcome of different assumptions on current density in the methods. In [15] a single basis function with piece-wise sinusoidal distribution along the dipole’s length, and uniform distribution along the dipole’s width were used, and the current’s direction was strictly along the dipole’s length. To calculate the Q-factor with our expressions, we used an EFIE current solution for the model with 239 RWG basis functions and a voltage gap excitation. Our model allows a non-uniform current distribution along the dipole’s width as well as arbitrary current direction within the dipole’s surface.

Fig. 2b shows the Q-factor as a function of the electric length $kl$ for an array of center-fed strip dipoles. The dipole element has length $l$ and width $w=l/40$, the periodicity is $p=1.2l$ in both directions of the array plane. Here a phasing angle $(\theta_{0},\phi_{0})=(0,0)$ is assumed, corresponding to a $z$-directed radiation. The voltage gap excitation is used. The electric $Q_{\rm e}$ and magnetic $Q_{\rm m}$ Q-factors in (63) are shown by dashed red and blue curves respectively, where a model with 239 RWG functions was used. The equivalent-tuned-bandwidth $Q_{\rm B}$, shown by magenta dotted line with markers, is calculated, for reflection threshold level $\Gamma_{0}=-10$ dB, from the voltage-gap input impedance, solved by CST solver. The $Q_{\rm Z}$ from the input-impedance solutions by CST and by our in-house MoM code are shown with solid green and dotted black curves respectively. We observe a fine agreement between all the curves at the frequencies below the first grating lobe, which occurs at $kl\simeq 5.2$. Above this electrical length, all methods capture the grating-lobe peaks of the Q-factor with a good agreement, however differences are observed between the peaks. The stored-energy Q-factor $\max(Q_{\rm e},Q_{\rm m})$ estimates the -10 dB bandwidth as expressed by $Q_{\rm B}$ better than $Q_{\rm Z}$ in the first half of that interval up till $kl\simeq 6.3$; above this value $Q_{\rm Z}$ appears to be closer to $Q_{\rm B}$. Note that the array is sufficiently wideband between the first two grating lobes ($Q\leq 5$), and the Q-factor model might not be reliable there [1]. At the grating lobe frequencies we observe that there are singularities of Q-factor. These singularities are related to the periodic resonance of the array. At low frequencies, we see that the stored electric energy dominates, which is typical for dipole-type antennas. The minimum of Q is located at $kl\simeq 3$, where the electric and the magnetic stored energies are equal. This electrical length corresponds to $l\approx 0.48\lambda$, where $\lambda$ is the wavelength. This agrees with the length of a single dipole required to achieve its first self-resonance [34].

VII-B Rectangular loop array

The next example is a rectangular loop (also called strip folded dipole) array. Fig. 3 depicts the Q-factors of such array with the strip width $w=l/10$ and periodicity $p=1.2l$. Each loop is fed at the middle of one of its shorter sides, and a phasing angle $(\theta_{0},\phi_{0})=(0,0)$ is considered. The loop model for the in-house MoM code consists of 798 RWG basis functions. The colours of the Q-factor curves are the same as in Fig. 2b in this and all the following examples.

A good agreement between the Q-factors is observed for frequencies below the grating lobe at $kl\simeq 5.2$. The magnetic stored energy dominates at the low frequencies, which is common for the loop-type antennas. The minimal value of the Q-factor below the grating lobe occurs at $kl\simeq 2.9$, where the magnetic and the electric stored energies are equal. For this electric length, the current density distribution is depicted in Fig. 3b (upper), where the normalized magnitude of current density is shown by a colourmap (from pale yellow for low magnitude to dark red for high magnitude), and the instantaneous direction of the current is represented by the blue arrows. We see that the length of the loop is approximately equal to one wavelength, which indicates a self-resonance of the loop. The loop at this electric length can be seen as two bent half-wavelength dipoles, fed in phase. Such a configuration has a maximum of its radiation pattern in the broadside direction [34]. Thus the surface waves in the array plane are low.

As a contrast, we observe an additional peak of the Q-factor at $kl\simeq 4.6$, below the grating lobe. For this electric length, the current density distribution is depicted in Fig. 3b (lower), and here two wavelengths fit in the loop. The loop can be seen as two symmetric pairs of half-wavelength dipoles, with a phase shift $180^{\circ}$ between dipoles in each pair. Such a configuration has a null of the radiation pattern in the broadside direction, and a maximum in the direction along the array plane. In this regime, the surface waves dominate over the radiated waves, which yields a high peak value of the Q-factor.

A remark should be made regarding the nature of the grating lobe peaks in this and all the other examples; numerically those peaks can be made arbitrarily large in amplitude by refinement of the frequency sample grid in the neighbourhood of those peaks. This indicates a singular behavior of the peaks, and also explains the differences in the peak amplitudes in numerical results by different methods.

Above the grating lobe we note that $Q_{\rm B}$ is below the prediction of $\max(Q_{\rm e},Q_{\rm m})$ and $Q_{\rm Z}$ and that the array is rather wideband $Q_{\rm B}\simeq 1-2$ in this region.

VII-C Capped dipole array

The Q-factor comparison for an array of capped dipoles is presented in Fig. 4 as a function of the electric size $kl$ of array elements. The element parameters, as illustrated by an inset, are $p=1.2l$, and $h=l/5$, and the strip width $w=l/10$. To compute the Q-factor, we used an element model with 931 RWG elements. The curves agree well in the whole range of $kl$ below the first grating lobe, excluding the region around $kl\simeq 5$. At this electrical size, the approximate length of the dipole’s section between the ends A and B (see the inset in Fig. 4) is one wavelength $\lambda$. The length of the section between points A and C (or B and C) is $\sim 1.5\lambda$. Thus, a multiple resonance of the structure occurs, which induces a dip of the impedance Q-factor $Q_{\rm Z}$. This effect is discussed for a finite case in [35, 36, 30], where there is a deviation between $Q_{\rm B}$ and $Q_{\rm Z}$. We observe a good agreement of all the methods as well above the first grating lobe, excluding a region around $kl\simeq 7.8$, where there is a slight difference between all the methods. As a contrast, recall Fig. 2b, where the methods give different results in the range of $kl$ above the grating lobe.

VII-D Array of mutually coupled dipoles

The Q-factors for an array of mutually coupled dipoles are shown in Fig. 5. The unit cell consists of two dipole elements, a parasitic dipole in the plane $z=0$ and a driven element in the plane $z=l/2$. The driven element is fed in its center. Each dipole is comprised of 239 RWG elements and the lengths and widths of the dipoles are the same. All the curves agree well in the whole range of $kl$ except of the interval $kl=[2.5,3.5]$. There, due to the multiple resonance of the structure, the $Q_{\rm Z}$ has a sharp dip [35], [36], [30]. We observe that the equivalent tuned impedance Q-factor, with increase of the reflection coefficient threshold $\Gamma_{0}$ converges towards $Q_{\rm Z}$. Indeed, the lower $\Gamma_{0}$ is, the more it highlights the single resonance caused by the matching element. We observe that for this particular array, the proposed stored-energy Q-factor predicts bandwidth better for less severe thresholds $\Gamma_{0}\geq-10$dB, while the impedance-based $Q_{\rm Z}$ is more accurate when $\Gamma_{0}<-10$dB.

VII-E Conical spiral array

A conical Archimedian spiral array is here used as an example of a non-flat structure having electric current density with nonzero $z$-component. The parametric equations, governing the Archimedian spiral’s center line is

[TABLE]

Here we let $\alpha=0.3$, the strip width $w=d/20$ and the period $p=1.2d$. The spiral geometry and the computed Q-factors for $N=0.25$ are depicted in Fig. 7. We observe a fine agreement between the methods below the grating lobe. The minimum of Q-factor is at electric size $kd=2.3$, which corresponds to the spiral’s length equal approximately to half wavelength. Above the grating lobe, all the methods give slightly different Q-factor values.

VII-F Discussion

In all presented examples the proposed stored-energy-based Q-factor agrees well overall at the frequencies below the first grating lobe. When the Q-factor is low ($Q\leq 5$), there exist cases with disagreement between the proposed Q-factor and the fractional bandwidth, see region around $kl\simeq 3$ in Fig. 5. In such regions, a sharp dip of Q-factor occurs, which typically related to multi-resonance behaviour [35, 36, 30].

VIII Conclusions

In this paper, the Q-factor expression is derived for lossless two-dimensionally-periodic arrays in free space. The proposed expression is given in terms of the electric currents only, and thus it accounts for an exact shape of an array element. The stored energies and radiated power, that enter the Q-factor formula, are expressed in the similar manner to the EFIE, typically used in MoM solvers, and hence the numerical implementation requires marginal modifications of an existing MoM code.

The numerical investigation indicates, that for narrow-band ($Q>5$) arrays the expression accurately predicts tuned fractional bandwidth, and agrees with input-impedance-$Q$ proposed in [1]. The proposed representation permits multiple propagating Floquet modes, and our simulations in this regime (i.e. at frequencies above the first grating lobe) demonstrate reasonable accuracy of bandwidth prediction. The accuracy of such prediction above the first gratinglobe appears to be better for dipole-type structures than for loop-type arrays, which is a phenomenon to consider in the future. In this region antenna elements tend to become electrically large with multiple resonances and the Q-factor description ceases to be valid.

The main advantages of the proposed Q-factor representation are: (i) it takes into account the exact shape of an array element and its current distribution, (ii) mathematically it has the same form as Q-factor expressions of finite radiators in [10, 12], which is suitable for establishing fundamental bounds as done for finite case in e.g. [18, 22, 23], (iii) it does not require a specified input port and thus applicable to larger class of problems, (iv) numerical computation of stored-energies and radiated-power kernels is not more difficult than computation of Z-matrix in MoM.

To conclude, the proposed expression is fully consistent with and generalizes the previous work [15], by extending the set of geometries and constraints for the Q-factor to include arbitrary-shaped metal inclusions on a periodic grid with arbitrary number of propagating modes.

Acknowledgment

The authors thank S. Amaya Maldonado for his earlier work on the in-house periodic MoM solver. The authors would like to acknowledge the support of the Swedish foundation for strategic research (SSF) under the grant ‘Complex analysis and convex optimization for electromagnetic design’ and the support of Vinnova center under the ChaseOn/iAA project.

Appendix A Vanishing of $\mathop{\mathrm{Im}}\int_{\delta U_{\rm top}}{\bm{E}}\times{\bm{H}}^{*}\cdot\hat{{\bm{n}}}\mathrm{d}S$

To show that $\mathop{\mathrm{Im}}\int_{\delta U_{\rm top}}{\bm{E}}\times{\bm{H}}^{*}\cdot\hat{{\bm{n}}}\mathrm{d}S$ vanishes we write out explicitly the electric ${\bm{E}}$ and magnetic ${\bm{H}}$ fields in terms of currents. The electric field (11) contains the scalar potential and vector potential terms. We start with evaluating the scalar-potential term at ${\bm{r}}_{2}\in\delta U_{\rm top}$ which is located at $z_{2}=d/2$. We notice that the evanescent-modes terms $(m,n)\in\mathbb{Z}^{2}\setminus\mathcal{P}$ are decaying exponentially with growth of $d$ as ${\bm{r}}_{1}\in\delta U_{\rm top}$, and the scalar-potential term is thus

[TABLE]

Here, $k_{\rm z}=\min_{(m,n)\in\mathbb{Z}^{2}\setminus\mathcal{P}}|k_{{\rm z}mn}|$ is the smallest longitudinal wavenumber among evanescent modes, and ${\bm{k}}_{mn}={\bm{k}}_{{\rm t}mn}+\hat{{\bm{z}}}k_{{\rm z}mn}$. Using the vector identity for a divergence of vector and scalar functions product, we obtain

[TABLE]

where $\mathbfcal{J}_{mn}=\int_{\Omega}\mathrm{e}^{\mathrm{j}{\bm{k}}_{mn}\cdot{\bm{r}}_{2}}{\bm{J}}({\bm{r}}_{2})\mathrm{d}v_{2}$. We express the vector potential (13) component in a similar form

[TABLE]

The magnetic field vector is treated in a similar manner, employing the vector identity for a curl of vector and scalar functions product

[TABLE]

The orthogonality of the Floquet modes, instrumental in the next step, follows from the integral over the top surface $\delta U_{\rm top}$ of the unit cell

[TABLE]

Note that this result is independent of the $z$-coordinate and is valid for any surface, which is a crossection of $U$, parallel to the $xy$-plane.

Orthogonality of the modes allows us to reduce the order of summation in the following integral, where the electric and magnetic fields are substituted using expressions (68), (69), (70)

[TABLE]

where the vector triple product identity was used in the last transition. It is evident that the expression is purely real, and thus $\mathop{\mathrm{Im}}\int_{\delta U_{\rm top}}{\bm{E}}\times{\bm{H}}^{*}\cdot\hat{{\bm{n}}}\mathrm{d}S=0$. Similarly $\mathop{\mathrm{Im}}\int_{\delta U_{\rm bottom}}{\bm{E}}\times{\bm{H}}^{*}\cdot\hat{{\bm{n}}}\mathrm{d}S=0$ for the bottom surface of $U_{d}$ as $d\to\infty$

Appendix B Vanishing of $\mathop{\mathrm{Re}}\int_{U}\{\nabla\cdot(\phi\nabla\phi^{*})\}\mathrm{d}v$

We apply the divergence theorem to the divergence term:

[TABLE]

We note that $\{(\nabla\phi({\bm{r}}),\phi({\bm{r}})\}\mathrm{e}^{\mathrm{j}{\bm{k}}_{\rm t00}\cdot{\bm{r}}})$ are periodic with lattice vectors $a\hat{{\bm{x}}}$ and $b\hat{{\bm{y}}}$, and thus $\phi\nabla\phi^{*}$ is periodic. Hence the integration over $\delta U_{\rm sides}$ at the right-hand side of (73) vanishes.

At $\delta U_{\rm top}$, the normal component of the gradient factor (68) is

[TABLE]

Similarly the scalar potential (27) is

[TABLE]

Combination of the last two equations with the orthogonality of modes (71) gives

[TABLE]

which vanish as $d\to\infty$. Similarly the integral over $U_{\rm bottom}$ vanishes, and thus the integral over the whole $\delta U$ vanishes.

Appendix C Order of integration and summation

In the first line of expression (31), the integration and summation order is as follows

[TABLE]

where $G_{mn}$ is a spectral term of the periodic Green’s function (9). For efficient and convenient evaluation of such expression, it is desirable to move the integration over the unit cell $U$ under both the integral over $\Omega\times\Omega$ and the infinite sum. Performing unit-cell-volume integration term-by-term significantly simplifies the expression, however interchange of the orders of operations (integration and infinite summation) may in general affect the result, and thus must be justified. We here consider (77) in the case where the frequency is in the open set defined by removing the grating-lobe frequencies, that is $|k_{zmn}|>0$. At the grating-lobe frequencies it is well known that the Floquet-mode representation of the Green’s function diverges.

The first step is to show that we are able to interchange the infinite summation and integration over $U\times\Omega\times\Omega$. We start with showing that the interchange is valid for the integration over $[U\times\Omega\times\Omega]_{\delta}=\{({\bm{r}},{\bm{r}}_{1},{\bm{r}}_{2}):|z-z_{1}|\geq\delta,|z-z_{2}|\geq\delta\}$. To connect with the original problem, we note that $\lim_{\delta\to 0}\int_{[U\times\Omega\times\Omega]_{\delta}}=\int_{U\times\Omega\times\Omega}$. As the summations over $(m,n)$ and $(p,q)$ are independent of each other, it suffices to show $\int_{U\times\Omega^{2}}\underset{(m,n)\in\mathbb{Z}^{2}\setminus\mathcal{P}}{\sum}{G}^{*}_{pq}{G}_{mn}{\bm{J}}_{1}\cdot{\bm{J}}^{*}_{2}\mathrm{d}v_{1}v_{2}v=\underset{(m,n)\in\mathbb{Z}^{2}\setminus\mathcal{P}}{\sum}\int_{U\times\Omega^{2}}{G}^{*}_{pq}{G}_{mn}{\bm{J}}_{1}\cdot{\bm{J}}^{*}_{2}\mathrm{d}v_{1}v_{2}v$ for fixed $(p,q)$. Thus, we consider closer the convergence of the sum in $(m,n)$ by focusing in each term on the factors, dependent on $m,n$.

Define the sequence of functions

[TABLE]

Here, we introduced the multi-index ${\bm{n}}=(m,n)$, $|{\bm{n}}|^{2}=m^{2}+n^{2}$. Note that there always exists $n_{0}$ such that the longitudinal wavenumber can be estimated as $|k_{{\rm z}mn}|\geq\alpha|{\bm{n}}|$ with some $\alpha>0$ whenever $|{\bm{n}}|\geq n_{0}$. For the purposes of the proof, it is sufficient to consider $N\geq n_{0}$ such that all ${\bm{n}}$ with $|{\bm{n}}|\geq n_{0}$ are not in $\mathcal{P}$.

We show that $\{f_{N}\}$ is uniformly convergent on $F=\{{\bm{r}}_{0}\in\mathbb{R}^{3}:x_{0}\in[0,a],y_{0}\in[0,b],z_{0}\in\mathbb{R},|z_{0}|\geq\delta\}$ by estimating

[TABLE]

where $M^{\prime}=\min(M,N)\geq n_{0}$ and $N^{\prime}=\max(M,N)$. This result is independent of ${\bm{\rho}}_{0}$ and $z_{0}$, and can be made arbitrarily small by choosing sufficiently large $M^{\prime}$. Thus, by the Cauchy criterion for uniform convergence, see e.g. [37, Thm. 7.8], $\{f_{N}\}$ converges uniformly on $F$. The uniform convergence allows us to interchange the integration and summation $\int_{A}f_{\infty}\mathrm{d}v=\lim_{N\to\infty}\int_{A}f_{N}\mathrm{d}v$ [37, Thm. 7.16]. Applying this twice yields

[TABLE]

Next, we show that for each term in the sum above, the order of integration is interchangeable, that is

[TABLE]

By Fubini’s theorem [38, Thm. 2.16.4], for any $h$ measurable function on $X\times Y$, if either $\int\int|h|\mathrm{d}\mu\mathrm{d}\nu<\infty$ or $\int\int|h|\mathrm{d}\nu\mathrm{d}\mu<\infty$, the result of an integral of $h$ over $X\times Y$ is independent of the order in which integrals are taken.

To see that it applies to our case, we first note that the function

[TABLE]

is bounded by a constant, independent of ${\bm{r}}_{1},{\bm{r}}_{2}$. Thus

[TABLE]

is finite for absolutely integrable currents ${\bm{J}}\in L_{1}$. Such currents are measurable since they are in $L_{1}$, and $G_{mn}G_{pq}^{*}$ is continuous and hence measurable on $[\Omega\times\Omega\times U]_{\delta}$. Thus, ${G}_{mn}({\bm{r}},{\bm{r}}_{1}){G}_{pq}^{*}({\bm{r}},{\bm{r}}_{2}){\bm{J}}({\bm{r}}_{1})\cdot{\bm{J}}^{*}({\bm{r}}_{2})$ is measurable, since products of measurable functions are measurable [37, Thm 11.18]. This satisfies the conditions of Fubini theorem, and hence (81) holds true.

By applying the convergence argument of (79), we can interchange $\int_{[\Omega\times\Omega\times U]_{\delta}}$ and summation. Finally, by taking limits in $\delta$, we obtain

[TABLE]

The next step is to reduce the kernel $g$, i.e. the integral over $U$ and the sum over $G_{mn}G_{pq}^{*}$ in  (84) as compactly as possible. We observe that $g$ in (32) is equivalent with

[TABLE]

Again, the interchange of the orders of operations (integration and infinite summation) may in general affect the result, and thus must be justified. We first show that the interchange is valid for integration over $U_{\delta}=\{{\bm{r}}=(x,y,z)\in\mathbb{R}^{3}:x\in[0,a],y\in[0,b],z\in\mathbb{R}\setminus(B_{1}\cup B_{2})\}$, where $B_{i}=\{z\in\mathbb{R}:|z-z_{i}|>\delta\}$, and subsequently demonstrate that the assertion holds also in the limit $\delta\to 0$.

From the uniform convergence argument (79) on $U_{\delta}$ it follows that

[TABLE]

We integrate over the cross-section, utilizing the orthogonality of modes (71) to reduce the sum

[TABLE]

It remains to show that the limit and the infinite sum are interchangeable. Define the sequence of functions

[TABLE]

By [37, Thm. 7.11], $\lim_{\delta\to c}\lim_{N\to\infty}h_{N}(\delta)=\lim_{N\to\infty}\lim_{\delta\to c}h_{N}(\delta)$ holds if $\{h_{N}\}$ is uniformly convergent on a set $E$ and $\lim_{\delta\to c}h_{N}(\delta)$ exists, given $c$ is a limit point of $E$.

For the case $z_{1}\neq z_{2}$ take $E=(0,|z_{1}-z_{2}|/2)$, and $c=0$. A straightforward integration yields

[TABLE]

It is evident that the limit at $\delta\to 0$ exists:

[TABLE]

To demonstrate the uniform convergence, we use Cauchy criterion

[TABLE]

Here, $M^{\prime}=\min(M,N)$ and $N^{\prime}=\max(M,N)$. This expression is independent of $\delta$ and can be made arbitrarily small by choosing sufficiently large $M^{\prime}$, and thus $\{h_{N}(\delta)\}$ converges uniformly. This, combined with existence of $\lim_{\delta\to 0}h_{N}(\delta)$ implies [37] that $\lim_{\delta\to 0}\lim_{N\to\infty}h_{N}(\delta)=\lim_{N\to\infty}\lim_{\delta\to 0}h_{N}(\delta)$ for $z_{1}\neq z_{2}$.

For the case $z_{1}=z_{2}$, we note that $B_{1}=B_{2}$ and take $E=(0,\infty)$. The sequence (88) after the direct integration is

[TABLE]

We note the existence of the limit $\lim_{\delta\to 0}h_{N}(\delta)=h_{N}(0)$ for any $N$, and use Cauchy criterion for showing the uniform convergence

[TABLE]

We note that this expression is independent of $\delta$ and can be made arbitrarily small by choosing appropriately large $M^{\prime}$. Thus the sequence is uniformly convergent, and the interchange of limits is valid for $z_{1}=z_{2}$ as well. This completes the proof.

This allows us to perform integration in the kernel expression (85) term by term, which yields

[TABLE]