# T3DRIS: Advancing Conformal RIS Design through
In-depth Analysis of Mutual Coupling Effects

Placido Mursia, Francesco Devoti, Marco Rossanese, Vincenzo Sciancalepore,
Gabriele Gradoni, Marco Di Renzo, Xavier Costa-Pérez P. Mursia, F. Devoti, M.
Rossanese and V. Sciancalepore, are with NEC Laboratories Europe, 69115
Heidelberg, Germany (email: name.surname@neclab.eu).
G. Gradoni is with University of Surrey, Stag Hill, University Campus,
Guildford GU2 7XH, UK (email: g.gradoni@surrey.ac.uk).
M. Di Renzo is with Université Paris-Saclay, CNRS, CentraleSupélec,
Laboratoire des Signaux et Systèmes, 91190 Gif-sur-Yvette, France (email:
marco.di-renzo@universite-paris-saclay.fr).
X. Costa-Pérez is with i2cat, ICREA, 08034 Barcelona, Spain and NEC
Laboratories Europe (email: xavier.costa@neclab.eu).
Email of the corresponding author<EMAIL_ADDRESS>
###### Abstract

This paper presents a theoretical and mathematical framework for the design of
a conformal reconfigurable intelligent surface (RIS) that adapts to non-planar
geometries, which is a critical advancement for the deployment of RIS on non-
planar and irregular surfaces as envisioned in smart radio environments.
Previous research focused mainly on the optimization of RISs assuming a
predetermined shape, while neglecting the intricate interplay between shape
optimization, phase optimization, and mutual coupling effects. Our
contribution, the T3DRIS framework, addresses this fundamental problem by
integrating the configuration and shape optimization of RISs into a unified
model and design framework, thus facilitating the application of RIS
technology to a wider spectrum of environmental objects. The mathematical core
of T3DRIS is rooted in optimizing the 3D deployment of the unit cells and
tuning circuits, aiming at maximizing the communication performance. Through
rigorous full-wave simulations and a comprehensive set of numerical analyses,
we validate the proposed approach and demonstrate its superior performance and
applicability over contemporary designs. This study—the first of its
kind—paves the way for a new direction in RIS research, emphasizing the
importance of a theoretical and mathematical perspective in tackling the
challenges of conformal RISs.

###### Index Terms:

Reconfigurable intelligent surfaces, smart radio environments, conformal
metasurfaces, mutual coupling, optimization.

##  I Introduction

###  I-A Background & Motivation

In the landscape of modern wireless communication systems, the advent of
reconfigurable intelligent surface (RIS) marks a significant milestone [1],
introducing the capability to dynamically manipulate the propagation of
electromagnetic (EM) waves. This property, together with its flexibility in
deployment and reconfiguration, low implementation cost and energy
consumption, positions RIS as a primary technology candidate for the upcoming
6th generation of wireless networks (6G) [2], opening up the possibility of
optimizing the communication channel to best serve different use cases that 6G
will support [3, 4]. Conventional implementations have predominantly explored
planar RIS architectures, which, while effective, encounter limitations in
adaptability and versatility, particularly when confronted with the complexity
of non-planar deployment scenarios [5].

Fig. 1 shows an interesting experiment: a wave from a signal source (on the
right) hits three different surfaces made of dipoles spaced evenly apart. When
the wave hits the flat surface (sub-figure ($a$)), it acts like a mirror,
bouncing the wave back in a tight, focused beam. However, when the surface is
curved (sub-figures ($b$) and ($c$)), an insightful behavior can be observed:
several beams are spreading out over a larger area. This prompts an
exploration into how the design of an EM surface shape can influence the way
it interacts with incoming waves, giving rise to the burgeoning field of
conformal or three dimensional (3D) RIS,111We use the terms 3D RIS and
conformal RIS interchangeably throughout the paper. sculpted to the unique
topography of their hosting environments [6].

Notwithstanding the progress in advanced printing techniques that facilitate
the fabrication of irregular electronic boards, the precision placement and
orientation of unit cells on 3D surfaces present formidable technical
challenges [7]. These include not only the physical constraints of inter-
element spacing, alignment, and orientation in 3D, but also the intricate
mathematical modeling necessary for effective dynamic control of such systems
[8, 9]. The fidelity of models based on infinite planar arrays falters when
confronted with the finite size and curvature of real-world structures,
leading to discrepancies between theoretical predictions and empirical
observations. Therefore, the optimization of the 3D RIS layout and its
corresponding configuration demands sophisticated algorithms tuned to complex
environmental dynamics and EM interactions.

Figure 1: Reflection properties of curved surfaces ($a$) planar, ($b$)
cylindrical and ($c$) super-cylindrical).

Recently, impedance-based models have been proposed to advance the
optimization and efficient design of RISs [10, 11]. It has been observed that
spacing the RIS elements at inter-distances smaller than $\lambda/2$ results
in a more compact layout, while still allowing for effective beamforming gains
despite the additional mutual coupling among the RIS elements, i.e.,
interactions between the elements of the antenna array. When dealing with non-
planar RISs, exploiting the mutual coupling at the design stage may benefit
the RIS configuration, leading to significant performance gains beyond state-
of-the-art (SotA) solutions [12, 13].

Through a comprehensive investigation of the fundamental properties and
performance characteristics of conformal RISs, we demonstrate that these 3D
structures can be optimized in their design phase and offer unprecedented
capabilities for manipulating EM waves. In particular, we propose a novel
solution, namely T3DRIS (Tailored 3D RIS, pronounced “tedris”), wherein the
arrangement of RIS elements is optimized while leveraging the mutual coupling
to maximize performance and keeping the RIS conformal to the support surface.
To the best of our knowledge, T3DRIS is the first attempt at optimizing the
layout of 3D RIS. Hence, it constitutes a breakthrough in the design of
intelligent surfaces, holding great promise for revolutionizing wireless
communication systems and advancing several applications for 6G.

###  I-B Contributions

The paper’s contributions can be elaborated upon as follows:

C1. Demonstration of 3D RIS Feasibility for Communication: The paper first and
foremost establishes the practicality of designing 3D RIS specifically for
communication purposes. This is a crucial step forward as it moves beyond the
traditional flat-panel designs, exploring the potential of RIS in three-
dimensional forms. This advancement is critical because 3D RIS can potentially
conform to various surfaces and shapes, offering more flexibility and
adaptability in various environments compared to their 2D counterparts.

C2. Mutual Coupling-Aware Model for 3D Conformal RIS: We delve into the
complex phenomenon of mutual coupling in 3D RIS with structures going beyond
planar geometry, a factor often overlooked in previous studies, which can
significantly impact performance. By developing a model that accounts for
these interactions, especially in RISs with non-uniform spacing between
elements, our paper provides a more accurate and realistic approach to future
RIS design. Moreover, it introduces significant research opportunities in
accurately modeling arbitrary impedance boundary conditions in conformal
structures.

C3. Framework for Joint Optimization of 3D Shape and Configuration: Building
on the mutual coupling-aware RIS channel model [10], this paper introduces a
novel framework for jointly optimizing both the physical shape of the RIS (3D
shape) and its operational configuration. This holistic approach ensures that
the design of the RIS is not only physically feasible but also operationally
effective, taking into account factors like signal strength, direction, and
interference. Apart from its immediate relevance to RIS-assisted wireless
networks, the introduced framework, which focuses on shaping an RIS by
determining the precise positioning of its elements, is directly applicable to
the development of three other cutting-edge technologies: smart skins [14],
fluid antennas [15], and movable antenna systems [16]. Each of these
technologies fundamentally relies on optimizing the positions of the radiating
components in antenna arrays. Consequently, this framework is expected to play
a significant role in shaping the design of future wireless systems.

C4. Validation through Simulation and Comparative Analysis: Lastly, we
validate our proposed approach using a commercial full-wave EM simulator,
which is crucial for testing the theoretical models in practical scenarios.
Additionally, we conduct an extensive numerical simulation study, comparing
our method against the most relevant solutions in the field. This
comprehensive validation process not only demonstrates the effectiveness of
our approach, but also highlights its significance and potential impact in the
context of existing solutions.

In summary, this paper represents a substantial step forward in the
development of conformal RIS for wireless applications. Its comprehensive
approach, from theoretical modeling to practical validation, addresses key
challenges in the field and opens up new possibilities for the design and
implementation of advanced smart surface technologies.

Notation. We denote matrices and vectors in bold text, while each of their
element is indicated in Roman font with a subscript.
$(\cdot)^{\operatorname{\mathrm{T}}}$, $(\cdot)^{\operatorname{\mathrm{H}}}$,
and $(\cdot)^{*}$ stand for vector or matrix transposition, Hermitian
transposition, and the complex conjugate operator, respectively. $\|\cdot\|$
denotes the L$2$-norm in the case of a vector and the spectral norm in the
case of a matrix, while the Frobenius norm is denoted by
$\|\cdot\|_{\mathrm{F}}$. The operators $\operatorname{\mathrm{tr}}(\cdot)$,
$\mathrm{diag}(\cdot)$, and $\left\lceil\cdot\right\rceil$ indicate the trace
of a matrix, a diagonal matrix, and the ceiling function, respectively. The
gradient operator with respect to the $N\times 1$ vector
$\operatorname{\mathbf{x}}$ is
$\nabla_{\operatorname{\mathbf{x}}}=[\frac{\partial}{\partial
x_{1}},\dots,\frac{\partial}{\partial x_{N}}]^{\operatorname{\mathrm{T}}}$.
$\operatorname{\mbox{$\mathbb{R}$}}$ and $\operatorname{\mbox{$\mathbb{C}$}}$
denote the sets of real and complex numbers, respectively. Lastly, $\cdot$ and
$\times$ represent dot and cross product, respectively.

Outline of the Paper. Section II gives a comprehensive analysis of the EM
implications arising from the conformal structure of RISs, together with their
effect on signal propagation. Section III details the system model. The
formulation of the joint shape and configuration optimization problem for 3D
RIS is outlined in Section IV, together with our proposed solution, T3DRIS.
The performance evaluation is provided in Section V, followed by an
exploration of related works in Section VI. Finally, Section VII closes the
paper with concluding remarks summarizing key insights and contributions.

##  II Fundamentals of Electromagnetic Wave Manipulation in 3D RIS

We outline four critical questions that address the issue of 3D RIS design
from an EM standpoint, which are detailed in the following.

_Q1. How do curved RIS affect signal propagation?_

The answer to this preliminary question has already been provided by [17]: The
finite size and non-zero curvature negatively affect the predictions that can
be made based on infinite planar array models. Also, there exists a
disagreement with measurements, based on spatial filters, such as frequency
selective surface (FSS) structures. Specifically, the limitations in modeling
accuracy are exacerbated in the case of RISs due to the presence of additional
re-radiation fields that depend on the non-uniform reflective state of the
unit cells. Floquet’s theory does not apply to wave propagation in curved
surfaces, as the periodicity is broken in the geometrical arrangement of the
elements, even if they are equispaced pairwise. Hence, full-wave simulation
methods based on unit cells terminated by periodic boundary conditions are
therefore ruled out.

The only viable option is therefore the mathematical modeling of the whole
surface. A comprehensive analysis of modeling methodologies for non-periodic
arrays [17] suggests that circuital methods, e.g., based on the method of
moments (MoM) and equivalent inductor-capacitor (LC) circuits (resonant
circuits), constitute a prominent solution in reducing the computation time
while providing accurate performance predictions. In the following sections,
we leverage a discrete model for RISs based on the canonical minimum
scattering approximation (CMSA) [10]. Specifically, transmit and receive
antenna arrays, as well as the RISs are modeled by leveraging the discrete
dipole approximation (DDA) [18], and dipole-to-dipole interactions (the EM
mutual coupling) are modeled by impedance matrices.

_Q2. Does mutual coupling play a role in the transition from planar to curved
radiating surfaces?_

In the context of non-periodic surfaces, the mutual coupling between dipoles
plays a major role in establishing optimized wave transformations at high
efficiency. The MoM provides a way to extend the impedance model beyond DDA as
it treats arbitrary surface current density and devises the radiated field
accounting for the influence of nearby unit cells, thus capturing mutual
coupling [19]. On the one hand, increasing the ability to control the
wavefronts requires the design of RISs with a large number of unit cells per
wavelength [11]. On the other hand, non-planar arrays are non-periodic and
have a finite size, whence they are also configured as geometrically non-
uniform [17].

