Characterization of Power Absorption Response of Periodic 3D Structures to Partially Coherent Fields

# Characterization of Power Absorption Response of Periodic 3D Structures to Partially Coherent Fields

Denis Tihon Stafford Withington Cambridge University, Cavendish Laboratory, J.J. Thomson Avenue, CB3 OHE Cambridge, UK Christopher N. Thomas Cambridge University, Cavendish Laboratory, J.J. Thomson Avenue, CB3 OHE Cambridge, UK Christophe Craeye Université catholique de Louvain, ICTEAM Institute, Place du Levant 3, 1348 Louvain-la-Neuve, Belgium
###### Abstract

In many applications of absorbing structures it is important to understand their spatial response to incident fields, for example in thermal solar panels, bolometric imaging and controlling radiative heat transfer. In practice, the illuminating field often originates from thermal sources and is only spatially partially coherent when reaching the absorbing device. In this paper, we present a method to fully characterize the way a structure can absorb such partially coherent fields. The method is presented for any 3D material and accounts for the partial coherence and partial polarization of the incident light. This characterization can be achieved numerically using simulation results or experimentally using the Energy Absorption Interferometry (EAI) that has been described previously in the literature. The absorbing structure is characterized through a set of absorbing functions, onto which any partially coherent field can be projected. This set is compact for any structure of finite extent and the absorbing function discrete for periodic structures.

\journal

josaa \datesCompiled July 14, 2019 \ociscodes(000.4430) Numerical approximation and analysis; (010.5620) Radiative transfer; (030.4070) Modes; (050.1755) Computational electromagnetic methods; (300.1030) Absorption. \doihttp://dx.doi.org/10.1364/ao.XX.XXXXXX

## 1 Introduction

Electromagnetic metamaterials with a range of behaviours not found in nature can be realised from 3D periodic structures where some or all of the dimensions of the unit-cell are smaller than the wavelength of the exciting field. Metamaterials engineered in this manner are expected to provide a new means of subwavelength imaging [1, 2], improved broadband absorbers [3, 4] and to allow controlled radiative heat transfer between surfaces [5], among many other uses. In many of the envisaged applications, particularly those at higher frequencies (THz), the incident radiation is thermal in origin and so is only expected to be partially spatially coherent and in a partial state of polarization on reaching the structure. Examples include absorbers for thermal solar panels and bolometric detectors, as well as coatings for controlled radiative cooling. Being able to model and characterize the response of 3D structures illuminated by partially coherent fields is therefore crucial for optimizing the performance of the metamaterials used. However, it is an area that has received little attention compared with understanding the response to fully coherent fields and frequency-dependent behaviour.

In a previous paper [6], a way of characterizing the absorption of fields in any state of coherence by a structure has been developed. A formal link between theoretical results and experimental data that can be collected using Energy Absorption Interferometry (EAI) [7] has been derived. However, for the sake of clarity, only the 2D case and scalar wave function have been considered in [6], with an extension to structures periodic along one direction. In this paper, we present a formalism for power absorption by arbitrary 3D structures that can take into account the partial coherence and polarization of any stationary incident fields, with a special treatment for the case of doubly periodic structures. We prove that it is possible to deal with the partially coherent nature of the excitation, either by decomposing the incident fields into a sum of fully coherent fields [8, 9] or by computing the set of natural absorption modes of the structure, which is independent from the excitation.

The opportunities of EAI are very large. In fact, the method is able to recover the natural dynamical modes of any energy absorbing structure. In the case of microwave and optical photon-counting detectors for quantum communications it is essential to avoid, or at least terminate carefully, electromagnetic modes that can only couple noise and stray light into the detector. In the case of energy harvesting components, antenna arrays and absorbers, it is essential to maximize the number of degrees of freedom available for collecting power. The same considerations apply to near-field energy and information transfer between separated or overlapping volumes. In a recent paper [10], the method has been extended to recovering the collective excitations of quantum-mechanical systems.

