# De-risking Carbon Capture and Sequestration with Explainable CO2 Leakage
Detection in Time-lapse Seismic Monitoring Images

Huseyin Tuna Erdinc,1 Abhinav Prakash Gahlot,2 Ziyi Yin,3
Mathias Louboutin,2 Felix J. Herrmann1,2,3

###### Abstract

With the growing global deployment of carbon capture and sequestration
technology to combat climate change, monitoring and detection of potential CO2
leakage through existing or storage induced faults are critical to the safe
and long-term viability of the technology. Recent work on time-lapse seismic
monitoring of CO2 storage has shown promising results in its ability to
monitor the growth of the CO2 plume from surface recorded seismic data.
However, due to the low sensitivity of seismic imaging to CO2 concentration,
additional developments are required to efficiently interpret the seismic
images for leakage. In this work, we introduce a binary classification of
time-lapse seismic images to delineate CO2 plumes (leakage) using state-of-
the-art deep learning models. Additionally, we localize the leakage region of
CO2 plumes by leveraging Class Activation Mapping methods.

## Introduction

According to the International Energy Agency and the International Panel on
Climate Change report (IPCC 2018), there is a need for a 50 percent reduction
of greenhouse gas emissions by 2050 to avoid an increase of 1.5 degrees
Celsius of Earth’s average temperature. This can only be achieved by reduced
dependence on fossil fuels, use of renewable sources of energy and large-scale
global deployment of carbon reduction technologies such as carbon capture and
sequestration (CCS). This technology consists of collection, transportation,
and injection of CO2 into an appropriate geologic storage reservoir for
extended time periods (tens of years). Especially, unlike other solutions, CCS
is considered a relatively low-cost, long-term and imminent solution. However,
potential CO2 leakage from the underground reservoirs due to pre-existing or
pressure-induced faults poses risks (Ringrose 2020). Thus, it is necessary to
de-risk CCS projects by monitoring CO2 plumes in order to accurately detect
and predict potential leakages as early as possible.

Time-lapse seismic monitoring has been introduced as a reliable technology to
monitor the CO2 dynamics in the Earth’s subsurface during carbon sequestration
(Lumley 2001) and is already in use at existing storage sites (Arts et al.
2008; Chadwick et al. 2010; Ringrose et al. 2013; Furre et al. 2017). In
essence, sequential (i.e once every 6 months/year/…) seismic datasets, called
vintages, are collected in the field over an area covering the storage
reservoir. Then, each seismic dataset is inverted to obtain high fidelity
images of the subsurface over time (Arts et al. 2008; Ayeni and Biondi 2010;
Yin, Louboutin, and Herrmann 2021). The evolution of the CO2 reservoir can
finally be visualized by subtracting the seismic images between different
points in time. However, due to the inherently weak and noisy amplitudes of
the CO2 reservoir’s response in those seismic images, detecting the presence
of potential irregularities, such as in CO2 plumes, corresponding to a leakage
is a challenging problem. To tackle this difficulty, we propose a machine
learning based detection method based on standard binary classification.

Recently, numerous methods leveraging machine learning have been introduced
for the detection of CO2 leakage based on a simple artificial neural network
(ANN) (Li et al. 2018), and a combination of convolutional neural networks
(CNN) and Long Short-Term Memory (LSTM) networks (Zhou et al. 2019). While
leading to accurate predictions, these methods usually rely solely on the
field recorded data rather than the subsurface seismic images. Besides,
practical considerations such as repeatability (the ability to record the data
in the exact same way every year) hinders their applicability to real world
cases. On the other hand, as we rely on visualizing the CO2 plumes in the
seismic image, we can take advantage of advanced seismic imaging techniques
designed for non-repeated seismic acquisition such as the joint recovery model
(JRM) (Oghenekohwo and Herrmann 2017a; Wason, Oghenekohwo, and Herrmann 2017;
Yin, Louboutin, and Herrmann 2021). Additionally, this imaging technique has
demonstrated higher fidelity imaging than sequential seismic imaging allowing
for easier detection of CO2 leakage.

We will show in the following sections that we can efficiently and accurately
detect CO2 from realistic seismic images recovered by JRM on synthetic but
representative models of the Earth subsurface. We demonstrate our method using
different state-of-the-art deep learning models in a transfer learning setting
to classify CO2 plume seismic images with regular (no-leakage) CO2 plume or
with CO2 leakage. As CO2 leakage detection needs trustworthiness, we further
unravel the decisions made by our models and utilize Class Activation Mapping
(CAM) methods (Zhou et al. 2015) to identify and visualize seismic image areas
crucial for model classification results. We show that the CAM result
accurately focuses on the leakage portion of the CO2 plume and reservoir,
validating that our network detects leakage based on state of the CO2
reservoir over time.

Our main contributions are organized as follows. First, we introduce the
classification models used for leakage detection and the CAM methods for
visualizing the area of interest in the seismic image. Second, we demonstrate
the accuracy of our models and qualitatively examine the results of our CAM
methods on a realistic synthetic set of CO2 plume images.

## Methodology

In order to speed up the training process and to compensate for the
overfitting that may occur with modest sized datasets, we rely on transfer
learning (Yosinski et al. 2014) using pre-trained state-of-the-art models as a
starting point. In particular, we consider four modern architectures known to
achieve high accuracy on standard dataset such as ImageNet-1k (Russakovsky et
al. 2014). The models used are VGG (Simonyan and Zisserman 2014), ResNet (He
et al. 2015), Vision Transformer (ViT) (Dosovitskiy et al. 2020), and Swin
Transformer (Swin) (Liu et al. 2021), all pre-trained on the standardized
ImageNet-1k dataset.

VGG: is a convolutional neural network (CNN) model that achieved significant
success in The ImageNet Large Scale Visual Recognition Challenge (ILSVRC)
competition in 2014 (Simonyan and Zisserman 2014). VGG consists of sequences
of convolution and maxpool layers. In our numerical experiments, the VGG16
variant with 16 trainable layers is used.

ResNet: is a CNN architecture with residual connections proposed to solve the
vanishing gradient problem in very deep networks (He et al. 2015). ResNet
consists of residual blocks and each residual block has convolution layers and
shortcut connections performing identity mapping. In our numerical
experiments, the ResNet34 variant with 34 trainable layers is used.

ViT: is an architecture based on transformer which is used in the field of
Natural Language Processing (NLP) (Vaswani et al. 2017). Internally, the
transformer learns a relationship between input token pairs, and uses 16x16
patches of images as input tokens (Dosovitskiy et al. 2020). In our numerical
experiments, the tiny ViT variant is used allowing lower memory and
computational imprint.

Swin: is a special type of ViT that represents image patches hierarchically by
starting from small-sized patches and gradually increasing the size through
merging to achieve scale-invariance property (Liu et al. 2021). Compared to
ViT, Swin transformer has superior (linear) computational efficiency by
computing self-attention within certain patches of a window. In our numerical
experiments, tiny Swin variant is used allowing lower memory and computational
imprint.

### CAM Methods

Deep learning models for classification are notoriously treated as “black
boxes” as they do not expose their internal knowledge or operations to its
users and do not provide interpretable results. To solve this problem, CAM
based saliency maps (heatmaps) were introduced to highlight the most class-
discriminative regions of to-be-classified input images (Zhou et al. 2015).
Since CO2 leakage requires high fidelity, transparent and interpretable
models, we use CAM to further make our model results explainable and highlight
the regions of the seismic image that are most relevant to the classification
results. In our study, we considered two CAM methods. First, Grad-CAM
(Selvaraju et al. 2019), a gradient-based CAM method considered as the state-
of-the-art in terms of explainability of neural networks for classification.
This CAM method extracts gradients from a specific layer of a model and
computes the weighted average of that specific layer’s activations. Second, we
consider Score-CAM (Wang et al. 2019), a perturbation based CAM method. Score-
CAM also computes the weighted average of activations of a user-specified
layer but, unlike Grad-CAM, Score-CAM relies on propagating (forward pass
through the network) a masked input image where the mask is obtained via
upsampling the activations of the user-defined layer. This CAM method provides
high accuracy and interpretable heatmaps and alleviates potential noise and
spread present in the gradient used for the Grad-CAM heatmaps.

Figure 1: Workflow for CO2 Leakage Monitoring

## Numerical Case Study

To generate the training dataset of CO2 plume evolution, we used five 2D
vertical slices extracted from the 3D Compass velocity model (E. Jones et al.
2012) shown in Fig. 1(a). This model is a synthetic but realistic model
representative of the complex geology of the southeast of the North Sea. The
dimension of each model (slice) used in our work is $2131$ X $4062$ m2. We
used FwiFlow (Li et al. 2020), to simulate the CO2 flow dynamics and JUDI
(Witte et al. 2019a) to model the seismic data and compute the seismic images
of the subsurface.

Hyperparameters | VGG16 | ResNet34 | ViT | Swin
---|---|---|---|---
Batch Size | 8 | 8 | 8 | 8
Learning Rate |  $5$x$10^{-5}$ |  $6$x$10^{-3}$ |  $4$x$10^{-3}$ | $10^{-3}$
Exp Decay Rate($\gamma$) | $0.95$ | $0.92$ | $0.98$ | $0.98$
Table 1: Training hyperparameters for the four models. All models were trained
with the same number of epochs and optimizer.

### Time-lapse reservoir and seismic simulation

We consider a realistic two well setting with a fixed injection well injecting
CO2 and a production well extracting brine from subsurface storage reservoir.
Injection of supercritical CO2 into saline aquifers is an example of multi-
phase flow in porous media. While we understand more complicated geothermal,
geochemical and geomechanical process may eventually be considered to model
the CO2 dynamics, we follow the two-phase immiscible incompressible flow
physics, which in its leading order describes the process of supercritical CO2
displacing brine in the pore space of the rock. The system is governed by
conservation of mass and Darcy’s law. We refer to the existing literature (Li
et al. 2020; Wen, Tang, and Benson 2021) (Li et al. 2020) for more details
about this physical system.

Using empirical relation and the Kozeny-Carman equation(Costa 2006), the
acoustic properties (velocity and density) from the Compass model were
converted into permeability and porosity (Fig. 1(b)) to simulate the multi-
phase flow (CO2 and brine in porous media) in the reservoir. We used
FwiFlow.jl (Li et al. 2020) to solve multi-phase flow equations based on the
finite volume method. We simulated the CO2 flow for a duration varying between
$7$ to $12$ years (Fig. 1(c)). The reservoir was initially filled with saline
water and we injected compressed CO2 at the rate of $1$ MT/day into the
reservoir for all simulations. In order to mimic CO2 leakage, we then created
a fracture at a random location along the top seal of the reservoir when the
pressure induced by the CO2 injection reaches a threshold of $15$ MPa. We then
converted back these simulated CO2 saturation snapshots over time into wave
properties with the patchy saturation model (Avseth, Mukerji, and Mavko 2010)
to obtain time-lapse subsurface models (Fig. 1(d)). We used this model because
at higher pressure condition, local fluid flow slows down resulting in an
acoustic velocity trend which follows patchy saturation
(Li2018ExperimentalSA).