Furthermore, RIS structures introduce non-uniformities on the local reflection
coefficient. Thus, the design of 3D RISs requires novel EM modeling methods
due to their non-uniform and non-planar structure. In general, further degrees
of freedom in wave control can be gained by intentionally introducing non-
uniform, and even disordered, arrangements of elements already within the
curved shell of a RIS. A notable example is constituted by hyper-uniform
disordered antenna arrays [20], where the side lobes can be effectively
suppressed thanks to the presence of non-uniformity. However, spurious lobes
can be generated by curved surfaces due to _creepy waves_ , which are those
surface waves that are converted into free-space radiation by propagating
along the surface curvature, and for which judicious physical modeling and
predictive design are imperative. Non-uniform element arrangements are hard to
model, but the mutual impedance model introduced in [10] for planar RIS
structures can inherently tackle non-uniformity.

_Q3. How does cell and surface curvature affect the radiated wave fields?_

Suppose a unit cell has a curvature in the order of a wavelength, the
impinging field induces an equivalent surface current density whose flow would
be perturbed by the non-planar structure. This is captured by an impedance
boundary condition that can be expressed in the appropriate local frame of
reference. For the sake of simplicity, let us consider a cylindrical structure
with coordinates $(\rho,\phi,z)$ representing the curvature radius, azimuth
angle, and elevation, respectively.222Note that the provided analysis can be
generalized for any surface geometry. Therefore, we can write the surface
impedance boundary condition as the following

$\textbf{E}(\rho,\phi,z)=\overline{\overline{\operatorname{\mathbf{Z}}}}_{s}(\phi,z)\cdot\left[\hat{\operatorname{\boldsymbol{\rho}}}\times\Delta\textbf{H}(\rho,\phi,z)\right],$
(1)

which relates the total (boundary) electric field $\textbf{E}(\rho,\phi,z)$ to
the magnetic fields at the left $\textbf{H}^{\textrm{l}}(\rho,\phi,z)$ and at
the right $\textbf{H}^{\textrm{r}}(\rho,\phi,z)$ of the equivalent impedance
sheet modelling the meta-surface, whence

$\Delta\textbf{H}(\rho,\phi,z)=\textbf{H}^{\textrm{l}}(\rho,\phi,z)-\textbf{H}^{\textrm{r}}(\rho,\phi,z),$
(2)

and where $\overline{\overline{\operatorname{\mathbf{Z}}}}_{s}$ is the non-
uniform (homogeneous) surface impedance that achieves a prescribed scattering
profile [21] and $\hat{\operatorname{\boldsymbol{\rho}}}$ radial versor
perpendicular to the surface. The local scattered field is thus related to the
induced current through
$\textbf{J}_{s}(\rho,\phi,z)=\hat{\operatorname{\boldsymbol{\rho}}}\,\times\,\Delta\textbf{H}(\rho,\phi,z)$.
Using this and decomposing the total field into incident and scattered
components, i.e.,
$\textbf{E}(\rho,\phi,z)=\textbf{E}^{\textrm{inc}}(\rho,\phi,z)+\textbf{E}^{\textrm{sca}}(\rho,\phi,z)$,
yields

$\textbf{E}^{\textrm{sca}}(\rho,\phi,z)=-\textbf{E}^{\textrm{inc}}(\rho,\phi,z)+\overline{\overline{\operatorname{\mathbf{Z}}}}_{s}(\phi,z)\cdot\textbf{J}_{s}(\rho,\phi,z).$
(3)

Consequently, the MoM can be used to devise the radiated fields of curved
surfaces through the reaction integrals [22, Eqs. (4) and (5)], which are used
to devise the elements of the finite size impedance matrix of the meshed
surface current [19] as follows

$\textbf{V}^{\textrm{MS}}=\textbf{Z}^{\textrm{MS}}\,\textbf{I}^{\textrm{MS}},$
(4)

where the surface port voltages $\textbf{V}^{\textrm{MS}}$ and currents
$\textbf{I}^{\textrm{MS}}$ are related by the impedance matrix
$\textbf{Z}^{\textrm{MS}}$ that is directly related to the (circuital) active
impedance matrix. Upon expansion of the surface current density onto the basis
functions $\textbf{u}_{n}^{\textrm{MS}}$ defined across a spatial mesh

$\textbf{J}_{s}=\sum_{n}I^{\textrm{MS}}_{n}\textbf{u}_{n}^{\textrm{MS}},$ (5)

where $I^{\textrm{MS}}_{n}$ is an element of $\textbf{I}^{\textrm{MS}}$ that
represents the associated current phasor, and the entries of
$\textbf{Z}^{\textrm{MS}}$ configure as a two-part impedance

$Z_{ij}^{\textrm{MS}}(n,m)=Z_{ij,G}^{\textrm{MS}}(n,m)+Z_{ij,Z_{s}}^{\textrm{MS}}(n,m),$
(6)

with individual parts given by an overlapping integral across the surface
space, given as

$Z_{ij,G}^{\textrm{MS}}(n,m)=\eta\,\int
u_{i}^{\textrm{MS}}(n)\,G_{ij}(n,m)\,u_{j}^{\textrm{MS}}(m)\,\,dS,$ (7)

where $G_{ij}(n,m)$ is the component $ij$ of the free-space dyadic Green
function propagating the fields from current element $n$ to current element
$m$, and with

$Z_{ij,Z_{s}}^{\textrm{MS}}(n,m)=\eta\,\int
u_{i}^{\textrm{MS}}(n)\,\overline{\overline{Z}}_{s,ij}(n,m)\,u_{j}^{\textrm{MS}}(m)\,\,dS,$
(8)

with free-space wave impedance $\eta=377\,\Omega$, and where the surface
impedance boundary condition couples current element $n$ with current element
$m$ through the component $ij$ of the impedance tensor in (3). A detailed
example for a circular modulated surface has been devised in [23] by using
entire-domain basis functions. In reflect-array and meta-surfaces, this
current follows the geometrical shape of the unit cell and is influenced by
the dielectric substrate, ground plane (in case the surface is reflective),
and embedded circuital element (if we prescribe unit cell reconfiguration).
This is captured by a modified expression of Green’s function that includes
multiple reflections between the unit cell interface and the ground plane
through the (lossy) dielectric substrate. While this is a computationally
intensive task, a fast algorithm for general - including curved - multilayered
media came to fruition over the last year [24]. Again, the MoM provides the
most general model to calculate the finite size impedance matrix from
arbitrary radiation structures. The case of planar RIS has been covered in
[25], where it is found that the end-to-end (E2E) impedance model holds also
in the case of real-life surface planar structures. As in [25], also in the
case of curved surfaces, the basis functions can be related to a passive
element, thus representing a mesh element of the surface boundary impedance,
or a port element, which represents the terminals loaded by the (RIS) control
circuitry. Therefore, in Section III we introduce the RIS self impedance
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}$
described as the port impedance $\operatorname{\mathbf{Z}}_{ij}^{\textrm{MS}}$
in (6).

_Q4. How can the impedance model be valid for real-life surface prototypes?_

It is worth pointing out that the E2E impedance model of RIS-aided wireless
channels can be used for arbitrary real-life structures [26, 27]. Besides the
ability of the MoM to tackle arbitrary RIS unit cells beyond the DDA described
above, there are further practical aspects that are captured by impedance-
based channel models just by extending the MoM formulation. For instance, the
presence of a supporting structure (walls, ground plane, etc.) can be
accounted for by using the multi-layer Green’s functions [24] in the reaction
integrals of the MoM in [25]. Inherently, the surface current densities
supported by arbitrary radiating structures usually entail characteristic
modes. In the case of RIS structures, two families of modes are included in
characteristic modes: i) structural modes, describing the radiation from the
base composite structure; ii) re-radiation modes, capturing the unit cell
reconfiguration and controlled by the tuning circuitry. We expect that the
radiation from curved RIS structures includes also a contribution from creepy
waves. We stress that the MoM predicts the full-wave radiation from
arbitrarily curved surfaces, but it remains a challenge to understand how to
control or even suppress the spurious radiation from creepy waves.

Hereafter, we assume that the RIS elements are electrically small and the
radius of curvature _does not bend the geometry of the unit cell_ and hence
does not change its resonance frequency. This is an important aspect for the
high radius of curvature of planar structure with patch-like unit cells, as
pointed out in [28], and will be investigated in future developments of the
present work. Furthermore, in order to make the treatment general and
accessible to communication-oriented studies, no substrate is assumed and a
dense arrangement of freestanding dipoles is considered, including planar
simultaneous transmitting and receiving RISs and omnisurfaces that have been
recently proposed in the literature to control the 3D space [29, 30].
Interestingly, the local planar structure approximation does not hold, as
there might not be a convenient system of reference where the array is
uniform. Consequently, approaches based on the vector-Legendre transformation
may become invalid [31]. Based on such considerations, in the next sections,
we introduce a new model, formulate an optimization problem, and describe a
corresponding optimization algorithm to design mutually coupled non-planar
RISs.

##  III System model

We consider an RIS-aided network, wherein a single-antenna base station (BS)
communicates with a single-antenna user equipment (UE) via an $N$-element RIS.
The center-point of the RIS is located at the origin, while the BS and the UE
are located in positions
$\operatorname{\mathbf{p}}_{\scriptscriptstyle\mathrm{BS}}\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times
1}$ and
$\operatorname{\mathbf{p}}_{\scriptscriptstyle\mathrm{UE}}\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times
1}$ with respect to the origin (the center-point of the RIS), respectively.
The selected scenario allows us to better focus on the 3D RIS layout
optimization problem, without losing generality. More complex scenarios with
multiple-antenna BS and multiple UEs under mobility conditions are left for
future work. We adopt the EM-consistent and mutual coupling aware channel
model in [10], according to which each RIS element is modeled as a thin wire
dipole driven by a tunable impedance. Specifically, the end-to-end channel is
given by

$\displaystyle
H=Y_{0}[Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}-\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{\operatorname{\mathrm{T}}}(\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}+\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}})^{-1}\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}]\in\operatorname{\mbox{$\mathbb{C}$}}$
(9)

where
$Y_{0}=Z_{\scriptscriptstyle\mathrm{L}}(Z_{\scriptscriptstyle\mathrm{L}}+Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{R}})^{-1}(Z_{\scriptscriptstyle\mathrm{G}}+Z_{\scriptscriptstyle\mathrm{T}\scriptscriptstyle\mathrm{T}})^{-1}\in\operatorname{\mbox{$\mathbb{C}$}}$
is a fixed constant, the subscripts $\mathrm{R},~{}\mathrm{T},~{}\mathrm{S}$
stand for the UE, the BS and the RIS, respectively, and $Z_{L}$ and $Z_{G}$
are the load impedances of the UE and the BS, respectively. Moreover,
$Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{R}}$ and
$Z_{\scriptscriptstyle\mathrm{T}\scriptscriptstyle\mathrm{T}}$ denote the self
impedances at the UE and BS, which are assumed to be known.
$Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}\in\operatorname{\mbox{$\mathbb{C}$}}$
represents the direct channel between the BS and the UE,
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
1}$ stands for the channel between the RIS and the UE, and
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
1}$ denotes the channel between the BS and the RIS. Lastly,
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
N}$ denotes the matrix of self and mutual impedances among the elements of the
RIS (derived as the port impedance
$\operatorname{\mathbf{Z}}_{ij}^{\textrm{MS}}$ in (6)), and
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
N}$ is the matrix of tunable RIS impedances.

The received signal at the UE is thus given by

$\displaystyle y=Hx+n\in\operatorname{\mbox{$\mathbb{C}$}}$ (10)

where $x$ is the transmitted symbol, $\operatorname{\mathbb{E}}[|x|^{2}]=P$
with $P$ being the total transmit power, and $n$ is the noise at the receiver
that is distributed as $\mathcal{CN}(0,\sigma_{n}^{2})$.

The RIS tunable impedance matrix is modeled as

$\displaystyle\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}=\mathrm{diag}\\{R_{0}+jb_{n}\\}_{n=1}^{N}$
(11)

where $R_{0}\geq 0$ is a constant that depends on the losses of the RIS
internal circuits, and $b_{n}\in\operatorname{\mbox{$\mathbb{R}$}}$ is the
tunable reactance at the $n$-th RIS element, i.e., the configuration of the
$n$-th element.

The entries of the matrix
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
N}$ can be computed by utilizing the analytical framework introduced in
[32].333Full-wave simulations can be used to compute the $Z$ parameters
required by our optimization approach, with no assumption on the element type
or scattering conditions [22]. Nonetheless, dipoles and CMSA allow us to
approach the problem analytically [10]. Likewise, the same framework can be
utilized to compute the vectors
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}}$
and
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}$,
as well as $Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}$.
Specifically, the self and mutual impedance between any pair of thin wire
dipoles $q$ and $p$ of length $\lambda/2$ and with position
$\operatorname{\mathbf{q}}_{q}\in\operatorname{\mbox{$\mathbb{R}$}}^{3}$ and
$\operatorname{\mathbf{q}}_{p}\in\operatorname{\mbox{$\mathbb{R}$}}^{3}$,
respectively, can be formulated as (see Appendix A for the details)