This paper is organized as follows. In Section 2, we consider the problem of absorption of partially coherent fields by arbitrary structures and the treatment of the structure is decoupled from the treatment of the fields. In Section 3, the decomposition of the partially coherent fields into a sum of coherent modes is presented. In Section 4, the Energy Absorption Interferometry (EAI) technique is described from a theoretical point of view and, in Section 5, we derive a relation between the data collected and the characterization of the structure studied. In Section 6, the case of structures periodic in two directions is studied and a simplified formulation is proposed. In Section 7, we discuss the convergence of the infinite summations and integrations that are involved in the formulation. The computation of the absorption modes is presented and commented on in Section 8. A summary of the whole method is presented in Section 9 and last, in Section 10, numerical results are presented and commented upon.

## 2 Power absorbed by 3D structures

Let us assume a 3D structure with any distribution of real permittivity (), permeability () and resistivity (). The resistivity accounts for all the loss mechanisms, including dielectric losses. For clarity, magnetic losses are not considered hereafter. However, the extension of the proposed method to such losses is straightforward. We will assume that the incident fields are temporally stationary, as is this generally the case for fields of thermal origin. For a stationary field, the cross-correlation function vanishes over time [11], which is equivalent to saying that the different frequency components of the field are fully incoherent. As a result, the total power absorbed is simply the sum of the power absorbed at each of the different frequencies present, which we may consider separately. time dependence of the fields and currents will be implicitly assumed hereafter, with the angular frequency for which the structure is characterized.

Let denote the incident electric field over the reference plane located above the structure. The basis formed using the , and axes is orthonormal and right-handed. For later convenience, we decompose the incident field into vector components resulting from Transverse-Electric (TE) and Transverse-Magnetic (TM) field modes:

 ~Ei(kt)=(~Ei,TE(kt)~Ei,TM(kt)), (1)

with or the amplitude of the plane wave whose real transverse wave-vector is and whose electric or magnetic fields are transverse to the plane , respectively. The polarization vectors of these plane waves are

 ETE(kt) =1kt(−ky^x+kx^y) (2a) ETM(kt) =kzktk0(kx^x+ky^y+k2t/kz^z), (2b)