Based on these models, we then simulated the baseline seismic survey
corresponding to the initial stage (before the injection of CO2) and the
monitor seismic survey corresponding to the final stage at the end of the
reservoir simulation (Fig. 1(e)). As mentioned in the introduction, it is very
difficult to exactly replicate the baseline and monitor surveys. In order to
mimic the realistic scenario in the field, the baseline and monitor datasets
were simulated using different acquisition geometries (position of the
measurements). Finally, we recovered the time-lapse seismic images using JRM
(Oghenekohwo and Herrmann 2017b; Wason, Oghenekohwo, and Herrmann 2017; Yin,
Louboutin, and Herrmann 2021) to alleviate potential noise and inaccuracies in
the seismic images in the case of non-replicated time-lapse surveys. These
recovered images along with the label (leakage/no-leakage) serve as the input
to the classification network. We generated a total of $1000$ leakage and
$870$ no-leakage scenarios, and computed the baseline, monitor and difference
images with the JRM method in each case.

Model | Accuracy | Precision | Recall | F1 | ROC-AUC
---|---|---|---|---|---
VGG16 | 0.920$\pm$ 0.089 | 0.941$\pm$ 0.133 | 0.921$\pm$ 0.081 | 0.927$\pm$ 0.075 | 0.920$\pm$ 0.076
ResNet34 | 0.948$\pm$ 0.020 | 0.982$\pm$ 0.028 | 0.928$\pm$ 0.044 | 0.948$\pm$ 0.040 | 0.967$\pm$ 0.019
ViT | 0.857$\pm$ 0.018 | 0.910$\pm$ 0.102 | 0.820$\pm$ 0.098 | 0.859$\pm$ 0.036 | 0.923$\pm$ 0.023
Swin | 0.836$\pm$ 0.036 | 0.881$\pm$ 0.108 | 0.818 $\pm$ 0.078 | 0.841$\pm$ 0.076 | 0.909$\pm$ 0.007
|  |  |  |  |
Table 2: Comparison of performance (for precision and recall, positives
represent leakage whereas negatives are no leakage) on the test dataset for
our four neural networks. The highest performance for each metric is
highlighted in bold. Figure 2: Grad-CAM and Score-CAM saliency maps overlayed
on the corresponding input seismic image containing a CO2 plume from leakage.
The CO2 plume can be seen on the seismic image as the high amplitude event at
1.3km depth and 1.8km in X.

### Training

The seismic difference images (difference between baseline and monitor
recovery results) were converted to 224x224 gray-scale images with bi-linear
interpolation and transformed into three channel images where each channel is
a copy of the actual gray-scale image. For the classification, the image
dataset was randomly split into an 80% training set and 20% test set. The
training set was then further divided into two parts, one for model parameter
updating (training) and another for hyperparameter tuning (validation). The
training hyperparameters from this second part are summarized in Table 1. For
training, we replaced the last fully connected layers (classification layers)
of each model with a new fully connected layer. We then trained the network
(Fig. 1(g)) in two steps. First, we only trained the last classification
layer, by freezing all the other layers, for 100 epochs. Since most of the
layers are fixed and do not need gradient updates, this first stage is
extremely cheap and computationally efficient. Second, we further trained the
full model for an additional 30 epochs to allow fine-tuning of all layers for
our specific classification task. Following standard practices in
classification settings, we used a binary cross-entropy loss function and the
Adam optimizer (Kingma and Ba 2015) for all models. Finally, after the
training (Fig. 1 (h)), we implemented the CAM based methods (Fig. 1 (i)). We
used the last convolutional layer activations for the CNN models, and the
activations preceding the last attention layer for the transformer-based
models.

## Analysis

We show on Table 2, different performance metrics on our testing dataset,
after training our four networks, with means and confidence intervals after 15
different runs. In detail, we show standard metrics such as accuracy,
precision, and recall. Additionally, we also show F1 score (Chinchor 1992),
that combine recall and accuracy, and area under curve of receiver operating
characteristic (ROC-AUC) (Bradley 1997) to further evaluate the classification
performance of models. We observe in Table 2 that the CNN models outperform
the transformer variants in all the metrics by a significant margin and that
ResNet34 achieves the best performance in all the measures of evaluation. This
result is consistent with the literature, hinting that despite being very
accurate on a specific task, transformers do not generalize well with our
modest sized dataset (Dosovitskiy et al. 2020). Additionally, we observe that
all models lead to better precision compared to recall (more false negatives
than false positives). This discrepancy can be attributed to the fact that
certain leakage images have very small CO2 leakage areas (up to a single
pixel) in the seismic images and are consequently very difficult to detect.

Second, we show in Fig. 2 the CAM results of each model on a single seismic
image from our test dataset. The high amplitude area shows the regions of the
seismic images that are most important to the classifier. As expected, those
heatmaps provide an explainable representation of the classification as the
high amplitudes align with the CO2 leakage part of the seismic image. We
observe that for the CNN, the saliency maps are well centered on the CO2
leakage portion despite being very coarse. Because of this coarseness, both
Grad-CAM and score-CAM provide similar results. On the other hand,
transformer-based networks lead to more focused saliency maps that target the
location of the CO2 leakage extremely well. We observe in that case, the
Score-CAM leads to reduction of aliases and noise compared to the Grad-CAM
results. This can be linked to the potential presence of noise in the
gradients of the transformers as the networks are very deep (Wang et al.
2019).

## Conclusion

We have introduced an interpretable deep-learning method for CO2 leakage
detection with very high accuracy on a synthetic but realistic model of a CO2
sequestration reservoir. First, we showed through four state-of-the-art models
that we can detect potential CO2 leakage from the recovered time-lapse seismic
images. Second, we demonstrated that CAM provides an interpretable and
accurate visualization of the CO2 plume in case of leakage. Additionally, we
showed that transformer-based models (ViT, Swin) led to more focused CAM and
that Score-CAM provided cleaner and therefore more explainable heatmaps. On
the other hand, we found that standard CNNs led to better classification
results and therefore better leakage detection. In particular, ResNet model
performed best and achieved a very high score above $90$% in every evaluation
metric. Future work will focus on improving the classification network to
achieve higher accuracy in leakage detection and on refining the heatmaps for
better explainability.

## Acknowledgments

This research was carried out with the support of Georgia Research Alliance
and partners of the ML4Seismic Center. The authors thank Philipp A. Witte at
Microsoft for the constructive discussion.

