# Contact-less Material Probing with Distributed Sensors: Joint Sensing and
Communication Optimization

Ali Kariminezhad, Member, IEEE, Soheil Gherekhloo, Member, IEEE,
and Aydin Sezgin, Senior Member, IEEE
A. Kariminezhad is with the Elektronische Fahrwerksysteme GmbH, Germany,
(Email: ali.kariminezhad@efs-auto.de), S. Gherekhloo is with the Chassis
Systems Control, Robert Bosch GmbH, Germany, (Email:
<EMAIL_ADDRESS>and A. Sezgin is with the Institute of Digital
Communication Systems, Ruhr-Universität Bochum (RUB), Germany, (Email:
aydin.sezgin@rub.de).

###### Abstract

The utilization of RF signals to probe material properties of objects is of
huge interest both in academia as well as industry. To this end, a setup is
investigated, in which a transmitter equipped with a two-dimensional multi-
antenna array dispatches a signal, which hits objects in the environment and
the reflections from the objects are captured by distributed sensors. The
received signal at those sensors are then amplified and forwarded to a
multiple antenna fusion center, which performs space-time post-processing in
order to optimize the information extraction. In this process, optimal design
of power allocation per object alongside sensors amplifications is of crucial
importance. Here, the power allocation and sensors amplifications is jointly
optimized, given maximum-ratio combining (MRC) at the fusion center. We
formulate this challenge as a sum-power minimization under per-object SINR
constraints, a sum-power constraint at the transmitter and individual power
constraints at the sensors. Moreover, the advantage of deploying zero-forcing
(ZF) and minimum mean-squared error (MMSE) at the fusion center is discussed.
Asymptotic analysis is also provided for the case that large number of sensors
are deployed in the sensing environment.

##  I Introduction

The material properties of objects can be obtained by capturing the response
of the object to an stimulating impulse. This process can be fulfilled by
classical contact probing which demands manual and exhaustive workload.
However, contact-less radio-frequency (RF) sensing is an alternative solution
for an autonomous sensing process. This advantage of wireless sensing does not
come for free and it requires sophisticated signal processing tasks both at
the transmission and reception sides. These processing tasks include
localization, channel estimation, synchronization and optimal resource
allocation.
The objects in the sensing environment are categorized to be either passive or
active signal sources. Identifying active sources requires passive sensing
systems which measure the intended quantity. For robustness of such systems,
multiple distributed sensors can be exploited to make the observation on the
sensing channel which is then forwarded to the fusion center for joint
processing. In a passive sensing system, the authors in [1] provide an
analytical solution for power allocation at the sensors under sensor sum-power
constraint. They study a single source scenario, where the mean squared of the
estimation error is minimized for unbiased estimators. This solution is not
tractable practically due to unrealizable power trading assumption among
distributed sensors. Consequently, the authors in [2] investigate a similar
problem albeit under individual power constraints. Passive sensing systems are
investigated in [3, 4, 5, 6] from different perspectives including
localization, scheduling and energy harvesting. Contrary to active sources,
passive sources require active sensing systems. These systems utilize
stimulation- and/or reflection-based approaches, where the information is
extracted from the backscattered channel. The authors in [7] study the optimal
waveform design in a collocated multi-input multi-output (MIMO) system with
single extended object. In that work, mutual information between the transmit
and received signals is maximized. Moreover, they study the waveform design
for minimizing the mean-squared of the channel estimation error. Moreover, the
authors in [8] study similar problem assuming multiple extended objects. For
non-collocated transceivers, the authors in [9] study the optimal waveform
design that maximizes the so-called Bhattacharyya distance. That work mainly
focuses on single-object detection in a single-clutter environment.

...123...Fusion Center Figure 1: RF-based stimulation of the objects in the
sensing environment with amplify-and-forward sensors. The two top objects are
the objects and the bottom one is the clutter.

In this paper, an active sensing system is exploited for material
characterization purposes. This characterization can be fulfilled by
estimating the second-order moment of the materials, assuming that their
impulse response is random and follows a Gaussian distribution. We address
this active system by utilizing a two-dimensional multi-antenna transmitter
which pre-processes the signal given the position of the objects in the
sensing environment. This design allows three-dimensional beamforming which
has been shown to be beneficial [10, 11, 12, 13, 14]. Here, we assume that the
position of the objects are known, however [15] and the references therein
investigate the localization problem thoroughly.

Having this information, the dispatched signal power at the objects surface is
maximized by transmitting in the direction of their steering vectors. Then,
the objects act differently to this incident signal by having various
emissions in a particular spectrum. Sensing the response of the objects in a
particular spectrum can help classifying the materials. For instance, in
photo-acoustic imaging, the object is stimulated at higher frequency spectrum,
however the response is captured at ultrasonic frequency range [16, 17, 18].
However, the reflection by the objects over the same transmit signal spectrum
can also be helpful for identification purposes. Therefore, at the same
frequency spectrum, the reflected signal from the objects are sensed,
amplified at multiple sensors, and then forwarded, i.e., amplify and forward
(AF), to the fusion center. The fusion center is equipped with multiple
single-antenna baseband units with high capacity links. The receive antennas
at the fusion center observe the noisy version of the forwarded signals from
the sensors, Fig. 1. Then, the fusion center performs post-processing for
object response detection. This detection is aimed to guarantee a certain
quality of service (QoS), which is quantized by signal-to-interference-plus-
noise power ratio (SINR) in this paper. To this end, available resources need
to be allocated optimally at the transmitter, sensors and fusion center. The
signals at the sensors are amplified and transmitted to the fusion center in
distinct time instants. Hence, the sensors operate in time-division multiple-
access (TDMA) perfectly (without any collision between their observations at
the fusion center). The observed signal at the multi-antenna fusion center
buffers the received signals and performs space-time processing.

We exploit maximum-ratio transmission (MRT) at the two-dimensional multi-
antenna transmitter. Moreover, various reception schemes are investigated for
performance comparison. These schemes are zero-forcing (ZF), maximum-ratio
combining (MRC) and minimum mean-squared error (MMSE). For MRC, the transmit
power allocation, and sensor amplification coefficients are optimized jointly.
Given per-object SINR constraints, sum-power at the multi-antenna transmitter
and sensors is minimized. This problem turns out to be a signomial problem
(SP), which is solved in this paper iteratively. As a benchmark, we compare
the performance of MRC and other receiver types with maximum amplification at
the sensors. Moreover, we propose a separate optimization algorithm for sum-
power minimization problem exploiting a MSE receiver. Asymptotic analysis is
also considered for the case that the fusion center is equipped with a massive
antenna array.

###  I-A Contribution

Conventionally, sensors are exploited for sensing quantities of interest,
e.g., temperature, which demand active sources. In this paper, for the first
time, we introduce an active sensing system for material characterization
purposes. In this system the objects are passive elements, hence a transmitter
is required for triggering purposes. In this system, the optimal power
allocation at the transmitter and optimal amplification at the sensors are
studied. The optimization problems, turn out to be non-convex problems. We
propose efficient algorithms to obtain good sub-optimal solutions in
polynomial time.

###  I-B Organization

The system model and the related assumptions are presented in section II. In
section III and section IV we discuss the applicable pre- and post-processing
schemes exploited in the problem, respectively. The optimization problems are
elaborated in section V. We provide asymptotic analysis in section VI. The
procedure of obtaining reliable second-order moment for material
characterization purposes in discussed in section VII. We provide the
numerical results in section VIII. Finally, we conclude the paper in section
IX

###  I-C Notation

Throughout the paper, we represent vectors in boldface lower-case letters
while the matrices are expressed in boldface upper-case. ${\bf{a}}^{H}$,
${\bf{a}}^{T}$, $\|\mathbf{a}\|_{1}$ and $\|\mathbf{a}\|_{2}$ represent
hermitian, transpose, $l_{1}$ and $l_{2}$ norms of vector $\mathbf{a}$,
respectively. The inner product of vectors $\mathbf{a}$ and $\mathbf{b}$ is
represented by $\langle\mathbf{a},\mathbf{b}\rangle$. The mutual information
between two random variables $x$ and $y$ is represented by $I(x;y)$. Moreover
the identity matrix is depicted by $\mathbf{I}$. The argument that
maximizes/minimizes a function $f(\mathbf{x})$ is represented by
$\mathbf{x}^{\star}$. The real component of a complex-valued variable is
denoted by $\Re\\{.\\}$.

Time slot 1Time slot 2Time slot
K$\mathbf{v}_{1}$$\mathbf{v}_{2}$$\mathbf{v}_{N}$$z_{1}$$z_{2}$$z_{N}$ Figure
2: Receiver structure at the fusion center. The received signal at time slot
$k$, consists of the transmit signal from sensor $k$. The received signal over
$K$ time slots are processed simultaneously for the $j$th object information
extraction, i.e., $l_{j},\ \forall j\in\mathcal{N}^{t}$. The communication
channel is assumed to be constant over time.

##  II System Model

We consider a sensing environment with $N^{t}$ objects of interest and $N^{c}$
clutters, i.e., $N=N^{t}+N^{c}$ objects in total. The RF signal from a single
multi-antenna transmitter with $M\times M^{{}^{\prime}}$ planner antenna array
stimulates the objects in the sensing environment, where there is a line-of-
sight (LoS) between the transmitter and objects. We assume that, the transmit
antennas are equi-distantly positioned in two-dimensional Cartesian basis
dimensions (uniform antenna array). Here, we consider the antenna at the
center of coordinates as the reference antenna. By assuming half-wavelength
distance between horizontal and vertical antennas, we obtain the steering
vector (including path loss) corresponding to object $i$ as

$\displaystyle\mathbf{a}_{i}=\sqrt{\nu_{i}}\big{[}1,$
$\displaystyle\cdots,e^{j\pi(m\sin\theta\sin\phi+m^{{}^{\prime}}\cos\phi)},\cdots$
$\displaystyle,e^{j\pi((M-1)\sin\theta\sin\phi+(M^{{}^{\prime}}-1)\cos\phi)}\big{]}^{T}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{C}}^{MM^{{}^{\prime}}},$
(1)

where $m\in\\{0,\cdots,M-1\\}$ and
$m^{{}^{\prime}}\in\\{0,\cdots,M^{{}^{\prime}}-1\\}$. Notice that the azimuth
and elevation angles of object $i$ are represented by $\theta_{i}$ and
$\phi_{i}$, respectively. Moreover, $\nu_{i}\leq 1,\ \forall i\in\mathcal{N}$
accounts for the path loss towards the $j$th object which is a function of
application environment, distance of the objects from the transmitter and
transmit signal frequency spectrum [19]. In this paper, we consider the effect
of distance on the path loss only. Therefore, $\nu_{j}=r_{j}^{-\gamma}$, where
$r_{j}$ is the distance of the $j$th object from the transmitter and $\gamma$
is the path loss coefficient. Hence, the objects located at the radius of 1
meter (unit distance) are considered as the reference objects (do not suffer
from path loss). This coefficient is experimentally observed to be in the
range of $2-5$. Here, we adopt $\gamma=2$, which is the path loss coefficient
in free space with line of sight (LoS) path.

...Tx$\mathbf{s}$$l_{1}$...$l_{N^{t}}$$l_{N^{t}+1}$...$l_{N^{t}+N^{c}}$$\mathbf{a}_{1}^{H}$$\mathbf{a}_{N^{t}}^{H}$$\mathbf{a}_{N^{t}+1}^{H}$$\mathbf{a}_{N^{t}+N^{c}}^{H}$...$y_{K}$$y_{1}$$g_{11}$$g_{1K}$$g_{N^{t}1}$$g_{N^{t}K}$$g_{(N^{t}+1)1}$$g_{(N^{t}+1)K}$$g_{(N^{t}+N^{c})1}$$g_{(N^{t}+N^{c})K}$$\mathbf{f}_{1}$$\mathbf{f}_{K}$...RX
Figure 3: Detailed illustration of the system.

Now, having the steering vectors of the objects, the transmit signal is formed
as

$\displaystyle\mathbf{s}=\sum_{j=1}^{N_{t}}\mathbf{u}_{j}\sqrt{p_{j}},$ (2)

where
$\mathbf{u}_{j}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{C}}^{MM^{{}^{\prime}}}$
specifies the beam direction towards the $j$th object and $p_{j}$ is the
transmit power allocated to the $j$th object. Notice that, the transmit
symbols are assumed to be known at the fusion center. Now, the reflected (or
radiated) signal from the objects (not necessarily in the same frequency
range) are given by

$\displaystyle x_{i}=\mathbf{a}^{H}_{i}\mathbf{s}l_{i},\ \forall
i\in\mathcal{N},$ (3)

where $\mathcal{N}=\\{\mathcal{N}^{t}\cup\mathcal{N}^{c}\\}$. Moreover,
$\mathcal{N}^{t}$ and its complement $\mathcal{N}^{c}$ are the sets of objects
and clutters, respectively, i.e., $\mathcal{N}^{t}=\\{1,...,N^{t}\\}$ and
$\mathcal{N}^{c}=\\{N^{t}+1,...,N^{t}+N^{c}\\}$. Moreover, the response of the
object $i$ to the incident RF signal is represented by $l_{i}$. Here, we
assume that
$l_{i}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{C}},\
\forall i\in\mathcal{N},$ have independent and identical zero-mean Gaussian
distribution. Therefore, estimating the second-order moment
$\mathbb{E}\\{|l_{i}|^{2}\\},\ \forall i$, helps classifying the objects. Now
the signal $x_{i},\ \forall i$, is sensed by $K$ distributed sensors. As it is
shown is Fig. 3, the received signal at the $k$th sensor is given by

$\displaystyle y_{k}=\sum_{i=1}^{N}g_{ik}x_{i}+n_{k},\ \forall
k\in\mathcal{K}=\\{1,...,K\\},$ (4)

where $g_{ik}\in\mathbb{C}$ is the channel from the object $i$ towards the
sensor $k$ and $n_{k}\in\mathbb{C}$ is the additive receiver noise at sensor
$k$. Here, we assume that the noise at all sensors $n_{k},\ \forall
k\in\mathcal{K},$ follow identical and independent zero-mean Gaussian
distribution. Moreover, the noise variance at sensor $k$ is given by
$\sigma^{2}_{k}$. The signal is amplified at the sensors and forwarded to the
fusion center in different time slots. One can easily write the signal
received by the fusion center with $R$ antennas at time slot $k\in\mathcal{K}$
as follow

$\displaystyle\mathbf{y}_{\text{fc},k}=\sqrt{\alpha}_{k}\mathbf{f}_{k}\left(\sum_{i=1}^{N}g_{ik}x_{i}+n_{k}\right)+\mathbf{n}_{\text{fc},k},$
(5)