with the wave vector of the plane wave, is the magnitude of the transverse wave-vector, the free-space wave number and the propagation constant of the plane-wave along the direction. The plus sign has to be chosen for and the minus sign otherwise (i.e. . and are the free space permittivity and permeability respectively.

The transformation from the spatial to the spectral representations takes the form

 ~Ei,TE(kt) =∬∞−∞Ei(x,y)⋅ETE(kt)exp(jkt⋅r)dr (3a) ~Ei,TM(kt) =∬∞−∞Ei(x,y)⋅ETM(kt)exp(jkt⋅r)dr, (3b)

with the position where the incident fields are estimated and the complex conjugate of . The integration is carried out over the reference plane . The dot product used throughout this paper is defined by . In the rest of this paper, integration from to should be assumed unless other limits are stated. The reverse operation, which consists of computing the spatial incident fields distribution in the plane using their plane wave spectrum, is given by

 Ei(x,y)=14π2∬(~Ei,TE(kt)ETE(kt)+~Ei,TM(kt)ETM(kt))exp(−jkt⋅r)dkt. (4)

The incident fields are partially coherent, which means that the total fields within the structure correspond to stochastic functions [11]. The mean power density dissipated over time at point can be computed using

 P(r)=12ρ(r)⟨|K(r)|2⟩, (5)

where angle brackets denote an ensemble averaging, which in this case is equivalent to a time-averaging due to the stationarity of the stochastic fields. The vector corresponds to the current density in . The currents being a linear function of the incident fields, they can be decomposed according to the contribution of each incident plane wave:

 K(r)=14π2∬~K––––o(r|kt)⋅~Ei(kt)dkt, (6)

with a tensor describing the currents induced in by incident plane waves of transverse wave-vector and both possible polarizations. If one is interested in the total power absorbed by the structure, substituting (6) into (5) and integrating over the whole structure, we obtain:

 P=132π4∭Ωρ(r)⟨(∬~K––––o(r|kt)⋅~Ei(kt)dkt)∗⋅(∬~K––––o(r|k′t)⋅~Ei(k′t)dk′t)⟩dr, (7)

with the integration over the structure. Rearranging the integration and averaging orders, and defining a double-dot operation according to

 a∗⋅A––––⋅b=(a∗b):A––––, (8)

we obtain:

 P=116π4∬∬⟨~E∗i(kt)~Ei(k′t)⟩:Po––––––(kt,k′t)dktdk′t, (9)

with

 (10)

the cross-spectral power density function, where and are used to denote the transpose and Hermitian transpose of dyadic respectively. Equation (10) is the extension of (6) in [6]. The dyadic notation arises from the treatment of the polarization of the incident fields.

In Equation (9), two different contributions can be highlighted: , which only depends on the structure, and the ensemble average of the incident fields, which is independent of the structure. Both factors can be treated separately because the excitation is decoupled from the structure properties and geometry. It is worth noting that if we are only interested in the power absorbed in a given region of the structure (e.g. detectors), then the domain of integration of Equation (10) can be limited to that region.

## 3 Coherent modes representation of the fields

In this section, we will focus on the partially coherent incident fields and expand them into a set of coherent modes. We start from , the cross spectral density dyadic [8], as defined by

 W––––––(r,r′)=⟨E∗i(r)Ei(r′)⟩, (11)

which describes the state of spatial coherence of the incident fields.

It can be proven that the behaviour of the function is governed by the wave equations with respect to the coordinates and . Therefore, knowing the value of that function on a reference plane is sufficient to know the value of that function everywhere in space. However, propagating that function through an optical system or within the absorbing structure can be computationally intensive [9]. An elegant way of easing the problem is to decompose the stochastic fields into a sum of non-stochastic (i.e. fully deterministic) modes, whose relative amplitudes are completely uncorrelated [8, 9]. Each of these modes can then be treated using classical coherent-field solvers. The corresponding expansion of reads

 W––––––(r,r′)=∑pλpE∗p(r)Ep(r′), (12)

with

 Ei(r)=∑papEp(r) (13)
 ⟨a∗paq⟩=λpδpq, (14)

where is the Kronecker delta. The propagation of the correlation function then corresponds to the independent propagation of each mode. The modes and their coefficients can be found solving the eigenvalue equation [8]

 (15)

This concept can be transferred to the spectral domain noticing that the relation between the spatial and the spectral domains is linear. Applying a Fourier transform to (11) and (12) with respect to both and coordinates, using the orthogonality of the Fourier kernel over the infinite plane and using definition (3), it can be shown that

 ~W––––––(kt,k′t)=∑pλp~Wp––––––––(kt,k′t) (16)
 ~Wp––––––––(kt,k′t)=~E∗p(kt)~Ep(k′t) (17)

with

 ~W––––––(kt,k′t)= (18) (⟨~E∗i,TE(kt)~Ei,TE(k′t)⟩⟨~E∗i,TE(kt)~Ei,TM(k′t)⟩⟨~E∗i,TM(kt)~Ei,TE(k′t)⟩⟨~E∗i,TM(kt)~Ei,TM(k′t)⟩).

Note that the phase factor of the Fourier transform of with respect to in both (11) and (12) is instead of due to the complex conjugation. This result can be introduced in Equation (9) to yield:

 P=116π4∑pλp∬∬~Wp––––––––(kt,k′t):Po––––––(kt,k′t)dktdk′t. (19)

## 4 Energy Absorption Interferometry

Energy Absorption Interferometry (EAI) is an experimental technique used to characterize the way a structure can absorb electromagnetic fields [7]. The Structure Under Test (SUT) is illuminated with two sources whose positions and relative phases can be controlled. Recording the powers absorbed for different positions and phase difference between sources, it is possible to predict the power absorbed by the structure when illuminated by any incident fields.

EAI can be implemented experimentally in many different ways, and over a wide range of wavelengths [10], depending on the way the SUT is illuminated (e.g: near-field probes, far field radiators) and the way the power absorbed by the sample is recorded. The power absorption can be measured directly when the SUT corresponds to photovoltaics and electromagnetic detectors, while, for other types of structures, the absorption can be recorded through the local temperature rise using bolometric calorimeters. A first demonstration has been reported at submillimeter-wavelengths [12].

Consider two horizontal point sources corresponding to oscillating currents and such that

 J1(r) =p1δ(r1) (20a) J2(r) =p2δ(r2), (20b)

where is the Dirac delta function in 3D and where is the component of the current density in the direction . The positions of the two sources, and respectively, both lie in the plane . The relative phase between the sources can be controlled and their positions and orientations changed. The SUT is located in the half-space . This experimental arrangement is illustrated in Figure 1.

We can define the matrix that provides the fields generated on the plane by a point source located at :

 ~Ei(kt)≡S––––(kt|r)⋅p. (21)

As a translation of the source in the spatial domain induces a linear phase slope in the angular spectrum domain, (21) simplifies to

 ~Ei(kt)=exp(jkt⋅r)So––––––(kt)⋅p. (22)

is the plane-wave decomposition of the fields over produced by a source located at . The different dyadic components are given by

 So,11=ηkyk02kzktexp(−jkzh) (23)
 So,12=−ηkxk02kzktexp(−jkzh) (24)
 So,21=−ηkx2ktexp(−jkzh) (25)
 So,22=−ηky2ktexp(−jkzh), (26)

with the free space impedance. Substituting (22) into (6), the currents induced in the structure by a point source in become

 K(r)≡\leavevmode\nobreak K–––––p(r|r1)⋅p1=14π2∬exp(jkt⋅r1)~K0––––––(r|kt)⋅So––––––(kt)⋅p1\leavevmode\nobreak dkt. (27)

If the system is excited by two different point sources and located at positions and respectively, whose relative phase shift is , starting from Equation (5) and using (27), the total power absorbed is

 (28)

Distributing the product, separating the integrals and using the same compact notation as in Equation (9), it is possible to extend Equation (16) of [6] to 3D geometries and polarized excitations:

 P=12(α11+α22+2α12cos(ϕ−β12)), (29)

with

 αpq\leavevmode\nobreak ejβpq=∭Ω(ρ(r)K–––––p†(r|rp)⋅K–––––p(r|rq))dr:(p∗ppq) (30)

and . The correlation function defined in [6] can then be extended to 3D according to

 (31)

Using EAI experiments, the magnitude and phase of can be obtained varying the relative phase of the sources and measuring the oscillating pattern in the absorption.

## 5 Relation between EAI and spectral correlation

In this section, we will provide a link between the data that can be gathered using EAI, which corresponds to the spatial correlation function , and the cross-spectral power density function as defined in (10). Starting from the definition of the correlation function given in (31), we can use (27) to obtain

 C––––(r1,r2)=1(2π)4∬∬exp(−jkt⋅r1)So––––––†(kt)⋅∭Ωρ(r)\leavevmode\nobreak ~K––––o†(r|kt)⋅~K––––o(r|k′t)dr⋅So––––––(k′t)exp(jk′t⋅r2)dk′tdkt. (32)

We simplify this expression using (10) and the change of variables

 r+=r2+r12r−=r1−r2k+t=kt+k′t2k−t=kt−k′t, (33)

to yield

 C––––(r1,r2)=18π4∬∬So––––––†(kt)⋅Po––––––(kt,k′t)⋅So––––––(k′t)×exp(−jk−t⋅r+)exp(−jk+t⋅r−)dk+tdk−t. (34)

As in [6], we can notice that Equation (34) corresponds to a double inverse Fourier transform. Therefore, using the relation between a function and its Fourier transform, we finally obtain

 (35)

with

 C––––Δ(r+,r−)=C––––(r++r−2,r+−r−2). (36)

To summarize the results, if the spatial spectrum of the sources is known (through the matrix ), one has to apply a Fourier transform on the measurements to obtain the cross correlation function in the spectral domain.

## 6 Dissipation in periodic structures

The computation of the total power absorbed by a structure using formula (9) requires the evaluation of several integrals. The first one is a volume integral over the structure. The two others are 2D integrals over the plane wave spectrum of the incident light. In case the structure studied is periodic, these expressions can be simplified. Only the case of 2D periodicity of period and along the and directions will be considered hereafter. The extension to any other type of periodicity is straightforward.

Let us qualify as quasiperiodic an excitation which shares the same periodicity as the structure within a linear phase shift. Using the Floquet theorem, it can be shown that the currents induced on a periodic structure by a quasiperiodic excitation share the same quasiperiodicity as the excitation. This means that the product , which appears in Equation (10), is periodic with periods and in the and directions, respectively, and can be decomposed into a Fourier series

 (37)

with

 Qmn––––––––––(z|kt,k′t)=1ab∬a,b0,0ρ(r)\leavevmode\nobreak ~K––––o†(r|kt)⋅~K––––o(r|k′t)×exp(j(k′t−kt)⋅r)exp(j2kmnt⋅r)dy\leavevmode\nobreak dx, (38)
 (39)

These expressions can be used in Equation (19). Changing the order of summation and integration and noticing that , the mean power dissipated in each unit cell (19) can be expressed as

 (40)
 Hmn––––––––––(kt)=12×∭Ω0ρ(r)–~K––o†(r|kt−kmnt)⋅~K––––o(r|kt+kmnt)dr, (41)

with the intersection of the volume with the volume of one unit cell of the structure.

As mentioned in [6] for the 2D case, both and tend to vanish for large spatial frequencies or large indices. Therefore the summation and the spectral integration can be truncated. The upper values to consider depend on the distance that the fields have to travel before reaching the structure, and the depth of penetration of the fields inside the structure. Thus, they depend on the geometry of the structure studied.

takes discrete values, so is also discrete with respect to the difference of spatial frequencies. As underlined in Equation (35), a Fourier relation links to the data collected by the EAI, which is invariant w.r.t. a shift by an integer number of cells. Therefore, is periodic with respect to . Consequently, to fully characterize a periodic structure, one of the two sources only has to scan one unit-cell, while the other one still has to scan the whole structure, decreasing the total amount of data needed.

A final observation concerns the computation of the power absorbed by an infinite periodic structure excited by a non-periodic source. The non-periodic source can always be decomposed into a sum (or integral) of quasiperiodic sources [13], each quasiperiodic component being periodic within a different phase factor. Then, the total power absorbed by the structure corresponds to the sum of the powers absorbed by the structure for each quasiperiodic component of the source taken separately. Indeed, the interference between different quasiperiodic components of the source will vanish when summing the power dissipated by each unit cell of the infinite structure. However, if one is interested in the power dissipated within one unit cell, the local interference between the different quasiperiodic components of the source does not vanish anymore. Thus, equation (19) must be used instead of Equation (40). Adapting the phase of the correlation between the different quasiperiodic components of the source, it is possible to compute the power dissipated in any unit-cell (see (10)).

## 7 Evolution of Hmn with height

In most practical cases, the source of the incident fields (e.g: thermal source or optical system) is not touching the absorbing structure. The incident fields then have to propagate to reach the structure. Consider an incident field in . In , the fields will change according to the factor

 ~Ei(h′|kt)=~Ei(h|kt)exp(−jkz(h−h′)). (42)

Using Equation (9), we then obtain the dissipated power as

 P=132π4∬∬⟨~E∗i(kt)~Ei(k′t)⟩×exp(−j(kz+k′z)h):Po––––––(kt,k′t)dktdk′t, (43)

where

 k′z=±√k20−k′2t, (44)

the sign of being subject to the same convention as previously defined. We can then introduce the notation , which corresponds to the cross-correlation function for absorption if the fields are emitted in

 (45)

Similarly, for the periodic case, we can introduce

 (46)

with

 kmnz=±√k20−|kt+2kmnt|2. (47)

The sign of following the conventions already defined for in (2). The main consequence of Equations (45) and (46) is that the amplitude of the correlation function decreases exponentially with respect to the distance between the sources and the structure in the evanescent part of the spectrum. It means that, for a predefined accuracy, the farther the sources, the smaller the domain on which the and functions need to be actually tabulated.

Moreover, if we consider that , we have

 kz=−√k20−k2t≃−√−k2t=−jkt. (48)

Therefore, for high indices, we have

 ∣∣Hmn––––––––––(z|kt+kmnt)∣∣=∣∣exp(−j(kz+kmnz)z)Hmn––––––––––(kt+kmnt)∣∣≃∣∣exp(−jkzz)exp(−|2kmnt+kt|z)Hmn––––––––––(kt+kmnt)∣∣≤∣∣exp(−(2kmnt−2k0)z)Hmn––––––––––(kt+kmnt)∣∣ (49)

This means that once a certain distance between the sources and the structure is exceeded, the high order correlation functions then decrease exponentially with respect to . Consequently, the further the sources are from the absorbing device, the smaller the number of correlations functions required to characterize the absorber’s behaviour, and the smaller the domain over which these functions must be tabulated. Since the results from EAI and the correlation function are related by a Fourier transform, it also means that the further the sources are from the structure, the coarser the spatial sampling of the data can be.

## 8 Computation of the absorption modes

As suggested in [7] and implemented in [14] for a 2D geometry, this mathematical framework can be used to compute the modes through which a structure can absorb power. The total power dissipated by the structure corresponds to the sum of the powers dissipated by each of these modes separately.

Equation (9) can be restated within the theory of the linear operators of Hilbert spaces. We consider the functional space made of all the possible incident field distributions and define the scalar product between two fields distributions and as follows:

 [Ep|Eq]=∬E†p(kt)⋅Eq(kt)\leavevmode\nobreak dkt. (50)

An unconventional ‘bra-ket’ notation has been used to avoid confusion with , which corresponds to the ensemble average of the fields. For any physical fields distribution, this scalar product is sesquilinear, positive-definite and exhibits a conjugate symmetry (i.e. ). In the framework of this Hilbert space, we can define the operator such that

 |^P~E]=116π4∬Po––––––(kt,k′t)⋅~E(k′t)\leavevmode\nobreak dk′t. (51)

 (52)