## References

  * Arts et al. (2008) Arts, R. J.; Chadwick, A.; Eiken, O.; Thibeau, S.; and Nooner, S. L. 2008.  Ten years’ experience of monitoring CO2 injection in the Utsira Sand at Sleipner, offshore Norway.  _First Break_ , 26.
  * Avseth, Mukerji, and Mavko (2010) Avseth, P.; Mukerji, T.; and Mavko, G. 2010.  _Quantitative seismic interpretation: Applying rock physics tools to reduce interpretation risk_.  Cambridge university press.
  * Ayeni and Biondi (2010) Ayeni, G.; and Biondi, B. 2010.  Target-oriented joint least-squares migration/inversion of time-lapse seismic data sets.  _Geophysics_ , 75.
  * Bradley (1997) Bradley, A. P. 1997.  The use of the area under the ROC curve in the evaluation of machine learning algorithms.  _Pattern Recognition_ , 30(7): 1145–1159.
  * Chadwick et al. (2010) Chadwick, A.; Williams, G.; Delepine, N.; Clochard, V.; Labat, K.; Sturton, S.; Buddensiek, M.-L.; Dillen, M.; Nickel, M.; Lima, A. L.; Arts, R.; Neele, F.; and Rossi, G. 2010.  Quantitative analysis of time-lapse seismic monitoring data at the Sleipner CO2 storage operation.  _The Leading Edge_ , 29(2): 170–177.
  * Chinchor (1992) Chinchor, N. 1992.  MUC-4 Evaluation Metrics.  In _Proceedings of the 4th Conference on Message Understanding_ , MUC4 ’92, 22–29. USA: Association for Computational Linguistics.  ISBN 1558602739.
  * Costa (2006) Costa, A. 2006.  Permeability-porosity relationship: A reexamination of the Kozeny-Carman equation based on a fractal pore-space geometry assumption.  _Geophysical Research Letters_ , 33(2).
  * Dosovitskiy et al. (2020) Dosovitskiy, A.; Beyer, L.; Kolesnikov, A.; Weissenborn, D.; Zhai, X.; Unterthiner, T.; Dehghani, M.; Minderer, M.; Heigold, G.; Gelly, S.; Uszkoreit, J.; and Houlsby, N. 2020.  An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale.  _arXiv_.
  * E. Jones et al. (2012) E. Jones, C.; A. Edgar, J.; I. Selvage, J.; and Crook, H. 2012.  Building Complex Synthetic Models to Evaluate Acquisition Geometries and Velocity Inversion Technologies.  _European Association of Geoscientists & Engineers_, cp-293-00580.
  * Furre et al. (2017) Furre, A.-K.; Eiken, O.; Alnes, H.; Vevatne, J. N.; and Kiær, A. F. 2017.  20 Years of Monitoring CO2-injection at Sleipner.  _Energy Procedia_ , 114: 3916–3926.  13th International Conference on Greenhouse Gas Control Technologies, GHGT-13, 14-18 November 2016, Lausanne, Switzerland.
  * He et al. (2015) He, K.; Zhang, X.; Ren, S.; and Sun, J. 2015.  Deep Residual Learning for Image Recognition.  _arXiv_.
  * IPCC (2018) IPCC. 2018.  Global warming of 1.5°C. An IPCC Special Report on the impacts of global warming of 1.5°C above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty.  _In Press_.
  * Kingma and Ba (2015) Kingma, D. P.; and Ba, J. 2015.  Adam: A Method for Stochastic Optimization.  In _ICLR (Poster)_.
  * Li et al. (2018) Li, B.; Zhou, F.; Li, H.; Duguid, A.; Que, L.; Xue, Y.; and Tan, Y. 2018.  Prediction of CO2 leakage risk for wells in carbon sequestration fields with an optimal artificial neural network.  _International Journal of Greenhouse Gas Control_ , 68: 276–286.
  * Li et al. (2020) Li, D.; Xu, K.; Harris, J. M.; and Darve, E. 2020.  Coupled Time-Lapse Full-Waveform Inversion for Subsurface Flow Problems Using Intrusive Automatic Differentiation.  _Water Resources Research_ , 56(8): e2019WR027032.  E2019WR027032 10.1029/2019WR027032.
  * Liu et al. (2021) Liu, Z.; Lin, Y.; Cao, Y.; Hu, H.; Wei, Y.; Zhang, Z.; Lin, S.; and Guo, B. 2021\.  Swin Transformer: Hierarchical Vision Transformer using Shifted Windows.  _arXiv_.
  * Lumley (2001) Lumley, D. E. 2001.  Time-lapse seismic reservoir monitoring.  _GEOPHYSICS_ , 66(1): 50–53.
  * Oghenekohwo and Herrmann (2017a) Oghenekohwo, F.; and Herrmann, F. J. 2017a.  Highly repeatable time-lapse seismic with distributed Compressive Sensing–-mitigating effects of calibration errors.  _The Leading Edge_ , 36(8): 688–694.  (The Leading Edge).
  * Oghenekohwo and Herrmann (2017b) Oghenekohwo, F.; and Herrmann, F. J. 2017b.  Improved time-lapse data repeatability with randomized sampling and distributed compressive sensing.  In _EAGE Annual Conference Proceedings_.  (EAGE, Paris).
  * Ringrose (2020) Ringrose, P. 2020.  _How to store CO2 underground: Insights from early-mover CCS Projects_ , volume 129.  Springer.
  * Ringrose et al. (2013) Ringrose, P.; Mathieson, A.; Wright, I.; Selama, F.; Hansen, O.; Bissell, R.; Saoula, N.; and Midgley, J. 2013.  The In Salah CO2 Storage Project: Lessons Learned and Knowledge Transfer.  _Energy Procedia_ , 37: 6226–6236.  GHGT-11 Proceedings of the 11th International Conference on Greenhouse Gas Control Technologies, 18-22 November 2012, Kyoto, Japan.
  * Russakovsky et al. (2014) Russakovsky, O.; Deng, J.; Su, H.; Krause, J.; Satheesh, S.; Ma, S.; Huang, Z.; Karpathy, A.; Khosla, A.; Bernstein, M.; Berg, A. C.; and Fei-Fei, L. 2014.  ImageNet Large Scale Visual Recognition Challenge.
  * Selvaraju et al. (2019) Selvaraju, R. R.; Cogswell, M.; Das, A.; Vedantam, R.; Parikh, D.; and Batra, D. 2019.  Grad-CAM: Visual Explanations from Deep Networks via Gradient-Based Localization.  _International Journal of Computer Vision_ , 128(2): 336–359.
  * Simonyan and Zisserman (2014) Simonyan, K.; and Zisserman, A. 2014.  Very Deep Convolutional Networks for Large-Scale Image Recognition.  _arXiv_.
  * Vaswani et al. (2017) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, L.; and Polosukhin, I. 2017.  Attention Is All You Need.  _arXiv_.
  * Wang et al. (2019) Wang, H.; Wang, Z.; Du, M.; Yang, F.; Zhang, Z.; Ding, S.; Mardziel, P.; and Hu, X. 2019.  Score-CAM: Score-Weighted Visual Explanations for Convolutional Neural Networks.  _arXiv_.
  * Wason, Oghenekohwo, and Herrmann (2017) Wason, H.; Oghenekohwo, F.; and Herrmann, F. J. 2017.  Low-cost time-lapse seismic with distributed compressive sensing–-Part 2: impact on repeatability.  _Geophysics_ , 82(3): P15–P30.  (Geophysics).
  * Wen, Tang, and Benson (2021) Wen, G.; Tang, M.; and Benson, S. M. 2021.  Towards a predictor for CO2 plume migration using deep neural networks.  _International Journal of Greenhouse Gas Control_ , 105: 103223.
  * Witte et al. (2019a) Witte, P. A.; Louboutin, M.; Kukreja, N.; Luporini, F.; Lange, M.; Gorman, G. J.; and Herrmann, F. J. 2019a.  A large-scale framework for symbolic implementations of seismic inversion algorithms in Julia.  _Geophysics_ , 84(3): F57–F71.  (Geophysics).
  * Yin, Louboutin, and Herrmann (2021) Yin, Z.; Louboutin, M.; and Herrmann, F. J. 2021.  Compressive time-lapse seismic monitoring of carbon storage and sequestration with the joint recovery model.  In _SEG Technical Program Expanded Abstracts_ , 3434–3438.  (IMAGE, Denver).
  * Yosinski et al. (2014) Yosinski, J.; Clune, J.; Bengio, Y.; and Lipson, H. 2014.  How transferable are features in deep neural networks?  _arXiv_.
  * Zhou et al. (2015) Zhou, B.; Khosla, A.; Lapedriza, A.; Oliva, A.; and Torralba, A. 2015.  Learning Deep Features for Discriminative Localization.  _arXiv_.
  * Zhou et al. (2019) Zhou, Z.; Lin, Y.; Zhang, Z.; Wu, Y.; Wang, Z.; Dilmore, R.; and Guthrie, G. 2019\.  A data-driven CO2 leakage detection using seismic data and spatial-temporal densely connected convolutional neural networks.  _International Journal of Greenhouse Gas Control_ , 90: 102790.


# Generalised Impedance Model of Wireless Links Assisted by Reconfigurable
Intelligent Surfaces

Keisuke Konno, Qiang Chen
Department of Communications Engineering
Graduate School of Engineering
Tohoku University
Japan
<EMAIL_ADDRESS>
&Sergio Terranova, Gabriele Gradoni
School of Mathematical Sciences
University of Nottingham
United Kingdom
<EMAIL_ADDRESS>

###### Abstract

We devise an end-to-end communication channel model that describes the
performance of RIS-assisted MIMO wireless links. The model borrows the
impedance (interaction) matrix formalism from the Method of Moments and
provides a physics-based communication model. In configurations where the
transmit and receive antenna arrays are distant from the RIS beyond a
wavelength, a reduced model provides accurate results for arbitrary RIS unit
cell geometry. Importantly, the simplified model configures as a cascaded
channel transfer matrix whose mathematical structure is compliant with widely
accepted, but less accurate, system level RIS models. A numerical validation
of the communication model is presented for the design of binary RIS
structures with scatterers of canonical geometry. Attained results are
consistent with path-loss models: For obstructed line-of-sight between
transmitter and receiver, the channel capacity of the (optimised) RIS-assisted
link scales as $R^{-2}$, with $R$ RIS-receiver distance at fixed transmitter
position. Our results shows that the applicability of communication models
based on mutual impedance matrices is not restricted to canonical minimum
scattering RIS unit cells.

NOTE: This work has been submitted to the IEEE TAP for possible publication.
Copyright may be transferred without notice, after which this version may no
longer be accessible.

##  1 Introduction

Reconfigurable intelligent surfaces (RIS) provides a viable technology for
molding the flow of electromagnetic wave energy within arbitrary propagation
environments [1]. Wave control capabilities demonstrated by the RIS include
anomalous reflection, beam-forming, focusing, and defocusing [2]. Their
operation has been studied in free-space, both near-field and far-field [3,
4], as well as in indoor and outdoor environments [5]. The modelling of RIS-
assisted wireless links have became a topic of paramount importance [6], and
several authors have proposed physics-based communications models [7]. The
scientific literature classifies first-principle electromagnetics model of the
RIS in two families: Continuous and Discrete models [8]. Continuous models
treat meta-surfaces as homogeneous sheets of currents whole effective
impedance and admittance are matched by design to perform anomalous
reflection. Discrete models treat individual pixels: In the first models
inspired from antenna arrays, the effective scattered field results from the
superposition of waves re-radiated upon specular reflection of individual RIS
unit cells. However, it has been shown that effective models based on re-
radiated fields: i) do not correctly capture the mutual coupling between unit
cells within finite size RIS structures [9, 10]; ii) do rely on local
periodicity which provides an approximate behaviour of the field re-radiated
by finite-size macro cells [11]; iii) do not accurately capture the effect of
loading circuitry on the re-radiated field, as well as parasitic effects
within the unit cell [12]. More recently, discrete models of the RIS have been
put forward using both impedance [13] and scattering [14] matrix formalisms.
Those models are capable of treating a finite-size RIS, capturing joint
amplitude and phase variability of the RIS unit cell response, and providing
an end-to-end (E2E) representation of the wireless link, i.e., include the
interaction of the RIS with the transmit and receive arrays associated to the
base-station (BS) to user-equipment (UE) respectively. A specific impedance-
based channel model has been formulated for arrays of dipoles [15], and
exploited to analyse the benefit of the knowledge of mutual coupling between
dipoles in channel gain optimisation [19]. While the model works in both near-
and far-field, the RIS unit cell has been assumed to be a minimal scattering
linear antenna, for which the surface current density is known in closed form
and the array active impedance can be computed efficiently including mutual
coupling [24]. The structural scattering of the open-circuit [21] baseline RIS
structure needs to be considered and included in impedance-based models. We
extend the model in [15] to arbitrary arrays of scatterers adopted for RIS
using the Method-of-Moment (MoM) formalism. The MoM code that we use has been
recently developed in [27], which has used the methodology originally proposed
in the theory of loaded scatterers [21], where port segments have been
separated from passive segments in the MoM discretisation. A MoM study has
been recently proposed for the path loss of a RIS-assisted SISO system, where
freestanding loaded loops have been assumed as RIS unit cells [28]. In our MoM
model, the ground plane has been included and the formulation is amenable to
accounting for single- as well as multi-layer substrates by using standard
methodologies [43]. Furthermore, our model provides an end-to-end
communication model that can be conveniently used for optimisation and
performance analysis of wireless links, and it naturally extends the model in
[15] beyond the canonical minimum scattering (CMS) approximation. This
extension has the benefit of making important work on CMS structures
practical: This work has been achieved in the literature lately, also
considering wide-band response of linear antennas [16], and including e.g.,
multi-path fading generated by clusters of dipoles in the impedance [17] and
the Green function [18] formalisms. Previous studies have discussed the
inclusion of circuital elements within periodic meta-surfaces with sub-
wavelength unit cells in MoM [25], There, a simple development of the MoM
interaction matrix allowed the computation of the scattered field for
arbitrary distribution of discrete loads. This boils down to add the impedance
matrix of the equivalent circuit of the tuneable elements to the MoM
interaction impedance matrix of port basis segments. In RIS design, it is
worth noticing that the inner coupling between loading circuitry and planar
metallic structure might generate parasitic effects, which are usually
captured by a more complicated (and structure dependent) equivalent circuit
[26]. In the far-field of a single RIS unit cell, our model reduces to a
cascaded formula that has the same mathematical structure of system levels
MIMO communication models [15]. However, the impedance matrix are augmented
with the self- and mutual- interactions involving scattering segments besides
port segments. Importantly, while system level models are not constructed in
physics-based grounds and are prone to efficient optimisation, our
electromagnetic compliant model preserves accuracy for arbitrary RIS unit cell
geometry. A numerical analysis of the the simplified communication model has
been carried out for two canonical unit cell geometries, and validated by the
(related) exact MoM code. The performance of a compact RIS-assisted MIMO link
are also assessed in terms of channel capacity.