where $\alpha_{k}$, $\mathbf{f}_{k}\in\mathbb{C}^{R}$, and
$\mathbf{n}_{\text{fc},k}\in\mathbb{C}\mathcal{N}(\mathbf{0}_{R},\sigma_{\text{f}c}^{2}\mathbf{I}_{R})$
denote the amplification factor of the $k$th sensor, the communication channel
from the $k$th sensor to the fusion center, and the received noise at the
fusion center in time slot $k$, respectively. At the fusion center, all
received signals over the $K$ time slots are stacked in the vector
$[\mathbf{y}_{\text{fc},1}^{T},\ \cdots,\mathbf{y}_{\text{fc},K}^{T}]^{T}$ and
post-processed using $\mathbf{V}=[\mathbf{v}_{1},\
...,\mathbf{v}_{N^{t}}]\in\mathbb{C}^{KR\times N^{t}}$ (see Fig 2). Doing
this, we obtain

$\displaystyle\mathbf{z}=\mathbf{V}^{H}\left(\sum_{i=1}^{N}\sqrt{\delta_{i}}l_{i}\mathbf{w}_{i}+\mathbf{n}^{{}^{\prime}}\right),$
(6)

where $\delta_{i}=|\mathbf{a}^{H}_{i}\mathbf{s}|^{2}$ represents the power of
the signal received at object $i$. Moreover, the vectors $\mathbf{w}_{i}$ and
$\mathbf{n}^{{}^{\prime}}$ are the equivalent channel and the noise vector,
respectively. They are defined as

$\displaystyle\mathbf{w}_{i}$
$\displaystyle=\begin{bmatrix}\sqrt{\alpha_{1}}g_{i1}\mathbf{f}^{T}_{1}&...&\sqrt{\alpha_{K}}g_{iK}\mathbf{f}^{T}_{K}\end{bmatrix}^{T},$
(7) $\displaystyle\mathbf{n}^{{}^{\prime}}$
$\displaystyle=\underbrace{\begin{bmatrix}\sqrt{\alpha_{1}}n_{1}\mathbf{f}^{T}_{1}&...&\sqrt{\alpha_{K}}n_{K}\mathbf{f}^{T}_{K}\end{bmatrix}^{T}}_{\mathbf{n}_{\text{s}}}+\mathbf{n}_{\text{fc}},$
(8)

where,
$\mathbf{n}_{\text{fc}}=[\mathbf{n}_{\text{fc},1}^{T},\cdots,\mathbf{n}_{\text{fc},K}^{T}]^{T}\sim\mathbb{C}\mathcal{N}(\mathbf{0}_{KR},\sigma_{\text{fc}}^{2}\mathbf{I}_{KR})$.
By introducing the vectors as in (7) and (8), the system under investigation
is simplified to a special form of a multiple access channel (MAC) as depicted
in Fig. 4. Here, the observations in time can be represented as the
observations in space by considering virtual antenna array. Hence, a virtual
array of $KR$ antennas can be assumed over a single observation instant.
Consequently, this channel is referred as the virtual multiple access channel
with interference (VMACI) throughout this paper. This channel consists of a
multiple access channel interfered by multiple clutters whose transmit powers
are functions of the allocated power to the MAC transmitters. However,
compared to the MAC, the VMACI has an extra degrees of freedom due to the
included amplification factors addressed by the parameter $\alpha_{i}$ in the
channel vectors. Here, we observe the dilemma, on one side, the sensor
amplifications improve the channel condition for objects
$\mathbf{w}_{i}(\boldsymbol{\alpha}),\ \forall i\in\mathcal{N}^{t}$ and
deteriorate the channel condition for clutters, i.e.,
$\mathbf{w}_{j}(\boldsymbol{\alpha}),\ \forall j\in\mathcal{N}^{c}$. On the
other side, the transmit power needs to be allocated to the objects, so that
less amount of power is received at the clutters surface. This is due to the
dependency of their reflected power on the reflected power from the objects.
Notice that the information symbol from the $i$th virtual user in VMACI, i.e.,
$m_{i}$, has the same differential entropy as the reflection coefficient from
the $i$th object, i.e., $l_{i}$, (i.e., $h(m_{i})=h(l_{i}),\ \forall i$). We
define the power of the virtual symbols $m_{i}$ by

$\displaystyle e_{i}=\mathbb{E}\\{|m_{i}|^{2}\\}=\delta_{i}Q_{i},\ \forall
i\in\\{\mathcal{N}\\},$ (9)

where the second-order moment of the reflection coefficient from the $i$th
object is represented by $Q_{i}=\mathbb{E}\\{|l_{i}|^{2}\\},\ \forall i$. The
estimation quality of $Q_{i}$ is a function of the beamforming vectors at the
transmitter, i.e., $\mathbf{u}_{j},\ \forall j\in\mathcal{N}^{t}$, which
controls the amount of received power at the objects surface. Moreover, sensor
amplification coefficients and post-processing matrix at the fusion center are
aimed to be optimized for performance improvement. In what follows, we assume
that only the steering vectors correspond with the objects are known at the
transmitter, i.e., $\mathbf{a}_{i},\ \forall i\in\mathcal{N}^{t}$. Whereas,
the following knowledge is given at the fusion center.

  1. 1.

$\mathbf{a}_{i}\ \forall i\in\mathcal{N}$,

  2. 2.

$g_{ik},\ \forall i\in\mathcal{N},k\in\mathcal{K}$,

  3. 3.

$\mathbf{f}_{k},\ \forall k\in\mathcal{K}$.

The channel between the sensors and the fusion center, i.e., $\mathbf{f}_{k},\
\forall k\in\mathcal{K}$, can be obtained priori in a channel training phase.
Moreover, we assume that the line-of-sight (LoS) is the dominant path between
the objects (objects and clutters) and the sensors. Hence, given the position
of the objects and the sensors, the LoS path of the channels $g_{ik},\ \forall
i\in\mathcal{N},k\in\mathcal{K}$ can be estimated and the non-LoS is ignored.
In the next section we discuss pre processing methods for three-dimensional
beamforming at the multi-antenna transmitter.

##  III Pre Processing

Given the steering vectors of the objects at the transmitter, maximum ratio
transmission (MRT) is the optimal scheme. In MRT, the transmit directions
toward the $j$th object ($j\in\mathcal{N}^{t}$) is adjusted to the
corresponding steering vectors, i.e., $\mathbf{a}_{j}$. Hence,

$\displaystyle\mathbf{u}_{j}=\frac{\mathbf{a}_{j}}{\|\mathbf{a}_{j}\|_{2}}=\frac{1}{\sqrt{MM^{{}^{\prime}}\nu_{j}}}\mathbf{a}_{j}.$
(10)

Utilizing this filter at the transmitter the received signal power at the
$j$th object (either a object or a clutter) is written as

$\displaystyle\delta_{j}=$
$\displaystyle\nu_{j}\left(MM^{{}^{\prime}}p_{j}+\sum_{\begin{subarray}{c}i=1\\\
i\neq
j\end{subarray}}^{N^{t}}p_{i}\mathbf{a}^{H}_{j}\mathbf{u}_{i}\mathbf{u}^{H}_{i}\mathbf{a}_{j}\right)\quad
j\in\mathcal{N}^{t},$ (11) $\displaystyle\delta_{j}=$
$\displaystyle{\nu_{j}\left(\sum_{i=1}^{N^{t}}p_{i}\mathbf{a}^{H}_{j}\mathbf{u}_{i}\mathbf{u}^{H}_{i}\mathbf{a}_{j}\right)}\quad\forall
j\in\mathcal{N}^{c}$ (12)

where $MM^{{}^{\prime}}$ is the antenna gain for the $j$th object. In the next
section, we elaborate the post processing schemes that are exploited at the
fusion center.

##  IV Post Processing

As can be noticed from (6), the post-processed signal includes both desired
and interference components. This can be seen by

$\displaystyle
z_{j}=\mathbf{v}_{j}^{H}\left(\underbrace{\sqrt{\delta_{j}}l_{j}\mathbf{w}_{j}}_{\text{desired}}+\underbrace{\sum_{\begin{subarray}{c}i=1\\\
i\neq
j\end{subarray}}^{N}\sqrt{\delta_{i}}l_{i}\mathbf{w}_{i}}_{\text{interference}}+\mathbf{n}^{{}^{\prime}}\right),\
\forall j\in\mathcal{N}^{t},$ (13)

where the $j$th column of the post-processing matrix $\mathbf{V}$ is denoted
by $\mathbf{v}_{j}\in\mathbb{C}^{KR}$. Now, the post-processing filters, power
allocation per object and signal amplification at the sensors need to be
designed to guarantee a certain threshold is differentiating the objects
$l_{j},\forall j\in\mathcal{N}^{t}$. The better this differentiating quality,
the more robust classification for material characterization purposes can be.
Intuitively, increasing the amount of information about the objects in the
received signal can guarantee an enhanced differentiation level. Using the
mutual information as the information measure, we can write

$\displaystyle I(z_{j};l_{j})=\log_{2}(1+\rho_{j}),\quad\forall
j\in\mathcal{N}^{t}$ (14)

where $\rho_{j}$ is the SINR corresponding with the $j$th object. Obviously,
mutual information is a monotonically increasing function in $\rho_{j}$. The
SINR corresponding to object $j$ is given by

$\displaystyle\rho_{j}$
$\displaystyle=\frac{\delta_{j}Q_{j}\mathbf{v}^{H}_{j}\mathbf{w}_{j}\mathbf{w}^{H}_{j}\mathbf{v}_{j}}{\Sigma_{j_{\text{int}}}+\Sigma_{j_{n_{\text{s}}}}+\Sigma_{j_{n_{\text{fc}}}}},\quad\forall
j\in\mathcal{N}^{t}$ (15)

where the interference and noise variances are

$\displaystyle\Sigma_{j_{\text{int}}}$
$\displaystyle=\mathbf{v}^{H}_{j}\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{N}\delta_{i}Q_{i}\mathbf{w}_{i}\mathbf{w}^{H}_{i}\mathbf{v}_{j}$
(16) $\displaystyle\Sigma_{j_{n_{\text{s}}}}$
$\displaystyle=\mathbf{v}^{H}_{j}\mathbf{A}_{n}\mathbf{v}_{j},$ (17)
$\displaystyle\Sigma_{j_{n_{\text{fc}}}}$
$\displaystyle=\sigma^{2}_{\text{fc}}||\mathbf{v}_{j}||^{2}_{2},$ (18)

respectively. The equivalent sensor noise covariance matrix observed in
decoding the information of the $j$th object is denoted by
$\mathbf{A}_{n}\in\mathbb{C}^{KR\times KR}$ which is a block diagonal matrix
with the $k$th block represented by
$\alpha_{k}\sigma^{2}_{n}\mathbf{f}_{k}\mathbf{f}^{H}_{k}$.
Having MRT at the transmitter, we consider the following signal combining
strategies,

  1. A.

maximum-ratio combining (MRC): maximizes signal-to-noise ratio (SNR). This is
SINR-optimal at low interference regime.

  2. B.

zero-forcing (ZF): maximizes signal-to-interference ratio (SIR). This is SINR-
optimal at high interference regime.

  3. C.

minimum mean-squared error (MMSE): Signal-to-interference-plus-noise ratio
(SINR) optimal.

$\text{T}_{\text{1}}$$m_{1}$...$\text{T}_{\text{N}}$$m_{N}$$\text{C}_{\text{1}}$$m_{N+1}$...$\text{C}_{\text{N}^{{}^{\prime}}}$$m_{N+N^{{}^{\prime}}}$Fusion
Center...$\mathbf{w}_{1}(\boldsymbol{\alpha})$...$\mathbf{w}_{N}(\boldsymbol{\alpha})$$\mathbf{w}_{N+1}(\boldsymbol{\alpha})$...$\mathbf{w}_{N+N^{{}^{\prime}}}(\boldsymbol{\alpha})$
Figure 4: Virtual multiple access channel with interference (VMACI). Compared
to conventional MAC, this channel has extra degrees-of-freedom imposed by the
designable channels.

In the rest of the paper, these receivers are studied in details and efficient
optimization algorithms are proposed for transmit power and sensor
amplification factor optimization.

###  IV-A Maximum-Ratio Combining

Assuming MRC at the receiver, the following signal combining vector maximizes
SNR,

$\displaystyle\mathbf{v}^{(\text{MR})}_{j}=\mathbf{w}_{j}\quad\forall
j\in\mathcal{N}^{t},$ (19)

which is less complex for practical implementations, however does not consider
the destructive effect of interference in the signal combining phase.
Utilizing MRC, we will minimize the sum transmit power plus sum power
amplification at the sensors jointly.

###  IV-B Zero-Forcing

Here, we enforce the interference to zero while decoding the signal of the
$j$th object. This can be done in space-time by

$\displaystyle\mathbf{v}^{(\text{ZF})}_{j}$
$\displaystyle=\text{null}\\{\mathbf{w}_{1},...,\mathbf{w}_{j-1},\mathbf{w}_{j+1},...,\mathbf{w}_{N}\\}.$
(20)

The ZF combining vector spans the null-space of the interference terms in
(20), however the optimal combining vector in this null-space for the $j$th
object is the $j$th column of

$\displaystyle\mathbf{V}^{(\text{ZF})^{\star}}$
$\displaystyle=\mathbf{W}\left(\mathbf{W}^{H}\mathbf{W}\right)^{-1},$ (21)

where
$\mathbf{W}=\begin{bmatrix}\mathbf{w}_{1},...,\mathbf{w}_{N}\end{bmatrix}$.

###  IV-C Minimum Mean-Squared Error (MMSE)

Optimal linear MMSE receiver can be considered by post processing in the
direction that minimizes the mean of the squared error for the $j$th object.
Mathematically

$\displaystyle\mathbf{v}^{(\text{MMSE})}_{j}=\text{arg}\min_{\mathbf{v}_{j}}\
\mathbb{E}\\{\|l_{j}-z_{j}(\mathbf{v}_{j})\|^{2}_{2}\\},$ (22)

which yields the following solution

$\displaystyle\mathbf{v}^{(\text{MMSE})}_{j}=\left(\sum_{i=1}^{N}\delta_{i}Q_{i}\mathbf{w}_{i}\mathbf{w}^{H}_{i}+\boldsymbol{\Sigma}_{n}\right)^{-1}\mathbf{w}_{j}$
(23)