Using the definition of in (10), it can be shown that the operator is Hermitian, so that its eigenvalues are real and the corresponding eigenvectors are orthogonal. From a physical point of view, one of the main consequences is that the power dissipated by any incident field is

 P=∑αΛα∣∣[Eα|~Ei]∣∣2, (53)

with the magnitude of . Therefore, the power is always real-valued, as required by physical law. Moreover, from a mathematical point of view, it means that a simple orthogonal projection of the incident fields on the different modes is required to obtain the total power absorbed. Therefore, once all the modes with a non-negligible eigenvalue are known, the exact projection of the incident fields into these modes does not require any knowledge about the other modes. This is convenient in the sense that only the incident fields are required, and not the total fields (incident and reflected). Last, for finite structures, it can be shown that the operator is Hilbert-Schmidt and therefore compact. It means that the number of modes whose contribution to absorption is not negligible is finite. Moreover, in the case of structures that are small with respect to the wavelength of operation, this number is generally small.

From a numerical point of view, the operator has to be discretized to compute these modes through an eigenvector decomposition. It can also be proven that, in the case of a periodic structure, each absorption mode is only composed of Floquet modes characterized by a common phase shift between consecutive unit cells. Therefore, the absorption modes of a periodic structure can be computed by simulating just one unit cell with a periodic solver and the appropriate periodic boundary conditions. This results in a naturally discretized operator, such that the numerical solution of the problem is straightforward.