The paper is organised as follows, in Sec. II we apply the method-of-moment
formulation to the communication system including the RIS, transmit and
receive arrays. We partition basis functions into port and passive segments.
We then devise an end-to-end communication model, and devise a reduced
cascaded model. In Sec. III we design two RIS structures based on loaded
scatterers, and we assess the model by numerical validations of the channel
capacity at variable receiver distance and direction. Finally, in Sec. IV we
draw a conclusion and future perspectives.

##  2 Formulation

###  2.1 Full-wave model

Figure 1: Equivalent circuit of RIS including transmitting/receiving antennas.

Here, an equivalent circuit of RIS including transmitting/receiving antennas
is modeled via method of moments (MoM) as shown in Fig. 1. The equivalent
circuit is expressed as a matrix equation as follows.

$({\mathbf{Z}}_{N\times N}+{\mathbf{Z}}_{N\times
N}^{L}){\mathbf{I}}_{N}={\mathbf{V}}_{N}$ (1)

Here, ${\mathbf{Z}}_{N\times N}$ is $N\times N$ impedance matrix,
${\mathbf{Z}}_{N\times N}^{L}$ is $N\times N$ load impedance matrix,
${\mathbf{I}}_{N}$ is $N$ dimensional current vector, and ${\mathbf{V}}_{N}$
is $N$ dimensional voltage vector, and $N$ is the total number of segments
(i.e. number of unknowns).

Equation (1) can be partitioned using block matrices/vectors as follows.

$\displaystyle\left[{\begin{array}[]{*{20}c}{\bf{V}}_{N_{T}}\\\ 0\\\ 0\\\
\end{array}}\right]$
$\displaystyle=\left[{\begin{array}[]{*{20}c}{{\bf{Z}}_{N_{T}\times
N_{T}}+{\bf{Z}}_{N_{T}\times N_{T}}^{L}}&{{\bf{Z}}_{N_{T}\times N_{R}}}\\\
{{\bf{Z}}_{N_{R}\times N_{T}}}&{{\bf{Z}}_{N_{R}\times
N_{R}}+{\bf{Z}}_{N_{R}\times N_{R}}^{L}}\\\ {{\bf{Z}}_{N_{S}\times
N_{T}}}&{{\bf{Z}}_{N_{S}\times N_{R}}}\\\ \end{array}}\right.$ (8)
$\displaystyle\left.{\begin{array}[]{*{20}c}{{\bf{Z}}_{N_{T}\times N_{S}}}\\\
{{\bf{Z}}_{N_{R}\times N_{S}}}\\\ {{\bf{Z}}_{N_{S}\times
N_{S}}+{\bf{Z}}_{N_{S}\times N_{S}}^{L}}\\\
\end{array}}\right]\left[{\begin{array}[]{*{20}c}{{\bf{I}}_{N_{T}}}\\\
{{\bf{I}}_{N_{R}}}\\\ {{\bf{I}}_{N_{S}}}\\\ \end{array}}\right]$ (15)

Here, $N_{T}$ is the number of segments of transmitting antennas, $N_{R}$ is
the number of segments of receiving antennas, and $N_{S}$ is the number of
segments of scatterers with variable load impedance (i.e.
$N_{T}+N_{R}+N_{S}=N$). ${\mathbf{Z}}_{i\times j}$ is $i\times j$ block
impedance matrix, where $i=N_{T},N_{R},N_{S}$ and $j=N_{T},N_{R},N_{S}$.
${\mathbf{Z}}_{N_{T}\times N_{T}}^{L}$ is a $N_{T}\times N_{T}$ diagonal
matrix corresponding to transmitting antennas and includes their load
impedance connected to the ports. ${\mathbf{Z}}_{N_{R}\times N_{R}}^{L}$ is a
$N_{R}\times N_{R}$ diagonal matrix corresponding to receiving antennas and
includes their load impedance connected to the ports.
${\mathbf{Z}}_{N_{S}\times N_{S}}^{L}$ is a $N_{S}\times N_{S}$ diagonal
matrix corresponding to scatterers with variable load impedance and includes
their load impedance connected to the ports. ${\mathbf{I}}_{N_{T}}$,
${\mathbf{I}}_{N_{R}}$, and ${\mathbf{I}}_{N_{S}}$ are $N_{T}$, $N_{R}$ and
$N_{S}$ dimensional block current vectors. ${\mathbf{V}}_{N_{T}}$ is a $N_{T}$
dimensional block voltage vector corresponding to the segments of the
transmitting antennas. A voltage-voltage MIMO channel matrix, compliant to
that used in wireless system optimization [19, 20], is derived from the
transfer matrix partitioning introduced in (15), which is based on [21], and
described in Appendix A. It configures as

$\displaystyle{\bf{V}}_{N_{R}}$ $\displaystyle={\bf{Z}}_{N_{R}\times
N_{R}}^{L}{\bf{I}}_{N_{R}}$ (16) $\displaystyle={\bf{H}}^{f}_{N_{R}\times
N_{T}}{\bf{V}}_{N_{T}}$ (17)

with

$\displaystyle{\bf{H}}^{f}_{N_{R}\times N_{T}}$
$\displaystyle={\bf{Z}}_{N_{R}\times N_{R}}^{L}({\bf{U}}_{N_{R}\times
N_{R}}-{\bf{Y}}_{N_{R}\times N_{S}}{\bf{Z}}_{N_{S}\times N_{S}}^{L}$
$\displaystyle\quad({\bf{U}}_{N_{S}\times N_{S}}+{\bf{Y}}_{N_{S}\times
N_{S}}{\bf{Z}}_{N_{S}\times N_{S}}^{L})^{-1}{\bf{Y}}_{N_{S}\times N_{R}}$
$\displaystyle\quad{\bf{Z}}_{N_{R}\times
N_{R}}^{L})^{-1}({\bf{U}}_{N_{R}\times N_{R}}+{\bf{Y}}_{N_{R}\times
N_{R}}{\bf{Z}}_{N_{R}\times N_{R}}^{L})^{-1}$
$\displaystyle\quad({\bf{Y}}_{N_{R}\times N_{T}}-{\bf{Y}}_{N_{R}\times
N_{S}}{\bf{Z}}_{N_{S}\times N_{S}}^{L}({\bf{U}}_{N_{S}\times N_{S}}$
$\displaystyle\quad+{\bf{Y}}_{N_{S}\times N_{S}}{\bf{Z}}_{N_{S}\times
N_{S}}^{L})^{-1}{\bf{Y}}_{N_{S}\times N_{T}})$ (18)

It is worth noticing that ${\bf{V}}_{N_{R}}$ is a receiving voltage vector and
entries of ${\bf{H}}^{f}_{N_{R}\times N_{T}}$ \- corresponding to the ports of
the transmitting/receiving antennas - are extracted to obtain complex-valued
channel gains [22]. The channel matrix obtained from (18) is a full-wave
expression as no approximation is introduced through its derivation. Finally,
we have verified that the channel matrix obtained from (18) perfectly agrees
with that obtained by full-wave analysis. Interestingly, the impedance
formulation in (18) does not configure as a cascaded model and it presents
inverse operators that cannot be expanded in Neumann series since their
spectral radius is not smaller than the unity. Hence, we support the
observation done in [23] concerning the mathematical structure of the
scattering formulation, whose spectral radius is smaller than the unity thus
allowing for a Born series type expansion, while the impedance formulation
does not allow for such an expansion.

###  2.2 Reduced model

Here, it is assumed that effect of mutual coupling among the transmitting
antennas, receiving antennas, and the RIS is negligible when their current
distribution is obtained (The assumption is reasonable when they are in far-
field region.). According to the assumption, current distribution of the
transmitting antennas are approximately obtained from Eq. (15) as follows.

$\displaystyle{\bf{V}}_{N_{T}}$ $\displaystyle=({\bf{Z}}_{N_{T}\times
N_{T}}+{\bf{Z}}_{N_{T}\times N_{T}}^{L}){\bf{I}}_{N_{T}}+{\bf{Z}}_{N_{T}\times
N_{R}}{\bf{I}}_{N_{R}}$ $\displaystyle\quad+{\bf{Z}}_{N_{T}\times
N_{S}}{\bf{I}}_{N_{R}}$ $\displaystyle\approx({\bf{Z}}_{N_{T}\times
N_{T}}+{\bf{Z}}_{N_{T}\times N_{T}}^{L}){\bf{I}}_{N_{T}}$
$\displaystyle{\bf{I}}_{N_{T}}$ $\displaystyle\approx({\bf{Z}}_{N_{T}\times
N_{T}}+{\bf{Z}}_{N_{T}\times N_{T}}^{L})^{-1}{\bf{V}}_{N_{T}}$ (19)

Next, current distribution of the RIS is obtained. According to the third
column of Eq. (15), current distribution of RIS is approximately obtained as
follows.

$\displaystyle 0$ $\displaystyle={\bf{Z}}_{N_{S}\times
N_{T}}{\bf{I}}_{N_{T}}+{\bf{Z}}_{N_{S}\times
N_{R}}{\bf{I}}_{N_{R}}+({\bf{Z}}_{N_{S}\times N_{S}}+{\bf{Z}}_{N_{S}\times
N_{S}}^{L}){\bf{I}}_{N_{S}}$ $\displaystyle\approx{\bf{Z}}_{N_{S}\times
N_{T}}{\bf{I}}_{N_{T}}+({\bf{Z}}_{N_{S}\times N_{S}}+{\bf{Z}}_{N_{S}\times
N_{S}}^{L}){\bf{I}}_{N_{S}}$ $\displaystyle{\bf{I}}_{N_{S}}$
$\displaystyle=-({\bf{Z}}_{N_{S}\times N_{S}}+{\bf{Z}}_{N_{S}\times
N_{S}}^{L})^{-1}{\bf{Z}}_{N_{s}\times N_{T}}{\bf{I}}_{N_{T}}$
$\displaystyle=-({\bf{Z}}_{N_{S}\times N_{S}}+{\bf{Z}}_{N_{S}\times
N_{S}}^{L})^{-1}{\bf{Z}}_{N_{s}\times N_{T}}$
$\displaystyle\quad({\bf{Z}}_{N_{T}\times N_{T}}+{\bf{Z}}_{N_{T}\times
N_{T}}^{L})^{-1}{\bf{V}}_{N_{T}}\quad(\because{\rm Eq}.(\ref{eq:Matrix-3}))$
(20)