$\displaystyle Z_{qp}=\frac{\eta}{8\pi}$
$\displaystyle\sum_{s_{o}=\\{-1,1\\}}e^{j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{qp}}\big{[}T_{0}(\operatorname{\boldsymbol{\delta}}_{qp}-2h\operatorname{\mathbf{e}}_{3},s_{0})+$
$\displaystyle
T_{0}(\operatorname{\boldsymbol{\delta}}_{qp}+2h\operatorname{\mathbf{e}}_{3},s_{0})-2T_{0}(\operatorname{\boldsymbol{\delta}}_{qp},s_{0})\big{]},$
(12)

where $\eta=377~{}\Omega$ is the free-space wave impedance, $\lambda$ is the
signal wavelength, $h$ is the half-length of the dipoles, and
$\operatorname{\boldsymbol{\delta}}_{qp}=\operatorname{\mathbf{q}}_{q}-\operatorname{\mathbf{q}}_{p}$.
Moreover, we have defined the function

$\displaystyle
T_{0}(\operatorname{\boldsymbol{\zeta}},s_{0})=E_{1}\Big{(}jk\big{(}\|\operatorname{\boldsymbol{\zeta}}\|+s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\zeta}}\big{)}\Big{)},$
(13)

where $s_{0}\in\\{-1,1\\}$, and $E_{1}(z)$ is the exponential integral defined
as $E_{1}(z)=\int_{z}^{\infty}\frac{e^{-t}}{t}dt$.

The expression in (III) is not applicable when the dipoles $q$ and $p$ are co-
linear [32]. In that case, the expression in [32, Eq. (9)] can be utilized via
numerical integration.

##  IV Problem formulation

We aim at maximizing the signal-to-noise ratio (SNR) at the UE by finding the
optimal values of the tunable impedances driving the RIS elements and the
optimal locations of the RIS elements, i.e., the RIS _configuration_ and the
RIS _shape_ , respectively. To this end, we denote the vector of RIS
configurations as
$\operatorname{\mathbf{b}}=[b_{1},\ldots,b_{N}]^{\operatorname{\mathrm{T}}}\in\operatorname{\mbox{$\mathbb{R}$}}^{N\times
1}$ and the set of RIS positions as
$\operatorname{\mathbf{Q}}=[\operatorname{\mathbf{q}}_{1},\ldots,\operatorname{\mathbf{q}}_{N}]\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times
N}$, which represent the optimization variables, accordingly. Note that in
(9),
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}$
and
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}$
depend on $\operatorname{\mathbf{Q}}$, and they can be computed by using
(III).

Therefore, we focus on finding the optimal $\operatorname{\mathbf{b}}$ and
$\operatorname{\mathbf{Q}}$ by solving the following optimization problem

###### Problem 1

$\displaystyle\displaystyle\max_{\operatorname{\mathbf{b}},\,\operatorname{\mathbf{Q}}}\,\,$
$\displaystyle|Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}\\!-\\!\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{\operatorname{\mathrm{T}}}(\operatorname{\mathbf{Q}})(\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}(\operatorname{\mathbf{Q}})\\!+\\!\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}(\operatorname{\mathbf{b}}))^{-1}\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}(\operatorname{\mathbf{Q}})|^{2}$
(14)
$\displaystyle\mathrm{s.t.}\quad\operatorname{\mathbf{Q}}\in\mathcal{Q}^{N},\quad\operatorname{\mathbf{b}}\in\mathcal{B}^{N},$
(15)

where $\mathcal{Q}$ and $\mathcal{B}$ represent the feasible sets of the RIS
elements positions and tunable reactances, respectively. Note also that we
assume $\mathcal{Q}$ to be a convex set (see Sections IV-B–IV-E). It is worth
specifying that the selection of the feasible set $\mathcal{Q}$ is dictated by
the deployment constraints of the surface. In other words, it ensures the
practicality of the obtained solution. However, constrained sets might lead to
suboptimal element arrangements from a pure communication standpoint.
Nevertheless, the selection of an unconstrained set yields an optimal RIS
layout for communication, which can be a valuable input for defining new
strategies for the design of conformal RISs.

Problem 1 is particularly challenging to solve given its non-convex
formulation, the inversion of the matrix
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}(\operatorname{\mathbf{Q}})+\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}(\operatorname{\mathbf{b}})$,
and the fact that, except for
$Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}$, the impedances
to be optimized depend on the locations of the RIS elements
$\operatorname{\mathbf{Q}}$. Moreover, we emphasize that SotA research works
focused only on the optimization of the RIS tunable reactances [33, 34, 35,
36], while no research work has tackled the joint optimization of the RIS
configuration and shape.

###  IV-A T3DRIS Algorithm

Hereafter, we introduce the proposed iterative algorithm to tackle Problem 1,
which is referred to as _T3DRIS_.

Let us denote the position of the RIS elements at the $i$-th iteration of
T3DRIS as $\operatorname{\mathbf{Q}}^{(i)}$, and the corresponding value of
the mutual impedances between the transmitter and the RIS, the self and mutual
impedances between the RIS elements, and between the latter and the receiver
as
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}=\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}(\operatorname{\mathbf{Q}}^{(i)})$,
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(i)}=\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}(\operatorname{\mathbf{Q}}^{(i)})$
and
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}=\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}(\operatorname{\mathbf{Q}}^{(i)})$.
Note that we depart from an initial condition in which the RIS elements are
located at the initial position $\operatorname{\mathbf{Q}}^{(0)}$, and the
corresponding self and mutual impedances are given by
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(0)}$,
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(0)}$
and
$\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(0)}$.444The
initial shape $\operatorname{\mathbf{Q}}^{(0)}$ can be arbitrarily chosen
without affecting the overall performance of the proposed algorithm, as shown
in Section V.

Therefore, we focus on optimizing the RIS configuration, which is specified by
(11). Specifically, at each iteration $i$ of the proposed algorithm, we
formulate the following optimization problem

###### Problem 2

$\displaystyle\max_{\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}}$
$\displaystyle|Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}-\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{\operatorname{\mathrm{T}}(i)}(\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(i)}+\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}})^{-1}\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}|^{2}$
(16)
$\displaystyle\mathrm{s.t.}\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}\\!=\\!\mathrm{diag}(R_{0}+jb_{k}),b_{k}\in\mathcal{B},\,k\\!=\\!1,\ldots,N,$
(17)

which can be solved by using [33, 34].

Once the optimized values of
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}$ are obtained, we
turn to the optimization of the positions of the RIS elements by formulating
the following optimization problem

###### Problem 3

$\displaystyle\max_{\operatorname{\mathbf{Q}}}~{}$
$\displaystyle|Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}-\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{\operatorname{\mathrm{T}}(i)}(\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(i)}+\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}})^{-1}\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}|^{2}$
$\displaystyle\mathrm{s.t.}~{}\operatorname{\mathbf{Q}}\in\mathcal{Q}^{N},$
(18)

which we tackle by using the projected gradient ascent method, as

$\displaystyle\operatorname{\mathbf{Q}}^{(i+1)}=P_{\mathcal{Q}}\Big{(}\operatorname{\mathbf{Q}}^{(i)}+\alpha_{i}\nabla
f(\operatorname{\mathbf{Q}}^{(i)})\Big{)},$ (19)

where $P_{\mathcal{Q}}$ is the projection operator onto the feasible set
$\mathcal{Q}$, $f(\operatorname{\mathbf{Q}}^{(i)})$ is the objective function
in Problem 3, which is expressed as a function of
$\operatorname{\mathbf{Q}}^{(i)}$, and
$\alpha_{i}>0\in\operatorname{\mbox{$\mathbb{R}$}}$ is the step size. We
define $P_{\mathcal{Q}}(\operatorname{\mathbf{Q}})$ as

$\displaystyle
P_{\mathcal{Q}}(\operatorname{\mathbf{Q}}_{0})=\displaystyle\arg\min_{\operatorname{\mathbf{Q}}\in\mathcal{Q}^{N}}\frac{1}{2}\|\operatorname{\mathbf{Q}}-\operatorname{\mathbf{Q}}_{0}\|^{2}_{\mathrm{F}}.$
(20)

Hence, the updating rule corresponds to solving the following problem

$\displaystyle\displaystyle\min_{\operatorname{\mathbf{Q}}^{(i+1)}\in\mathcal{Q}^{N}}\|\operatorname{\mathbf{Q}}^{(i+1)}-\operatorname{\mathbf{Q}}^{(i)}-\alpha_{i}\nabla
f(\operatorname{\mathbf{Q}}^{(i)})\|_{\mathrm{F}}^{2}.$ (21)

The gradient of $f(\operatorname{\mathbf{Q}}^{(i)})$ is challenging to compute
due to the expression of each $Z_{qp}$ in (III) and the fact that each
$Z_{qp}$ appears inside the inverse of
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}$.
Hence, we employ the Neumann series approximation to tackle this problem [35].
Specifically, at each iteration of the algorithm, we optimize only a small
variation of the position of the RIS elements as
$\operatorname{\mathbf{Q}}^{(i+1)}=\operatorname{\mathbf{Q}}^{(i)}+\Delta\operatorname{\mathbf{Q}}$,
where
$\Delta\operatorname{\mathbf{Q}}\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times
N}$ denotes the small variations. More precisely, we update the location of
one RIS element at a time. Hence, for each $k=1,\ldots,N$, we have
$\operatorname{\mathbf{q}}_{k}^{(i+1)}=\operatorname{\mathbf{q}}_{k}^{(i)}+\Delta\operatorname{\mathbf{q}}_{k}$.
As a result, at the $i$-th iteration of the algorithm, we obtain
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(i+1)}=\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(i)}+\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)},~{}\forall
k$, where the matrix
$\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
N}$ is non-zero only in $k$-th row and column, which correspond to the $k$-th
RIS element being updated. Moreover, due to
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}$
being symmetric by definition, also
$\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}=(\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)})^{{\operatorname{\mathrm{T}}}}$
applies. Similarly, we have

$\displaystyle{\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i+1)}=\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}+\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}}\quad\text{and}\quad\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i+1)}=\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}+\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)},$
(22)

where only the $k$-th entry of
$\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
1}$ and
$\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
1}$ is non-zero.

Let
$\operatorname{\mathbf{G}}^{(i)}=(\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(i)}+\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}})^{-1}$,
and using the Neumann series approximation, the inverse matrix in Problem 3 is
approximated as

$\displaystyle(\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(i)}+\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}})^{-1}\approx\operatorname{\mathbf{G}}^{(i)}-\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{G}}^{(i)}=\tilde{\operatorname{\mathbf{\Gamma}}}^{(i)},$
(23)

provided that
$\|\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}^{(k)}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\|\ll
1$ [35, 36]. Furthermore, note that
$\operatorname{\mathbf{G}}^{(i)}=(\operatorname{\mathbf{G}}^{(i)})^{\operatorname{\mathrm{T}}}$.

For ease of writing, we refer to the objective function in Problem 3 as

$\displaystyle f(\Delta\operatorname{\mathbf{q}}_{k})$
$\displaystyle=|h(\Delta\operatorname{\mathbf{q}}_{k})|^{2}$
$\displaystyle=|Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}-\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{\operatorname{\mathrm{T}}(i+1)}(\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(i+1)}+\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}})^{-1}\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i+1)}|^{2}.$
(24)

From (23), we can approximate the function
$h(\Delta\operatorname{\mathbf{q}}_{k})$ as

$\displaystyle h(\Delta\operatorname{\mathbf{q}}_{k})\approx\,$ $\displaystyle
Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}-(\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}+\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)})^{\operatorname{\mathrm{T}}}\tilde{\operatorname{\mathbf{\Gamma}}}^{(i)}(\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}+\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}).$
(25)

Moreover, we introduce the variables

$\displaystyle
Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}$
$\displaystyle=Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{T}}-\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{\operatorname{\mathrm{T}}(i)}\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}\in\operatorname{\mbox{$\mathbb{C}$}},$
(26)
$\displaystyle\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}$
$\displaystyle=\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
1},$ (27)
$\displaystyle\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}$
$\displaystyle=\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{z}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
1},$ (28)

and rewrite (25) in a more convenient form as

$\displaystyle h$ $\displaystyle(\Delta\operatorname{\mathbf{q}}_{k})\approx
Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}-(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}+(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}$
$\displaystyle+(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}-(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}$
$\displaystyle-(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{G}}^{(i)}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}+(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{G}}^{(i)}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}$
$\displaystyle+(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{G}}^{(i)}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}.$
(29)

Also, the matrix
$\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}$
can be expressed as