As mentioned previously, only the modes whose contribution to absorption is non-negligible need to be kept, resulting in a significant reduction in the amount of data that must be stored. The data required to fully characterize the absorption by non-periodic structures reduces from a 4D function to a compact set of 2D functions corresponding to the absorption modes for and . For periodic structures, a discrete set of discrete functions has to be stored for each possible phase shift between consecutive unit cells. That data can be reorganized into a discrete set of 2D functions, where each function is defined over the first Brillouin zone. Generally, both sets for periodic and non-periodic structures are small if the structure or its unit cell are not large with respect to the wavelength of the incident fields.

## 9 Absorption by a structure

Now that all necessary operators are available, we can address the problem of power absorption by a periodic structure illuminated by a spatially partially coherent field. However, we will first briefly summarize the results obtained so far.

In Section 3, we considered the treatment of spatially partially coherent fields by decomposition into coherent modes. The knowledge of the response of the structure to coherent fields is then sufficient to calculate the absorbed power. In Section 2, we provided a alternative scheme for characterizing power absorption by any structure, based on the computation of the 4D cross-spectral power density function . Then, in Section 6, we showed that, provided that the structure is periodic, a 2D compact set of 2D functions is sufficient to fully characterize the way in which the periodic structure can absorb fields.

Once this information has been retrieved using simulations or experiments, the collected data needs to be stored. The data can be reduced in size by the use of the so-called absorption modes (cf. Section 8).