After that, current distribution of the receiving antennas is approximately
obtained from the second column of Eq. (15) with the help of Eqs. (19) and
(20).

$\displaystyle 0$ $\displaystyle={\bf{Z}}_{N_{R}\times
N_{T}}{\bf{I}}_{N_{T}}+({\bf{Z}}_{N_{R}\times N_{R}}+{\bf{Z}}_{N_{R}\times
N_{R}}^{L}){\bf{I}}_{N_{R}}$ $\displaystyle\quad+{\bf{Z}}_{N_{R}\times
N_{S}}{\bf{I}}_{N_{S}}$ $\displaystyle{\bf{I}}_{N_{R}}$
$\displaystyle=-({\bf{Z}}_{N_{R}\times N_{R}}+{\bf{Z}}_{N_{R}\times
N_{R}}^{L})^{-1}$ $\displaystyle\quad({\bf{Z}}_{N_{R}\times
N_{T}}{\bf{I}}_{N_{T}}+{\bf{Z}}_{N_{R}\times N_{S}}{\bf{I}}_{N_{S}})$
$\displaystyle=-({\bf{Z}}_{N_{R}\times N_{R}}+{\bf{Z}}_{N_{R}\times
N_{R}}^{L})^{-1}({\bf{Z}}_{N_{R}\times N_{T}}({\bf{Z}}_{N_{T}\times N_{T}}$
$\displaystyle\quad+{\bf{Z}}_{N_{T}\times
N_{T}}^{L})^{-1}{\bf{V}}_{N_{T}}-{\bf{Z}}_{N_{R}\times
N_{S}}({\bf{Z}}_{N_{S}\times N_{S}}+{\bf{Z}}_{N_{S}\times N_{S}}^{L})^{-1}$
$\displaystyle{\bf{Z}}_{N_{s}\times N_{T}}({\bf{Z}}_{N_{T}\times
N_{T}}+{\bf{Z}}_{N_{T}\times N_{T}}^{L})^{-1}{\bf{V}}_{N_{T}})$
$\displaystyle=({\bf{Z}}_{N_{R}\times N_{R}}+{\bf{Z}}_{N_{R}\times
N_{R}}^{L})^{-1}(-{\bf{Z}}_{N_{R}\times N_{T}}+{\bf{Z}}_{N_{R}\times N_{S}}$
$\displaystyle\quad({\bf{Z}}_{N_{S}\times N_{S}}+{\bf{Z}}_{N_{S}\times
N_{S}}^{L})^{-1}{\bf{Z}}_{N_{s}\times N_{T}})$
$\displaystyle\quad\times({\bf{Z}}_{N_{T}\times N_{T}}+{\bf{Z}}_{N_{T}\times
N_{T}}^{L})^{-1}{\bf{V}}_{N_{T}}$ (21)

The receiving voltage vector ${\bf V}_{N_{R}}$ of the receiving antennas is
expressed as follows.

${\bf{V}}_{N_{R}}={\bf{Z}}_{N_{R}\times
N_{R}}^{L}{\bf{I}}_{N_{R}}={\bf{H}}_{N_{R}\times N_{T}}{\bf{V}}_{N_{T}}$ (22)

where the voltage-voltage MIMO channel matrix is obtained as

$\displaystyle{\bf{H}}_{N_{R}\times N_{T}}$
$\displaystyle={\bf{Z}}_{N_{R}\times N_{R}}^{L}({\bf{Z}}_{N_{R}\times
N_{R}}+{\bf{Z}}_{N_{R}\times N_{R}}^{L})^{-1}$
$\displaystyle\quad\times(-{\bf{Z}}_{N_{R}\times N_{T}}+{\bf{Z}}_{N_{R}\times
N_{S}}\,\,\Phi_{S}\,\,{\bf{Z}}_{N_{s}\times N_{T}})$
$\displaystyle\quad\times({\bf{Z}}_{N_{T}\times N_{T}}+{\bf{Z}}_{N_{T}\times
N_{T}}^{L})^{-1}{\bf{V}}_{N_{T}}$ (23)

and the RIS state matrix is given by the admittance matrix

$\Phi_{S}=({\bf{Z}}_{N_{S}\times N_{S}}+{\bf{Z}}_{N_{S}\times
N_{S}}^{L})^{-1}$ (24)

which is consistent with the communication model in [15] in the far-field
approximation, that has been used in sum-rate optimization [20]. Because the
block load impedance matrix ${\bf{Z}}_{N_{R}\times N_{R}}^{L}$ is a diagonal
matrix, its entries are zero except for the ports of the receiving antennas.
Finally, entries of ${\bf{H}}_{N_{R}\times N_{T}}$ corresponding to the ports
of the transmitting/receiving antennas are extracted and the channel matrix is
obtained.

###  2.3 Structural scattering

In NLOS between transmit and receive array, the cascaded communication model
in [15] vanishes if all the dipole ports are in open circuit, on account of
the minimum scattering assumption [29]. Importantly, the MoM based formulation
leads to a non-zero RIS state matrix when unit cells are unloaded. The reason
lies in the presence of passive segments that capture the structural
scattering of the bare meta-surface structure. This can be seen by considering
the active impedance of a single unit cell of the RIS. Suppose we have a
linear antenna discretised by three adjacent segments, of which only one
embeds port terminals. Adopting the MoM formalism, the RIS state matrix for
the three segments is found from

${\bf{Z}}_{N_{S}\times N_{S}}=\begin{pmatrix}Z&Z_{c}&Z_{c}\\\
Z_{c}&Z_{p}&Z_{c}\\\ Z_{c}&Z_{c}&Z\end{pmatrix}$ (25)

where, without loss of generality, mutual coupling impedance $Z_{c}$ has been
assumed uniform within the macro-pixel, with input impedance at the port
segment $Z_{p}$ and at the passive segments $Z$. The load impedance at the
port segment reads

${\bf{Z}}_{N_{S}\times N_{S}}^{L}=\begin{pmatrix}0&0&0\\\ 0&Z_{L}&0\\\
0&0&0\end{pmatrix}$ (26)

Inserting (25) and (26) in (24) yields the matrix inversion

$\Phi_{S}={\begin{pmatrix}Z&Z_{c}&Z_{c}\\\ Z_{c}&Z_{p}+Z_{L}&Z_{c}\\\
Z_{c}&Z_{c}&Z\end{pmatrix}}^{-1}$ (27)

which gives, in the open circuit ($oc$) case, i.e., $Z_{L}\rightarrow\infty$

$\Phi_{S}^{oc}=\begin{pmatrix}\frac{Z}{Z^{2}-Z_{c}^{2}}&0&\frac{Z_{c}}{Z^{2}-Z_{c}^{2}}\\\
0&0&0\\\
\frac{-Z_{c}}{Z^{2}-Z_{c}^{2}}&0&\frac{Z}{Z^{2}-Z_{c}^{2}}\end{pmatrix}$ (28)

that encodes the passive scattering of the planar structure. Note that
differently from [15], the RIS state admittance does not go to zero for
unloaded unit cells, which was a consequence of the CMS approximation.

##  3 Numerical results

###  3.1 Design of RIS

Figure 2: Geometry of RIS Figure 3: Unitcell model of infinite array of RIS
elements over PEC ground plane (Left: Planar dipole elements, Right:
Rectangular loop elements).

Here, it is demonstrated that our proposed expression works for calculating
channel capacity between transmitting and receiving antennas in the presence
of RIS. Geometry of the RIS is shown in Fig. 2. Here, either of planar dipole
element or planar loop element is used as RIS element. Phase of scattering
field of the RIS elements is tunable by variable load impedances. Variable
load impedances loaded with RIS elements are simply modeled as open or short
terminations, i.e. the RIS is a so-called single bit RIS. Transmitting and
receiving antennas are planar dipole array antennas. Origin of the coordinate
system is correspnding to the center of the RIS and a position vector pointing
the center of the $mn$th RIS element is ${\bf r}_{mn}$. A position vector
pointing the center of the transmitting antenna is ${\bf r}_{t}$ and that
pointing the center of the receiving antenna is ${\bf r}_{r}$. All of the
antenna elements and RIS elements are discretized into triangular segments and
their current distribution is expanded using Rao-Wilton-Glisson (RWG) basis
function [30].

In advance of designing the RIS, phase of reflection coefficient of the RIS
elements is obtained numerically. Here, method of moments based on periodic
Green’s function (PGF) is used for numerical analysis of the RIS elements in
an infinite periodic array [31], [32]. Poor convergence of the PGF is enhanced
by introducing Ewald transformation with the optimum splitting parameter
[33]-[35]. Singularity at a point where the source/observation points are is
overlapping is annihilated using L’Hospital rule [36], [37].

Figure 4: On/off states of RIS elements composed of dipole elements.  Figure
5: On/off states of RIS elements composed of loop elements.  Figure 6:
Scattering patterns of RIS composed of dipole elements.  Figure 7: Scattering
patterns of RIS composed of loop elements.

The geometry of unit cell model composed of planar dipole elements and planar
loop elements is shown in Fig. 3. Phase of reflection coefficient of the unit
cell model is obtained for plane wave of normal incidence and is shown in
Table 1. According to the Table 1, it is found that phase shift of 180 deg. is
approximately available by switching variable load impedances from on to off,
and vice versa. Therefore, it can be said that the planar dipole elements and
planar loop elements with variable load impedances are applicable to a single-
bit RIS.

The single-bit RIS is designed using the planar dipole elements and planar
loop elements. Based on the phase of reflection coefficient, on/off states of
the RIS elements are designed so that the scattering fields of the RIS
elements are superposed with the minimum error at specific main beam
direction. Here, dimensions of the RIS and transmitting antennas are shown in
TABLE 2. The RIS is designed so that its main beam is directed to
$\theta=0,30^{\circ}$. The RIS is analysed using the MoM with the layered
media Green’s function(LMGF) and infinite ground plane of the RIS is modeled
using the LMGF [38]. A direct wave component is extracted from the LMGF and
calculated in the spatial domain so that the convergence of the LMGF is
enhanced [39], [40]. Singularity at a point where the source/observation
points are is overlapping is annihilated using coordinate transformation and
analytic integral [41], [42]. So-called Sommerfeld integral in the LMGF is
performed efficiently via numerical interpolation by using Taylor expansions
[43], [44].