$\displaystyle\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}=\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{e}}_{k}^{\operatorname{\mathrm{T}}}+\operatorname{\mathbf{e}}_{k}(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)})^{\operatorname{\mathrm{T}}}-\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{e}}_{k}\operatorname{\mathbf{e}}_{k}^{\operatorname{\mathrm{T}}},$
(30)

where
$\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
1}$ is the only non-zero column of
$\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}$,
$\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\in\operatorname{\mbox{$\mathbb{C}$}}$
is the element in position $(k,k)$ of
$\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}$,
or equivalently the only non-zero element of
$\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}$,
and $\operatorname{\mathbf{e}}_{k}$ is the $k$-th column of the identity
matrix of size $N$. Hence, we obtain

$\displaystyle{(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)})^{\operatorname{\mathrm{T}}}}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}$
$\displaystyle={(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}$
(31)
$\displaystyle(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{G}}^{(i)}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}$
$\displaystyle=\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}(\operatorname{\mathbf{g}}_{k}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{g}}_{k}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}$
$\displaystyle=\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}(\bar{\operatorname{\mathbf{g}}}_{k}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)},$
(32)

where $\operatorname{\mathbf{g}}_{k}^{(i)}$ is the $k$-th column of
$\operatorname{\mathbf{G}}^{(i)}$, and we have defined
$\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}=g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}+g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}-g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\operatorname{\mathbf{e}}_{k}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
1}$ and
$\bar{\operatorname{\mathbf{g}}}_{k}^{(i)}=g_{k,k}^{(i)}\operatorname{\mathbf{g}}_{k}^{(i)}+g_{k,k}^{(i)}\operatorname{\mathbf{g}}_{k}^{(i)}-g_{k,k}^{(i)}g_{k,k}^{(i)}\operatorname{\mathbf{e}}_{k}$.
Similarly, we obtain

$\displaystyle\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}$
$\displaystyle={(\widetilde{\operatorname{\mathbf{G}}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)},$
(33)
$\displaystyle{(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\mathbf{\Delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\operatorname{\mathbf{G}}^{(i)}}$
$\displaystyle={(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)})^{\operatorname{\mathrm{T}}}\widetilde{\operatorname{\mathbf{G}}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}},$
(34)

where the matrices
$\widetilde{\operatorname{\mathbf{G}}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}$
and
$\widetilde{\operatorname{\mathbf{G}}}_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}}^{(i)}$
are constructed as

$\displaystyle\widetilde{\operatorname{\mathbf{G}}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}$
$\displaystyle=\begin{bmatrix}g_{1,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}{(\operatorname{\mathbf{g}}_{1}^{(i)})^{\operatorname{\mathrm{T}}}}+g_{11}^{(i)}(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}-g_{1,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}g_{11}^{(i)}\operatorname{\mathbf{e}}_{k}^{\operatorname{\mathrm{T}}}\\\
\vdots\\\
g_{N,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}{(\operatorname{\mathbf{g}}_{N}^{(i)})^{\operatorname{\mathrm{T}}}}+g_{NN}^{(i)}(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}-g_{N,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}g_{NN}^{(i)}\operatorname{\mathbf{e}}_{k}^{\operatorname{\mathrm{T}}}\end{bmatrix}^{\operatorname{\mathrm{T}}}$
$\displaystyle=\begin{bmatrix}\tilde{\operatorname{\mathbf{g}}}_{1,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}&\ldots&\tilde{\operatorname{\mathbf{g}}}_{N,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}\end{bmatrix}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
N},$ (35)
$\displaystyle\widetilde{\operatorname{\mathbf{G}}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}$
$\displaystyle=\begin{bmatrix}g_{1,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}{(\operatorname{\mathbf{g}}_{1}^{(i)})^{\operatorname{\mathrm{T}}}}+g_{11}^{(i)}{(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)})^{\operatorname{\mathrm{T}}}}-g_{1,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}g_{11}^{(i)}\operatorname{\mathbf{e}}_{k}^{\operatorname{\mathrm{T}}}\\\
\vdots\\\
g_{N,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}{(\operatorname{\mathbf{g}}_{N}^{(i)})^{\operatorname{\mathrm{T}}}}+g_{NN}^{(i)}{(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)})^{\operatorname{\mathrm{T}}}}-g_{N,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}g_{NN}^{(i)}\operatorname{\mathbf{e}}_{k}^{\operatorname{\mathrm{T}}}\end{bmatrix}^{\operatorname{\mathrm{T}}}$
$\displaystyle=\begin{bmatrix}\tilde{\operatorname{\mathbf{g}}}_{1,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}&\ldots&\tilde{\operatorname{\mathbf{g}}}_{N,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\end{bmatrix}\in\operatorname{\mbox{$\mathbb{C}$}}^{N\times
N}.$ (36)

Since the locations of the RIS elements are updated one at a time, we can
further simplify the approximate objective function of Problem 3 in (25) by
expressing it as a function of the columns of
$\widetilde{\operatorname{\mathbf{G}}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}$
and
$\widetilde{\operatorname{\mathbf{G}}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}$
as

$\displaystyle h$ $\displaystyle(\Delta\operatorname{\mathbf{q}}_{k})\approx
Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}-\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}+(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}$
$\displaystyle+\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}(\tilde{\operatorname{\mathbf{g}}}_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}-g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}-G_{k,k}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}$
$\displaystyle+(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)})^{\operatorname{\mathrm{T}}}\tilde{\operatorname{\mathbf{g}}}_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}+\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}(\bar{\operatorname{\mathbf{g}}}_{k}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}.$
(37)

Hence, after applying the Neumann series we can rewrite Problem 3 as the
following

###### Problem 4

$\displaystyle\max_{\Delta\operatorname{\mathbf{q}}_{k}}~{}$
$\displaystyle\big{\rvert}Z_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}-\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}+(\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}+\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}(\tilde{\operatorname{\mathbf{g}}}_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}$
$\displaystyle-
g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}-G_{k,k}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}$
$\displaystyle+(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)})^{\operatorname{\mathrm{T}}}\tilde{\operatorname{\mathbf{g}}}_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}+\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}(\bar{\operatorname{\mathbf{g}}}_{k}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}\big{\rvert}^{2},\,\forall
k$
$\displaystyle\mathrm{s.t.}\quad\Delta\operatorname{\mathbf{q}}_{k}\,:\,\|\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}^{(k)}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\|\ll
1,$ (38)

where the constraint is needed to make the Neumann series approximation
accurate by design. We are now in the position of computing the gradient of
the approximate objective function of Problem 3 in (19) with respect to the
(small) update on the position of the $k$-th RIS element. We remark that
$\operatorname{\mathbf{q}}_{k}\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times
1}$ is the $k$-th column of the matrix
$\operatorname{\mathbf{Q}}\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times N}$,
which contains the position of the RIS elements and
$\Delta\operatorname{\mathbf{q}}_{k}$ represents its update.

The gradient of $f(\Delta\operatorname{\mathbf{q}}_{k})$ is given by the
following

$\displaystyle\nabla f$
$\displaystyle(\Delta\operatorname{\mathbf{q}}_{k})=\nabla|h(\Delta\operatorname{\mathbf{q}}_{k})|^{2}=\nabla\big{(}h(\Delta\operatorname{\mathbf{q}}_{k})h^{*}(\Delta\operatorname{\mathbf{q}}_{k})\big{)}=$
$\displaystyle=\big{(}\nabla
h(\Delta\operatorname{\mathbf{q}}_{k})\big{)}h^{*}(\Delta\operatorname{\mathbf{q}}_{k})+h(\Delta\operatorname{\mathbf{q}}_{k})\big{(}\nabla^{*}h(\Delta\operatorname{\mathbf{q}}_{k})\big{)}.$
(39)

From the expression in (37), the gradient of
$h(\Delta\operatorname{\mathbf{q}}_{k})$ is approximated by

$\displaystyle\nabla h(\Delta\operatorname{\mathbf{q}}_{k})\approx$
$\displaystyle\nabla\delta_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}\bigg{[}(\tilde{\operatorname{\mathbf{g}}}_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\\!-\\!g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}\\!-\\!G_{k,k}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}\\!+\\!(\bar{\operatorname{\mathbf{g}}}_{k}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}\bigg{]}$
$\displaystyle+\nabla\delta_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}\bigg{[}(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)})^{\operatorname{\mathrm{T}}}\tilde{\operatorname{\mathbf{g}}}_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\\!+\\!(\bar{\operatorname{\mathbf{g}}}_{k}^{(i)})^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}\delta_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}\\!-\\!g_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\\!-\\!G_{k,k}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}\bigg{]}$
$\displaystyle+\nabla(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)})^{\operatorname{\mathrm{T}}}\bigg{[}\operatorname{\mathbf{g}}_{\scriptscriptstyle\mathrm{R}\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}\\!+\\!\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}\tilde{\operatorname{\mathbf{g}}}_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(i)}\\!\\!+\\!\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}\delta_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}\bar{\operatorname{\mathbf{g}}}_{k}^{(i)}\\!\\!+\\!\tilde{\operatorname{\mathbf{g}}}_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(i)}\delta_{k,\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}\bigg{]},$
(40)

while the gradient of $h^{*}(\Delta\operatorname{\mathbf{q}}_{k})$ is obtained
from (40) as follows

$\displaystyle\nabla h^{*}(\Delta\operatorname{\mathbf{q}}_{k})=\big{(}\nabla
h(\Delta\operatorname{\mathbf{q}}_{k})\big{)}^{*}.$ (41)

In addition, we note that the gradient of
$\delta_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}$,
$\delta_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}$ and
$\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)}$,
as well as their complex conjugates, boil down to the gradient of the single
impedances $Z_{qp}$. In particular, we have that

$\displaystyle\nabla\delta_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{R}}^{(k)}=\nabla
Z_{k\scriptscriptstyle\mathrm{R}},\quad\nabla\delta_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{T}}^{(k)}=\nabla
Z_{k\scriptscriptstyle\mathrm{T}},\quad\text{and}$
$\displaystyle\nabla(\operatorname{\boldsymbol{\delta}}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}^{(k)})^{\operatorname{\mathrm{T}}}=\begin{bmatrix}\nabla
Z_{1,k}^{(k)},\ldots,\nabla Z_{N,k}^{(k)}\end{bmatrix},$ (42)

where $k$ is the index of the $k$-th RIS element. Finally, the gradient of
$Z_{kl}$ for any two antennas $l$ and $k$, with respect to the antenna $k$ is
formulated for $\lambda/2$-length dipoles as

$\displaystyle\nabla Z_{kl}=\frac{\eta}{8\pi}\\!\\!$
$\displaystyle\sum_{s_{o}=\\{-1,1\\}}\\!\\!\Big{(}jks_{0}\operatorname{\mathbf{e}}_{3}e^{jks_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}(\operatorname{\mathbf{q}}_{k}-\operatorname{\mathbf{q}}_{l}))}g(\operatorname{\mathbf{q}}_{k}-\operatorname{\mathbf{q}}_{l},s_{0})$
$\displaystyle+e^{jks_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}(\operatorname{\mathbf{q}}_{k}-\operatorname{\mathbf{q}}_{l})}\nabla
g(\operatorname{\mathbf{q}}_{k}-\operatorname{\mathbf{q}}_{j},s_{0})\Big{)},$
(43)

where the function
$g(\operatorname{\mathbf{q}}_{k}-\operatorname{\mathbf{q}}_{l},s_{0})$ and its
gradient $\nabla
g(\operatorname{\mathbf{q}}_{k}-\operatorname{\mathbf{q}}_{l},s_{0})$ are
available in (71) and (73) of Appendix B, respectively. Moreover, note that
$\nabla Z_{kl}=-\nabla Z_{kl}$. Finally, the position of the $k$-th RIS
element is updated by solving the problem

###### Problem 5

$\displaystyle\displaystyle\min_{\operatorname{\mathbf{q}}_{k}^{(i+1)}}$
$\displaystyle\|\operatorname{\mathbf{q}}_{k}^{(i+1)}-\operatorname{\mathbf{q}}_{k}^{(i)}-\alpha_{i}\nabla
f(\Delta\operatorname{\mathbf{q}}_{k})\|^{2}$ (44)
$\displaystyle\mathrm{s.t.}\quad\operatorname{\mathbf{q}}_{k}^{(i+1)}\in\mathcal{Q}$
(45)
$\displaystyle\quad\quad\,\,\operatorname{\mathbf{q}}_{k}^{(i+1)}\,:\,\|\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}^{(k)}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\|\ll
1,$ (46)

where $\nabla f(\Delta\operatorname{\mathbf{q}}_{k})$ is obtained by
substituting (43) into (40) and (41), which are in turn substituted into (39).