where
$\boldsymbol{\Sigma}_{n}=\mathbb{E}\\{\mathbf{n}^{{}^{\prime}}\mathbf{n}^{{}^{\prime
H}}\\}$. Hence, we require the knowledge of the signal amplification factors
in the sensors (i.e., $\alpha_{k},\ \forall k$) and transmit power $p_{j}$,
which in turn need to be optimized. Hence, the design parameters become inter-
connected, which in turn renders the problem to be non-convex. In section V an
iterative procedure is proposed for this type of receiver.

In what follows, we formulate the sum-power-plus-sum-amplification
minimization problem.

##  V Optimization Problem

In this section, we formulate sum-power minimization problems for different
receiver types under object SINR constraints. This problem is a non-convex
problem and it includes transmit power and power consumption at the sensors as
the optimization parameters. Utilizing maximum-ratio combining receiver, we
formulate an approximate problem for resolving the non-convexity, which is
improved iteratively. Moreover, the purposed algorithm captures the
interaction in optimizing transmit power and sensor power consumption through
joint optimization. This joint optimization algorithm is only applicable if
the receiver applies suboptimal MRC. Zero-forcing and optimal linear receiver
MMSE, are the other types of receivers, for which the transmit power and
sensors power consumption are optimized separately. This separate optimization
is due to the difficulty in jointly optimizing the transmit power and sensor
amplification factor for these types of receivers.

###  V-A Maximum-Ratio Combining

By exploiting maximum-ratio transmission (MRT) at the transmitter and maximum-
ratio combining (MRC) at the fusion center, the sum transmit power plus sum
power amplification minimization problem is formulated as

$\displaystyle\min_{p_{j},\forall
j{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in\mathcal{N}^{t},\alpha_{k}},\forall
k{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in\\{1,2,\ldots,K\\}}}\quad$
$\displaystyle\sum_{j=1}^{N^{t}}p_{j}+\sum_{k=1}^{K}\alpha_{k}$ (24) subject
to $\displaystyle\rho_{j}\geq\psi_{j},$ (24a)
$\displaystyle\sum_{j=1}^{N^{t}}p_{j}\leq P_{\text{max}},$ (24b)
$\displaystyle\alpha_{k}\leq\alpha_{\text{max}},$ (24c)

where the $j$th object SINR demand is defined by $\psi_{j}$. Furthermore, the
sum transmit power is restricted by $P_{\text{max}}$ and the maximum
amplification power of each sensor is limited by $\alpha_{\text{max}}$ as in
(24b) and (24c), respectively. Evidently, the objective function is affine,
however SINR constraints in (24a) produce a non-convex set. The SINR
expression for the $j$th object is written as

$\displaystyle\rho_{j}=\frac{\Sigma_{j_{\text{des}}}}{\Sigma_{j_{\text{int}}}+\Sigma_{j_{n_{\text{s}}}}+\Sigma_{j_{n_{\text{fc}}}}},$
(25)

where

$\displaystyle\Sigma_{j_{\text{des}}}$
$\displaystyle=\delta_{j}Q_{j}\left(\sum_{k=1}^{K}\alpha_{k}|g_{jk}|^{2}\|\mathbf{f}_{k}\|^{2}_{2}\right)^{2},$
(26) $\displaystyle\Sigma_{j_{\text{int}}}$ $\displaystyle=\sum_{i\neq
j}\delta_{i}Q_{i}\left|\sum_{k=1}^{K}\alpha_{k}g_{jk}g^{*}_{ik}\|\mathbf{f}_{k}\|^{2}_{2}\right|^{2},$
(27) $\displaystyle\Sigma_{j_{n_{\text{s}}}}$
$\displaystyle=\sum_{k=1}^{K}\sigma^{2}_{n_{k}}\alpha^{2}_{k}|g_{jk}|^{2}\|\mathbf{f}_{k}\|^{4}_{2},$
(28) $\displaystyle\Sigma_{j_{n_{\text{fc}}}}$
$\displaystyle=\sigma^{2}_{\text{fc}}\sum_{k=1}^{K}\alpha_{k}|g_{jk}|^{2}\|\mathbf{f}_{k}\|^{2}_{2}.$
(29)

Notice that, $\delta_{j},\ \forall j\in\mathcal{N}$ are functions of the
transmit power as defined in (11) and (12). Here, we assume the thermal noise
variance at the fusion center and sensors are equal, i.e.,
$\sigma^{2}=\sigma_{\text{fc}}=\sigma_{n_{k}},\ \forall k\in\mathcal{K}$.
Having the SINR for the $j$th object, the constraint (24a) can be reformulated
as

$\displaystyle\frac{\psi_{j}\left(\Sigma_{j_{\text{int}}}+\Sigma_{j_{n_{\text{s}}}}+\Sigma_{j_{n_{\text{fc}}}}\right)}{\Sigma_{j_{\text{des}}}}\leq
1,\ \forall j\in\mathcal{N}^{t},$ (30)

where $\Sigma_{j_{\text{des}}}$, $\Sigma_{j_{n_{\text{s}}}}$ and
$\Sigma_{j_{n_{\text{fc}}}}$ are posynomials in the optimization parameters,
i.e., $p_{j},\ \forall j\in\mathcal{N}^{t}$ and $\alpha_{k},\ \forall
k\in\mathcal{K}$.

######  Lemma 1.

$\Sigma_{j_{\text{int}}}$ is a signomial function in general. It can be a
posynomial if the following constraint holds, ($\forall i\neq j\ \text{and}\
\forall k\neq l,\forall\psi\in\mathbb{Z}$)

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2\psi\pi-\frac{\pi}{2}\leq\left(\measuredangle
g_{jk}g^{*}_{jl}g^{*}_{ik}g_{il}\right)\leq 2\psi\pi+\frac{\pi}{2}.}$ (31)

###### Proof.

The proof is provided in Appendix X-A. ∎

Here, we consider the general case, where $\Sigma_{j_{\text{int}}}$ is a
signomial. Thus, one can write $\Sigma_{j_{\text{int}}}$ as the difference of
two posynomials as follows

$\displaystyle\Sigma_{j_{\text{int}}}=\Sigma^{(1)}_{\text{int}}-\Sigma^{(2)}_{\text{int}},$
(32)

where both terms $\Sigma^{(1)}_{\text{int}}$ and $\Sigma^{(2)}_{\text{int}}$
are positive. Hence, by plugging (32) into the inequality constraint (30), we
obtain

$\displaystyle\frac{\psi_{j}\left(\Sigma^{(1)}_{j_{\text{int}}}+\Sigma_{j_{n_{\text{s}}}}+\Sigma_{j_{n_{\text{fc}}}}\right)}{\Sigma_{j_{\text{des}}}+\psi_{j}\Sigma^{(2)}_{j_{\text{int}}}}\leq
1,\ \forall j\in\mathcal{N}^{t}.$ (33)

The left hand-side of the inequality constraint (33) is the division of
posynomials, which can not be converted to a convex function. Problem (24) is
a signomial program (SP) [20], which can be converted to a complementary
geometric program (GP). This program allows upperbound constraint on the
division of two posynomials. The denominator of (33) is approximated by a
monomial function (known as condensation method [21]) based on the following
lower-bound

$\displaystyle\sum_{k}c_{k}\mu_{k}\geq\prod_{k}\mu^{c_{k}}_{k},\ \text{where}\
c_{k}\geq 0,\ \sum_{k}c_{k}=1$ (34)

This states the relationship between arithmetic and geometric mean. Lower-
bound on the denominator of (33) operates as an upper-bound on the whole
expression. By defining $\hat{\mu}_{k}=c_{k}\mu_{k}$ we get

$\displaystyle\sum_{k}\hat{\mu}_{k}\geq\prod_{k}\left(\frac{\hat{\mu}_{k}}{c_{k}}\right)^{c_{k}}$
(35)

Algorithm 1 MRC, Joint Optimization

1:Determine feasible $\alpha^{(0)}_{k},\ \forall k$ and $p^{(0)}_{j},\ \forall
j$,

2:Calculate $\hat{\mu}^{(0)}_{jk},\ \forall j,k$ for the given
$\alpha^{(0)}_{k},\ \forall k$ and $p^{(0)}_{j},\ \forall j$,

3:Determine $\Sigma^{(0)}_{j_{D}},\ \forall j$ as in (37),

4:Determine $c^{(0)}_{jk},\ \forall j,k$ from (39),

5:Set $q=1$,

6:Define
$E^{(q)}=\sum_{j=1}^{N^{t}}p^{(q)}_{j}+\sum_{k=1}^{K}\alpha^{(q)}_{k}$,
$E^{(-1)}=0$,

7:while $E^{(q-1)}-E^{(q-2)}$ large do

8: Plug the lower-bound in (38) into the constraint (30),

9: Solve the following geometric program (GP),
$\displaystyle\min_{p^{(q)}_{j},\alpha^{(q)}_{k},\ \forall j,k}$
$\displaystyle\ E^{(q)}$ (36) s. t.
$\displaystyle\frac{\psi_{j}\left(\Sigma_{\text{int}}+\Sigma_{n_{k}}+\Sigma_{n_{\text{fc}}}\right)}{\tilde{\Sigma}_{D}(c^{(q-1)}_{jk})}\leq
1,\ \forall j,$ (36a) $\displaystyle\sum_{j=1}^{N^{t}}p^{(q)}_{j}\leq
P_{\text{max}},$ (36b) $\displaystyle\alpha^{(q)}_{k}\leq\alpha_{\text{max}},$
(36c)

10: Solutions are $p^{(q)\star}_{j},\ \forall j$ and $\alpha^{(q)\star}_{k},\
\forall k,$

11: Determine
$E^{(q)}=\sum_{j=1}^{N^{t}}p^{(q)\star}_{j}+\sum_{k=1}^{K}\alpha^{(q)\star}_{k}$,

12: Determine $\hat{\mu}^{(q)}_{jk},\ \forall j,k$, $\Sigma^{(q)}_{j_{D}},\
\forall j$ and $c^{(q)}_{jk},\ \forall j,k$,

13: $q=q+1$.

14:end while

Now, we utilize this inequality in the SINR constraints of the $j$th object in
(33). The denominator of (33) is rewritten as the summation of monomials by

$\displaystyle\Sigma_{j_{D}}=\Sigma_{j_{\text{des}}}+\psi_{j}\Sigma^{(2)}_{j_{\text{int}}}=\sum_{k=1}^{K+K^{{}^{\prime}}}\hat{\mu}_{jk}$
(37)

where $\hat{\mu}_{jk}$ are the individual monomials. Furthermore,
$K^{{}^{\prime}}$ is the number of monomials that the posynomial
$\Sigma^{(2)}_{j_{\text{int}}}$ consists of, which can be quantified according
to lemma 1. Then from (35) and (37) we obtain,

$\displaystyle\Sigma_{j_{D}}=\sum_{k=1}^{K+K^{{}^{\prime}}}\hat{\mu}_{jk}\geq\prod_{k=1}^{K+K^{{}^{\prime}}}\left(\frac{\hat{\mu}_{jk}}{c_{jk}}\right)^{c_{jk}}=\tilde{\Sigma}_{j_{D}}(c_{jk}),$
(38)

where $\tilde{\Sigma}_{D}$ is a function of $c_{jk}$, which needs to be
optimized to fulfill the inequality with equality. For that, $c_{jk}$ must be
a function of $\alpha_{k},\ \forall k$ as

$\displaystyle c^{\star}_{jk}=\frac{\hat{\mu}_{jk}}{\Sigma_{j_{D}}}.$ (39)

Due to the inter-dependency of the optimization parameters, we optimize
$\alpha_{k},\ \forall k$ and $c_{jk}$ successively in an iterative fashion.
That means, $c_{jk},\ \forall j,k$ is optimized for the current iteration
based on the solution of $\hat{\mu}_{jk}$ in the previous iteration. The
procedure is explained in Algorithm 1 elaborately. Notice that the lower-bound
in (38) is the approximation of $\Sigma_{D}$ around any feasible $\alpha_{k},\
\forall k$, though sub-optimal. Hence, by improving $\alpha_{k}\ \forall k$
and $p_{j}\ ,\forall j\in\mathcal{N}^{t}$ at each iteration, $c_{jk}$ is
calculated, which in turn is utilized for the next iteration. The convergence
of the algorithm is numerically illustrated in section VIII.

Algorithm 2 MMSE, Separate Optimization

1:Initialize $\boldsymbol{\alpha}^{(0)}=\alpha_{\text{max}}\mathbf{1}_{K}$

2:Initialize $\mathbf{p}^{(0)}=\frac{P_{\text{max}}}{N^{t}}\mathbf{1}_{N}$

3:Set $q=1$

4:Define $E^{(q)}=\|\boldsymbol{\alpha}^{(q)}\|_{1}+\|\mathbf{p}^{(q)}\|_{1}$,
$E^{-1}=0$

5:Define $\delta^{(q)}_{i}$ as a function of $\mathbf{p}^{(q)}$ from (11) and
(12)

6:Define $\mathbf{w}^{(q)}_{i}$ as a function of $\boldsymbol{\alpha}^{(q)}$
from (7)

7:Define
$\boldsymbol{\Sigma}^{(q)}_{i}=\mathbb{E}\\{\mathbf{n}^{(q)^{\prime}}\mathbf{n}^{(q)^{\prime
H}}\\}$ as a function of $\boldsymbol{\alpha}^{(q)}$ from (8)

8:Set
$\mathbf{v}^{(q)}_{j}=\left(\sum_{i=1}^{N}\delta^{(q)}_{i}Q_{i}\mathbf{w}^{(q)}_{i}\mathbf{w}^{{(q)}^{H}}_{i}+\boldsymbol{\Sigma}^{(q)}_{n}\right)^{-1}\mathbf{w}^{(q)}_{j}$

9:while $E^{(q-1)}-E^{(q-2)}$ large do

10: $\mathbf{p}^{(q)}=\text{power-
opt}(\boldsymbol{\alpha}^{(q-1)},\mathbf{v}^{(q-1)}_{j})$: linear program

11: Update $\mathbf{v}^{(q-1)}_{j},\ \forall j$

12: $\boldsymbol{\alpha}^{(q)}=\text{amp-
opt}(\mathbf{p}^{(q)},\mathbf{v}^{(q-1)}_{j})$: signomial program

13: $q=q+1$

14:end while

$1$$3$$5$$7$$9$$11$$12$$12.2$$12.4$$12.6$$12.8$$13$$13.2$$q$ (iteration
index)$\min\ \sum_{j=1}^{N}p_{j}+\sum_{k=1}^{K}\alpha_{k}$ [dB] Figure 5:
Maximum-ratio transmission (MRT) at the transmitter and maximum-ratio
combining (MRC) at the fusion center.

###  V-B Minimum Mean-Squared Error (MMSE)