On/off states of the RIS elements are shown in Figs. 4 and 5. It is found that
on/off states of the RIS elements are switched corresponding to the main beam
direction or position of the transmitting antenna. The typical phase mask of a
lens emerges upon RIS unit cell optimization in both dipole and loop based
reflecting sufaces. Normalised scattering patterns of the RISs composed of the
planar dipole elements and the planar loop elements are shown in Figs. 6 and
7, respectively. It is found that main beam of the RIS is directed to the
specific angle because on/off states of the RIS elements are switched so that
their scattering fields are with the minimum error at the direction.

Table 1: Phase of reflection coefficient of planar dipole and planar loop elements. | Planar dipole | Planar loop
---|---|---
Type of elements | ($l=0.5\lambda,w=0.01\lambda,$ | ($a=b=0.4\lambda,w=0.01\lambda,$
and their dimensions | $h=0.25\lambda,d_{x}=0.5\lambda,$ | $h=0.25\lambda,d_{x}=d_{y}=0.6\lambda)$
| $d_{y}=0.7\lambda)$ |
Phase of | On state | 6.4 | -56
reflection | (Short termination) |  |
coefficient [deg.] | Off state | 164.6 | 122.8
| (Open termination) |  |
Table 2: Dimensions of the RIS and transmitting/receiving antennas.  RIS | Tx antenna | Rx antenna
---|---|---
Planar dipole elements | Planar loop elements |  |
$l$ | $0.5\lambda$ | $a=b$ | $0.4\lambda$ | $l$ | $0.5\lambda$ | $l$ | $0.5\lambda$
$w$ | $0.01\lambda$ | $w$ | $0.01\lambda$ | $w$ | $0.01\lambda$ | $w$ | $0.01\lambda$
$d_{x}$ | $0.5\lambda$ | $d_{x}$ | $0.6\lambda$ | $d_{xt}$ | $0.5\lambda$ | $d_{xr}$ | $0.5\lambda$
$d_{y}$ | $0.7\lambda$ | $d_{y}$ | $0.6\lambda$ | $M_{xt}$ | $2$ | $M_{xr}$ | $2$
$h$ | $0.25\lambda$ | $h$ | $0.25\lambda$ | ${\bf r}_{t}$ | $(0,0,5\lambda)$ | ${\bf r}_{r}$ | $(R\cos\theta,0,R\sin\theta)$
$M_{x}=M_{y}$ | $11$ | $M_{x}=M_{y}$ | $11$ |  | $(0,0,10\lambda)$ |  | $\theta=0,30^{\circ}$

###  3.2 Channel capacity

Here, it is demonstrated that our proposed expression of the channel capacity
is applicable to a wireless communication system including RIS. Dimensions of
the transmitting/receiving antennas and the RIS are shown in Table 2. Channel
capacity is expressed as follows.

$C=\log_{2}\left|{{\bf{I}}+{\bf{HH}}^{\dagger}{\frac{\gamma}{M_{xt}}}}\right|$
(29)

Here, ${\bf{I}}$ is an unit matrix, ${\bf H}$ is a $M_{tr}\times M_{xr}$
channel matrix extracted from ${\bf{H}}_{N_{R}\times N_{T}}$, and $\gamma$ is
signal-to-noise ratio (i.e. Ratio between transmitting power and noise power
at receiving antennas. )

Simulated channel capacity is shown in Figs. 8 and 9 for the RIS composed of
the planar dipole array antennas. Also, simulated channel capacity is shown in
Figs. 10 and 11 for the RIS composed of the planar loop array antennas. It is
shown that the channel capacity obtained using the proposed expression agrees
well with that obtained using the full-wave analysis except for small error
when the receiving antennas are too close to the RIS (e.g. $R=1\lambda$). This
is predicted at large surface-receiver distances $R$ and in absence of LOS,
for which

${\bf{Z}}_{N_{R}\times N_{S}}\sim\frac{\textbf{a}_{N_{R}\times N_{S}}}{R},$
(30)

as explained in [45], and where $\textbf{a}_{N_{R}\times N_{S}}$ is the array
steering vector and upon substitution in (23) and (29) yields

$C\sim\left|{{\bf{HH}}^{\dagger}{\frac{\gamma}{M_{xt}}}}\right|\sim\frac{1}{R^{2}}.$
(31)

Importantly, this is consistent (and provides a link) with path-loss models
(see for example [46, 47]), and provides an assessment of the impedance model
for arbitrary RIS structures. The small error of the channel capacity obtained
using the proposed expression stems from the effect of the mutual coupling
among transmitting antennas, receiving antennas and the RIS that is neglected
in its formulation. However, it can be said that applicability of the proposed
expression is kept because the receiving antennas are usually far from the RIS
in practical system.

On the other hand, numerical results should be further discussed from the
physical point of view. The receiving antennas receive direct wave from the
transmitting antenna and scattering wave from the RIS simultaneously.
Therefore, the channel capacity reflects the result of the superposition of
the direct/scattering waves. For example, as shown in the channel capacity
corresponding to the RIS whose angle of main beam is directed to $\theta=0$,
gradient of the channel capacity between the transmitting/receiving antennas
is non-constant when the receiving antennas receive both of the
direct/scattering waves. The non-constant gradient of the channel capacity
comes from phase difference between the direct/scattering waves. On the other
hand, as shown in the channel capacity corresponding to the RIS whose angle of
main beam is directed to $\theta=30^{\circ}$, the gradient of the channel
capacity between the transmitting/receiving antennas is proportional to
$\frac{1}{R^{2}}$ when the $R$ becomes large. The $\frac{1}{R^{2}}$ gradient
of the channel capacity is found because direct wave from the transmitting
antennas is negligible and $R$ approaches to boundary of far-field region of
the RIS (e.g. $\frac{2D^{2}}{\lambda}=148\lambda$, where $D$ is diagonal line
of the RIS composed of the planar dipole elements). According to the numerical
results and discussions, it can be said that the proposed expression is
applicable to calculating the channel capacity between the
transmitting/receiving antennas in the presence of the RIS.

Figure 8: Channel capacity of RIS composed of planar dipole elements (${\bf
r}_{t}=(0,0,5\lambda)$).  Figure 9: Channel capacity of RIS composed of planar
dipole elements (${\bf r}_{t}=(0,0,10\lambda)$).  Figure 10: Channel capacity
of RIS composed of planar loop elements (${\bf r}_{t}=(0,0,5\lambda)$).
Figure 11: Channel capacity of RIS composed of planar loop elements (${\bf
r}_{t}=(0,0,10\lambda)$).

##  4 Conclusion

We have obtained a cascaded communication model for RIS-assisted wireless
links from a first-principle electromagnetic model. The formulation has been
achieved by partitioning the interaction matrices of the method-of-moments. We
have found that the model configures as a cascaded product of interaction
matrices in the far-field of the RIS unit cell. This extends communication
models based on impedance matrices that has been originally derived under the
canonical minimum scattering approximation. The system-level performance of
RIS-assisted end-to-end MIMO communications can be accurately simulated with
Method of Moments based channel matrices. Finite-size RIS unit cells with
arbitrary geometry and dimensions can be characterised and subsequently
optimised to achieve prescribed electromagnetic shaping functions. Future work
is required to devise computationally efficient methods for multilayer and
transmissive RIS structures.

APPENDIX A: Full-wave Expression of Channel Matrix
Here, full-wave expression of a channel matrix is derived in the same manner
with [27]. Starting from Eq. (15), following expression is obtained.

$\displaystyle\left[{\begin{array}[]{*{20}c}{\bf{V}}_{N_{T}}\\\
-{\bf{Z}}_{N_{R}\times N_{R}}^{L}{\bf{I}}_{N_{R}}\\\ -{\bf{Z}}_{N_{S}\times
N_{S}}^{L}{\bf{I}}_{N_{S}}\\\ \end{array}}\right]$
$\displaystyle=\left[{\begin{array}[]{*{20}c}{{\bf{Z}}_{N_{T}\times
N_{T}}+{\bf{Z}}_{N_{T}\times N_{T}}^{L}}\\\ {{\bf{Z}}_{N_{R}\times
N_{T}}}&{{\bf{Z}}_{N_{R}\times N_{R}}}\\\ {{\bf{Z}}_{N_{S}\times
N_{T}}}&{{\bf{Z}}_{N_{S}\times N_{R}}}\\\ \end{array}}\right.$ (38)
$\displaystyle\quad\left.{\begin{array}[]{*{20}c}{{\bf{Z}}_{N_{T}\times
N_{S}}}\\\ {{\bf{Z}}_{N_{R}\times N_{S}}}\\\ {{\bf{Z}}_{N_{S}\times N_{S}}}\\\
\end{array}}\right]\left[{\begin{array}[]{*{20}c}{{\bf{I}}_{N_{T}}}\\\
{{\bf{I}}_{N_{R}}}\\\ {{\bf{I}}_{N_{S}}}\\\ \end{array}}\right]$ (45)
$\displaystyle{\bf V}_{N}^{\prime}$ $\displaystyle={\bf Z}_{N\times
N}^{\prime}{\bf I}_{N}$ (46)

Here, ${\bf V}_{N}^{\prime}\equiv[{\bf{V}}_{N_{T}},-{\bf{Z}}_{N_{R}\times
N_{R}}^{L}{\bf{I}}_{N_{R}},-{\bf{Z}}_{N_{S}\times
N_{S}}^{L}{\bf{I}}_{N_{S}}]^{T}$ and ${\bf Z}_{N\times N}^{\prime}$ is an
$N\times N$ impedance matrix where the block load impedance matrix at the
transmitting antenna is embedded. Eq. (46) can be solved by multiplying
${\mathbf{Y}}_{N\times N}\equiv\\{{\bf Z}_{N\times N}^{\prime}\\}^{-1}$ with
its both sides and following expression is obtained.