In Problem 5, the second constraint is again needed to ensure convergence of
the Neumann series. In particular, this latter constraint is imposed by
operating a line search on the maximum value of $\alpha_{i}$, which ensures
that the update of the $k$-th position of the RIS element is small enough to
reach convergence of the Neumann series. Moreover, in order to guarantee
convergence of the overall algorithm, we also impose a monotonic increase of
the objective function in Problem 4, which is achieved by ensuring that
$f(\operatorname{\mathbf{q}}_{k}^{i+1})\geq
f(\operatorname{\mathbf{q}}_{k}^{i})$.

In the following, we describe several possible strategies to define the
feasible set $\mathcal{Q}$ for the positions of the RIS elements.

###  IV-B 3D sphere

As a first example, we consider a 3D sphere of given radius $R$. Hence,
$\mathcal{Q}$ is modeled as

$\displaystyle\mathcal{Q}=\\{\operatorname{\mathbf{q}}\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times
1}\,:\,\|\operatorname{\mathbf{q}}\|\leq R\\},$ (47)

where the elements of
$\operatorname{\mathbf{q}}_{k}=\begin{bmatrix}x_{k}\\!\\!&\\!\\!y_{k}\\!\\!&\\!\\!z_{k}\end{bmatrix}^{\operatorname{\mathrm{T}}}$
represent the Cartesian coordinates of the $k$-th RIS element.

Hence, in this case, the projection operation is achieved as

$\displaystyle
P_{\mathcal{Q}}(\operatorname{\mathbf{q}})=\begin{cases}\operatorname{\mathbf{q}}&\text{if
}\|\operatorname{\mathbf{q}}\|\leq R\\\
\frac{\operatorname{\mathbf{q}}}{\|\operatorname{\mathbf{q}}\|}R&\text{otherwise}\end{cases}.$
(48)

###  IV-C Planar array

In this case, the projection operation is achieved by fixing either the first
or second component of each $\operatorname{\mathbf{q}}_{k}$ in Cartesian
coordinates to a specified value. Indeed, the former case corresponds to the
RIS laying on the $y-z$ plane, while the latter is associated with the RIS
laying on the $x-z$ plane. Moreover, the remaining coordinates are given
within specified limits (e.g., in the case of an RIS on the $y-z$ plane) as

$\displaystyle y_{\mathrm{min}}\leq y_{k}\leq y_{\mathrm{max}},\quad
z_{\mathrm{min}}\leq z_{k}\leq z_{\mathrm{max}}\quad\forall k.$ (49)

Let $x$ denote the fixed coordinate of the planar array, therefore, we modify
Problem 5 as

###### Problem 6

$\displaystyle\begin{array}[]{cl}\displaystyle\min_{\operatorname{\mathbf{q}}_{k}^{(i+1)}}&\|\operatorname{\mathbf{q}}_{k}^{(i+1)}-\operatorname{\mathbf{q}}_{k}^{(i)}+\alpha_{i}\nabla
f(\Delta\operatorname{\mathbf{q}}_{k})\|^{2}\\\
\mathrm{s.t.}&x_{k}^{(i+1)}=x\quad\forall k\\\ &y_{\mathrm{min}}\leq
y_{k}^{(i+1)}\leq y_{\mathrm{max}}\quad\forall k\\\ &z_{\mathrm{min}}\leq
z_{k}^{(i+1)}\leq z_{\mathrm{max}}\quad\forall k\\\
&\operatorname{\mathbf{q}}_{k}^{(i+1)}\,:\,\|\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}^{(k)}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\|\ll
1.\end{array}$ (55)

###  IV-D Spherical array

In this case, we express each vector $\operatorname{\mathbf{q}}_{k}$ in a
spherical coordinate system as

$\displaystyle\operatorname{\mathbf{q}}_{k}=\begin{bmatrix}\rho_{k}&\theta_{k}&\phi_{k}\end{bmatrix}^{\operatorname{\mathrm{T}}}\forall
k,$ (56)

where $\rho_{k}$ represents the distance from the RIS reference point, while
$\theta_{k}$ and $\phi_{k}$ represent the azimuth and elevation angles,
respectively. Moreover, we have that
$\rho_{k}=\sqrt{x_{k}^{2}+y_{k}^{2}+z_{k}^{2}}$,
$\theta_{k}=\arctan\frac{y_{k}}{x_{k}}$, and
$\phi_{k}=\arccos\frac{z_{k}}{\rho_{k}}$. Accordingly, the gradient of the
objective function needs to be re-scaled as

$\displaystyle\nabla
f(\Delta\operatorname{\mathbf{q}}_{k})=\begin{bmatrix}\frac{d\,f}{d\,\rho_{k}}&\frac{1}{\rho_{k}\sin\phi_{k}}\frac{d\,f}{d\,\theta_{k}}&\frac{1}{\rho_{k}}\frac{d\,f}{d\,\phi_{k}}\end{bmatrix}^{\operatorname{\mathrm{T}}}.$
(57)

Let $R$ denote the chosen radius, we thus modify Problem 5 as follows

###### Problem 7

$\displaystyle\begin{array}[]{cl}\displaystyle\min_{\operatorname{\mathbf{q}}_{k}^{(i+1)}}&\|\operatorname{\mathbf{q}}_{k}^{(i+1)}-\operatorname{\mathbf{q}}_{k}^{(i)}+\alpha_{i}\nabla
f(\Delta\operatorname{\mathbf{q}}_{k})\|^{2}\\\
\mathrm{s.t.}&\rho_{k}^{(i+1)}=R\quad\forall k\\\
&\theta_{\mathrm{min}}\leq\theta_{k}^{(i+1)}\leq\theta_{\mathrm{max}}\\\
&\phi_{\mathrm{min}}\leq\phi_{k}^{(i+1)}\leq\phi_{\mathrm{max}}\\\
&\operatorname{\mathbf{q}}_{k}^{(i+1)}\,:\,\|\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}^{(k)}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\|\ll
1.\end{array}$ (63)

###  IV-E Cylindrical array

In this case, we modify the formulation in (56) as
$\rho_{k}=\sqrt{x_{k}^{2}+y_{k}^{2}}$,
$\theta_{k}=\arctan\frac{y_{k}}{x_{k}}$, and $\phi_{k}=z_{k}$, where $z_{k}$
represents the height of the location of the $k$-th RIS element. Similarly, as
in the spherical array case, the gradient of the objective function needs to
be re-scaled as

$\displaystyle\nabla
f(\Delta\operatorname{\mathbf{q}}_{k})=\begin{bmatrix}\frac{d\,f}{d\,\rho_{k}}&\frac{1}{\rho_{k}}\frac{d\,f}{d\,\theta_{k}}&\frac{d\,f}{d\,\phi_{k}}\end{bmatrix}^{\operatorname{\mathrm{T}}}.$
(64)

The complete formalization of the proposed optimization framework is provided
in Algorithm 1, while Fig. 2 provides a graphical view of its corresponding
block diagram. We remark here that appropriate modelling of the set of
feasible RIS element positions $\mathcal{Q}$ opens the possibility to employ
the proposed T3DRIS framework not only during the RIS design stage, but also
_a posteriori_ , i.e., for the optimization of existing antenna systems and
their inherently flexible geometry (e.g., smart skins [14], fluid antennas
[15], and movable antenna systems [16]).

Figure 2: Graphical description of the proposed T3DRIS approach.

###  IV-F Convergence analysis

Our proposed approach consists of an alternating optimization between the two
optimization variables, namely the RIS impedance matrix
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}$ and the 3D
position of each of its elements $\operatorname{\mathbf{Q}}$, which are
disjoint thanks to the structure of Problem 1. Hence, we can decouple our
optimization by alternating between the solution of Problem 2 and Problem 3.
The former is optimized via the approach in [34], which is shown to converge
to a stationary point of Problem 2 with a monotonically increasing behavior of
its objective function. Whereas, for a given
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}$ the optimization
of $\operatorname{\mathbf{Q}}$ is given in Problem 3. Firstly, our proposed
approach proceeds by applying the Neumann series approximation for each $k$-th
RIS element, thus obtaining Problem 4. Secondly, we apply a projected gradient
ascent [37], where the projection onto the feasible set $\mathcal{Q}$ is
detailed in Problem 5. At each iteration of the algorithm, we carefully select
the step size $\alpha_{i}$ such that the following conditions are met: i)
convergence of the Neumann series, which is guaranteed by the condition
$\|\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}^{(k)}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\|\ll
1$, ii) monotonic increase of the objective function in Problem 4, which is
achieved by ensuring that $f(\operatorname{\mathbf{q}}_{k}^{i+1})\geq
f(\operatorname{\mathbf{q}}_{k}^{i})$, with equality implying that
$\operatorname{\mathbf{q}}_{k}^{i}$ is a stationary point of Problem 4, and
iii) feasibility of the solution, which is ensured by considering only an
updated position of each $k$-th RIS element belonging to the feasible set
$\mathcal{Q}$. The aforementioned choice of the step size $\alpha_{i}$
together with the assumption on the convexity of $\mathcal{Q}$ guarantees
convergence to a stationary point of Problem 4 [37].

Therefore, starting from the solution of Problem 2 with its optimized
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}$, the projected
gradient ascent can never decrease the objective function thanks to the
aforementioned choice of $\alpha_{i}$. Lastly, we note that the objective
function of Problem 1 is bounded from above thanks to the physical nature of
the problem, i.e., the passive nature of the RIS, and since the location of
the RIS elements is bounded. Finally, we thus conclude that every limit point
of the sequence of solutions in
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}$ and $\mathcal{Q}$
obtained with the proposed approach is a stationary point of the original
Problem 1.

Algorithm 1 T3DRIS

1:Initialize $\epsilon$, $\mathrm{I}_{\mathrm{max}}$, $\mathrm{SNR}^{(0)}$,
$i=0$, and $\operatorname{\mathbf{Q}}^{(0)}$

2:Given $\operatorname{\mathbf{Q}}^{(0)}$, calculate all $Z_{qp}$ by using
(III)

3:while
$\big{(}\mathrm{SNR}^{(i)}-\mathrm{SNR}^{(i+1)}\big{)}/\mathrm{SNR}^{(i)}>\epsilon$
$\wedge$ $i<\mathrm{I}_{\mathrm{max}}$ do

4: For a fixed $\operatorname{\mathbf{Q}}^{(i)}$, optimize
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}$ using [34]

5: for $k=1,\ldots,N$ do

6: Calculate $\nabla f(\Delta\operatorname{\mathbf{q}}_{k})$ by substituting
(42) into (40) and (41), which are in turn substituted in (39)

7: Find the maximum value of $\alpha_{i}$ that ensures
$\|\operatorname{\mathbf{G}}^{(i)}\operatorname{\mathbf{\Delta}}^{(k)}_{\scriptscriptstyle\mathrm{S}\scriptscriptstyle\mathrm{S}}\|\ll
1$ and $f(\operatorname{\mathbf{q}}_{k}^{i+1})\geq
f(\operatorname{\mathbf{q}}_{k}^{i})$ via bisection

8: Set $\operatorname{\mathbf{q}}_{k}^{i+1}$ as the solution of Problem 5

9: Recalculate all $Z_{qp}$ by using (III)

10: end for

11: $i=i+1$

12:end while

13:Output: $\operatorname{\mathbf{Q}}$ and
$\operatorname{\mathbf{Z}}_{\scriptscriptstyle\mathrm{RIS}}$

###  IV-G Computational complexity

The overall computational complexity of the proposed T3DRIS approach is given
by the number of iterations necessary for convergence, i.e.,
$\mathrm{I}_{\mathrm{conv}}$, times the sum of the computational complexity of
the algorithm in [34], which is used for the optimization of the RIS
configuration, and the computational complexity of the 3D position update of
each RIS unit cell via iterative Neumann series and projected gradient ascent.
In our single-antenna UE setup, the former is given by
$\mathcal{O}(N^{4}+N^{3})$, where $N$ is the number of RIS elements [34].
Whereas, the complexity of the 3D RIS element positions is dictated
essentially by the calculation of the gradient $\nabla
f(\Delta\operatorname{\mathbf{q}}_{k})$, $k=1,\ldots,N$, whose complexity is
$\mathcal{O}(N)$, the bisection search for the maximum value of $\alpha_{i}$,
which requires $\gamma^{\scriptscriptstyle\mathrm{bis}}$ iterations, with
$\gamma^{\scriptscriptstyle\mathrm{bis}}\leq\left\lceil\log_{2}(\frac{1}{\epsilon^{\scriptscriptstyle\mathrm{bis}}})\right\rceil$
and $\epsilon^{\scriptscriptstyle\mathrm{bis}}$ the tolerance used to reach
convergence, and the solution of Problem 5, which essentially boils down to
the projection onto the feasible set $\mathcal{Q}$ and is again
$\mathcal{O}(N)$ in the considered scenarios. Lastly, all mutual impedances
$Z_{qp}$ are re-calculated for every update of the RIS unit cells. In
particular, for every update of each $\operatorname{\mathbf{q}}_{k}$, the
mutual impedance with the BS, UE and the remaining $N-1$ RIS unit cells need
to be re-calculated. The overall complexity of T3DRIS is thus given by
$\mathcal{O}(\mathrm{I}_{\mathrm{conv}}(N^{4}+N^{3}+N(\gamma^{\scriptscriptstyle\mathrm{bis}}+N+3)))$.