Expanding Equation (53), the total absorbed power for any coherent excitation becomes

 P=∑αΛα∬E†α(kt)⋅~Ei(kt)dkt∬~E†i(k′t)⋅Eα(k′t)dk′t. (54)

Rearranging this expression and noticing that a partially coherent excitation can be rephrased as a sum of coherent excitations that are orthogonal to each other, it can be shown that the total power absorbed by a structure when it is illuminated by partially coherent fields is

 P=∑αΛα∬∬E†α(kt)⋅~W––––––(kt,k′t)⋅Eα(k′t)dktdk′t. (55)

In compact notation, this becomes

 P=∑αΛα[Eα|^WEα]. (56)

Notice the resemblance between (56) and (52). Modes and are obtained using an eigenvalue decomposition of the operators and , respectively. One can see that, provided that this eigenvalue decomposition has been performed for at least one of the two operators, the partially coherent nature of the fields can be handled using either (19) or (55). If one wants to study the way in which given partially coherent fields can be absorbed by any structure, the eigenvalue decomposition of the incident fields, which is independent from the structure, can be performed once and the power absorbed by each structure computed using (19), while if one is interested in the way a given structure can absorb any partially coherent incident fields, it is more interesting to compute the absorption modes of the structure, which are independent from the incident fields, and apply (55) to compute the power absorbed for each excitation.