$\displaystyle{\mathbf{I}}_{N}$ $\displaystyle={\mathbf{Y}}_{N\times
N}{\mathbf{V}}_{N}^{\prime}$
$\displaystyle\left[{\begin{array}[]{*{20}c}{{\mathbf{I}}_{N_{T}}}\\\
{{\mathbf{I}}_{N_{R}}}\\\ {{\mathbf{I}}_{N_{S}}}\\\ \end{array}}\right]$
$\displaystyle=\left[{\begin{array}[]{*{20}c}{{\mathbf{Y}}_{N_{T}\times
N_{T}}}&{{\mathbf{Y}}_{N_{T}\times N_{R}}}\\\ {{\mathbf{Y}}_{N_{R}\times
N_{T}}}&{{\mathbf{Y}}_{N_{R}\times N_{R}}}\\\ {{\mathbf{Y}}_{N_{S}\times
N_{T}}}&{{\mathbf{Y}}_{N_{S}\times N_{R}}}\\\ \end{array}}\right.$ (53)
$\displaystyle\quad\left.{\begin{array}[]{*{20}c}{{\mathbf{Y}}_{N_{T}\times
N_{S}}}\\\ {{\mathbf{Y}}_{N_{R}\times N_{S}}}\\\ {{\mathbf{Y}}_{N_{S}\times
N_{S}}}\\\
\end{array}}\right]\left[{\begin{array}[]{*{20}c}{{\mathbf{V}}_{N_{T}}}\\\
{-{\mathbf{Z}}_{N_{R}\times N_{R}}^{L}{\mathbf{I}}_{N_{R}}}\\\
{-{\mathbf{Z}}_{N_{S}\times N_{S}}^{L}{\mathbf{I}}_{N_{S}}}\\\
\end{array}}\right]$ (60)

According to Eq. (60), following three block matrix equations are obtained.

$\displaystyle{}{\mathbf{I}}_{N_{T}}={\mathbf{Y}}_{N_{T}\times
N_{T}}{\mathbf{V}}_{N_{T}}-{\mathbf{Y}}_{N_{T}\times
N_{R}}{\mathbf{Z}}_{N_{R}\times N_{R}}^{L}{\mathbf{I}}_{N_{R}}$
$\displaystyle\quad\quad\quad\quad-{\mathbf{Y}}_{N_{T}\times
N_{S}}{\mathbf{Z}}_{N_{S}\times N_{S}}^{L}{\mathbf{I}}_{N_{S}}$ (61)
$\displaystyle{\mathbf{I}}_{N_{R}}={\mathbf{Y}}_{N_{R}\times
N_{T}}{\mathbf{V}}_{N_{T}}-{\mathbf{Y}}_{N_{R}\times
N_{R}}{\mathbf{Z}}_{N_{R}\times N_{R}}^{L}{\mathbf{I}}_{N_{R}}$
$\displaystyle\quad\quad\quad\quad-{\mathbf{Y}}_{N_{R}\times
N_{S}}{\mathbf{Z}}_{N_{S}\times N_{S}}^{L}{\mathbf{I}}_{N_{S}}$ (62)
$\displaystyle{\mathbf{I}}_{N_{S}}={\mathbf{Y}}_{N_{S}\times
N_{T}}{\mathbf{V}}_{N_{T}}-{\mathbf{Y}}_{N_{S}\times
N_{R}}{\mathbf{Z}}_{N_{R}\times N_{R}}^{L}{\mathbf{I}}_{N_{R}}$
$\displaystyle\quad\quad\quad\quad-{\mathbf{Y}}_{N_{S}\times
N_{S}}{\mathbf{Z}}_{N_{S}\times N_{S}}^{L}{\mathbf{I}}_{N_{S}}$ (63)

(63) can be transformed as follows.

$\displaystyle{\bf{I}}_{N_{S}}$ $\displaystyle=({\bf{U}}_{N_{S}\times
N_{S}}+{\bf{Y}}_{N_{S}\times N_{S}}{\bf{Z}}_{N_{S}\times N_{S}}^{L})^{-1}$
$\displaystyle\quad({\bf{Y}}_{N_{S}\times
N_{T}}{\bf{V}}_{N_{T}}-{\bf{Y}}_{N_{S}\times N_{R}}{\bf{Z}}_{N_{R}\times
N_{R}}^{L}{\bf{I}}_{N_{R}})$ (64)

Here, ${\mathbf{U}}_{N_{S}\times N_{S}}$ is $N_{S}\times N_{S}$ unit matrix.
The next, Eq. (62) can be transformed as follows.

$\displaystyle{\bf{I}}_{N_{R}}$ $\displaystyle=({\bf{U}}_{N_{R}\times
N_{R}}+{\bf{Y}}_{N_{R}\times N_{R}}{\bf{Z}}_{N_{R}\times N_{R}}^{L})^{-1}$
$\displaystyle\quad({\bf{Y}}_{N_{R}\times
N_{T}}{\bf{V}}_{N_{T}}-{\bf{Y}}_{N_{R}\times N_{S}}{\bf{Z}}_{N_{S}\times
N_{S}}^{L}{\bf{I}}_{N_{S}})$ (65)

Equation (64) is substituted into Eq. (65) and following expression is
obtained.

$\displaystyle{\bf{I}}_{N_{R}}$ $\displaystyle=({\bf{U}}_{N_{R}\times
N_{R}}+{\bf{Y}}_{N_{R}\times N_{R}}{\bf{Z}}_{N_{R}\times
N_{R}}^{L})^{-1}({\bf{Y}}_{N_{R}\times N_{T}}{\bf{V}}_{N_{T}}$
$\displaystyle\quad-{\bf{Y}}_{N_{R}\times N_{S}}{\bf{Z}}_{N_{S}\times
N_{S}}^{L}({\bf{U}}_{N_{S}\times N_{S}}+{\bf{Y}}_{N_{S}\times
N_{S}}{\bf{Z}}_{N_{S}\times N_{S}}^{L})^{-1}$
$\displaystyle\quad({\bf{Y}}_{N_{S}\times
N_{T}}{\bf{V}}_{N_{T}}-{\bf{Y}}_{N_{S}\times N_{R}}{\bf{Z}}_{N_{R}\times
N_{R}}^{L}{\bf{I}}_{N_{R}}))$ (66) $\displaystyle=({\bf{U}}_{N_{R}\times
N_{R}}+{\bf{Y}}_{N_{R}\times N_{R}}{\bf{Z}}_{N_{R}\times
N_{R}}^{L})^{-1}\left\\{\right.({\bf{Y}}_{N_{R}\times N_{T}}$
$\displaystyle\quad-{\bf{Y}}_{N_{R}\times N_{S}}{\bf{Z}}_{N_{S}\times
N_{S}}^{L}({\bf{U}}_{N_{S}\times N_{S}}+{\bf{Y}}_{N_{S}\times
N_{S}}{\bf{Z}}_{N_{S}\times N_{S}}^{L})^{-1}$
$\displaystyle\quad{\bf{Y}}_{N_{S}\times
N_{T}}){\bf{V}}_{N_{T}}+{\bf{Y}}_{N_{R}\times N_{S}}{\bf{Z}}_{N_{S}\times
N_{S}}^{L}({\bf{U}}_{N_{S}\times N_{S}}$
$\displaystyle\quad+{\bf{Y}}_{N_{S}\times N_{S}}{\bf{Z}}_{N_{S}\times
N_{S}}^{L})^{-1}{\bf{Y}}_{N_{S}\times N_{R}}{\bf{Z}}_{N_{R}\times
N_{R}}^{L}{\bf{I}}_{N_{R}}\left.\right\\}\\!\\!\\!\\!$ (67)

Equation (67) can be solved for ${\bf{I}}_{N_{R}}$ and ${\bf{I}}_{N_{R}}$ is
obtained as follows.

$\displaystyle{\bf{I}}_{N_{R}}$ $\displaystyle=({\bf{U}}_{N_{R}\times
N_{R}}-{\bf{Y}}_{N_{R}\times N_{S}}{\bf{Z}}_{N_{S}\times
N_{S}}^{L}({\bf{U}}_{N_{S}\times N_{S}}+{\bf{Y}}_{N_{S}\times N_{S}}$
$\displaystyle\quad{\bf{Z}}_{N_{S}\times N_{S}}^{L})^{-1}{\bf{Y}}_{N_{S}\times
N_{R}}{\bf{Z}}_{N_{R}\times N_{R}}^{L})^{-1}({\bf{U}}_{N_{R}\times N_{R}}$
$\displaystyle\quad+{\bf{Y}}_{N_{R}\times N_{R}}{\bf{Z}}_{N_{R}\times
N_{R}}^{L})^{-1}({\bf{Y}}_{N_{R}\times N_{T}}-{\bf{Y}}_{N_{R}\times N_{S}}$
$\displaystyle\quad{\bf{Z}}_{N_{S}\times N_{S}}^{L}({\bf{U}}_{N_{S}\times
N_{S}}+{\bf{Y}}_{N_{S}\times N_{S}}{\bf{Z}}_{N_{S}\times N_{S}}^{L})^{-1}$
$\displaystyle\quad{\bf{Y}}_{N_{S}\times N_{T}}){\bf{V}}_{N_{T}}$ (68)

Finally, ${\bf{Z}}_{N_{R}\times N_{R}}^{L}$ is multiplied by both sides of Eq.
(68), following expressions are obtained. Here, ${\bf{V}}_{N_{R}}$ and
${\bf{H}}^{f}_{N_{R}\times N_{T}}$ are defined as in (17) and (18),
respectively.

## Acknowledgement

KK and QC are supported by the Ministry of Internal Affairs and Communications
of Japan, grant number JPJ000254. ST and GG are supported in part by EU H2020
RISE-6G project, grant number 101017011. GG is also supported by the Royal
Society Industry Fellowship, grant number INF\R2\192066, and by the EPSRC,
grant numbers EP/V048937/1 and EP/X038491/1.