We remark here that, while the computational complexity of T3DRIS may increase
significantly with the number of RIS elements, the proposed approach is
intended to be applied offline during the network planning phase. Indeed, the
optimization of the RIS shape needs to be carried out before its fabrication
and is not intended to be changed dynamically.

(a) ULA, $\lambda/16$ inter-el. spacing.

(b) UPA, $\lambda/2$ inter-el. spacing.

(c) Cylindrical, $\lambda/2$ inter-el. spacing and radius $0.1$ m.

(d) Spherical, $\lambda/2$ inter-el. spacing and radius $0.1$ m.

Figure 3: Four different RIS structures (black color) with their corresponding
non-planar element arrangements, optimally designed with T3DRIS in the
unconstrained case (green color), whose performance in terms of SNR versus
algorithm iterations is shown on the right-hand side.

##  V Performance evaluation

We consider a scenario wherein the center-point of the RIS is located at the
origin, while the single-antenna BS and UE are located at $[1.3,0,0]~{}$m and
$[0.98,0.56,-0.65]~{}$m, respectively, which correspond to a distance of $1.3$
m and azimuth of elevation angles of $30^{\circ}$, $120^{\circ}$ from the
perspective of the RIS, respectively, while the transmit signal from the BS is
incident from the broadside of the RIS. We consider a total of $N=100$
$\lambda/2$-length dipoles forming the RIS, unless otherwise stated.555In
conventional RIS designs, the distance between the element centers is kept at
$\lambda/2$ while the size of the elements themselves, such as patch antennas,
must be slightly smaller to avoid superposition or contact. The initial
positions $\operatorname{\mathbf{Q}}^{(0)}$ of the dipoles follow a regular
conformal structure. In particular, we consider four specific structures, $i$)
uniform linear array (ULA) with $\lambda/16$ inter-element spacing666The setup
at $\lambda/16$ aims at keeping the overall dimension of the ULA proportional
with the other initial structures considered, which results in a comparable
average pathloss among RIS elements and the BS/UE in all the experiments.,
$ii$) uniform planar array (UPA) with $\lambda/2$ inter-element spacing,
$iii$) cylindrical RIS with $\lambda/2$ inter-element spacing and radius $0.1$
m, and $iv$) spherical RIS with $\lambda/2$ inter-element spacing and radius
$0.1$ m, as shown in Fig. 3. Moreover, for each of the specific initial RIS
structures, we distinguish two cases: _i_) unconstrained and _ii)_ constrained
optimization. In the former case, the feasible set for the positions of the
RIS elements is a 3D sphere of radius $R=0.05$ m centered at the origin,
whereas in the latter case, the optimized RIS shape must preserve the general
geometrical shape of the initial position, which are chosen having in mind
practical constraints on RIS fabrication (e.g., cylinder of a given radius, a
plane with given maximum height and width, etc.). In this case, the limits on
the Cartesian/polar coordinates are given as $1.5$ times the ones
corresponding to the initial RIS shape. The operating frequency is set to $30$
GHz, i.e., $\lambda=1$ cm, while we set $R_{0}=0.2$ and $Y_{0}=1$ without loss
of generality. The transmit power at the BS and noise power are set as $P=10$
dBm and $\sigma_{n}^{2}=-80$ dBm, respectively. In all the experiments, the
maximum number of iterations and the convergence threshold are set to
$I_{\scriptscriptstyle\mathrm{max}}=2000$ and $\epsilon=10^{-6}$,
respectively. Lastly, the set of feasible RIS impedances is given by
$\mathcal{B}=[-5000,188]~{}\Omega$. The relevant simulation parameters are
summarized in Table I.

TABLE I: Simulation parameters.

Parameter | Value | Parameter | Value
---|---|---|---
$\operatorname{\mathbf{p}}_{\scriptscriptstyle\mathrm{BS}}$ |  $[1.3,~{}0,~{}0]$ m | $\operatorname{\mathbf{p}}_{\scriptscriptstyle\mathrm{UE}}$ |  $[0.98,~{}0.56,~{}-0.65]$ m
$\mathrm{R}$ |  $0.05$ m | $N$ | $100$
$\rm{R}_{0}$ | $0.2~{}\Omega$ | $P$ |  $10$ dBm
$I_{\scriptscriptstyle\mathrm{max}}$ | $2000$ | $\epsilon$ | $10^{-6}$
$\lambda$ |  $1$ cm | $Y_{0}$ | $1$
$\sigma_{n}^{2}$ |  $-80$ dBm | $\mathcal{B}$ | $[-5000,188]~{}\Omega$

(a) ULA, $\lambda/16$ inter-el. spacing.

(b) UPA, $\lambda/2$ inter-el. spacing.

(c) cylindrical, $\lambda/2$ inter-el. spacing and radius $0.1$ m.

(d) Spherical, $\lambda/2$ inter-el. spacing and radius $0.1$ m.

Figure 4: Four different RIS structures (black color) with their corresponding
non-planar element arrangements, optimally designed with T3DRIS in the
constrained case (green color), whose performance in terms of SNR versus
algorithm iterations is shown on the right-hand side.

###  V-A Numerical results

####  V-A1 Unconstrained case

Fig. 3 shows how the initial RIS structures can be optimally transformed
within a 3D space thanks to the application of T3DRIS. For the four considered
case studies, we can draw the following conclusions: $i$) the optimized RIS
elements spread out at different heights along the $y$\- and $z$-axis
following regular _striped_ patterns, and $ii$) the RIS elements tend to move
closer to the BS and the UE. Indeed, moving the positions of the RIS elements
towards the BS and UE causes a reduction in the path loss between them, which
ultimately increases SNR. Since our focus is mostly on the effects of the RIS
shape, we chose the radius of the 3D sphere wherein the RIS elements can be
freely moved small enough to limit such effects.

Interestingly, in the case of the ULA, the initial RIS structure is not able
to steer energy at different elevation angles, hence the spread over the
$z$-axis is more remarkable. It is noticeable that the optimized positions of
the dipoles significantly increase the resulting SNR, hence empirically
demonstrating that the optimization of the RIS structure is a viable approach
to enhance communication performance.

As for the UPA geometry, we see that optimizing the position and spacing of
the dipoles leads to a noticeable SNR improvement. It is also noticeable that,
similarly to the previous example, the spread of the elements in space
exhibits a regular and periodic pattern, with stripes that are spaced at
regular intervals in the $x$-$y$ plane. Whereas, the cylindrical and spherical
initial geometries tend to result in an optimized surface wherein the elements
are concentrated along equally-spaced curved arches on the $x$-$y$ plane,
which are more pronounced in the spherical case.

Overall, in all the considered scenarios we notice that the optimized RIS
structure follows a certain periodic structure, which mimics the periodicity
of the incoming and reflected plane waves.

Figure 5: Bar plot of the maximum SNR obtained with the conventional and
optimized RIS shapes.

####  V-A2 Constrained case

In Fig. 4, we demonstrate the effectiveness of our proposed T3DRIS approach in
optimizing the element arrangements over a specific constrained shape. As for
the unconstrained case, we notice that in all four considered scenarios, the
optimized RIS element arrangements follow period patterns with non-uniform
inter-element spacings. Moreover, we notice that compared with the
unconstrained case, the SNR improvement is less pronounced, which is due to a
reduction in the degrees of freedom in the set of feasible RIS element
positions.

In Fig. 5, we show the initial SNR for each RIS structure (labeled as _SotA_)
and compare it against the SNR obtained after executing T3DRIS for both the
unconstrained and constrained cases. Note that the initial SNR value
corresponds to the performance of conventional RIS devices, i.e., whose shape
is not optimized for the given application scenario. Notably, the highest gain
compared to the initial shape is obtained for the UPA geometry: this
exacerbates the need to design metamorphic RISs as compared to planar
surfaces. In Fig. 6, we show the distribution of the normalized inter-element
spacings at the $90$th percentile, i.e., $\lambda/d_{ij}$ with $d_{qp}$ the
distance between elements $q$ and $p$, obtained with T3DRIS as compared to the
initial shape.

###  V-B Full-wave electromagnetic simulations

(a) ULA.

(b) UPA.

(c) Cylindrical.

(d) Spherical.

Figure 6: Distribution of the normalized inter-element spacings ($90$th
percentile), i.e., $\lambda/d_{ij}$ with $d_{qp}$ the distance between the RIS
elements $q$ and $p$.

In this section, we provide full-wave EM simulations to validate our proposed
approach against existing methods, using CST Studio Suite [38]. In agreement
with the proposed analytical model for conformal RISs, the full-wave
simulations are implemented by adopting a thin wire dipole as a unit cell. As
for the numerical results in Section V-A, the choice of freestanding dipoles
is motivated by the: $i)$ availability of closed-form mathematical expressions
for the mutual impedance (see (III)), which facilitates the development of
optimization algorithms, and allows for a deep understanding of the role of
mutual coupling in the re-radiation pattern [32, 39]; $ii)$ freedom of
arrangement through curved shells, i.e., as detailed in Section II, the
resonant frequency of the unit cell does not depend on the radius of curvature
of the RIS, similar to the case of planar structures [28]; In this regard, we
compare our proposed approach against the method in [8]. Here, the authors
design the optimized RIS configuration for a given cylindrical shape, by
assuming a conventional geometric channel model based on steering vectors and
ideal unitary reflection coefficients at the unit cells. However, such kind of
models inherently assume that the impinging signal at the unit cell can be
reflected with arbitrary delay and maximum efficiency, which in general is
feasible only if the unit cell surface area is much larger than the signal
wavelength. Moreover, the mutual coupling across the RIS board is not
accounted for. Note that the relationship between reflection coefficients and
load impedances at the RIS is not trivial and is still an open topic [40, 41].
Therefore, in order to operate a fair comparison with our impedance-based
model, we transform the RIS phase shifts obtained with the method in [8]
$\\{\theta_{k}\\}_{k=1}^{N}$ into load impedances $\\{Z_{k}\\}_{k=1}^{N}$ as
[42]

$\displaystyle e^{j\theta_{k}}=\frac{Z_{k}-Z_{0}}{Z_{k}+Z_{0}}\quad\forall k,$
(65)

where $Z_{0}$ is the characteristic (self) impedance of the RIS unit cell,
which can be obtained with (III). In (65), the solution for each $Z_{k}$ leads
to a complex number with both real and imaginary part. In order to have
coherence with the RIS load impedance model in (11), we fix the real part to
the parameter $R_{0}$ and retain the imaginary part obtained by inverting
(65).

(a)  (b)  (c)  (b) SotA algorithm in [8]. (c) T3DRIS -
constrained case. (d) T3DRIS -
unconstrained case.

Figure 7: Heatmap of the full-wave simulation of the radiated power in azimuth
and elevation at a cylindrical RIS of radius $0.1$ m originated by the
conventional method [8], and the optimized RIS shape and configuration via the
proposed T3DRIS approach, and corresponding horizontal cuts at elevation angle
of $120^{\circ}$.

Fig. 7 shows the $2$D heatmap of the beampattern expressed as the power
irradiated along the azimuth and elevation directions for the cylindrical RIS
structure in Section V-A, i.e., with radius $0.1$ m, which is optimized
according to [8] and by using the proposed T3DRIS approach in both the
unconstrained and constrained cases. We note that T3DRIS achieves a gain of
$5$ dBi and $3.5$ dBi for the two cases, respectively, in terms of radiated
power towards the location of the UE, which is highlighted with a red arrow,
as compared to the scheme proposed in [8]. Also, we show the polar plots of
horizontal cuts at elevation angle of $120^{\circ}$, thus demonstrating that
T3DRIS provides a narrower beamwidth, expressed in terms of the half-power
beamwidth (HPBW), and thus a higher directivity, thanks to the employed EM-
compliant and mutual coupling aware RIS model, and the joint optimization of
the RIS shape and configuration. It is worth noting that the pointing
direction is slightly different for the three considered cases. This is
explained by considering that after the application of the proposed T3DRIS
approach, the center of mass of the RIS shape is changed, and therefore the
required pointing direction towards the UE.