## 10 Numerical results

To illustrate the technique, we computed the functions for a periodic structure mimicking the typical design of room-temperature bolometers [16]. The structure consists of a squared lattice of gold nanospheres used as an absorbing layer deposited on a suspended membrane. The membrane is 200 nm thick and the radius of the spheres is 100 nm. The spacing between two consecutive spheres is 200 nm. We positioned the spheres 5 nm above the surface of the layer to avoid some numerical issues. However, this distance being smaller that of the wavelength, it has a negligible impact on the final results. The free-space wavelength of the incident fields is 750 nm, which corresponds to a maximum of absorption for a 100 nm free-standing single golden sphere [15]. Plasmonic effects of gold at that frequency are taken into account by using a negative value for the real part of the permittivity, while losses mechanisms are taken into account adding a negative imaginary part to the permittivity. At this wavelength, the complex relative permittivity of gold is then [17] and the permittivity of the suspended membrane is . The structure can be seen on Figure 2. The plane coincides with the top of the spheres.

The coherent solver used to perform the full-wave simulation of the structure is based on the Poggio-Miller-Chan-Harrington-Wu-Tsay (PMCHWT) formulation applied to the Method of Moments [19, 20, 21]. The periodicity of the structure is handled through the use of the periodic Green’s function [22]. The implementation of the code is described in [23]. The sphere and the interfaces between the layer and the air are discretized using 690 Rao-Wilton-Glisson (RWG) [24] and 800 rooftop basis functions, respectively.

First, we computed using (41). The results can be seen in Figure 3 for different and different entries of the matrix. The first two graphs having purely real values, their phase is not displayed. The log-magnitude of the field is shown for better readability.

First we consider Fig. 3(a). Two circles centered on are visible. The inner one corresponds to the limit of the visibility angle. Inside that circle, modes are propagative in free space while outside modes are evanescent. The function has to vanish when approaching this circle from the inside, as the power carried by the incident plane wave is vanishing for grazing incidence. The second circle corresponds to a leaky surface mode of the structure (see for example [18] and references therein). Most of the other noticeable inhomogeneities correspond to the replicas of these circles into other Brillouin zones [25].

Looking at Fig. 3(b), we can observe similar circles. Interestingly, the wave vector of the surface leaky wave is slightly different for TM or TE excitations ( vs. , respectively), corresponding to the fact that there exist two different surface leaky waves. The central square is due to the excitation of this TM leaky wave through its harmonics, while the others inhomogeneities in the visible region are due to the excitation of the TE leaky wave.