## References

  * [1] E. C. Strinati et al., “Wireless Environment as a Service Enabled by Reconfigurable Intelligent Surfaces: The RISE-6G Perspective,” 2021 Joint European Conference on Networks and Communications & 6G Summit (EuCNC/6G Summit), Porto, Portugal, 2021, pp. 562-567, doi: 10.1109/EuCNC/6GSummit51104.2021.9482474.
  * [2] M. Di Renzo et al., ”Smart Radio Environments Empowered by Reconfigurable Intelligent Surfaces: How It Works, State of Research, and The Road Ahead,” in IEEE Journal on Selected Areas in Communications, vol. 38, no. 11, pp. 2450-2525, Nov. 2020, doi: 10.1109/JSAC.2020.3007211.
  * [3] H. Zhang, N. Shlezinger, F. Guidi, D. Dardari, Y.C. and Eldar, “6G wireless communications: From far-field beam steering to near-field beam focusing,” arXiv preprint arXiv:2203.13035, 2022.
  * [4] G. Bartoli, A. Abrardo, N. Decarli, D. Dardari, and M. Di Renzo, “Spatial Multiplexing in Near Field MIMO Channels with Reconfigurable Intelligent Surfaces.” arXiv preprint arXiv:2212.11057, 2023.
  * [5] X. Pei et al., “RIS-Aided Wireless Communications: Prototyping, Adaptive Beamforming, and Indoor/Outdoor Field Trials,” in IEEE Transactions on Communications, vol. 69, no. 12, pp. 8627-8640, Dec. 2021, doi: 10.1109/TCOMM.2021.3116151.
  * [6] G. Gradoni, M. Di Renzo, A. Diaz-Rubio, S. Tretyakov, C. Caloz, Z. Peng, A. Alu, G. Lerosey, M. Fink, V. Galdi, and T.J., Cui, “Smart radio environments.” arXiv preprint arXiv:2111.08676, 2021.
  * [7] M. Di Renzo, A. Ahmed, A. Zappone, V. Galdi, G. Gradoni, M. Moccia, and G. Castaldi, G., “Digital reconfigurable intelligent surfaces: On the impact of realistic reradiation models,” arXiv preprint arXiv:2205.09799, 2022.
  * [8] M. Di Renzo, F. H. Danufane, S. Tretyakov, ”Communication Models for Reconfigurable Intelligent Surfaces: From Surface Electromagnetics to Wireless Networks Optimization,” Proceedings of the IEEE, Vol. 110, No. 9, September 2022.
  * [9] F. H. Danufane, M. D. Renzo, J. de Rosny and S. Tretyakov, ”On the Path-Loss of Reconfigurable Intelligent Surfaces: An Approach Based on Green’s Theorem Applied to Vector Fields,” in IEEE Transactions on Communications, vol. 69, no. 8, pp. 5573-5592, Aug. 2021.
  * [10] R. J. Williams, P. Ramırez-Espinosa, J. Yuanand, E. de Carvalho, “Electromagnetic Based Communication Model for Dynamic Metasurface Antennas,” IEEE Transactions on Wireless Communications, vol. 21, no. 10, pp. 8616-8630, Oct. 2022.
  * [11] V. Degli-Esposti, E. M. Vitucci, M. D. Renzo and S. A. Tretyakov, ”Reradiation and Scattering From a Reconfigurable Intelligent Surface: A General Macroscopic Model,” in IEEE Transactions on Antennas and Propagation, vol. 70, no. 10, pp. 8691-8706, Oct. 2022.
  * [12] S. V. Hum, M. Okoniewski and R. J. Davies, ”Modeling and Design of Electronically Tunable Reflectarrays,” in IEEE Transactions on Antennas and Propagation, vol. 55, no. 8, pp. 2200-2210, Aug. 2007.
  * [13] R. J. Williams, P. Ramírez-Espinosa, J. Yuan and E. de Carvalho, “Electromagnetic Based Communication Model for Dynamic Metasurface Antennas,” IEEE Transactions on Wireless Communications, vol. 21, no. 10, pp. 8616-8630, Oct. 2022, doi: 10.1109/TWC.2022.3167783.
  * [14] S. Shen, B. Clerckx and R. Murch, R., “Modeling and architecture design of reconfigurable intelligent surfaces using scattering parameter network analysis.” IEEE Transactions on Wireless Communications, vol. 21, no. 2, pp.1229-1243, 2022.
  * [15] 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 Communications Letters, vol. 10, no. 5, pp. 938-942, May 2021, doi: 10.1109/LWC.2021.3050826.
  * [16] M. Akrout, F. Bellili, A. Mezghani, and J.A. Nossek, 2023. “Physically Consistent Models for Intelligent Reflective Surface-assisted Communications under Mutual Coupling and Element Size Constraint,” arXiv preprint arXiv:2302.11130.
  * [17] P. Mursia, S. Phang, V. Sciancalepore, G. Gradoni, and M. Di Renzo, 2023. “Modeling and optimization of reconfigurable intelligent surfaces in propagation environments with scattering objects,” arXiv preprint arXiv:2302.01739.
  * [18] R. Faqiri, C. Saigre-Tardif, G. C. Alexandropoulos, N. Shlezinger, M. F. Imani and P. del Hougne, “PhysFad: Physics-Based End-to-End Channel Modeling of RIS-Parametrized Environments With Adjustable Fading,” IEEE Transactions on Wireless Communications, vol. 22, no. 1, pp. 580-595, Jan. 2023.
  * [19] X. Qian and M. D. Renzo, “Mutual Coupling and Unit Cell Aware Optimization for Reconfigurable Intelligent Surfaces,” IEEE Wireless Communications Letters, vol. 10, no. 6, pp. 1183-1187, June 2021, doi: 10.1109/LWC.2021.3061449.
  * [20] A. Abrardo, D. Dardari, M. Di Renzo and X. Qian, “MIMO Interference Channels Assisted by Reconfigurable Intelligent Surfaces: Mutual Coupling Aware Sum-Rate Optimization Based on a Mutual Impedance Channel Model,” IEEE Wireless Communications Letters, vol. 10, no. 12, pp. 2624-2628, Dec. 2021, doi: 10.1109/LWC.2021.3109017.
  * [21] J. Mautz and R. Harrington, “Modal analysis of loaded N-port scatterers,” in IEEE Transactions on Antennas and Propagation, vol. 21, no. 2, pp. 188-199, March 1973, doi: 10.1109/TAP.1973.1140431.
  * [22] M. Najafi, V. Jamali, R. Schober and H. V. Poor, “Physics-Based Modeling and Scalable Optimization of Large Intelligent Reflecting Surfaces,” IEEE Transactions on Communications, vol. 69, no. 4, pp. 2673-2691, April 2021, doi: 10.1109/TCOMM.2020.3047098.
  * [23] A. Rabault, L.L. Magoarou, J. Sol, G.C. Alexandropoulos, N. Shlezinger, H.V. Poor, and P. del Hougne, 2023. “On the tacit linearity assumption in common cascaded models of RIS-parametrized wireless channels.” arXiv preprint arXiv:2302.04993.
  * [24] M. Di Renzo, V. Galdi, and G. Castaldi, “Modeling the Mutual Coupling of Reconfigurable Metaurfaces,” arXiv preprint arXiv:2210.08619, 2022.
  * [25] S. Hubert, A. Vallecchi, A. Schuchinsky and C. Craeye, “Modeling of impedance-loaded sub-wavelength metasurfaces,” 2017 International Conference on Electromagnetics in Advanced Applications (ICEAA), Verona, Italy, 2017, pp. 1041-1044, doi: 10.1109/ICEAA.2017.8065438.
  * [26] S. V. Hum, M. Okoniewski and R. J. Davies, “Modeling and Design of Electronically Tunable Reflectarrays,” IEEE Transactions on Antennas and Propagation, vol. 55, no. 8, pp. 2200-2210, Aug. 2007, doi: 10.1109/TAP.2007.902002.
  * [27] K. Konno, Q. Yuan, Q. Chen, K. Yokokawa, J. Goto, and T. Fukawasa, “Efficient Method of Moments for Numerical Analysis of Antennas with Variable Load Impedance,” IEEE Trans. Antennas Propag., vol. 68, no. 12, pp. 8233-8237, Dec. 2020.
  * [28] D. Sarkar, and Y. Antar, “Electromagnetic Insights into Path Loss Modelling of IRS-Assisted SISO Links: Method-of-Moment Based Analysis,” Frontiers in Communications and Networks, p.59, 2021.
  * [29] W. K. Kahn and H. Kurss, “Minimum-scattering antennas,” IEEE Transactions on Antennas and Propagation, vol. 13, no. 5, pp. 671-675, September 1965, doi: 10.1109/TAP.1965.1138529.
  * [30] S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antennas Propag., vol. AP-30, no. 3, pp. 409-418, May 1982.
  * [31] M. G. Floquet, “Sur les équations différentielles linéaires à coefficients périodiques,” Annale École Normale Siiperieur, pp. 47-88, 1883.
  * [32] B. A. Munk, Frequency Selective Surfaces: Theory and Design, New York, NY, USA: Wiley, 2000.
  * [33] P. P. Ewald, “Dispersion und Doppelbrechung von Elektronengittern (Kristallen),” Dissertation, München, 1912, also Ann. Phys. 49, p. 1, 1916.
  * [34] P. P. Ewald, “Die Berechnung optischer und elektrostatischer Gitterpotentiale,” Ann. Phys. 64, pp. 253–287, 1921.
  * [35] K. E. Jordan, G. R. Richter, and P. Sheng, “An efficient numerical evaluation of the Green’s function for the Helmholtz operator on periodic structures,” J. Comp. Phys., vol. 63, pp. 222–235, 1986.
  * [36] I. Stevanoviæ  and J. R. Mosig, “Periodic Green’s function for skewed 3-D lattices using the Ewald transformation,” Microw. Opt. Tech. Lett., vol. 49, no. 6, pp. 1353-1357, Jun. 2007.
  * [37] J. Su, X.-W. Xu, M. He, and K. Zhang, “Integral-equation analysis of frequency selective surfaces using Ewald transformation and lattice symmetry” Progress In Electromagnetics Research, vol. 121, pp. 249-269, 2011.
  * [38] W. C. Chew, Waves and Fields in Inhomogeneous Media, IEEE Press, NY 1995.
  * [39] W. C. Chew, J. L. Xiong, and M. A. Saville, “A matrix-friendly formulation of layered medium Green’s function,” IEEE Antennas Wireless Propag. Lett., vol. 5, pp. 490-494, 2006.
  * [40] Y. P. Chen, W. C. Chew, and L. Jiang, “A new Green’s function formulation for modeling homogeneous objects in layered medium,” IEEE Trans. Antennas Propag., vol. 60, no. 10, pp. 4766-4776, Oct. 2012.
  * [41] R. D. Graglia, “On the numerical integration of the linear shape functions times the 3-D Green’s function or its gradient on a plane triangle,” IEEE Trans. Antennas Propag., vol. 41, no. 10, pp. 1448-1455, Oct. 1993.
  * [42] D. R. Wilton, S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak, and C. M. Butler, “Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains,” IEEE Trans. Antennas Propag., vol. 32, no. 3, pp. 276-281, March 1984.
  * [43] 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, 2017.
  * [44] C. A. Balanis, Advanced Engineering Electromagnetics, 2nd ed., John Willey & Sons, pp. 967-979.
  * [45] M. Akrout, V. Shyianov, F. Bellili, A. Mezghani, and R.W. Heath, 2022. “Super-Wideband Massive MIMO,” arXiv preprint arXiv:2208.01556.
  * [46] M. Di Renzo, F. Habibi Danufane, X. Xi, J. de Rosny and S. Tretyakov, ”Analytical Modeling of the Path-Loss for Reconfigurable Intelligent Surfaces – Anomalous Mirror or Scatterer ?,” 2020 IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications (SPAWC), Atlanta, GA, USA, pp. 1-5, 2020.
  * [47] J. C. B. Garcia, A. Sibille and M. Kamoun, ”Reconfigurable Intelligent Surfaces: Bridging the Gap Between Scattering and Reflection,” in IEEE Journal on Selected Areas in Communications, vol. 38, no. 11, pp. 2538-2547, Nov. 2020.