In this section, we investigate the performance of MMSE in minimizing the sum
of transmit power and sensor amplification. Here, we proceed with the
optimization of transmit power allocation and sensor amplification separately.
Notice that, sum-power minimization requires a power balancing at the
transmitter and sensors. Hence, optimal power allocation at the transmitter
becomes a function of sensors amplification solution. Moreover, MMSE receiver
vector depends on both transmit power and sensor amplification solutions. This
inter-dependency is resolved by iterative optimization. In this iterative
method, by initial guess for $\mathbf{p}=[p_{1},...,p_{N}]$ and
$\boldsymbol{\alpha}=[\alpha_{1},...,\alpha_{K}]$, the initial MMSE receiver
beamforming for the $j$th object $\mathbf{v}_{j},\ \forall j$ is obtained from
(23). Then, for a given initial sensor amplification and MMSE beamforming, the
transmit power is minimized to satisfy per object SINR constraints. Therefore,
the transmit power minimization problem is written as

$\displaystyle\text{power-opt}=\arg\min_{\mathbf{p}}\quad$
$\displaystyle\|\mathbf{p}\|_{1}$ (40) subject to
$\displaystyle\rho_{j}\geq\psi_{j},\ \forall j$ (40a)
$\displaystyle\|\mathbf{p}\|_{1}\leq P_{\text{max}},$ (40b)

where $\rho_{j},\ \forall j$ is the expression in (15), given
$\boldsymbol{\alpha}$ and $\mathbf{v}_{j},\ \forall j$. This optimization
problem is a linear program which can be solved efficiently by simplex method.
Now, the MMSE beamforming vector is updated by the initial guess of the sensor
amplification and the power allocation solution from (40). With this updated
MMSE receiver, the sensor amplification minimization problem is formulated as

$\displaystyle\text{amp-opt}=\arg\min_{\boldsymbol{\alpha}}\quad$
$\displaystyle\|\boldsymbol{\alpha}\|_{1}$ (41) subject to
$\displaystyle\rho_{j}\geq\psi_{j},\ \forall j$ (41a)
$\displaystyle\alpha_{k}\leq\alpha_{\text{max}},\ \forall k.$ (41b)

This problem can be formulated as a signomial problem (SP) and be solved
iteratively similar to problem (24). The MMSE receiver is now updated again
given the solution from (40) and (41). Hence, this procedure requires an inner
iteration for optimizing over sensor amplification, and an outer iteration for
updating power allocation and MMSE receiver. This procedure is explained in
Algorithm 2 in details. Notice that, this separate optimization algorithm can
be utilized for the case that the fusion center exploits ZF. In that case,
step 8 in Algorithm 2 is replaced by the expression in (20).

$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$8$$10$$12$$14$$16$$18$$20$$22$$24$$\psi$
(SINR demands)$\min\ \sum_{j=1}^{N}p_{j}+\sum_{k=1}^{K}\alpha_{k}$ [dB]MRC,
$r_{3}=2$MRC, $r_{3}=0.5$ZFMMSE, $r_{3}=2$MMSE, $r_{3}=0.5$ (a) opt. Tx-power,
max. amp.

$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$-8$$-4$$0$$4$$8$$12$$16$$20$$24$$\psi$ (SINR
demands)$\min\ \sum_{j=1}^{N}p_{j}+\sum_{k=1}^{K}\alpha_{k}$ [dB]MRC,
$r_{3}=2$MRC, $r_{3}=0.5$ZFMMSE, $r_{3}=2$MMSE, $r_{3}=0.5$ (b) opt. Tx-power,
opt. amp

Figure 6: (a) optimization of $p_{j},\ \forall j$, with sensor maximum
amplification, i.e., $\alpha_{k}=\alpha_{m}ax,\ \forall k$. (b) Successive
optimization of $p_{j},\ \forall j$ and $\alpha_{k},\ \forall k$.

##  VI Asymptotic Analysis:
fusion center with massive antenna array or Massive number of sensors?

Suppose that the system operator has the option to increase the number
antennas at the fusion center and/or to increase the number of sensors. As a
trade-off for this complexity burden, it is of crucial importance to analyze
the performance improvement. Recall that, the sensors observations are
amplified and forwarded to the fusion center in different time instants
(TDMA). Consider the channel matrix of the VMACI as

$\displaystyle\mathbf{W}=[\mathbf{w}_{1},\cdots,\mathbf{w}_{N}]$ (42)

where $\mathbf{w}_{i},\ \forall i\in\mathcal{N}$ are defined in (7). Notice
that, the information of the $i$th object, i.e., $l_{i}$, is received at the
fusion center over the channel vector $\mathbf{w}_{i}$.

######  Lemma 2.

The coherence value of matrix $\mathbf{W}$ approaches to zero, as $K\gg N$
holds and the number of sensors go to infinity, i.e., $K\rightarrow\infty$.

###### Proof.

The proof is provided in Appendix X-B. ∎

######  Lemma 3.

Maximum-ratio combining is the optimal receiver if the coherence value of
matrix $\mathbf{W}$ is zero, i.e., $\mu=0$.

###### Proof.

The proof is provided in Appendix X-C. ∎

From Lemma 3, by exploiting MRC at the fusion center, the approximate per-
object SNR is maximized as

$\displaystyle\digamma_{j}(\mathbf{p})=\max_{\boldsymbol{\alpha}}\quad$
$\displaystyle\frac{\Sigma_{j_{\text{des}}}}{\Sigma_{j_{n_{\text{s}}}}+\Sigma_{j_{n_{\text{fc}}}}},$
(43) subject to $\displaystyle\alpha_{k}\leq\alpha_{\text{max}},\ \forall k.$
(43a)

where $\Sigma_{j_{\text{des}}}$, $\Sigma_{j_{n_{\text{s}}}}$ and
$\Sigma_{j_{n_{\text{fc}}}}$ are functions of $\mathbf{p}$ and
$\boldsymbol{\alpha}$ as in (26), (28) and (29).

######  Lemma 4.

The objective function in problem (43) is maximized when the constraint holds
with equality, i.e., $\alpha_{k}=\alpha_{\text{max}},\ \forall k$

###### Proof.

The objective function in problem (43) is monotonically increasing in
$\alpha_{k},\forall k$ given that $\alpha_{k}>0,\ \forall k$. The proof can be
found in the Appendix X-D. ∎

From lemma 4, the achievable SNR per-object $\digamma_{j}(\mathbf{p}),\
\forall j$ are only functions of power allocation at the transmitter. This
power allocation problem is a linear program, which can be solved efficiently.

##  VII Material Characterization

So far, we have discussed the design of system parameters with allocating the
available resources optimally. In this task, we assumed that the second-order
moment of the objects are given. However, for robust classification purposes,
the second-order moments of the objects need to be estimated. Notice that the
estimation quality is a function of the INR after post-processing. Therefore,
in this section, we propose an iterative forward-backward
transmission/reception scheme. In this scheme, after per-object post-
processing, the estimate of the second-order moments of the objects are
provided to optimize the transmit power and sensor amplification. The
procedure is elaborated in Algorithm. 3.

Algorithm 3 Recursive Method for Waveform Optimization and object Second-order
Moment Acquisition

1:Initialize $Q^{(0)}_{j}=1,\ \forall j\in\mathcal{N}^{t}$,

2:Set $q=1$,

3:while $Q^{(q)}_{j}-Q^{(q-1)}_{j}$ do large

4: Determine $\mathbf{p}^{(q)}$ and $\boldsymbol{\alpha}^{(q)}$, from
Algorithm 1 or 2,

5: Estimate $Q^{(q)}_{j}$ from $z_{j}$ in (13),

6: $q=q+1$,

7:end while

$\displaystyle\frac{\partial\ln(\text{SNR}_{j})}{\partial\alpha_{k}}$
$\displaystyle=\frac{2b_{jk}}{\sum_{l=1}^{K}\alpha_{l}b_{jl}}-\frac{2\alpha_{k}b_{kj}\|\mathbf{f}_{k}\|^{2}+b_{jk}}{\sum_{l=1}^{K}{\alpha^{2}_{l}b_{jl}\|\mathbf{f}_{l}\|^{2}}+\sum_{l=1}^{K}\alpha_{l}b_{jl}}$
(44)
$\displaystyle=\frac{2b_{jk}\left(\sum_{l=1}^{K}{\alpha^{2}_{l}b_{jl}\|\mathbf{f}_{l}\|^{2}}+\sum_{l=1}^{K}\alpha_{l}b_{jl}\right)-\left(2\alpha_{k}b_{kj}\|\mathbf{f}_{k}\|^{2}+b_{jk}\right)\sum_{l=1}^{K}\alpha_{l}b_{jl}}{\left(\sum_{l=1}^{K}\alpha_{l}b_{jl}\right)\left(\sum_{l=1}^{K}{\alpha^{2}_{l}b_{jl}\|\mathbf{f}_{l}\|^{2}}+\sum_{l=1}^{K}\alpha_{l}b_{jl}\right)}$
(45)

##  VIII Numerical Results

In this section, we provide the simulation results for a two-object single-
clutter environment, i.e., $N^{t}=2,N^{c}=1$. The number of antennas at the
transmitter is assumed to be 4, i.e., $M=2$ and $M^{{}^{\prime}}=2$ (two
antennas per dimension). Moreover, the distance between the antennas is 1cm.
The fusion center is equipped with 10 antennas, i.e., $R=10$. We assume 3
objects are placed in space with following azimuth and elevation

$\displaystyle\boldsymbol{\theta}=[20\ 45\ 70],\quad\boldsymbol{\phi}=[40\ 30\
85].$ (46)

where, the first two objects are the objects of interest and the last object
is a clutter. The two objects are assumed to be in unit distance, i.e.,
$r_{1}=r_{2}=1$, and the clutter is located once at $r_{3}=0.5$ and once at
$r_{3}=2$. The sensors amplification factor is assumed to be equal to 2. The
noise variance at the sensors and fusion center are assumed to be equal to
$0.5$, i.e., $\sigma^{2}_{\text{fc}}=\sigma^{2}_{n_{k}}=0.5,\ \forall k$.
Considering MRC at the receiver, the sum-power minimization problem is solved
iteratively under per-object SINR constraints. This problem is a signomial
program, which is turned to a geometric program according to Algorithm 1 and
solved iteratively until convergence. The convergence of the algorithm is
depicted in Fig. 5 for a per-object SINR constraint equal to 1, where we
observe the fast convergence. Assuming maximum amplification at the sensors,
the transmit power minimization problem is also a signomial problem, which is
solved in a similar manner as elaborate in Algorithm 1. The minimum sum-power
consumption for this case (maximum amplification) is compared to the case with
optimal amplification in Fig. 6(a) and 6(b). In these figures, we observe that
MRC is optimal when the SINR demands are sufficiently low. Notice that, at
sufficiently low SINR demands, joint optimization of the transmit power and
the sensor amplification is crucial. Higher SINR demands can not be satisfied
by MRC. Since, ZF outperforms MRC as the interference increases, it is
efficient to zero force the interference. Intuitively, by zero-forcing
processing at the fusion center, interference-free signaling dimensions
becomes less than the number of available dimensions $KR$. This is due to
reserving $N-1$ dimension for null steering. This leaves us with $KR-N+1$
signaling dimensions. Therefore, comparing ZF and MRC we notice the trade-off
between sacrificing some dimensions in expense of obtaining interference-free
dimensions, and utilizing all dimensions. Finally, we observe that in expense
of extra complexity, the MMSE receiver outperforms MRC and ZF at the fusion
center. As can be seen from Fig. 6(a) and Fig. 6(b), the performance of MMSE
receiver approaches the performance of ZF at sufficiently high SINR demands,
only when the clutter is located at larger distance compared to the objects.
Having no clutters in the sensing environment, i.e.,
$\mathcal{N}^{c}=\emptyset$, ZF is the optimal receiver at sufficiently high
SINR demands.

##  IX Conclusion

In this paper, we have introduced an active sensing system for material
characterization purposes. In this system, several sensors are deployed in the
sensing environment to forward their observation to the fusion center for
complex computation. These observations include the response of the objects to
the incident signal. For optimal system design, the sum-transmit power-plus-
sum sensor amplifications is minimized under per-object SINR constraints. The
optimization problems are non-convex. We provide efficient algorithms to
obtain good sub-optimal solutions in polynomial time.

##  X Appendix

###  X-A Proof of Lemma 1

The expression in (27) is rewritten as

$\displaystyle\sum_{i\neq
j}\delta_{i}Q_{i}\left(\underbrace{\sum_{k=1}^{K}\alpha_{k}\|\mathbf{f}_{k}\|^{2}g_{jk}g^{*}_{ik}\sum_{l=1}^{K}\alpha_{l}\|\mathbf{f}_{l}\|^{2}g^{*}_{jl}g_{il}}_{\Gamma(\boldsymbol{\alpha)}}\right),$
(47)

where we define the expression in the braces as $\Gamma(\boldsymbol{\alpha})$,
with $\boldsymbol{\alpha}=[\alpha_{1},\cdots,\alpha_{K}]$. Notice that,
$\Gamma(\boldsymbol{\alpha})$ is the summation of $K^{2}$ monomial functions.
The monomials corresponding with $k=l$ have real positive values, since
$g_{jk}g^{*}_{ik}g_{jl}^{*}g_{il}=|g_{jk}|^{2}|g_{ik}|^{2}\geq 0$. The
monomials corresponding with $k\neq l$ are

$\displaystyle\alpha_{k}\alpha_{l}\|\mathbf{f}_{k}\|^{2}\|\mathbf{f}_{l}\|^{2}\left(g_{jk}g^{*}_{ik}g_{jl}^{*}g_{il}+g^{*}_{jk}g_{ik}g_{jl}g_{il}^{*}\right),$
$\displaystyle=2\alpha_{k}\alpha_{l}\|\mathbf{f}_{k}\|^{2}\|\mathbf{f}_{l}\|^{2}\Re\\{g_{jk}g^{*}_{ik}g_{jl}^{*}g_{il}\\},$
(48)

which do not necessarily yield a positive value. Hence,
$\Sigma_{j_{\text{int}}}$ is a signomial function in $p_{j}\ \forall
j\in\mathcal{N}^{t}$ and $\alpha_{k},\ k\in\mathcal{K}$. The expression
$\Sigma_{j_{int}}$ is posynomial, if (48) is positive. This condition can be
written as

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2\psi\pi-\frac{\pi}{2}\leq\measuredangle\left(g_{jk}g^{*}_{jl}g^{*}_{ik}g_{il}\right)\leq
2\psi\pi+\frac{\pi}{2}\quad\forall\psi\in\mathbb{Z}}.$ (49)

###  X-B Proof of Lemma 2

The coherence value of matrix $\mathbf{W}$ is defined as