The graph (c) in Fig. 3 clearly exhibits the features associated with the visibility limit and the two leaky surface waves. In addition, some straight lines for which the coupling between TE and TM modes vanishes appear due to the symmetries of the structure. Finally, in graph (e), we can see that the circles are centered on . The reciprocal lattice vectors for a periodicity of nm corresponding to , it is coherent with the expected results.

We also computed the natural absorption modes of the structure and their corresponding eigenvalues for periodic excitations with no phase shift between consecutive unit cells. The eigenvalues were computed for different distances between the reference plane and the absorbing structure. The results can be seen in Figure 4. We observe that the closer the sources are to the structure, the higher is the number of different modes through which the structure can absorb energy. It is consistent if we consider that evanescent modes can transport energy if there is some reflection on the surface of the structure, but that the amplitude of these modes will decrease over the distance separating the reference plane and the structure. Further, we see that, as the distance of the sources from the observer increases, the far-field solution is eventually recovered.

As a final illustration of the technique, we consider the following example. The structure simulated here above is illuminated by an incident plane wave which is partially polarized into TE mode. This plane wave is split into two parts of equal amplitude and one of the two beams goes through a TE polarizer. The incident fields are not perfectly single-frequency and are therefore characterized by a finite coherence length, which is supposed to be much larger than the typical dimensions of one unit-cell of the structure under study. The two beams follow paths with different lengths before reaching the structure, so that the amplitude of the two waves is not perfectly coherent. The non-polarized wave reaching the structure is characterized by a transverse wave vector and the polarized beam is characterized by the transverse wave vector . These two wave vectors correspond to an identical phase shift between consecutive unit cells, so that the power absorbed due to one mode depends on the amplitude of the other mode. We want to compute the way in which this excitation can be absorbed by the structure.

First, the fields exciting the structure can be characterized using the spectral cross-correlation function . The excitation consisting of only two plane waves, we have

 ~W––––––(k(1)t,k(1)t) =(0.5+A000.5−A) (57a) ~W––––––(k(1)t,k(2)t) =(B(0.5+A)000) (57b) ~W––––––(k(2)t,k(2)t) =(0.5+A000), (57c)

with a real factor that describes the degree of polarization of the incident wave and a complex factor that depends on the coherence length of the excitation and the path difference of the two beams. corresponds to a perfectly polarized TE wave while corresponds to a perfectly unpolarized wave. If the two beams travel over the same distance, they arrive in phase and . But if the lengths the two beams travel are different, becomes complex and its amplitude tends to decrease, until the difference between paths is much bigger than the coherence length, in which case . Using (18), it can be shown that .

We can decompose this correlation matrix into a sum of properly weighted perfectly coherent modes (cf. Equation (12)), which provides

 ~E1(k(1)t)=(01)\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ~E1 (k(2)t)=(00) λ1=0.5−A (58a) ~E2(k(1)t)=1√2(10)\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ~E2 (k(2)t)=1√2(B/|B|0) λ2=(1+|B|)(0.5+A) (58b) ~E3(k(1)t)=1√2(10)\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ~E3 (k(2)t)=1√2(−B/|B|0) λ3=(1−|B|)(0.5+A) (58c)

From a physical point of view, we can see that the partially coherent excitation can be decomposed into

1. a TM plane wave of transverse wave vector ,

2. two TE plane waves of equal amplitude and transverse vectors and ,

3. two TE plane waves of opposite amplitude and transverse vectors and .

The amplitude of each of these modes depends on the degree of coherence of the incident fields. The total power dissipated inside the structure will correspond to the sum of the powers absorbed by each of these modes separately.

Alternatively, we can compute the natural modes of absorption of the structure from the data computed previously (see Fig. 3). Since we are dealing with a far-field excitation, only four terms have to be taken into account in Equation (53), whose associated fields distribution are a linear combination of the four propagating Floquet modes associated to the phase shift between consecutive unit cells.

In matrix form, the eigenmodes and the associated eigenvalues are

 (k(2)t)=1√2(0−1) Λ1=5.75 (59a) (k(2)