###  V-C Analysis of the convergence of T3DRIS

In Fig. 8, on the left-hand side, we empirically show that T3DRIS is
convergent in a finite number of iterations in the considered reference
scenarios. We note that the constrained case requires in general fewer
iterations, since the feasible space is greatly reduced compared to the
unconstrained case. Moreover, on the right-hand side, we show the simulation
time versus the number of RIS elements $N$ for the four considered cases, with
color codes given on the left-hand side, thus validating the analysis provided
in Sections IV-F and IV-G.

##  VI Related work

Non-conventionally shaped antennas and arrays are becoming increasingly
popular due to their multifunctional capabilities and deployments on
supporting surfaces going far beyond traditional flat structures [43]. A
mechanical metasurface constructed from filamentary metal traces that can
self-evolve to morph into a wide range of 3D target shapes is presented in
[44]. In [45], a radiating structure of combined loop dipoles is placed on top
of a flexible Polychlorierte Biphenyle (PCB) substrate, enabling the antenna
to be bent and to be adapted to the geometry of the surface of deployment. A
flexible patch antenna made by conductive fabric encapsulated in a flexible
polymer substrate for wearable devices is presented in [46]. The use of
conductive fabric is proposed in [47] and [48], where the resulting radiating
elements can be directly embedded into clothes. Moreover, the performance of
frequency selective curved antenna arrays is assessed in [49]. In [50], a
cylindrical array is proposed to focus the beams in the near field. A flexible
antenna array capable of self-adaptation to the surface of deployment is
proposed in [51]. Conformal antenna arrays are seen as an aerodynamically
efficient solution to deploy large arrays on air and space-borne applications
in [52, 53]. In [54] a cylindrical array of antennas is proposed to obtain
$360^{\circ}$ beamforming.

Figure 8: T3DRIS convergence and simulation time versus the number of RIS
elements $N$.

Recently, non-planar structures are investigated in the RIS realm. A case
study focused on conformal RIS for vehicle-to-vehicle (V2V) communication is
considered in [55], wherein the performance of an RIS to be mounted on the
side of a car is studied. In particular, the RIS elements are distributed over
a cylindrical surface to be mounted on the side of vehicles and to be
configured with a predefined (static) reflection angle. A metamorphic smart
surface with shape-changing properties is proposed in [9]. Such surfaces can
be deployed on objects whose shape changes over time, e.g., on curtains and
blinds. Interestingly, it is also demonstrated that non-planar surfaces can
match the performance of UPAs with up to $68\%$ fewer antenna elements. In
[12], the shape of a retro-reflective passive array is mechanically changed to
modify the reflection pattern and embed radar-readable information in the
reflected signal. In [13], the geometry of each unit cell is changed to obtain
a frequency-tunable smart surface system. This is achieved by employing
antenna elements designed as a rollable thin plastic film with flexible copper
strips whose length is tuned to change the surface operating frequency for
multi-network coverage extension applications. The concept of non-planar RISs
in wearable devices to enable low-power sensing through the modulation of the
reflected signal towards the BS is envisioned in [3].

The interesting properties of non-regularly (conformal) shaped arrays come
with a non-negligible mutual coupling among the configurable elements, which
has to be taken into account for their optimization and synthesis [56]. To
overcome such a problem, the characterization of the admittance matrix of
antenna arrays serves as a valuable optimization proxy that incorporates both
the mutual coupling and the configuration of the array elements. Hence, it
provides a viable way to synthesize the desired response [56, 57]. In [35],
the admittance matrix is leveraged to model the E2E RIS-aided channel
response. In the absence of mutual coupling, this approach leads to a closed-
form expression for the optimal RIS configuration. When mutual coupling is
present, iterative optimization algorithms are derived to account for its
influence. In [36], it is demonstrated that accounting for the mutual coupling
during the optimization stage results in better system performance. This idea
is pushed to the next stage in [33, 34], wherein the admittance is used to
take into account the presence of non-controllable passive scatterers (i.e.,
the multipath).

While conformal RISs are gaining momentum, related studies are limited in
terms of arbitrarily shaped structures, and overlook the crucial aspect of
accurately modeling the mutual coupling among the elements. To address this
fundamental research gap, we present a pioneering optimization framework that
enables the joint optimization of the RIS response and spatial distribution of
its elements, while adhering to the constraints imposed by the shape of the
deployment support and exploiting the mutual coupling by design. To the best
of our knowledge, this is the first work addressing these fundamental aspects
of RIS design in the context of applications and deployments in future smart
radio environments.

##  VII Conclusions

This paper has investigated the complex reflective properties of conformal
surfaces. Our findings unequivocally demonstrate that non-planar RIS designs
and ad-hoc configurations surpass conventional planar approaches in terms of
versatility and adaptability, thus heralding a paradigm shift in RIS
optimization. Our findings demonstrate that such a counter-intuitive
integration of non-planar geometries enables the dynamic control over EM waves
in 3D, opening up new avenues for innovation and optimization in wireless
communications. In particular, we leverage a recently proposed mutual coupling
model for RISs to propose T3DRIS: an iterative algorithm that jointly
optimizes the configuration and locations of the RIS elements in 3D, subject
to realistic deployment and implementation constraints. Our empirical and
numerical evaluation using full-wave simulations show the superiority of
T3DRIS with respect to SotA solutions. The insights gained from this paper
serve as a _foundation for further exploration and advancements in the field
of smart surface technologies_.

##  Appendix A Self/mutual impedance between RIS elements

We start from the closed-form expression of the mutual impedance between any
two antennas $p$ and $q$ available in [32], by assuming that all the dipole
elements have the same half-length $h$

$\displaystyle Z_{qp}\\!=\\!\frac{\eta
c}{8\pi}\\!\\!\sum_{s_{o}=\\{-1,1\\}}s_{0}e^{js_{0}\frac{2\pi}{\lambda}h}\bigg{[}I_{qp}(h,s_{0})+I_{qp}(-h,s_{0})$
$\displaystyle-2\cos(2\pi h/\lambda)I_{qp}(0,s_{0})\bigg{]},$ (66)

where $c=\frac{1}{\sin(2\pi h/\lambda)}$, and the function $I_{qp}(\xi,s_{0})$
is given by

$\displaystyle I_{qp}(\xi,s_{0})$
$\displaystyle=\mathcal{I}(-s_{0},\operatorname{\boldsymbol{\delta}}_{qp}-\xi\operatorname{\mathbf{e}}_{3},-h,0)$
$\displaystyle+\mathcal{I}(s_{0},\operatorname{\boldsymbol{\delta}}_{qp}-\xi\operatorname{\mathbf{e}}_{3},0,h),$
(67)

where
$\operatorname{\boldsymbol{\delta}}_{qp}=\operatorname{\mathbf{q}}_{q}-\operatorname{\mathbf{q}}_{p}$,
with
$\operatorname{\mathbf{q}}_{q}\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times
1}$ and
$\operatorname{\mathbf{q}}_{p}\in\operatorname{\mbox{$\mathbb{R}$}}^{3\times
1}$ denoting the positions of the two considered dipoles. Moreover, the
function $\mathcal{I}(s_{0},\operatorname{\boldsymbol{\xi}},L,U)$ is defined
as

$\displaystyle\mathcal{I}(s_{0},\operatorname{\boldsymbol{\xi}},L,U)=s_{0}e^{j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\xi}}}\big{[}$
$\displaystyle
T_{0}(\operatorname{\boldsymbol{\xi}}+L\operatorname{\mathbf{e}}_{3},s_{0})$
$\displaystyle-
T_{0}(\operatorname{\boldsymbol{\xi}}+U\operatorname{\mathbf{e}}_{3},s_{0})\big{]},$
(68)

with $T_{0}(\operatorname{\boldsymbol{\zeta}},s_{0})$ given in (13). When
setting $h=\lambda/4$, we obtain $c=1$, $\cos(2\pi h/\lambda)=0$ and
$e^{j\frac{2\pi}{\lambda}s_{0}h}=s_{0}j$, which allows us to simplify the
expression in (66) as

$\displaystyle Z_{qp}=\frac{j\eta}{8\pi}\sum_{s_{o}=\\{-1,+1\\}}\big{[}$
$\displaystyle I_{qp}(h,s_{0})+I_{qp}(-h,s_{0})\big{]}.$ (69)

By substituting (67) and (68) into (69), we obtain (III).

##  Appendix B Gradient of $Z_{qp}$ with respect to the position of the RIS
elements

We compute the expression of the gradient of $Z_{qp}$ with respect to the
position $\operatorname{\mathbf{q}}_{q}$ and $\operatorname{\mathbf{q}}_{p}$
of the dipoles $q$ and $p$, respectively. For simplicity, we rewrite (III) as

$\displaystyle
Z_{qp}=\frac{\eta}{8\pi}\sum_{s_{o}=\\{-1,+1\\}}e^{j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}}g(\operatorname{\boldsymbol{\delta}}_{qp},s_{0}).$
(70)

with
$\operatorname{\boldsymbol{\delta}}_{qp}=\operatorname{\mathbf{q}}_{q}-\operatorname{\mathbf{q}}_{p}$,
and

$\displaystyle g(\operatorname{\boldsymbol{\delta}}_{qp},s_{0})=$
$\displaystyle
T_{0}(\operatorname{\boldsymbol{\delta}}_{qp}-2h\operatorname{\mathbf{e}}_{3},s_{0})+T_{0}(\operatorname{\boldsymbol{\delta}}_{qp}+2h\operatorname{\mathbf{e}}_{3},s_{0})$
$\displaystyle-2T_{0}(\operatorname{\boldsymbol{\delta}}_{qp},s_{0}).$ (71)

The gradient of (70) with respect to $\operatorname{\boldsymbol{\delta}}_{qp}$
is

$\displaystyle\nabla
Z_{qp}=\frac{\eta}{8\pi}\\!\\!\\!\sum_{s_{o}=\\{-1,1\\}}\\!\\!\\!\Big{(}j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}e^{j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{qp}}g(\operatorname{\boldsymbol{\delta}}_{qp},s_{0})$
$\displaystyle+e^{j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\delta}}_{qp}}\nabla
g(\operatorname{\boldsymbol{\delta}}_{qp},s_{0})\Big{)},$ (72)

with

$\displaystyle\nabla g(\operatorname{\boldsymbol{\delta}}_{qp},s_{0})=$
$\displaystyle\nabla
T_{0}(\operatorname{\boldsymbol{\delta}}_{qp}-2h\operatorname{\mathbf{e}}_{3},s_{0})+\nabla
T_{0}(\operatorname{\boldsymbol{\delta}}_{qp}+2h\operatorname{\mathbf{e}}_{3},s_{0})$
$\displaystyle-2\nabla T_{0}(\operatorname{\boldsymbol{\delta}}_{qp},s_{0}).$
(73)

Since $\frac{d\,E_{1}(x)}{d\,x}=-\frac{e^{-x}}{x}$, with
$E_{1}(x)=\int_{x}^{\infty}\frac{e^{-t}}{t}dt$ the exponential integral, we
have that

$\displaystyle\nabla
T_{0}(\operatorname{\boldsymbol{\zeta}},s_{0})=-\Big{(}\frac{\operatorname{\boldsymbol{\zeta}}}{\|\operatorname{\boldsymbol{\zeta}}\|}+s_{0}\operatorname{\mathbf{e}}_{3}\Big{)}\frac{e^{-j\frac{2\pi}{\lambda}(\|\operatorname{\boldsymbol{\zeta}}\|+s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\zeta}})}}{\|\operatorname{\boldsymbol{\zeta}}\|+s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}\operatorname{\boldsymbol{\zeta}}}.$
(74)

Moreover, we have that
$\nabla\operatorname{\boldsymbol{\delta}}_{qp}^{\operatorname{\mathrm{T}}}=\operatorname{\mathbf{I}}_{3}$,
and
$\nabla\operatorname{\boldsymbol{\delta}}_{qp}^{\operatorname{\mathrm{T}}}=-\operatorname{\mathbf{I}}_{3}$.
We are now in the position of obtaining $\nabla Z_{qp}$ by substituting
(73)-(74) in (72), which gives

$\displaystyle\nabla Z_{qp}\\!\\!=$
$\displaystyle\frac{\eta}{8\pi}\\!\\!\\!\sum_{s_{o}=\\{-1,1\\}}\\!\\!\\!\\!\\!\Big{(}j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}e^{j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}(\operatorname{\mathbf{q}}_{q}-\operatorname{\mathbf{q}}_{p})}g(\operatorname{\mathbf{q}}_{q}-\operatorname{\mathbf{q}}_{p},s_{0})$
$\displaystyle+e^{j\frac{2\pi}{\lambda}s_{0}\operatorname{\mathbf{e}}_{3}^{\operatorname{\mathrm{T}}}(\operatorname{\mathbf{q}}_{q}-\operatorname{\mathbf{q}}_{p})}\nabla
g(\operatorname{\mathbf{q}}_{q}-\operatorname{\mathbf{q}}_{p},s_{0})\Big{)}.$
(75)