$\displaystyle\mu=\max_{i,j\in\mathcal{N},i\neq j}\
\frac{|\langle\mathbf{w}_{i},\mathbf{w}_{j}\rangle|}{\|\mathbf{w}_{i}\|_{2}\|\mathbf{w}_{j}\|_{2}}.$
(50)

Notice that,
$|\langle\mathbf{w}_{i},\mathbf{w}_{j}\rangle|=|\mathbf{w}^{H}_{i}\mathbf{w}_{j}|,\
\forall i,j\in\mathcal{N},\ i\neq j$. Hence, we prove that the inner product
of the columns of matrix $\mathbf{W}$ approaches to zero as
$K\rightarrow\infty$ and $K\gg N$. To do this, we write

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{|\langle\mathbf{w}_{i},\mathbf{w}_{j}\rangle|}{\|\mathbf{w}_{i}\|_{2}\|\mathbf{w}_{j}\|_{2}}=\frac{1}{KR\sqrt{\gamma_{i}\gamma_{j}}}\sum_{l=1}^{KR}w^{*}_{il}w_{jl},}$
(51)

where $\gamma_{i}=\frac{\mathbf{w}^{H}_{i}\mathbf{w}_{i}}{KR}>0$. Then, from
the law of large numbers we obtain

$\displaystyle\lim_{KR\rightarrow\infty}\frac{1}{KR{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sqrt{\gamma_{i}\gamma_{j}}}}\sum_{l=1}^{KR}w^{*}_{il}w_{jl}=\frac{\mathbb{E}\\{w_{i}^{*}w_{j}\\}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sqrt{\gamma_{i}\gamma_{j}}}},$
(52)

where $w_{i}$ is a random variable that represents the realizations of the
vector $\mathbf{w}_{i}$. These random variables are fully-correlated for the
case that $K=1$, since $\mathbf{w}_{i}=\mathbf{w}_{j}$. This does not imply
any limitation on the number of antennas at the fusion center, i.e., $R$ can
be infinitely large. However by increasing the number of sensors, the
correlation between $w_{i}$ and $w_{j}$ decreases. Hence, for
$K\rightarrow\infty$,

$\displaystyle\mathbb{E}\\{w_{i}^{*}w_{j}\\}=0.$ (53)

which results in zero coherence value for the matrix $\mathbf{W}$, if the
number of columns ($N$) is limited, i.e., $N\ll K$.

###  X-C Proof of Lemma 3

If $\mu=0$, the inner product of the columns of $\mathbf{W}$ approaches to
zero. Thus, by using MRC ($\mathbf{v}_{j}=\mathbf{w}_{j},\ \forall
j\in\mathcal{N}^{t}$) at the fusion center, we obtain

$\displaystyle\mathbf{v}^{H}_{j}\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{N}\delta_{i}Q_{i}\mathbf{w}_{i}\mathbf{w}^{H}_{i}\mathbf{v}_{j}=0,\quad\forall
j\in\mathcal{N}^{t}.$ (54)

This shows that the MRC suppresses the interference, while the SNR is
maximized. Hence, MRC is the optimal receiver if the coherence value of matrix
$\mathbf{W}$ is zero.

###  X-D Proof of Lemma 4

In order to prove Lemma 4, we show that
$\text{SNR}_{j}=\frac{\Sigma_{j_{\text{des}}}}{\Sigma_{j_{n_{\text{s}}}}+\Sigma_{j_{n_{\text{fc}}}}}$
is monotonically increasing in $\alpha_{k}>0,\ \forall k$. To this end, we
compute the first derivative of $\ln(\text{SNR}_{j})$ and show that it is
strictly positive. We define $b_{jk}=|g_{jk}|^{2}\|\mathbf{f}_{k}\|^{2}$.
Then, using the assumption of the system model, i.e.,
$\sigma_{\text{fc}}=\sigma_{nk}$, the derivative of $\ln(\text{SNR})$ is given
at the top of the page in (44). One can easily rewrite (44) as (45). Since
$b_{jk}>0$, for $\alpha_{k}>0$, the denominator in (45) is always positive.
Now, we need to show that the numerator is strictly positive. The numerator of
(45) can be written as

$\displaystyle
2b_{jk}\sum_{l=1}^{K}{\alpha^{2}_{l}b_{jl}\|\mathbf{f}_{l}\|^{2}}+\left(b_{jk}-2\alpha_{k}b_{jk}\|\mathbf{f}_{k}\|^{2}\right)\sum_{l=1}^{K}\alpha_{l}b_{jl}$
(55)

Since $b_{jk}\sum_{l=1}^{K}\alpha_{l}b_{jl}>0$, (55) is lower bounded by
$2\zeta(\boldsymbol{\alpha})$, where

$\displaystyle\zeta(\boldsymbol{\alpha})=b_{jk}\sum_{l=1}^{K}{\alpha^{2}_{l}b_{jl}\|\mathbf{f}_{l}\|^{2}}-\alpha_{k}b_{jk}\|\mathbf{f}_{k}\|^{2}\sum_{l=1}^{K}\alpha_{l}b_{jl}.$
(56)

Now, we get the lower-bound of $\zeta(\boldsymbol{\alpha})$ by determining the
extremum of the function w.r.t. $\alpha_{l},\ \forall
l{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\neq
k}$, and showing that the single extrema is a minimum. Hence,

$\displaystyle\frac{\partial\zeta(\boldsymbol{\alpha})}{\partial\alpha_{l}}$
$\displaystyle=2{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}b_{jk}}\alpha_{l}b_{jl}\|\mathbf{f}_{l}\|^{2}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}b_{jk}}\alpha_{k}b_{jl}\|\mathbf{f}_{k}\|^{2}=0$
(57)
$\displaystyle\Rightarrow\alpha_{l}=\frac{\alpha_{k}\|\mathbf{f}_{k}\|^{2}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}\|\mathbf{f}_{l}\|^{2}}\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{
if }l\neq k}$ (58)
$\displaystyle\frac{\partial^{2}\zeta(\boldsymbol{\alpha})}{\partial\alpha_{l}^{2}}=2b_{jl}\|\mathbf{f}_{l}\|^{2}>0$
(59)

Notice that $\frac{\partial\zeta(\boldsymbol{\alpha})}{\partial\alpha_{k}}=0$.
Therefore, we obtain the minimum of $\zeta(\boldsymbol{\alpha})$ by replacing
$\alpha_{l}=\frac{\alpha_{k}\|\mathbf{f}_{k}\|^{2}}{2\|\mathbf{f}_{l}\|^{2}},\
\forall
l{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\neq
k}$ into (56). Doing this, we get the minimum value for
$\zeta(\boldsymbol{\alpha})$ which is equal to zero, i.e.,
$\zeta^{\star}(\boldsymbol{\alpha})=0$. Therefore,
$\zeta(\boldsymbol{\alpha})$ is lower-bounded by zero and as a result (55) is
strictly positive, and thereby the expression in (45) is strictly positive.
Hence, we conclude that the per-object SNR is monotonically increasing in
$\alpha_{k},\ \forall k$ and $\alpha^{\star}_{k}=\alpha_{\text{max}},\ \forall
k$.

## References

  * [1] G. Alirezaei, O. Taghizadeh, and R. Mathar, “Optimum Power Allocation With Sensitivity Analysis for Passive Radar Applications,” _IEEE Sensors Journal_ , vol. 14, no. 11, pp. 3800–3809, Nov 2014.
  * [2] T. Wang, R. C. de Lamare, and A. Schmeink, “Joint receiver design and power allocation strategies for multihop wireless sensor networks,” in _2011 8th International Symposium on Wireless Communication Systems_ , Nov 2011, pp. 86–90.
  * [3] H. Godrich, A. P. Petropulu, and H. V. Poor, “Power Allocation Strategies for Target Localization in Distributed Multiple-Radar Architectures,” _IEEE Transactions on Signal Processing_ , vol. 59, no. 7, pp. 3226–3240, July 2011.
  * [4] Y. Shen, W. Dai, and M. Z. Win, “Optimal power allocation for active and passive localization,” in _2012 IEEE Global Communications Conference (GLOBECOM)_ , Dec 2012, pp. 3713–3718.
  * [5] L. Guo, Z. Ning, Q. Song, L. Zhang, and A. Jamalipour, “A QoS-Oriented High-Efficiency Resource Allocation Scheme in Wireless Multimedia Sensor Networks,” _IEEE Sensors Journal_ , vol. 17, no. 5, pp. 1538–1548, March 2017.
  * [6] D. Zhang, Z. Chen, M. K. Awad, N. Zhang, H. Zhou, and X. S. Shen, “Utility-Optimal Resource Management and Allocation Algorithm for Energy Harvesting Cognitive Radio Sensor Networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 34, no. 12, pp. 3552–3565, Dec 2016.
  * [7] Y. Yang and R. S. Blum, “MIMO radar waveform design based on mutual information and minimum mean-square error estimation,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 43, no. 1, pp. 330–343, January 2007.
  * [8] A. Leshem, O. Naparstek, and A. Nehorai, “Information theoretic adaptive radar waveform design for multiple extended targets,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 1, no. 1, pp. 42–55, June 2007.
  * [9] S. Jeong, S. Khalili, O. Simeone, A. Haimovich, and J. Kang, “Multistatic cloud radar systems: joint sensing and communication design,” _Transactions on Emerging Telecommunications Technologies_ , vol. 27, no. 5, pp. 716–730, 2016, ett.3020.
  * [10] S. M. Razavizadeh, M. Ahn, and I. Lee, “Three-dimensional beamforming: A new enabling technology for 5g wireless networks,” _IEEE Signal Processing Magazine_ , vol. 31, no. 6, pp. 94–101, Nov 2014.
  * [11] Y. H. Nam, B. L. Ng, K. Sayana, Y. Li, J. Zhang, Y. Kim, and J. Lee, “Full-dimension MIMO (FD-MIMO) for next generation cellular technology,” _IEEE Communications Magazine_ , vol. 51, no. 6, pp. 172–179, June 2013.
  * [12] J. Koppenborg, H. Halbauer, S. Saur, and C. Hoek, “3D beamforming trials with an active antenna array,” in _2012 International ITG Workshop on Smart Antennas (WSA)_ , March 2012, pp. 110–114.
  * [13] M. Karaman, I. O. Wygant, Ã. Oralkan, and B. T. Khuri-Yakub, “Minimally Redundant 2-D Array Designs for 3-D Medical Ultrasound Imaging,” _IEEE Transactions on Medical Imaging_ , vol. 28, no. 7, pp. 1051–1061, July 2009.
  * [14] A. C. Dhanantwari, S. Stergiopoulos, L. Song, C. Parodi, F. Bertora, P. Pellegretti, and A. Questa, “An efficient 3D beamformer implementation for real-time 4D ultrasound systems deploying planar array probes,” in _IEEE Ultrasonics Symposium, 2004_ , vol. 2, Aug 2004, pp. 1421–1424 Vol.2.
  * [15] S. Gezici, Z. Tian, G. B. Giannakis, H. Kobayashi, A. F. Molisch, H. V. Poor, and Z. Sahinoglu, “Localization via ultra-wideband radios: a look at positioning aspects for future sensor networks,” _IEEE Signal Processing Magazine_ , vol. 22, no. 4, pp. 70–84, July 2005.
  * [16] K. J. Francis, P. Rajalakshmi, and S. Channappayya, “Wavelet domain frequency interpolation for photo-acoustic tomography,” in _2014 International Conference on Medical Imaging, m-Health and Emerging Communication Systems (MedCom)_ , Nov 2014, pp. 6–9.
  * [17] M. Qu, M. Mehrmohammadi, R. Truby, K. Homan, and S. Emelianov, “Magneto-photo-acoustic imaging using dual-contrast agent,” in _2010 IEEE International Ultrasonics Symposium_ , Oct 2010, pp. 511–514.
  * [18] H. A. A. Naam, M. O. Idrees, A. Awad, O. S. Abdalsalam, and F. Mohamed, “Non invasive blood glucose measurement based on Photo-Acoustic Spectroscopy,” in _2015 International Conference on Computing, Control, Networking, Electronics and Embedded Systems Engineering (ICCNEEE)_ , Sept 2015, pp. 1–4.
  * [19] H. Bertoni, _Radio Propagation for Modern Wireless Systems_ , ser. Wireless communications. Prentice Hall PTR, 2000.
  * [20] S. Boyd, S.-J. Kim, L. Vandenberghe, and A. Hassibi, “A tutorial on geometric programming,” _Optimization and Engineering_ , vol. 8, no. 1, pp. 67–127, Apr 2007.
  * [21] M. Chiang, “Geometric Programming for Communication Systems,” _Commun. Inf. Theory_ , vol. 2, no. 1/2, pp. 1–154, Jul. 2005. [Online]. Available: http://dx.doi.org/10.1516/0100000005


# Reflection, transmission and surface susceptibility tensor in two-
dimensional materials

Luca Dell’Anna1 Yu He1,2 Michele Merano1<EMAIL_ADDRESS>1
Dipartimento di Fisica e Astronomia G. Galilei, Universit$\grave{a}$ degli
studi di Padova, via Marzolo 8, 35131 Padova, Italy
2 College of Physics, Sichuan University, Chengdu 610064, China.

###### Abstract

In a recent experiment, the out-of-plane surface susceptibility of a single-
layer two-dimensional atom crystal in the visible spectrum has been measured.
This susceptibility gives a measurable contribution to the reflectivity of
two-dimensional materials. Here we provide a complete theoretical description
of the reflective properties, considering incoming $s$ and $p$ polarized plane
waves at any angle of incidence on the crystal, computing local, reflected and
transmitted electromagnetic fields. We finally connect the microscopic
polarizability to both the in-plane and the out-of-plane macroscopic surface
susceptibilities.

††preprint: APS/123-QED

##  I INTRODUCTION

Two-dimensional (2D) crystals are the thinnest materials that can be produced
[1, 2]. They are composed of only one atomic plane like graphene or Boron-
Nitride [3, 4], or one molecular plane like transition metal dichalcogenides
[5, 6]. This planar configuration make them highly anisotropic in the
direction perpendicular to the crystal plane. For instance, they can have in-
plane macroscopic dimensions [7, 8, 9, 10] while keeping an out-of-plane
thickness of the order of one atom. This particular geometry suggests that
also their optical response should be anisotropic, at least in the vertical
direction to the crystal plane. Ab-initio calculations indeed predict this
anisotropy for a 2D crystal [11, 12, 13]. Notwithstanding the enormous
progresses in thin films and 2D materials optical characterization [14, 15,
16, 17, 18, 19, 20, 21], it turns out that a measurement of the out-of-plane
optical constants of a 2D crystal is still a difficult task [22].

These materials are usually deposited on a substrate whose optical response is
added to that of the monolayer. This is enough to hide the contribution coming
from the out-of-plane optical constants of the single-layer crystal [22]. In
thin-film optics, the sensitivity of a measurement to vertical anisotropy is
dependent on the path length of the light through the film, which is extremely
limited when we deal with these atomically thin crystals [22]. As a
consequence, optical experiments of monolayers, deposited on some substrate,
have access only to the in-plane optical constants of the samples that are
studied [6, 14, 15, 16, 19, 23, 20, 21].

Only recently an experiment successfully measured the out-of-plane optical
constants of 2D crystals, namely graphene and monolayer $\rm MoS_{2}$ [24].
The substrate contribution was eliminated by a complete immersion of these
single-layer crystals in a transparent polymer. Experimental data say that the
out-of-plane surface susceptibility ($\chi_{\bot}$) is a measurable and finite
quantity different from the in-plane surface susceptibility ($\chi_{\|}$),
while the out-of-plane surface conductivity was zero within the experimental
errors.

One might expect that a 2D crystal with atomic thickness along the vertical
direction has a vanishing $\chi_{\bot}$, being $\chi_{\bot}$ a macroscopic
quantity. Indeed some of the most cited papers about optics in 2D crystals
assume a null $\chi_{\bot}$ [25, 26, 27, 28]. All the other papers describe
the linear optical response of a monolayer assuming isotropy even in the
vertical direction [23, 22, 6]. The aim of this work is to prove that
$\chi_{\bot}$ is indeed finite and different from $\chi_{\|}$ in a rigorous
analytical way.

We will start with a microscopic description of a 2D crystal, as done in Ref.
[29], and then will connect it to the macroscopic one, by means of a procedure
similar to that used in Ref. [30]. The results in Refs. [29, 30] were limited
to normal incidence, therefore, they had no access to $\chi_{\bot}$.

##  II MICROSCOPIC THEORY

Figure 1:  $s$\- ($p$-) polarized electric (magnetic) field incident with an
angle $\theta$ on a single-layer two-dimensional atomic crystal, modeled by
dipoles placed on a regular lattice in the $x$-$y$ plane. The two half-spaces,
separated by the crystal, are denoted by the labels $1$ and $2$.

We consider an insulating 2D crystal, in the vacuum, composed by atoms with
isotropic polarizability $\alpha$ placed on a 2D Bravais lattice (Fig.1). A
plane wave propagating in a generic direction specified by the unit vector
$\hat{s}=\cos\phi\sin\theta\,\hat{x}-\sin\phi\sin{\theta}\,\hat{y}+\cos{\theta}\,\hat{z}$
which is defined by the polar angles $\theta$ (the incident angle) and $\phi$
(the azimuthal angle). Where $\hat{x}$, $\hat{y}$, $\hat{z}$ are unit vectors
along the $x$, $y$ and $z$ directions. The incident electric field varies
periodically in time $\vec{E}_{i}(t)=\vec{E}_{i}\,e^{i\omega t}$ where
$\omega$ is the frequency from which one can define $k=\omega/c=2\pi/\lambda$,
where $\lambda$ is the wavelength of the incident field. The charge
distribution in the crystal is distorted by the electric field, which
generates oscillations of the electric dipoles placed at the Bravais lattice
sites. The total electric field which acts on a single dipole is the local
field $\vec{E}_{loc}$. It is the superposition of the incident field
$\vec{E}_{i}$ and the contributions arising from all the other dipoles.
Because of translation invariance, we make the hypothesis that the local field
at any lattice point has the same magnitude and have the same frequency of the
incident field but with a different phase.

###  II.1 Computation of the local field

We compute the local field acting on the dipole at the origin of our reference
system (Fig.1). The contributions to this local field are given by the fields
coming from the oscillating dipoles placed at the points labelled by the
integers $m$ and $n$, which are given by [31, 32, 33]

$\displaystyle\vec{E}_{mn}=\,$ $\displaystyle\frac{{\cal
P}_{mn}}{4\pi\epsilon_{0}r_{mn}^{3}}\left[3(\tilde{\vec{p}}\cdot\hat{{r}}_{mn})\hat{{r}}_{mn}-\tilde{\vec{p}}\right.$
$\displaystyle\left.-\frac{1}{c^{2}}(\hat{{r}}_{mn}\times\ddot{\vec{p}})\times\hat{{r}}_{mn}\right]$
(1)

where
$\tilde{\vec{p}}=\vec{p}(t-\frac{r_{mn}}{c})+\frac{r_{mn}}{c}\dot{\vec{p}}(t-\frac{r_{mn}}{c})$.
$\vec{p}$ is the induced dipole moment at each lattice point, $\epsilon_{0}$
is the vacuum permittivity, $c$ is the velocity of light in vacuum, and
$\vec{{\emph{r}}}_{mn}=r_{mn}(\cos\vartheta_{mn}\hat{{\emph{x}}}+\sin\vartheta_{mn}\hat{{\emph{y}}})$
while $\hat{{\emph{r}}}_{mn}=\vec{{\emph{r}}}_{mn}/r_{mn}$. The phase in the
prefactor is given by

${\cal P}_{mn}=e^{i(\omega t-k\,r_{mn})}\,e^{-ik\delta_{mn}}$ (2)

with a phase shift $\delta_{mn}$ whose modulus is the distance of the point
$(r_{mn}\cos\vartheta_{mn},r_{mn}\sin\vartheta_{mn},0)$ from the plane wave
crossing the origin of the reference frame with equation
$\hat{s}\cdot\vec{r}=0$. As a result, in quite general terms, the phase shift
$\delta_{mn}$ can be written as

$\delta_{mn}=r_{mn}\sin\theta\cos(\vartheta_{mn}+\phi)$ (3)

The components of $\vec{{\emph{E}}}_{mn}$ are, therefore,

$\displaystyle E_{mn,x}=\frac{e^{i\omega
t}}{4\pi\epsilon_{0}\,r_{mn}^{3}}\,e^{-ik\,(r_{mn}+\delta_{mn})}\left\\{p_{x}\left[\left(3\cos^{2}\vartheta_{mn}-1\right)(1+ik\,r_{mn})+k^{2}r_{mn}^{2}\sin^{2}\vartheta_{mn}\right]\right.$
(4)
$\displaystyle\phantom{------------------}+\left.p_{y}\cos\vartheta_{mn}\sin\vartheta_{mn}\left[3(1+k\,r_{mn})-k^{2}r_{mn}^{2}\right]\right\\}$
$\displaystyle E_{mn,y}=\frac{e^{i\omega
t}}{4\pi\epsilon_{0}\,r_{mn}^{3}}\,e^{-ik\,(r_{mn}+\delta_{mn})}\left\\{p_{y}\left[\left(3\sin^{2}\vartheta_{mn}-1\right)(1+ik\,r_{mn})+k^{2}r_{mn}^{2}\cos^{2}\vartheta_{mn}\right]\right.$
(5)
$\displaystyle\phantom{------------------}+\left.p_{x}\cos\vartheta_{mn}\sin\vartheta_{mn}\left[3(1+k\,r_{mn})-k^{2}r_{mn}^{2}\right]\right\\}$
$\displaystyle E_{mn,z}=\frac{e^{i\omega
t}}{4\pi\epsilon_{0}\,r_{mn}^{3}}\,e^{-ik\,(r_{mn}+\delta_{mn})}p_{z}\left[k^{2}r_{mn}^{2}-1-ik\,r_{mn}\right]$
(6)

Without loss of generality we will take $\phi=\pi/2$ (see Fig. 1) so that the
phase shift reduces to

$\delta_{mn}=-r_{mn}\sin\theta\sin\vartheta_{mn}$ (7)

Summing over $\vartheta_{mn}$, the last terms in Eqs. (4), (5) are zero, since
those terms are odd under $\vartheta_{mn}\rightarrow\vartheta_{mn}+\pi$ for
our choice of the reference frame. As a result the sum over the contributions
to the local filed coming from all the dipoles is parallel to the dipole
sitting at the origin.
For the square lattice $r_{mn}=a\sqrt{m^{2}+n^{2}}$, where $a$ is the lattice
parameter, and $\tan\vartheta_{mn}=n/m$, while for the triangular lattice
$r_{mn}=a\sqrt{(m+n/2)^{2}+3n^{2}/4}$ and
$\tan\vartheta_{mn}=\sqrt{3}\,n/(2m+n)$. After summing over the dipoles, and
using

$\vec{{\emph{p}}}=\alpha\epsilon_{0}\,\vec{{\emph{E}}}_{loc}$ (8)

we get numerical evidences that

$\displaystyle\sum^{\prime}_{m,n}E_{mn,x}=\frac{\alpha}{4\pi
a^{3}}\left(C_{0}+ika\frac{C_{1}}{\cos\theta}\right)E_{loc,x}$ (9)
$\displaystyle\sum^{\prime}_{m,n}E_{mn,y}=\frac{\alpha}{4\pi
a^{3}}\Big{(}C_{0}+ika\,C_{1}\cos\theta\Big{)}E_{loc,y}$ (11)
$\displaystyle\sum^{\prime}_{m,n}E_{mn,z}=\frac{\alpha}{4\pi
a^{3}}\Big{(}-2C_{0}+ika\,C_{1}f(\theta)\Big{)}E_{loc,z}$

where the prime on the summation symbol indicates that the origin, $n=m=0$, is
excluded from the sum. $C_{0}$ is the static result already reported in our
previous paper [29] ($C_{0}\approx 4.517$ for the square lattice and
$C_{0}\approx 5.517$ for the triangular lattice) and $C_{1}=-2\pi Na^{2}$,
where $N$ is the density of dipoles ($N=1/a^{2}$ for the square lattice,
$N=2/\sqrt{3}a^{2}$ for the triangular lattice).
We notice that the $\theta$-dependence vanishes in the static limit, namely in
the long-wavelength limit, $k\rightarrow 0$, and appears in the three field
components with a $1/\cos\theta$, $\cos\theta$ and a combination of them

$f(\theta)=\frac{1}{\cos\theta}-\cos\theta=\sin\theta\tan\theta.$ (12)

As shown in Fig. 2, the $\theta$-dependence in Eq. (11) is the same for either
square and triangular lattices. Analogous checks have been performed to verify
$1/\cos\theta$ and $\cos\theta$ in Eqs. (9) and (11) respectively.

Figure 2: $f(\theta)$ as a function of the angle of incidence $\theta$. The
square dots refer to the square lattice and the triangular dots refer to
triangle lattice. The dots are numerical results, summing over almost $10^{4}$
sites, while the line is given by Eq. (12).

The local field at one site is due to the incident field and the contributions
from the other oscillating dipoles. We has, therefore, to solve the following
equation

$\displaystyle\vec{E}_{loc}=\vec{E}_{i}+\sum_{m,n}^{\prime}\vec{E}_{mn}$ (13)

We can now easily find the local field $\vec{E}_{loc}$ from Eqs. (9), (11),
(11) getting

$\displaystyle E_{loc,x}=\frac{E_{ix}}{1-\frac{C_{0}\alpha}{4\pi
a^{3}}+ik\frac{\alpha N}{2}\frac{1}{\cos\theta}}$ (16) $\displaystyle
E_{loc,y}=\frac{E_{iy}}{1-\frac{C_{0}\alpha}{4\pi a^{3}}+ik\frac{\alpha
N}{2}\cos\theta}$ $\displaystyle
E_{loc,z}=\frac{E_{iz}}{1+\frac{C_{0}\alpha}{2\pi a^{3}}+ik\frac{\alpha
N}{2}\sin\theta\tan\theta}$

These equations show that the local field, and therefore also the macroscopic
polarization vector, as we will see, are connected to the incident field via a
diagonal tensor.

##  III MACROSCOPIC THEORY

In this section, we will find the expression for the reflected and transmitted
fields and the relation in between the microscopic polarizability $\alpha$ and
the macroscopic surface susceptibility of the crystal. Because of momentum
conservation, the reflected field will propagate along the unit vector
$\hat{s}_{r}=-\sin{\theta}\ \hat{y}-\cos{\theta}\ \hat{z}$, and the
transmitted field along the same unit vector $\hat{s}$ of the incident field.

The first macroscopic quantity that we can introduce is the surface
polarization $\vec{P}$

$\vec{P}=N\vec{p}=\epsilon_{0}\alpha N\vec{E}_{loc}$ (17)

using Eq. (8) and it is related to the incident field from Eqs. (16)-(16). The
surface susceptibility is a diagonal tensor that, similarly to that connecting
in between $\vec{E}_{loc}$ and $\vec{E}_{i}$, it connects $\vec{P}$ with the
macroscopic field $\vec{E}$,

$\begin{pmatrix}P_{x}\\\ P_{y}\\\
P_{z}\end{pmatrix}=\epsilon_{0}\begin{pmatrix}\chi_{xx}&0&0\\\
0&\chi_{yy}&0\\\ 0&0&\chi_{zz}\end{pmatrix}\begin{pmatrix}E_{x}\\\ E_{y}\\\
E_{z}\end{pmatrix}$ (18)

In order to find the expressions for $\chi_{xx}$, $\chi_{yy}$, $\chi_{zz}$ and
$E_{x}$ $E_{y}$, $E_{z}$, we will consider the case of $s$ and $p$ polarized
incident waves.

Figure 3:  For an $s$ polarized wave, the total macroscopic electromagnetic
field is the sum of the incident field plus the radiation reaction field due
to the surface electric current associated to the surface polarization
$P_{x}$. For a $p$ polarized wave, the total macroscopic electromagnetic field
is the sum of the incident field plus the radiation reaction fields due to the
surface electric current associated to $P_{y}$ and the surface magnetic
current (see text) associated to $P_{z}$. The radiation reaction fields
propagate along $\hat{s_{r}}$ (reflected field) and $\hat{s}$ (transmitted
field).

###  III.1 s-wave

For an $s$ polarized incident wave the electric field oscillates along the
$x$-direction (see Fig.1). We choose $\vec{E}_{i}=E_{i}\,\hat{x}$,
$\vec{E}_{r}=E_{r}\,\hat{x}$ and $\vec{E}_{t}=E_{t}\,\hat{x}$ for the
incident, reflected and transmitted beams. Equations (16)-(18) imply that only
the components $E_{loc,x}$, $P_{x}$ and $E_{x}$ are non-null in this case.
This allows us to fix $\chi_{xx}$. The polarization vector varies in time
since $\vec{E}_{loc}\propto e^{i\omega t}$, therefore, along the crystal
plane, it gives rise to an in-plane surface electric current $J_{P_{x}}$,
which is given by