## References

  * [1] M. Di Renzo, A. Zappone _et al._ , “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2450–2525, 2020.
  * [2] W. Long, R. Chen _et al._ , “A promising technology for 6g wireless networks: Intelligent reflecting surface,” _J. Commun. Informat. Netw._ , vol. 6, no. 1, pp. 1–16, 2021.
  * [3] M. Di Renzo, M. Debbah _et al._ , “Smart radio environments empowered by reconfigurable ai meta-surfaces: An idea whose time has come,” _EURASIP J. Wireless Commun. Netw._ , vol. 2019, no. 1, pp. 1–20, 2019\.
  * [4] S. Basharat, S. A. Hassan _et al._ , “Reconfigurable intelligent surfaces: Potentials, applications, and challenges for 6g wireless networks,” _IEEE Wireless Commun._ , vol. 28, no. 6, pp. 184–191, 2021.
  * [5] Z. Cui and S. Pollin, “Impact of reconfigurable intelligent surface geometry on communication performance,” _IEEE Wireless Commun. Lett._ , vol. 12, no. 5, pp. 898–902, 2023.
  * [6] Y. Lin, S. Jin _et al._ , “Conformal irs-empowered mimo-ofdm: Channel estimation and environment mapping,” _IEEE Trans. Commun._ , vol. 70, no. 7, pp. 4884–4899, 2022.
  * [7] J. Budhu, L. Szymanski, and A. Grbic, “Design of planar and conformal, passive, lossless metasurfaces that beamform,” _IEEE J. Microw._ , vol. 2, no. 3, pp. 401–418, 2022.
  * [8] M. Mizmizi, R. A. Ayoubi _et al._ , “Conformal metasurfaces: A novel solution for vehicular communications,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 4, pp. 2804–2817, 2023.
  * [9] R. I. Zelaya, R. Ma, and W. Hu, “Towards 6G and Beyond: Smarten Everything with Metamorphic Surfaces,” in _Proc. 20th ACM HotNets Workshop_ , 2021, pp. 155–162.
  * [10] G. Gradoni and M. Di Renzo, “End-to-End Mutual Coupling Aware Communication Model for Reconfigurable Intelligent Surfaces: An Electromagnetic-Compliant Approach Based on Mutual Impedances,” _IEEE Wireless Commun. Lett._ , vol. 2337, no. c, pp. 1–5, 2021.
  * [11] M. Di Renzo, F. H. Danufane, and S. Tretyakov, “Communication models for reconfigurable intelligent surfaces: From surface electromagnetics to wireless networks optimization,” _Proc. of the IEEE_ , vol. 110, no. 9, pp. 1164–1209, 2022.
  * [12] J. Nolan, K. Qian, and X. Zhang, “RoS: Passive smart surface for roadside-to-vehicle communication,” in _Proc. ACM SIGCOMM_ , 2021, pp. 165–178.
  * [13] R. Ma, R. I. Zelaya, and W. Hu, “Softly, Deftly, Scrolls Unfurl Their Splendor: Rolling Flexible Surfaces for Wideband Wireless,” in _Proc. ACM Mobicom_ , 2023, pp. 1–15.
  * [14] G. Oliveri, P. Rocca _et al._ , “Holographic smart EM skins for advanced beam power shaping in next generation wireless environments,” _IEEE J. Multiscale Multiphysics Comput. Techniques_ , vol. 6, pp. 171–182, 2021.
  * [15] K.-K. Wong, W. K. New _et al._ , “Fluid antenna system—Part I: Preliminaries,” _IEEE Commun. Lett._ , vol. 27, no. 8, pp. 1919–1923, 2023\.
  * [16] L. Zhu, W. Ma, and R. Zhang, “Modeling and performance analysis for movable antenna enabled wireless communications,” arXiv, 2210.05325, 2022.
  * [17] F. Mattiello, G. Leone, and G. Ruvio, “Analysis and characterization of finite-size curved frequency selective surfaces,” _Studies in Engineering and Technology_ , vol. 2, no. 1, 2015.
  * [18] P. C. Chaumet, “The discrete dipole approximation: A review,” _Mathematics_ , vol. 10, no. 17, p. 3049, 2022.
  * [19] W. Gibson, _The Method of Moments in Electromagnetics_. CRC Press, 2007.
  * [20] H. Zhang, Q. Cheng _et al._ , “Hyperuniform disordered distribution metasurface for scattering reduction,” _Applied Physics Lett._ , vol. 118, no. 10, 2021.
  * [21] F. Yang and Y. Rahmat-Samii, _Surface Electromagnetics: With Applications in Antenna, Microwave, and Optical Engineering_. Cambridge University Press, 2019.
  * [22] D. Barbarić, M. Bosiljevac, and Z. Šipuš, “Analysis of curved metasurfaces based on method of moments,” in _Proc. European Conf. Antennas Propag._ , 2020, pp. 1–5.
  * [23] M. Bodehou, D. González-Ovejero _et al._ , “Method of moments simulation of modulated metasurface antennas with a set of orthogonal entire-domain basis functions,” _IEEE Trans. Antennas Propag._ , vol. 67, no. 2, pp. 1119–1130, 2019.
  * [24] K. Konno, Q. Chen, and R. J. Burkholder, “Fast computation of layered media green’s function via recursive taylor expansion,” _IEEE Antennas and Wireless Propag. Lett._ , vol. 16, pp. 1048–1051, 2016.
  * [25] K. Konno, S. Terranova _et al._ , “Generalised impedance model of wireless links assisted by reconfigurable intelligent surfaces,” 2023.
  * [26] P. Mursia, T. Mazloum _et al._ , “Empirical validation of the impedance-based ris channel model in an indoor scattering environment,” in _Proc. European Conf. Antennas Propag._ , 2023.
  * [27] P. Zheng, R. Wang _et al._ , “Mutual coupling in ris-aided communication: Model training and experimental validation,” 2024.
  * [28] T. Dimitrijevic, J. Jokovic _et al._ , “Tlm mesh impact on textile antenna modelling,” in _Proc. Microwave Mediterranean Symposium_ , 2022, pp. 1–4.
  * [29] Y. Liu, X. Mu _et al._ , “STAR: Simultaneous transmission and reflection for 360° coverage by intelligent surfaces,” _IEEE Wireless Commun._ , vol. 28, no. 6, pp. 102–109, 2021.
  * [30] H. Zhang, S. Zeng _et al._ , “Intelligent omni-surfaces for full-dimensional wireless communications: Principles, technology, and implementation,” _IEEE Commun. Mag._ , vol. 60, no. 2, pp. 39–45, 2022\.
  * [31] Z. Sipus, M. Bosiljevac, and S. Skokic, “Analysis of curved frequency selective surfaces,” in _Proc. European Conf. Antennas Propag._ , 2007, pp. 1–5.
  * [32] M. Di Renzo, V. Galdi, and G. Castaldi, “Modeling the mutual coupling of reconfigurable metasurfaces,” in _Proc. European Conf. Antennas Propag._ , 2023, pp. 1–4.
  * [33] P. Mursia, S. Phang _et al._ , “SARIS: Scattering aware reconfigurable intelligent surface model and optimization for complex propagation channels,” _IEEE Wireless Commun. Lett._ , vol. 12, no. 11, pp. 1–1, 2023\.
  * [34] H. E. Hassani, X. Qian _et al._ , “Optimization of RIS-Aided MIMO – A Mutually Coupled Loaded Wire Dipole Model,” 2023. [Online]. Available: http://arxiv.org/abs/2306.09480
  * [35] X. Qian and M. Di Renzo, “Mutual coupling and unit cell aware optimization for reconfigurable intelligent surfaces,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 6, pp. 10–14, 2020.
  * [36] A. Abrardo, D. Dardari _et al._ , “MIMO Interference Channels Assisted by Reconfigurable Intelligent Surfaces: Mutual Coupling Aware Sum-Rate Optimization Based on a Mutual Impedance Channel Model,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 12, pp. 2624–2628, 2021.
  * [37] D. P. Bertsekas, _Nonlinear programming_. Athena Scientific, 1999.
  * [38] Dassault Systèmes, “CST Studio Suite: Electromagnetic field simulation software,” https://www.3ds.com/products-services/simulia/products/cst-studio-suite/, 2022\.
  * [39] J. A. Russer, D. Semmler _et al._ , “Analysis of intelligent surface-aided mimo communication systems,” in _Proc. 26th Int. ITG WSA & SCC_, 2023, pp. 1–5.
  * [40] M. Nerini, S. Shen _et al._ , “A universal framework for multiport network analysis of reconfigurable intelligent surfaces,” 2023.
  * [41] A. Abrardo, A. Toccafondi, and M. D. Renzo, “Design of reconfigurable intelligent surfaces by using s-parameter multiport network theory – optimization and full-wave validation,” 2023.
  * [42] D. M. Pozar, _Microwave Engineering_. Wiley Press, 2011.
  * [43] B. Mohamadzade, R. B. Simorangkir _et al._ , “Recent developments and state of the art in flexible and conformal reconfigurable antennas,” _Electronics_ , vol. 9, no. 9, p. 1375, 2020.
  * [44] Y. Bai, H. Wang _et al._ , “A dynamically reprogrammable metasurface with self-evolving shape morphing,” _Nature_ , vol. 609, pp. 701–708, 2022.
  * [45] S.-J. Ha, Y.-B. Jung _et al._ , “Reconfigurable beam-steering antenna using dipole and loop combined structure for wearable applications,” _ETRI journal_ , vol. 34, no. 1, pp. 1–8, 2012.
  * [46] R. B. Simorangkir, Y. Yang _et al._ , “A method to realize robust flexible electronically tunable antennas using polymer-embedded conductive fabric,” _IEEE Trans. Antennas Propag._ , vol. 66, no. 1, pp. 50–58, 2017.
  * [47] A. Di Natale and E. Di Giampaolo, “A reconfigurable all-textile wearable uwb antenna,” _Progress In Electromagnetics Research C_ , vol. 103, pp. 31–43, 2020.
  * [48] S. M. Salleh, M. Jusoh _et al._ , “Textile antenna with simultaneous frequency and polarization reconfiguration for wban,” _IEEE Access_ , vol. 6, pp. 7350–7358, 2017.
  * [49] N. Yuan, X.-C. Nie _et al._ , “Accurate analysis of conformal antenna arrays with finite and curved frequency selective surfaces,” _J. Electromagn. Waves Appl._ , vol. 21, no. 13, pp. 1745–1760, 2007.
  * [50] Y. F. Wu and Y. J. Cheng, “Proactive conformal antenna array for near-field beam focusing and steering based on curved substrate integrated waveguide,” _IEEE Trans. Antennas Propag._ , vol. 67, no. 4, pp. 2354–2363, 2019.
  * [51] B. D. Braaten, S. Roy _et al._ , “A self-adapting flexible (selflex) antenna array for changing conformal surface applications,” _IEEE Trans. Antennas Propag._ , vol. 61, no. 2, pp. 655–665, 2012.
  * [52] T. E. Morton and K. M. Pasala, “Pattern synthesis of conformal arrays for airborne vehicles,” in _Proc. IEEE Aerospace Conf._ , vol. 2. IEEE, 2004, pp. 1030–1039.
  * [53] B. A. Nunna and V. K. Kothapudi, “Design and analysis of x-band conformal antenna array for spaceborne synthetic aperture radar applications,” in _Proc. Second Int. Conf. Computational Electron. Wireless Commun._ Springer, 2023, pp. 193–204.
  * [54] S. Vellucci, M. Longhi _et al._ , “Phase-gradient huygens metasurface coatings for dynamic beamforming in linear antennas,” 2023.
  * [55] D. Tagliaferri, M. Mizmizi _et al._ , “Conformal intelligent reflecting surfaces for 6g v2v communications,” in _Proc. 1st Int. Conf. 6G Netw._ IEEE, 2022, pp. 1–8.
  * [56] R. Coe and A. Ishimaru, “Optimum scattering from an array of half-wave dipoles,” _IEEE Trans. Antennas Propag._ , vol. 18, no. 2, pp. 224–230, 1970.
  * [57] J. Mautz and R. Harrington, “Modal analysis of loaded n-port scatterers,” _IEEE Trans. Antennas Propag._ , vol. 21, no. 2, pp. 188–199, 1973.