$J_{P_{x}}=\partial P_{x}/\partial t=i\omega P_{x}=i\omega\epsilon_{0}\alpha
NE_{loc,x}$ (19)

This surface current generates a macroscopic electromagnetic field that
propagates in two different directions. The first goes along the unit vector
$\hat{s}_{r}$ and it corresponds to the reflected field. The second one goes
along the unit vector $\hat{s}$ and, in superposition with the incident field,
it gives rise to the transmitted field (see Fig. 3). The total macroscopic
field is given by the incident field plus the macroscopic field generated by
the dipoles. It must obey the following boundary conditions [34, 35, 36, 25,
26, 37]

$\displaystyle\hat{z}\wedge\left(\vec{E}_{2}-\vec{E}_{1}\right)=0$ (20)
$\displaystyle\hat{z}\wedge\left(\vec{H}_{2}-\vec{H}_{1}\right)=J_{P_{x}}\,\hat{x}$
(21)

where the subscript 1 (2) refers to the limit of the macroscopic field when it
approaches the crystal from above (below). These equations show that the
component of the macroscopic electric field generated by the dipoles is
continuous across the crystal surface whereas the component of the macroscopic
magnetic field due to the dipoles is discontinuous. Using

$\eta\vec{H}=\hat{s}\wedge\vec{E}$ (22)

where $\eta$ is the wave impedance of the vacuum, and since
$\hat{z}\wedge\hat{s}\wedge\vec{E}=-\vec{E}\cos\theta$, we can write Eqs.
(20), (21) and (18) in terms of the electric field components

$\displaystyle E_{i}+E_{r}$ $\displaystyle=E_{t}$ (23) $\displaystyle
E_{i}-E_{r}$ $\displaystyle=E_{t}+\frac{\eta J_{P_{x}}}{\cos\theta}$ (24)
$\displaystyle E_{t}$ $\displaystyle=\frac{P_{x}}{\epsilon_{0}\chi_{xx}}$ (25)

We have three equations with three unknown variables: $E_{r}$, $E_{t}$, and
$\chi_{xx}$. Defining the reflection and the transmission coefficients as
$r_{s}=E_{r}/E_{i}$ and $t_{s}=E_{t}/E_{i}$, respectively, we obtain

$\displaystyle r_{s}$
$\displaystyle=\frac{-ik\chi_{xx}}{ik\chi_{xx}+2\cos\theta}$ (26)
$\displaystyle t_{s}$ $\displaystyle=1+r_{s}$ (27) $\displaystyle\chi_{xx}$
$\displaystyle=\frac{4\pi a^{3}N\alpha}{4\pi a^{3}-C_{0}\alpha}$ (28)

The Fresnel coefficients for the $s$ polarized wave are the same already
obtained in Ref. [28]. The surface susceptibility is the same obtained for
normal incidence [29, 30], showing that non-local effects are excluded in our
theory.

Our microscopic theory allows us also to compute the radiation reaction
electric field acting on the dipole at the origin of our reference frame [30],
getting

$\displaystyle E_{R_{x}}$ $\displaystyle=\frac{ikaC_{1}\alpha}{4\pi
a^{3}\cos\theta}\,E_{loc,x}=-\frac{\eta J_{P_{x}}}{2\cos\theta}$ (29)

This is clearly a macroscopic quantity as well as $J_{P_{x}}$. We now show
that it corresponds to the macroscopic field generated by the oscillating
dipoles. Using Eq. (29) in Eq. (24) we can rewrite Eqs. (23)-(25) getting

$\displaystyle E_{r}$ $\displaystyle=E_{R_{x}}$ $\displaystyle E_{t}$
$\displaystyle=E_{i}+E_{R_{x}}$ $\displaystyle\chi_{xx}$
$\displaystyle=\frac{P_{x}}{\epsilon_{0}(E_{i}+E_{R_{x}})}$ (30)

finding that $E_{R_{x}}$ is equal to $E_{r}$ and, in superposition with
$E_{i}$, it gives the transmitted field. For that reason we can identify
$E_{R_{x}}$ with the macroscopic electric field generated by the oscillating
dipoles that propagates along $\hat{s}$ and $\hat{s}_{r}$ directions (see Fig.
3), while, from the last equation, the full macroscopic electric field on the
crystal plane can be identified with the sum of the radiation reaction field
and the incident field, namely $E_{x}=E_{i}+E_{R_{x}}$.

###  III.2 p-wave

For a $p$ polarized wave the incident, reflected and transmitted magnetic
fields oscillate along $x$-direction (see Fig.1). We choose, therefore,
$\vec{H}_{i}=-H_{i}\,\hat{x}$, $\vec{H}_{r}=-H_{r}\,\hat{x}$ and
$\vec{H}_{t}=-H_{t}\,\hat{x}$. Since we computed the local electric field, it
is useful to first write explicitly the incident, reflected and transmitted
electric fields

$\displaystyle\vec{E}_{i}$
$\displaystyle=E_{i}\cos{\theta}\,\hat{y}+E_{i}\sin{\theta}\,\hat{z}=E_{iy}\,\hat{y}+E_{iz}\,\hat{z}$
$\displaystyle\vec{E}_{r}$
$\displaystyle=-E_{r}\cos{\theta}\,\hat{y}+E_{r}\sin{\theta}\,\hat{z}=E_{ry}\,\hat{y}+E_{rz}\,\hat{z}$
$\displaystyle\vec{E}_{t}$
$\displaystyle=E_{t}\cos{\theta}\,\hat{y}+E_{t}\sin{\theta}\,\hat{z}=E_{ty}\,\hat{y}+E_{tz}\,\hat{z}$
(31)

where $E_{i}=\eta H_{i}$, $E_{r}=\eta H_{r}$, $E_{t}=\eta H_{t}$. Equations
(16)-(18) imply that we can fix separately $\chi_{yy}$ and $\chi_{zz}$. In
this case we have one varying macroscopic surface polarization along the
crystal plane ($P_{y}$) and one perpendicular to the crystal plane ($P_{z}$).
These two components induce two surface currents $J_{P_{y}}$ and $J_{P_{z}}$
respectively. They generate two macroscopic fields propagating along the
$\hat{s}$ direction and two along the $\hat{s_{r}}$ direction. According to
the superposition principle the total macroscopic field is the sum of the
incident field and the macroscopic fields generated by the two currents
considered separately (see Fig. 3). The expression of the in-plane electric
current is

$J_{P_{y}}=\partial P_{y}/\partial t=i\omega P_{y}=i\omega\epsilon_{0}\alpha
NE_{loc,y}$ (32)

In the radiation zone the electromagnetic field, due to an electric dipole
oscillating in the $\hat{z}$ direction, is identical to an electromagnetic
field due to a magnetic dipole oscillating in the $-\hat{x}$ direction [32].
Hence, a polarization $P_{z}$ generates an out-of-plane electric surface
current equivalent to an in-plane magnetic surface current, responsible for
the discontinuity of the macroscopic electric field [34, 35, 36, 37]. Using
$P_{z}\propto E_{loc,z}\propto E_{i,z}$ and $\vec{E}_{i}\propto e^{i(\omega
t+ky\sin\theta)}$, we have

$J_{P_{z}}\,\hat{x}=-\frac{1}{\epsilon_{0}}\hat{z}\wedge\vec{\nabla}P_{z}=\frac{ik\,\sin{\theta}}{\epsilon_{0}}{P_{z}}\,\hat{x}$
(33)

We must solve two separate sets of boundary conditions for the macroscopic
field

$\begin{array}[]{ >{\displaystyle}l >{\displaystyle{}}l @{\hspace{1em}} | @{\hspace{1em}} >{\displaystyle}l }&\hat{z}\wedge\left(\vec{E}_{2}-\vec{E}_{1}\right)=0\hfil\hskip 10.00002pt\vrule\hskip 10.00002pt&\hat{z}\wedge\left(\vec{E}_{2}-\vec{E}_{1}\right)=J_{P_{z}}\,\hat{x}\\\ \vspace{1em}\displaystyle&\hat{z}\wedge\left(\vec{H}_{2}-\vec{H}_{1}\right)=J_{P_{y}}\,\hat{y}\hfil\hskip 10.00002pt\vrule\hskip 10.00002pt&\hat{z}\wedge\left(\vec{H}_{2}-\vec{H}_{1}\right)=0\\\ \end{array}$ (34)

In terms of the magnetic field components, Eqs. (34) and (18) become

$\begin{array}[]{ >{\displaystyle}r >{\displaystyle{}}l @{\hspace{0.1em}} | @{\hspace{0.2em}} >{\displaystyle}l }H_{i}-H_{r\|}&=H_{t\|}\hfil\hskip 1.00006pt\vrule\hskip 1.99997pt&H_{i}-H_{r\bot}=H_{t\bot}+\frac{J_{P_{z}}}{\eta\cos{\theta}}\\\ H_{i}+H_{r\|}&=H_{t\|}+J_{P_{y}}\hfil\hskip 1.00006pt\vrule\hskip 1.99997pt&H_{i}+H_{r\bot}=H_{t\bot}\\\ \eta\cos{\theta}\ H_{t\|}&=\frac{P_{y}}{\epsilon_{0}\chi_{yy}}\hfil\hskip 1.00006pt\vrule\hskip 1.99997pt&\eta\sin{\theta}\ H_{t\bot}=\frac{P_{z}}{\epsilon_{0}\chi_{zz}}\\\ \end{array}$ (35)

where $H_{r\|}$ ($H_{r\bot}$) is the contributions to the reflected field due
to the in-plane (out-of-plane) surface current and $H_{t\|}$ and $H_{t\bot}$
are the total macroscopic magnetic fields immediately below the crystal for
the two cases considered separately. We have, then, the following solutions

$\begin{array}[]{ >{\displaystyle}r >{\displaystyle{}}l @{\hspace{1em}} | @{\hspace{0.3em}} >{\displaystyle}l }\hfil&\chi_{yy}=\frac{4\pi a^{3}N\alpha}{4\pi a^{3}-C_{0}\alpha}\hfil\hskip 10.00002pt\vrule\hskip 3.00003pt&\chi_{zz}=\frac{2\pi a^{3}N\alpha}{2\pi a^{3}+C_{0}\alpha}\\\ \hfil&\frac{H_{r\|}}{H_{i}}=\frac{ik\chi_{yy}\cos{\theta}}{ik\chi_{yy}\cos{\theta}+2}\hfil\hskip 10.00002pt\vrule\hskip 3.00003pt&\frac{H_{r\bot}}{H_{i}}=\frac{-ik\chi_{zz}\sin{\theta}\tan{\theta}}{ik\chi_{zz}\sin{\theta}\tan{\theta}+2}\\\ \hfil&\frac{H_{t\|}}{H_{i}}=\frac{2}{ik\chi_{yy}\cos{\theta}+2}\hfil\hskip 10.00002pt\vrule\hskip 3.00003pt&\frac{H_{t\bot}}{H_{i}}=\frac{2}{ik\chi_{zz}\sin{\theta}\tan{\theta}+2}\\\ \end{array}$ (36)

This finding is a very interesting result because it shows that our crystal
has an isotropic in-plane surface susceptibility
($\chi_{\|}=\chi_{xx}=\chi_{yy}$) different from the out-of-plane
susceptibility ($\chi_{\bot}=\chi_{zz}$), namely we are dealing with an
uniaxial crystal. We can also evaluate the ratio $\chi_{\bot}/\chi_{\|}$.
Choosing the value of $a=0.26$ nm, for atom polarizabilities $\alpha$
typically varying in the range between 1 and 30 $\rm cm^{-3}$, the
$\chi_{\bot}/\chi_{\|}$ ratio varies in between 0.93 and 0.1 (0.94 and 0.18)
for a triangular lattice (for a square lattice). This values are of the same
order of magnitude of what has been observed [24].

The full reflected field from the crystal is $H_{r}=H_{r\|}+H_{r\bot}$ (see
Fig. 3), so that we can finally derive the reflection coefficient for 2D
crystals, $r_{p}={H_{r}}/{H_{i}}$, getting

$r_{p}=\frac{ik\chi_{yy}\cos{\theta}}{ik\chi_{yy}\cos{\theta}+2}-\frac{ik\chi_{zz}\sin{\theta}\tan{\theta}}{ik\chi_{zz}\sin{\theta}\tan{\theta}+2}$
(37)

We compute now the radiation-reaction electric field acting on the dipole at
the origin along the $y$ and $z$ directions

$\displaystyle E_{R_{y}}$ $\displaystyle=\frac{ikaC_{1}\alpha\cos\theta}{4\pi
a^{3}}E_{loc,y}=-\frac{\eta J_{y}\cos\theta}{2}$ (38)

and

$\displaystyle E_{R_{z}}$
$\displaystyle=\frac{ikaC_{1}\alpha(\sin\theta)^{2}}{4\pi
a^{3}\cos{\theta}}E_{loc_{z}}=-\frac{ik(\sin\theta)^{2}}{2\epsilon_{0}\cos{\theta}}P_{z}$
$\displaystyle=-\frac{\tan\theta}{2\epsilon_{0}}\frac{\partial}{\partial
y}P_{z}=-\frac{\tan\theta}{2}J_{P_{z}}.$ (39)

and we show that these macroscopic quantities correspond to the macroscopic
fields generated by the oscillating dipoles. We can rewrite Eqs. (35) in terms
of the electric fields

$\begin{array}[]{ >{\displaystyle}r >{\displaystyle{}}l @{\hspace{1em}} | @{\hspace{1em}} >{\displaystyle}l }\hfil&E_{iy}+E_{ry\|}=E_{ty\|}\hfil\hskip 10.00002pt\vrule\hskip 10.00002pt&E_{iz}+E_{rz\bot}=E_{tz\bot}\\\ \hfil&E_{iy}-E_{ry\|}=E_{ty\|}-2E_{R_{y}}\hfil\hskip 10.00002pt\vrule\hskip 10.00002pt&E_{iz}-E_{rz\bot}=E_{tz\bot}-2E_{R_{z}}\\\ \hfil&\epsilon_{0}\chi_{yy}E_{ty\|}=P_{y}\hfil\hskip 10.00002pt\vrule\hskip 10.00002pt&\epsilon_{0}\chi_{zz}E_{ty\bot}={P_{z}}\\\ \end{array}$ (40)

which give the following solutions

$\begin{array}[]{ >{\displaystyle}r >{\displaystyle{}}l @{\hspace{1em}} | @{\hspace{1em}} >{\displaystyle}l }\hfil&E_{ry\|}=E_{R_{y}}\hfil\hskip 10.00002pt\vrule\hskip 10.00002pt&E_{rz\bot}=E_{R_{z}}\\\ \hfil&E_{ty\|}=E_{iy}+E_{R_{y}}\hfil\hskip 10.00002pt\vrule\hskip 10.00002pt&E_{tz\bot}=E_{iz}+E_{R_{z}}\\\ \hfil&\chi_{yy}=\frac{P_{y}}{\epsilon_{0}(E_{iy}+E_{R_{y}})}\hfil\hskip 10.00002pt\vrule\hskip 10.00002pt&\chi_{zz}=\frac{P_{z}}{\epsilon_{0}(E_{iz}+E_{R_{z}})}\\\ \end{array}$ (41)

We can interpret these equations in the following way. The in-plane (out-of-
plane) electric surface current $J_{P_{y}}$ ($J_{P_{z}}$) generates a
macroscopic electromagnetic field propagating along the $\hat{s}$ and
$\hat{s}_{r}$ directions. We identify $E_{R_{y}}$ ($E_{R_{z}}$) as the
electric field component along the $y$-direction ($z$-direction) of this
macroscopic field. From Eqs. (III.2)), we can write the first two equations on
the left-hand-side of Eqs. (41) as

$\displaystyle H_{r\|}$
$\displaystyle=-\frac{E_{R_{y}}}{\eta\cos{\theta}}=-H_{R\|}$ $\displaystyle
H_{t\|}$
$\displaystyle=H_{i}+\frac{E_{R_{y}}}{\eta\cos{\theta}}=H_{i}+H_{R\|}$ (42)

and the first two equations on the right-hand-side of Eqs. (41) as

$\displaystyle H_{r\bot}$
$\displaystyle=\frac{E_{R_{z}}}{\eta\sin{\theta}}=H_{R\bot}$ $\displaystyle
H_{t\bot}$
$\displaystyle=H_{i}+\frac{E_{R_{z}}}{\eta\sin{\theta}}=H_{i}+H_{R\bot}$ (43)

The full reflected field (see Fig. 3) is, therefore,

$H_{r}=-H_{R\|}+H_{R\bot}$ (44)

while the full transmitted field is given by the superposition of the incident
field and the two macroscopic fields propagating in the $\hat{s}$ direction
(see Fig. 3), generated by the two currents,

$H_{t}=H_{i}+H_{R\|}+H_{R\bot}$ (45)

We can finally derive the transmission coefficient, defined by
$t_{p}=H_{t}/H_{i}$, getting

$t_{p}=1-\frac{ik\chi_{yy}\cos{\theta}}{ik\chi_{yy}\cos{\theta}+2}-\frac{ik\chi_{zz}\sin{\theta}\tan{\theta}}{ik\chi_{zz}\sin{\theta}\tan{\theta}+2}$
(46)

Finally the equations in the last line of Eqs. (41) provide the macroscopic
electric fields $E_{y}$ and $E_{z}$. They are related to the radiation-
reaction electric fields. Along the $y$-direction we have
$E_{y}=E_{i}\cos{\theta}+E_{R_{y}}$ while along the $z$-direction we have
$E_{z}=E_{i}\sin{\theta}+E_{R_{z}}$.

##  IV CONCLUSION

We have shown that a 2D crystal, microscopically composed by atoms with an
isotropic polarizability $\alpha$, has both a macroscopic $\chi_{\|}$ and a
macroscopic $\chi_{\bot}$ with $\chi_{\bot}\neq\chi_{\|}$ and $\chi_{\bot}\neq
0$, in contrast to some theoretical models, generally used to derive
reflection and transmission coefficients for these materials, which assume
either a null $\chi_{\bot}$ [26, 25, 27, 28] or isotropy [23, 6, 15, 16, 19].
Moreover we have confirmed that $\chi_{\bot}<\chi_{\|}$ and that the ratio
$\chi_{\bot}/\chi_{\|}$ is compatible with experimental observations.

Our microscopic theory provides both the local and the radiation-reaction
fields when an incoming electromagnetic plane wave shines on a 2D material.
The radiation-reaction fields coincide with the macroscopic fields scattered
by the crystal, and they contribute to the reflected and the transmitted
radiations. This analysis allows us, then, to derive the complete set of
Fresnel coefficients and their dependence on the angle of incidence.

In particular we describe the reflection of an incident $p$ wave, as obtained
by the superposition of the scattered electromagnetic fields due to an in-
plane electric and an in-plane magnetic surface currents. The first current is
related to a macroscopic surface polarization oscillating in the crystal
plane, the second one to a macroscopic surface polarization oscillating
perpendicularly to the crystal plane. Combining the boundary conditions for
these two currents we provide a simple expression for the macroscopic electric
fields along the $y$ and $z$ directions, which turn to be the sum of the
related components of the incident field with the radiation-reaction field
along the same directions. For an $s$ polarized wave, only an in-plane
electric surface current is involved, so the reflected field is given by the
scattered electromagnetic field due to only an in-plane current and the
macroscopic electric field is the sum of the incident field plus the radiation
reaction field.

###### Acknowledgements.

Y.H acknowledges financial support from China Scholarship Council. L.D was in
charge of the microscopic theory, M.M of the macroscopic theory, all the
authors discussed the physical model and wrote the paper.

## References

  * Novoselov  _et al._ [2004] K. S. Novoselov, A. K. Geim, S. V. Morozov, D.-e. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric field effect in atomically thin carbon films, Science 306, 666 (2004).
  * Novoselov  _et al._ [2005] K. Novoselov, D. Jiang, F. Schedin, T. Booth, V. Khotkevich, S. Morozov, and A. Geim, Two-dimensional atomic crystals, Proc. Natl. Acad. Sci. U.S.A. 102, 10451 (2005).
  * Gorbachev  _et al._ [2011] R. V. Gorbachev, I. Riaz, R. R. Nair, R. Jalil, L. Britnell, B. D. Belle, E. W. Hill, K. S. Novoselov, K. Watanabe, T. Taniguchi, A. K. Geim, and P. Blake, Hunting for monolayer boron nitride: Optical and raman signatures, Small 7, 465 (2011).
  * Wickramaratne  _et al._ [2018] D. Wickramaratne, L. Weston, and C. G. Van de Walle, Monolayer to bulk properties of hexagonal boron nitride, J. Phys. Chem. C 122, 25524 (2018).
  * Mak  _et al._ [2010] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Atomically Thin ${\mathrm{MoS}}_{2}$: A New Direct-Gap Semiconductor, Phys. Rev. Lett. 105, 136805 (2010).
  * Li  _et al._ [2014a] Y. Li, A. Chernikov, X. Zhang, A. Rigosi, H. M. Hill, A. M. van der Zande, D. A. Chenet, E.-M. Shih, J. Hone, and T. F. Heinz, Measurement of the optical dielectric function of monolayer transition-metal dichalcogenides: ${\mathrm{MoS}}_{2}$, $\mathrm{Mo}\mathrm{S}{\mathrm{e}}_{2}$, ${\mathrm{WS}}_{2}$, and $\mathrm{WS}{\mathrm{e}}_{2}$, Phys. Rev. B 90, 205422 (2014a).
  * Huang  _et al._ [2015] Y. Huang, E. Sutter, N. N. Shi, J. Zheng, T. Yang, D. Englund, H.-J. Gao, and P. Sutter, Reliable exfoliation of large-area high-quality flakes of graphene and other two-dimensional materials, ACS Nano 9, 10612 (2015).
  * Magda  _et al._ [2015] G. Z. Magda, J. Peto, G. Dobrik, C. Hwang, L. P. Biro, and L. Tapaszto, Exfoliation of large-area transition metal chalcogenide single layers, Sci. Rep. 5, 1 (2015).
  * Desai  _et al._ [2016] S. B. Desai, S. R. Madhvapathy, M. Amani, D. Kiriya, M. Hettick, M. Tosun, Y. Zhou, M. Dubey, J. W. Ager III, D. Chrzan, and A. Javey, Gold-mediated exfoliation of ultralarge optoelectronically-perfect monolayers, Adv. Mater. 28, 4053 (2016).
  * Huang  _et al._ [2020] Y. Huang, Y.-H. Pan, R. Yang, L.-H. Bao, L. Meng, H.-L. Luo, Y.-Q. Cai, G.-D. Liu, W.-J. Zhao, Z. Zhou, _et al._ , Universal mechanical exfoliation of large-area 2d crystals, Nat. Commun. 11, 1 (2020).
  * Trevisanutto  _et al._ [2010] P. E. Trevisanutto, M. Holzmann, M. Côté, and V. Olevano, Ab initio high-energy excitonic effects in graphite and graphene, Phys. Rev. B 81, 121405 (2010).
  * Guilhon  _et al._ [2019] I. Guilhon, M. Marques, L. Teles, M. Palummo, O. Pulci, S. Botti, and F. Bechstedt, Out-of-plane excitons in two-dimensional crystals, Phys. Rev. B 99, 161201 (2019).
  * Liu  _et al._ [2020] X. Liu, J. Yang, and W. Guo, Semiempirical van der Waals method for two-dimensional materials with incorporated dielectric functions, Phys. Rev. B 101, 045428 (2020).
  * Chang  _et al._ [2014] Y.-C. Chang, C.-H. Liu, C.-H. Liu, Z. Zhong, and T. B. Norris, Extracting the complex optical conductivity of mono-and bilayer graphene by ellipsometry, Appl. Phys. Lett. 104, 261909 (2014).
  * Liu  _et al._ [2014] H. L. Liu, C. C. Shen, S. H. Su, C. L. Hsu, M. Y. Li, and L. J. Li, Optical properties of monolayer transition metal dichalcogenides probed by spectroscopic ellipsometry, Appl. Phys. Lett. 105, 201905 (2014).
  * Li  _et al._ [2014b] W. Li, A. G. Birdwell, M. Amani, R. A. Burke, X. Ling, Y.-H. Lee, X. Liang, L. Peng, C. A. Richter, J. Kong, D. J. Gundlach, and N. V. Nguyen, Broadband optical properties of large-area monolayer cvd molybdenum disulfide, Phys. Rev. B 90, 195434 (2014b).
  * Stenzel  _et al._ [2015] O. Stenzel _et al._ ,  _The physics of thin film optical spectra_ (Springer, 2015).
  * Tompkins and Hilfiker [2015] H. G. Tompkins and J. N. Hilfiker,  _Spectroscopic Ellipsometry: Practical Application to Thin Film Characterization_ (Momentum Press, 2015).
  * Morozov and Kuno [2015] Y. V. Morozov and M. Kuno, Optical constants and dynamic conductivities of single layer mos2, mose2, and wse2, Appl. Phys. Lett. 107, 083103 (2015).
  * Jayaswal  _et al._ [2018] G. Jayaswal, Z. Dai, X. Zhang, M. Bagnarol, A. Martucci, and M. Merano, Measurement of the surface susceptibility and the surface conductivity of atomically thin mos2 by spectroscopic ellipsometry, Opt. Lett. 43, 703 (2018).
  * Elliott  _et al._ [2020] J. D. Elliott, Z. Xu, P. Umari, G. Jayaswal, M. Chen, X. Zhang, A. Martucci, M. Marsili, and M. Merano, Surface susceptibility and conductivity of ${\mathrm{mos}}_{2}$ and ${\mathrm{wse}}_{2}$ monolayers: A first-principles and ellipsometry characterization, Phys. Rev. B 101, 045414 (2020).
  * Nelson  _et al._ [2010] F. Nelson, V. Kamineni, T. Zhang, E. Comfort, J. Lee, and A. Diebold, Optical properties of large-area polycrystalline chemical vapor deposited graphene by spectroscopic ellipsometry, Appl. Phys. Lett. 97, 253110 (2010).
  * Blake  _et al._ [2007] P. Blake, E. Hill, A. Castro Neto, K. Novoselov, D. Jiang, R. Yang, T. Booth, and A. Geim, Making graphene visible, Appl. Phys. Lett. 91, 063124 (2007).
  * Xu  _et al._ [2021] Z. Xu, D. Ferraro, A. Zaltron, N. Galvanetto, A. Martucci, L. Sun, P. Yang, Y. Zhang, Y. Wang, Z. Liu, _et al._ , Optical detection of the susceptibility tensor in two-dimensional crystals, Commun. Phys. 4, 1 (2021).
  * Falkovsky and Pershoguba [2007] L. Falkovsky and S. Pershoguba, Optical far-infrared properties of a graphene monolayer and multilayer, Phys. Rev. B 76, 153410 (2007).
  * Hanson [2008] G. W. Hanson, Dyadic green’s functions and guided surface waves for a surface conductivity model of graphene, J. Appl. Phys. 103, 064302 (2008).
  * Zhan  _et al._ [2013] T. Zhan, X. Shi, Y. Dai, X. Liu, and J. Zi, Transfer matrix method for optics in graphene layers, J. Phys: Condens Matter 25, 215301 (2013).
  * Merano [2016] M. Merano, Fresnel coefficients of a two-dimensional atomic crystal, Phys. Rev. A 93, 013832 (2016).
  * Dell’Anna and Merano [2016] L. Dell’Anna and M. Merano, Clausius-mossotti lorentz-lorenz relations and retardation effects for two-dimensional crystals, Phys. Rev. A 93, 053808 (2016).
  * Merano [2017] M. Merano, Role of the radiation-reaction electric field in the optical response of two-dimensional crystals, Ann. Phys. 529, 1700062 (2017).
  * [31] R. P. Feynmann, Lectures on physics, 1st ed.(addison-wesley,boston, ma, 1964), pp. 21–9 [formula (21.26) has an error of sign].
  * Jackson [1999] J. D. Jackson, Classical electrodynamics, 3rd ed. (wiley, new york, 1998), p. 411. (1999).
  * Born and Wolf [2013] M. Born and E. Wolf,  _Principles of optics: electromagnetic theory of propagation, interference and diffraction of light_ (Elsevier, 2013).
  * Senior and Volakis [1987] T. B. Senior and J. L. Volakis, Sheet simulation of a thin dielectric layer, Radio Sci. 22, 1261 (1987).
  * Idemen [1990] M. Idemen, Universal boundary relations of the electromagnetic field, J. Phys.Soc. Jpn. 59, 71 (1990).
  * Kuester  _et al._ [2003] E. F. Kuester, M. A. Mohamed, M. Piket-May, and C. L. Holloway, Averaged transition conditions for electromagnetic fields at a metafilm, IEEE Trans. Antennas Propag. 51, 2641 (2003).
  * Majérus  _et al._ [2018] B. Majérus, E. Dremetsika, M. Lobet, L. Henrard, and P. Kockaert, Electrodynamics of two-dimensional materials: Role of anisotropy, Phys. Rev. B 98, 125419 (2018).