#  On Energy Efficiency and Fairness Maximization in RIS-Assisted MU-MISO
mmWave Communications  ††thanks: This publication has emanated from research
conducted with the financial support of Science Foundation Ireland under Grant
Number 13/RC/2077_P2.

Ahmed Magbool, Vaibhav Kumar, and Mark F. Flanagan
School of Electrical and Electronic Engineering, University College Dublin,
Belfield, Dublin 4, Ireland
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Reconfigurable intelligent surfaces (RISs) are considered to be a promising
solution to overcome the blockage issue in the millimeter-wave (mmWave) band.
Energy efficiency is an important performance metric in mmWave systems with a
large number of antennas. However, due to the severe path loss in mmWave
systems, resource allocation algorithms tend to allocate most of the resources
for the benefit of the users with higher channel gains. In this paper, we
propose a lexicographic-based approach to find the optimal power allocation,
RIS phase shift matrix, and analog precoders that maximize both energy
efficiency and user fairness. We solve the corresponding multi-objective
optimization problem in two stages. In the first stage, we maximize the energy
efficiency, and in the second stage we maximize the fairness subject to a
minimum energy efficiency constraint. We propose an alternating optimization
procedure to solve the optimization problem in each stage. The optimal power
allocation is found using Dinkelbach’s method in the first stage and using
convex optimization techniques in the second stage, the RIS phase shift matrix
is found using a gradient ascent algorithm, and the analog precoder is
determined using beam alignment. Numerical results show that the proposed
algorithm can achieve an excellent trade-off between the energy efficiency and
fairness by boosting the minimum weighted rate with a minor and controllable
reduction in the energy efficiency.

###### Keywords:

Reconfigurable intelligent surfaces, mmWave communications, energy efficiency,
user fairness, lexicographic approach, Dinkelbach’s method.

##  I Introduction

To cope with the ever-increasing demand for better connectivity and higher
data rates, the millimeter-wave (mmWave) band has been investigated as a
spectrum enabler for the existing fifth-generation (5G) and the upcoming
sixth-generation (6G) of cellular wireless communication networks [1]. To
address the high sensitivity to blockages in the mmWave band, the use of
reconfigurable intelligent surfaces (RISs) has been proposed to create a
programmable propagation environment by controlling the phase of the incident
signals using low-cost nearly-passive reflecting elements [2].

Recent research contributions draw more attention to the energy efficiency of
mmWave systems [3], where the spectral efficiency should be enhanced at a
minimal expense in additional power consumption. Two solutions for EE
maximization are proposed in [4]; one is gradient-based and the other uses a
sequential fractional programming based alternating approach. A genetic
approach is proposed in [5] based on a covariance matrix adaptation evolution
strategy and Dinkelbach’s method to maximize the EE in mmWave and terahertz
(THz) systems. The optimal beamformer design is obtained in [6] along with an
asymptotic channel capacity characterization under hardware impairments.
Furthermore, a long short-term memory (LSTM) deep neural network has been
utilized to boost the EE in [7].

Due to the high path loss in mmWave systems, users served by the same base
station (BS) can have significantly different channel strengths. When user
fairness is ignored, EE optimization algorithms tend to give users with very
small channel gains a very low quality of service (QoS), as these users make
only a minor contribution to the EE objective function. One solution to avoid
this is to incorporate QoS constraints into the optimization problem, e.g., by
forcing the rate of each user to be greater than a specific threshold.
However, determining these threshold rates can be problematic as it can make
the optimization problem infeasible. Moreover, the simulation results in [4]
show that setting high thresholds can severely impact the EE of the system. To
avoid this, a multi-objective function is formulated in [8] which combines the
sum rate and Jain’s fairness index; a successive convex approximation method
is proposed to solve the resulting non-convex optimization problem. The
minimum rate in the system is maximized in [9] whereas the minimum EE of the
individual users is targeted in [10].

In this paper, we target the maximization of both EE and fairness in RIS-
assisted multi-user multiple-input-single-output (MU-MISO) mmWave systems.
First, we formulate a multi-objective optimization which focuses on maximizing
the system’s EE as well as the fairness between users, which is expressed in
terms of proportional rates. The former is an important aspect of system
design while the latter is crucial to assuring the QoS of end users. Next, we
use a lexicographic method to divide the optimization problem into two stages.
In the first stage, we maximize the EE only, and in the second stage, we
maximize the fairness subject to an EE constraint derived from the solution
found in the first stage. In each stage, we use an alternating optimization
(AO) framework to find the optimal power allocation, RIS phase shift matrix
and analog precoders. We provide extensive numerical results to evaluate the
performance of the proposed method, and demonstrate its improved performance
with respect to the state-of-the-art EE and fairness maximization techniques
of [4] and [10].

Notations: Bold lowercase and uppercase letters denote vectors and matrices,
respectively. $|\cdot|$ and $\text{arg}(\cdot)$ are the magnitude and the
angle of a complex number, respectively (for a complex vector, they are
assumed to operate element-wise). $||\mathbf{\cdot}||_{2}$ is the Euclidean
norm. $\mathbf{(\cdot)}^{\text{T}}$ is matrix transpose. $\mathbf{I}_{a}$ is
an $a\times a$ identity matrix. $\boldsymbol{0}_{a\times b}$ and
$\boldsymbol{1}_{a\times b}$ are $a\times b$ matrices all of whose elements
are equal to zero and one, respectively. $\text{diag}(a_{1},\dots,a_{N})$ is a
diagonal matrix with the elements $a_{1},\dots,a_{N}$ on the main diagonal and
zeros elsewhere. $\mathbb{C}$ and $\mathbb{R}^{+}$ indicate the set of the
complex and non-negative real numbers, respectively. $j=\sqrt{-1}$ is the
imaginary unit. $\mathbb{E}[\cdot]$ stands for the expectation operator.
$\mathcal{CN}(0,b)$ represents a complex Gaussian random distribution with a
mean of $0$ and a variance of $b$, and $\mathcal{U}(a,b)$ represents a uniform
distribution between $a$ and $b$. $\nabla_{\mathbf{a}}f(\cdot)$ represents the
gradient of the real-valued function $f(\cdot)$ with respect to $\mathbf{a}$.

##  II System Model

Figure 1: System model for the proposed RIS-assisted MU-MIMO mmWave
communication system.

We consider a downlink mmWave system in which a base station (BS) equipped
with $M$ antennas communicates with $K$ single-antenna users as shown in Fig.
1. The direct links between the BS and the users are assumed to be blocked and
the communication is achieved via an RIS with $N$ reflecting elements. In the
following subsections, we describe the received signal model and channel
model.

We assume that the users are served using orthogonal resources (e.g.,
frequency-division multiple access (FDMA))111Orthogonal signaling limits the
rate of the users, but this assumption allows us to keep the rate expressions
tractable. Some other works proposed mitigating the multi-user interference
using zero-forcing beamforming (ZFBF). However, the performance of ZFBF is not
optimal if some RIS-user channels in the system are highly correlated, which
is typical in mmWave systems. ZFBF also requires an additional hardware
complexity compared to the constant-modulus analog beamformer. Interested
readers can refer to [3] and [11] for more details.. The signal received by
user $k$ is

$y_{k}=\mathbf{h}_{2,k}\boldsymbol{\Theta}\mathbf{H}_{1,k}\mathbf{v}_{k}\sqrt{p_{k}}s_{k}+n_{k},$
(1)

where $\mathbf{h}_{2,k}\in\mathbb{C}^{1\times N}$ is the channel between the
$k$-th user and the RIS,
$\boldsymbol{\Theta}=\text{diag}(\vartheta_{1},\dots,\vartheta_{N})\in\mathbb{C}^{N\times
N}$ is the RIS phase shift matrix with $\vartheta_{n}=e^{j\theta_{n}}$ and
$\theta_{n}\in[0,2\pi)$, $\mathbf{H}_{1,k}\in\mathbb{C}^{N\times M}$ is the
BS-RIS channel corresponding to the $k$-th user,
$\mathbf{v}_{k}\in\mathbb{C}^{M\times 1}$ is the analog precoder for user $k$
with each element having a magnitude of $1/\sqrt{M}$, $p_{k}\in\mathbb{R}^{+}$
is the power allocated to user $k$ with $\sum_{k=1}^{K}p_{k}\leq
P_{\text{max}}$ where $P_{\text{max}}$ is the total power budget, $s_{k}$ is
the symbol to be transmitted to the $k$-th user satisfying
$\mathbb{E}[|s_{k}|^{2}]=1$, and
$n_{k}\sim\mathcal{C}\mathcal{N}(0,\sigma^{2})$ is the complex additive white
Gaussian noise (AWGN).

##  III Problem Formulation

In this section, we formulate the multi-objective EE and fairness optimization
problem. We assume that channel state information is available at the
transmitter. The rate of the $k$-th user can be expressed as

$R_{k}(\mathbf{p},\boldsymbol{\Theta},\mathbf{v}_{k})=B\log_{2}\left(1+\frac{1}{\sigma^{2}}p_{k}|\mathbf{h}_{2,k}\boldsymbol{\Theta}\mathbf{H}_{1,k}\mathbf{v}_{k}|^{2}\right),$
(2)

where $B$ is the bandwidth and $\mathbf{p}=[p_{1},\dots,p_{k}]^{\text{T}}$.

On the other hand, the total power consumption of the system can be written as
[4]

$\mathcal{P}_{\text{tot}}(\mathbf{p})=P_{\text{BS}}+\xi\sum\nolimits_{k=1}^{K}p_{k}+NP_{\theta}+\sum\nolimits_{k=1}^{K}P_{\text{U},k},$
(3)

where $P_{\text{BS}}$ is the total hardware static power consumption at the
BS, $\xi$ is the power amplifier efficiency at the BS, $P_{\theta}$ is the
power consumption of each phase shifter at the RIS, and $P_{\text{U},k}$ is
the hardware static power dissipated by user $k$.

The energy efficiency (EE) is then defined as

$\eta_{\text{EE}}(\mathbf{p},\boldsymbol{\Theta},\mathbf{V})=\frac{R_{\text{sum}}(\mathbf{p},\boldsymbol{\Theta},\mathbf{V})}{\mathcal{P}_{\text{tot}}(\mathbf{p})},$
(4)

where
$R_{\text{sum}}(\mathbf{p},\boldsymbol{\Theta},\mathbf{V})=\sum_{k=1}^{K}R_{k}(\mathbf{p},\boldsymbol{\Theta},\mathbf{v}_{k})$,
and $\mathbf{V}=[\mathbf{v}_{1},\dots,\mathbf{v}_{k}]$.

Because of the high path loss in mmWave systems, some users may have
significantly lower channel gains than others. For such a situation,
maximizing (4) leads to allocating most of the resources in favor of the
dominant users and ignoring those with small channel gains. As a result,
achieving fairness among users in the system is also of high importance. This
can be done by maximizing the following objective function:

$F(\mathbf{p},\boldsymbol{\Theta},\mathbf{V})=\min_{k\in\\{1,\dots,K\\}}\frac{R_{k}(\mathbf{p},\boldsymbol{\Theta},\mathbf{v}_{k})}{w_{k}},$
(5)

where $w_{k}$ is a weight assigned to the $k$-th user by the service provider
and represents the relative QoS of that user.

The general EE and fairness optimization problem can then be formulated as a
multi-objective optimization problem as follows:

$\displaystyle\max_{\mathbf{p},\boldsymbol{\Theta},\mathbf{V}}\ \ \
[\eta_{\text{EE}}(\mathbf{p},\boldsymbol{\Theta},\mathbf{V}),F(\mathbf{p},\boldsymbol{\Theta},\mathbf{V})]$
(6a) s.t. $\displaystyle|\vartheta_{n}|=1,\ \ \forall n=1,\dots,N,$ (6b)
$\displaystyle\sum\nolimits_{k=1}^{K}p_{k}\leq P_{\text{max}},$ (6c)
$\displaystyle p_{k}\geq 0,\ \ \forall k=1,\dots,K,$ (6d)
$\displaystyle|v_{k,m}|=\frac{1}{\sqrt{M}},\forall k=1,\dots,K,\forall
m=1,\dots,M,$ (6e)

where $v_{k,m}$ is the $m$-th element of $\mathbf{v}_{k}$. This multi-
objective optimization problem is challenging due to the presence of two
objective functions that can have positive or negative correlation trends. To
address this issue, we adopt a lexicographic approach [12], which divides (6)
into two stages.

In the first stage, EE maximization is targeted by solving the following
optimization problem:

$\displaystyle\eta_{\text{EE}}^{*}=\max_{\mathbf{p},\boldsymbol{\Theta},\mathbf{V}}\eta_{\text{EE}}(\mathbf{p},\boldsymbol{\Theta},\mathbf{V})\
\text{s.t.}\
(\ref{ps_con_EEf}),(\ref{power_con_EEf}),(\ref{power_pos_EEf}),(\ref{CM}).$
(7)

Once $\eta_{\text{EE}}^{*}$ is obtained, the fairness can be maximized subject
to an EE constraint in the second stage via

$\displaystyle\max_{\mathbf{p},\boldsymbol{\Theta},\mathbf{V}}$ $\displaystyle
F(\mathbf{p},\boldsymbol{\Theta},\mathbf{V})$ (8a) s.t.
$\displaystyle(\ref{ps_con_EEf}),(\ref{power_con_EEf}),(\ref{power_pos_EEf}),\text{and
}(\ref{CM}),$ (8b)
$\displaystyle\eta_{\text{EE}}(\mathbf{p},\boldsymbol{\Theta},\mathbf{V})\geq\rho\eta_{\text{EE}}^{*},$
(8c)

where $\rho\in[0,1]$ is a design parameter. This optimization problem
indicates that the fairness should be maximized at the expense of reducing the
EE by a factor not greater than $(1-\rho)$.

The optimal solutions to the optimization problems (7) and (8) are difficult
to obtain directly due to the non-convexity of these problems with respect to
the optimization variables and the coupling of the optimization variables in
the objective functions. In the next section, we propose an alternating
optimization framework to solve (7) and (8).

##  IV Proposed Solution

In this section, we present the method of solving (7) and (8) using an AO
method for solving each optimization problem.

###  IV-A Stage 1: EE Maximization

1\. Power allocation: For fixed values of the RIS phase-shift matrix and
analog precoders, we target maximizing (4) with respect to the power
allocation vector $\mathbf{p}$ subject to the constraints (6c) and (6d).
Dinkelbach’s algorithm [13] can be used to efficiently solve this fractional
programming problem. Specifically, the objective function (4) can be
transformed into a concave function as 222In this manuscript, having some
arguments omitted in a multi-variable function indicates that their values are
being held fixed.

$\mathcal{G}(\mathbf{p},\omega^{(t-1)})=R_{\text{sum}}(\mathbf{p})-\omega^{(t-1)}\mathcal{P}_{\text{tot}}(\mathbf{p}),$
(9)

where $\omega^{(t)}\triangleq\eta_{\text{EE}}(\mathbf{p}^{(t)})$ and
$\mathbf{p}^{(t)}$ denotes the value of $\mathbf{p}$ in the $t$-th iteration.
The optimal power allocation can be obtained by finding the optimal value of
$\mathbf{p}$ that maximizes the concave function (9) subject to the
constraints (6c) and (6d) (this can be achieved using a standard algorithm,
such as the interior-point method or active-set method, via the CVX toolbox
[14]) with $\omega^{(t)}$ updated in each iteration until convergence. A
summary of the power allocation algorithm for EE maximization is given in
Algorithm 1.

1

Input: $\boldsymbol{\Theta}$, $\mathbf{V}$,
$\mathbf{p}^{(1)}\\!=\\!\tfrac{P_{\max}}{K}\boldsymbol{1}_{K\times 1}$,
$\omega^{(1)}\\!=\\!0$, $\tau\\!=\\!1$, $\epsilon\\!>\\!0$

2

Output: $\mathbf{p}=\mathbf{p}^{(\tau)}$

3

4repeat

5

6 Solve the concave optimization problem
$\mathbf{p}^{(\tau+1)}=\operatorname*{arg\,max}_{\mathbf{p}}\mathcal{G}(\mathbf{p},\omega^{(\tau)})\
\text{s.t.}\ (\ref{power_con_EEf}),(\ref{power_pos_EEf})$, where
$\mathcal{G}(\cdot)$ is given by (9);

7

8
$\omega^{(\tau+1)}=\eta_{\text{EE}}(\mathbf{p}^{(\tau+1)},\boldsymbol{\Theta},\mathbf{V})$;

9

10 $\tau\leftarrow\tau+1$;

11

12until  _$|\omega^{(\tau)}-\omega^{(\tau-1)}| <\epsilon$ _;

Algorithm 1 Power allocation based on Dinkelbach’s method for EE maximization.

2\. Phase shift matrix: Since the phase shift matrix appears only in the
numerator of (4), it is sufficient to maximize $R_{\text{sum}}$ subject to the
constraint (6b).

We use the following substitutions for simplicity of illustration:

$\mathbf{H}_{1,k}\mathbf{v}_{k}=[b_{1,k}e^{j\gamma_{1,k}},b_{2,k}e^{j\gamma_{2,k}},\dots,b_{N,k}e^{j\gamma_{N,k}}]^{\text{T}},$
(10)
$\mathbf{h}_{2,k}=[c_{1,k}e^{j\mu_{1,k}},c_{2,k}e^{j\mu_{2,k}},\dots,c_{N,k}e^{j\mu_{N,k}}],$
(11)

where $b_{n,k},c_{n,k}\in\mathbb{R}^{+}$, and
$\gamma_{n,k},\mu_{n,k}\in[0,2\pi)$. Using (10) and (11), we can express the
effective channel gain as

$\begin{split}h_{k}&=|\mathbf{h}_{2,k}\boldsymbol{\Theta}\mathbf{H}_{1,k}\mathbf{v}_{k}|^{2}\\\
&\begin{split}=\sum_{i=1}^{N}b^{2}_{i,k}c^{2}_{i,k}+&2\sum_{m=1}^{N-1}\sum_{\ell=m+1}^{N}\big{[}b_{m,k}c_{m,k}b_{\ell,k}c_{\ell,k}\\\
&\cos(\theta_{m}-\theta_{\ell}+\psi_{m,k}-\psi_{\ell,k})\big{]},\end{split}\end{split}$
(12)

where $\psi_{i,k}=\gamma_{i,k}+\mu_{i,k}$. Therefore

$\frac{\partial
h_{k}}{\partial\theta_{n}}(\boldsymbol{\theta})=\\!-2b_{n,k}c_{n,k}\sum_{i=1,i\neq
n}^{N}\\!\\!\\!\\!b_{i,k}c_{i,k}\sin(\theta_{n}-\theta_{i}+\psi_{n,k}-\psi_{i,k}),$
(13)

and

$\frac{\partial R_{k}}{\partial\theta_{n}}(\boldsymbol{\theta})=\frac{1}{\ln
2}\bigg{(}\frac{p_{k}}{\sigma^{2}+p_{k}h_{k}}\bigg{)}\frac{\partial
h_{k}}{\partial\theta_{n}}(\boldsymbol{\theta}),$ (14)

where $\boldsymbol{\theta}=[\theta_{1},\dots,\theta_{n}]^{\text{T}}$. Also, we
can obtain

$\frac{\partial
R_{\text{sum}}}{\partial\theta_{n}}(\boldsymbol{\theta})=\sum_{k=1}^{K}\frac{\partial
R_{k}}{\partial\theta_{n}}(\boldsymbol{\theta}),$ (15)

and the update of $\boldsymbol{\theta}$ should follow

$\boldsymbol{\theta}^{(t+1)}=\boldsymbol{\theta}^{(t)}+\alpha\nabla_{\boldsymbol{\theta}}R_{\text{sum}}(\boldsymbol{\theta}^{(t)}),$
(16)

where $\boldsymbol{\theta}^{(t)}$ denotes the value of $\boldsymbol{\theta}$
at the $t$-th iteration, $\alpha>0$ is the step size, and
$\nabla_{\boldsymbol{\theta}}R_{\text{sum}}(\boldsymbol{\theta})=[\frac{\partial
R_{\text{sum}}}{\partial\theta_{1}}(\boldsymbol{\theta}),\dots,\frac{\partial
R_{\text{sum}}}{\partial\theta_{N}}(\boldsymbol{\theta})]^{\text{T}}$.

3\. Analog precoder: The $k$-th user’s analog precoder optimization problem is

$\displaystyle\max_{\mathbf{v}_{k}}|\mathbf{h}_{2,k}\boldsymbol{\Theta}\mathbf{H}_{1,k}\mathbf{v}_{k}|^{2}\
\text{s.t.}\ (\ref{CM}).$ (17)

Let

$\mathbf{h}_{2,k}\boldsymbol{\Theta}\mathbf{H}_{1,k}=[s_{1}e^{j\beta_{1}},\dots,s_{M}e^{j\beta_{M}}],$
(18)

and

$\mathbf{v}_{k}=\frac{1}{\sqrt{M}}[e^{j\varphi_{1}},\dots,e^{j\varphi_{M}}],$
(19)

where $s_{n}\in\mathbb{R}^{+}$ and $\beta_{n},\varphi_{n}\in[0,2\pi)$. Then
the objective function in (17) can be written as

$|\mathbf{h}_{2,k}\boldsymbol{\Theta}\mathbf{H}_{1,k}\mathbf{v}_{k}|^{2}\\!=\\!\frac{1}{\sqrt{M}}\bigg{|}\\!\sum_{i=1}^{M}s_{i}e^{j(\varphi_{i}+\beta_{i})}\\!\bigg{|}^{2}\\!\leq\\!\frac{1}{\sqrt{M}}\bigg{|}\sum_{i=1}^{M}s_{i}\bigg{|}^{2},$
(20)

where the upper bound is obtained (using the triangle inequality) by setting
$\varphi_{m}=-\beta_{m}\ \forall m=1,\dots,M$.

###  IV-B Stage 2: Fairness Maximization with an EE Constraint

1\. Power allocation: Next, we target the optimization problem (8) with
respect to the power allocation vector $\mathbf{p}$. To solve this problem, we
introduce an intermediate variable to represent the minimum weighted rate as

$z=\min_{k\in\\{1,\dots,K\\}}\frac{R_{k}}{w_{k}},$ (21)

so that (8) becomes

$\displaystyle\max_{\mathbf{p},z}\ \ \ $ $\displaystyle z$ (22a) s.t.
$\displaystyle(\ref{power_con_EEf})\text{ and }(\ref{power_pos_EEf}),$ (22b)
$\displaystyle
R_{\text{sum}}(\mathbf{p})\geq\rho\eta_{\text{EE}}^{*}\mathcal{P}_{\text{tot}}(\mathbf{p}),$
(22c) $\displaystyle\frac{1}{w_{k}}R_{k}(\mathbf{p})\geq z,\ \ \forall
k=1,\dots,K,$ (22d)

which is a concave optimization problem amenable to solution via standard
solvers such as CVX.

2\. Phase shift matrix: As the next step in the AO approach, we consider the
optimization problem (8) with respect to the phase shift matrix. It is obvious
that the derivative of the minimum function in (5) is non-smooth and has
discontinuities at points where $\frac{R_{i}}{w_{i}}=\frac{R_{k}}{w_{k}}$,
where $i,k\in\\{1,\dots,K\\}$ and $i\neq k$. Therefore, the gradient ascent
method cannot be applied directly. To resolve this issue, we use the well-
known log-sum-exp (LSE) approximation to approximate $F$ as [15]

$F\approx\hat{F}=-\frac{1}{\zeta}\ln\left(\sum\nolimits_{k=1}^{K}e^{-\zeta
R_{k}/w_{k}}\right),$ (23)

where $\zeta>0$. This approximation is accurate for sufficiently large values
of $\zeta$.

The partial derivative of $\hat{F}$ with respect to $\theta_{n}$ can be
written as

$\frac{\partial\hat{F}}{\partial\theta_{n}}(\boldsymbol{\theta})=\frac{\sum_{k=1}^{K}\frac{1}{w_{k}}\frac{\partial
R_{k}}{\partial\theta_{n}}e^{-\zeta R_{k}/w_{k}}}{\sum_{k=1}^{K}e^{-\zeta
R_{k}/w_{k}}},$ (24)

which is a continuous function of $\boldsymbol{\theta}$. Therefore, the
expression (24) can be used (with a large value of $\zeta$) as an
approximation to the partial derivative $\frac{\partial
F}{\partial\theta_{n}}$.

The update of $\boldsymbol{\theta}$ should then follow

$\boldsymbol{\theta}^{(t+1)}=\boldsymbol{\theta}^{(t)}+\alpha\nabla_{\boldsymbol{\theta}}\hat{F}(\boldsymbol{\theta}^{(t)}),$
(25)

where
$\nabla_{\boldsymbol{\theta}}\hat{F}(\boldsymbol{\theta})=[\frac{\partial\hat{F}}{\partial\theta_{1}}(\boldsymbol{\theta}),\dots,\frac{\partial\hat{F}}{\partial\theta_{N}}(\boldsymbol{\theta})]^{\text{T}}$.
However, the constraint (8c) should be checked after each update of
$\boldsymbol{\theta}$; if it is violated by the new value of
$\boldsymbol{\theta}$, the update is not accepted by the algorithm.

3\. Analog precoder: Since it is required to maximize the effective channel
gain of the user with the minimum weighted rate, we can follow the same
procedure as in EE maximization (Staege 1) to obtain the analog precoder.

The two-stage optimization procedure is summarized in Algorithm 2.

1

Input: $\zeta\geq 0$, $\rho\in[0,1]$, $\alpha>0$,
$\boldsymbol{\Theta}^{(1)}=\mathbf{I}_{N}$,
$\mathbf{V}^{(1)}=\tfrac{1}{\sqrt{M}}\boldsymbol{1}_{M\times K}$,
$\eta_{\mathrm{EE}}^{(1)}=0$, $t=1$ $\epsilon>0$

2

Output: $\mathbf{p}=\mathbf{p}^{(t)}$,
$\boldsymbol{\Theta}=\boldsymbol{\Theta}^{(t)}$, $\mathbf{V}=\mathbf{V}^{(t)}$

3

/* Stage 1: EE maximization */

4

5repeat

6

7 Compute $\mathbf{p}^{(t+1)}$ using Algorithm 1;

8

9 $\boldsymbol{\theta}^{(1)}\leftarrow\boldsymbol{0}_{N\times 1}$, $i=1$;

10

11 repeat

12

13
$\boldsymbol{\theta}^{(i+1)}=\boldsymbol{\theta}^{(i)}+\alpha\nabla_{\boldsymbol{\theta}}R_{\mathrm{sum}}\big{(}\boldsymbol{\theta}^{(i)}\big{)}$
;

14

15 $i\leftarrow i+1$;

16

17 until  _$\|\boldsymbol{\theta}^{(i)}-\boldsymbol{\theta}^{(i-1)}\|_{2}
<\epsilon$_;

18
$\boldsymbol{\Theta}^{(t+1)}=\mathrm{diag}\big{(}\exp\big{(}j\boldsymbol{\theta}^{(i)}\big{)}\big{)}$;

19

20 for  _$k=1,2,\ldots,K$_ do

21

22
$\boldsymbol{\beta}_{k}=\text{arg}(\mathbf{h}_{2,k}\boldsymbol{\Theta}\mathbf{H}_{1,k})$;

23

24 $\mathbf{v}_{k}^{(t+1)}=\frac{1}{\sqrt{M}}e^{-j\boldsymbol{\beta}_{k}}$;

25

26 end for

27

28 $\mathbf{V}^{(t+1)}=[\mathbf{v}_{1}^{(t+1)},\dots,\mathbf{v}_{K}^{(t+1)}]$;

29

30 Compute
$\eta_{\text{EE}}^{(t+1)}(\mathbf{p}^{(t+1)},\\!\boldsymbol{\Theta}^{(t+1)},\\!\mathbf{V}^{(t+1)})$
using (4);

31

32 $t\leftarrow t+1$;

33

34until  _$|\eta_{\text{EE}}^{(t)}-\eta_{\text{EE}}^{(t-1)}| <\epsilon$ _;

35$\eta_{\text{EE}}^{*}=\eta_{\text{EE}}^{(t)}$;

36

/* Stage 2: Fairness maximization with EE constraint */

37

38$F^{(t)}=F(\mathbf{p}^{(t)},\boldsymbol{\Theta}^{(t)},\mathbf{V}^{(t})$;

39

40repeat

41

42 Solve the concave optimization problem
$\big{(}\mathbf{p}^{(t+1)},z^{(t+1)}\big{)}=\underset{\mathbf{p},z}{\operatorname*{arg\,max}}\
\ z\ \ \text{s.t.}\
(\ref{power_con_EEf}),\\!(\ref{power_pos_EEf}),\\!(\ref{conn1}),\\!(\ref{conn2})$;

43 $\boldsymbol{\theta}^{(1)}=\boldsymbol{0}_{N\times 1}$;

44

45 repeat

46

47
$\boldsymbol{\theta}^{(i+1)}=\boldsymbol{\theta}^{(i)}+\alpha\nabla_{\boldsymbol{\theta}}\hat{F}\big{(}\boldsymbol{\theta}^{(i)}\big{)}$;

48

49 Compute $\eta_{\text{EE}}(\boldsymbol{\theta}^{(i+1)})$ using (4);

50

51 if
_$\eta_{\text{EE}}(\boldsymbol{\theta}^{(i+1)})\geq\rho\eta_{\text{EE}}^{*}$_
then

52
$\boldsymbol{\theta}^{(i+1)}=\boldsymbol{\theta}^{(i)}-\alpha\nabla_{\boldsymbol{\theta}}\hat{F}(\boldsymbol{\theta}^{(i)})$;

53

54 break out of the loop ;

55

56 end if

57

58 $i\leftarrow i+1$;

59

60 until  _$\|\boldsymbol{\theta}^{(i)}-\boldsymbol{\theta}^{(i-1)}\|_{2}
<\epsilon$_;

61 $\boldsymbol{\Theta}^{(t+1)}=\text{diag}(e^{j\boldsymbol{\theta}^{(i)}})$;

62

63 Compute $\mathbf{V}^{(t+1)}$ as per lines 9-13;

64

65 Compute
$F^{(t+1)}(\\!\mathbf{p}^{(t+1)},\\!\boldsymbol{\Theta}^{(t+1)},\\!\boldsymbol{V}^{(t+1)}\\!)$
using (5);

66

67 $t\leftarrow t+1$;

68

69until  _$|F^{(t)}-F^{(t-1)}| <\epsilon$_;

Algorithm 2 Lexicographic approach for EE and fairness maximization.

$0$$3$$6$$9$$12$$15$$0$$40$$80$$120$Iteration numberAverage EE
[Mbits/s/Joule]EE maximization (Stage 1)Fairness maximization (Stage 2),
$\rho=0.85$Fairness maximization (Stage 2), $\rho=0.5$

Figure 2: Average EE vs. iteration number for EE maximization (Stage 1) and
fairness maximization (Stage 2) at $P_{\max}=25$ dBm.

$0$$3$$6$$9$$12$$15$$0$$5$$10$$15$$20$$26$Iteration numberAvg. min. weighted
rate [Mbits/s]Fairness maximization (Stage 2), $\rho=0.5$Fairness maximization
(Stage 2), $\rho=0.85$EE maximization (Stage 1)

Figure 3: Avg. min. weighted rate vs. iteration number for EE maximization
(Stage 1) and fairness maximization (Stage 2) at $P_{\max}=25$ dBm.

$5$$15$$25$$35$$45$$0$$40$$80$$120$$160$$200$$P_{\max}$ [dBm]Average EE
[Mbits/s/Joule]EE maximization [4]Proposed, $\rho=0.85$Proposed, $\rho=0.5$EE-
aware fairness maximization [10]

Figure 4: Average EE as the power budget varies for EE maximization [4], EE-
aware fairness maximization [10], and the proposed approach.

$5$$15$$25$$35$$45$$0$$0.2$$0.4$$0.6$$0.8$$1$$P_{\max}$ [dBm]Average Jain’s
fairness indexEE-aware fairness maximization [10]Proposed, $\rho=0.5$Proposed,
$\rho=0.85$EE maximization [4]

Figure 5: Average Jain’s fairness index as the power budget varies for EE
maximization [4], EE-aware fairness maximization [10], and the proposed
approach.

$4$$100$$200$$300$$400$$0$$100$$200$$300$$400$$450$Number of RIS
elementsAverage EE [Mbits/s/Joule]EE maximization [4]Proposed,
$\rho=0.85$Proposed, $\rho=0.5$EE-aware fairness maximization [10]

Figure 6: Average EE as the number of RIS elements increases for EE
maximization [4], EE-aware fairness maximization [10], and the proposed
approach.

$4$$100$$200$$300$$400$$0.2$$0.4$$0.6$$0.8$$1$Number of RIS elementsAverage
Jain’s fairness indexEE-aware fairness maximization [10]Proposed,
$\rho=0.5$Proposed, $\rho=0.85$EE maximization [4]

Figure 7: Average Jain’s fairness index as the number of RIS elements
increases for EE maximization [4], EE-aware fairness maximization [10], and
the proposed approach.

##  V Numerical Results

In this section, the performance of the proposed method is demonstrated and
compared with that of two benchmark schemes, the pure EE maximization scheme
of [4], and the EE-aware fairness maximization scheme of [10].

A square cell with a side of 150 m is considered in simulations with a BS
located at the point $(0,0,10)$ m having $M=16$ antennas that serve four users
with $w_{k}\sim\mathcal{U}(1,4)$. FDMA transmission is used with a center
frequency for the $k$-th user of $f_{k}=27.75+(2k-1)0.0625\,\mathrm{GHz}$
(each user has a bandwidth of 125 MHz). A spacing of $\mathbf{\lambda_{c}}/2$
is assumed between two adjacent antennas at the BS, where $\lambda_{c}$ is the
carrier wavelength for the system center frequency $f_{c}=28\,\mathrm{GHz}$.
The number of RIS elements is assumed to be $N=64$ with both vertical and
horizontal spacings of $\mathbf{\lambda_{c}}/2$ between adjacent elements. We
assume that the locations of the users are uniformly distributed in the
cuboidal region with opposite corners at points $(30,-75,0)$ m and
$(150,75,2)$ m. The total hardware static power at the BS is
$P_{\text{BS}}=9\,\mathrm{dBW}$, while that of the $k$-th UE is
$P_{\text{U},k}=10\,\mathrm{dBm}$. The power consumption of each RIS phase
shifter is assumed to be $P_{\theta}=1\,\mathrm{dBm}$ and the circuit
dissipated power coefficient at the BS $\xi$ is set to $1.2$ [4]. The minimum
function approximation parameter $\zeta$ is set to $50$ and the algorithmic
convergence threshold is selected to be $\epsilon=1\times 10^{-3}$.

Monte Carlo simulation is performed to obtain the results by averaging the
performance metric over 1000 combinations of random channels and user
locations. The widely used Saleh-Valenzuela (SV) channel model [16] is
adopted, where the number of paths is assumed to be $L=4$ [4, 5]. The path
gain of the LoS component is calculated using the free space path loss
formula, and a difference of 15 dBm is assumed between the LoS and NLoS path
gains [17]. The small-scale fading is modeled as a Gaussian random variable
with a 0 dB mean and a 1 dB variance.

Two metrics are used to compare the performance of the different methods; the
EE (defined in (4)) and Jain’s fairness index for weighted rates, which is
defined as [18]

$F_{J}=\frac{(\sum_{k=1}^{K}R_{k}/w_{k})^{2}}{K\sum_{k=1}^{K}(R_{k}/w_{k})^{2}}\in[1/K,1],$
(26)

where $F_{J}=1$ indicates that all users have the same weighted rate
$R_{k}/w_{k}$.

To examine the convergence of the objective functions (4) and (5) during
Algorithm 2, we plot the EE and the minimum weighted rate versus the iteration
number in Figs. 4 and 4, respectively, for the EE maximization problem (Stage
1) and the fairness maximization problem (Stage 2) for two values of $\rho$
($0.85$ and $0.5$) and $P_{\text{max}}=25\,\mathrm{dBm}$. It can be noted that
in all cases the proposed algorithm converges within a small number of
iterations. In addition, the minimum weighted rate produced by allocating the
resources to maximize the EE (Stage 1) is very low. However, when 15% and 50%
maximum reductions in the EE are allowed (corresponding to $\rho=0.85$ and
$\rho=0.5$, respectively), the minimum weighted rate can be increased
significantly. The actual reductions in the EE are $13$% and $42$% for
$\rho=0.85$ and $\rho=0.5$ respectively, which are less than maximum allowed
reductions $1-\rho$.

Figs. 4 and 7 show the EE and Jain’s fairness index of the proposed method
versus the maximum transmit power, respectively, for $\rho=0.85$ and
$\rho=0.5$. The proposed method is compared with the benchmark algorithms [4]
and [10]. It can be noticed that all curves saturate when $P_{\text{max}}$
exceeds 40 dBm as the additional power does not help in increasing the EE
objective function (4). On the one hand, the method in [4] achieves higher EE
than the other algorithms. However, it produces a very low fairness index,
especially at low power budgets. On the other hand, [10] achieves excellent
fairness even at low power budgets, but lower EE than the other methods. The
proposed method achieves a controllable trade-off between the aforementioned
extreme situations. If the value of $\rho$ is chosen carefully (e.g.,
$\rho=0.85$), a very good level of fairness can be attained with with a pre-
designed minor reduction in the EE.

Figs. 7 and 7 depict the EE and Jain’s fairness index for
$P_{\text{max}}=25\,\mathrm{dBm}$ as a function of the number of RIS elements,
respectively. The relationship between the EE and the number of RIS elements
follows a log-like trend. Increasing the number of RIS elements makes the
numerator of (4) increase logarithmically while the denominator increases
linearly. Since the power consumption of the RIS elements is small, the ratio
(4) saturates. In other terms, adding extra RIS elements can have a
significant contribution towards boosting both the EE and fairness. In the
case of fairness maximization, a small number of RIS elements (i.e., $N\leq
50$) is sufficient to achieve a Jain’s fairness score of 1. A larger number of
of RIS elements allows the proposed method to achieve a very good level of
fairness (i.e., higher than 0.9) with higher fairness for $\rho=0.5$ than for
$\rho=0.85$. This comes at the expense of a minor penalty in the EE as can be
seen in Fig. 7.

##  VI Conclusion

This paper presents a two-stage lexicographic approach to find the optimal
power allocation, RIS phase shift matrix, and analog precoders that maximize
both EE and fairness in an RIS-assisted mmWave MU-MISO system. First, the EE
is maximized, then the minimum weighted rate is maximized subject to a
constraint on the allowable reduction in the EE. The optimization problems in
both stages are solved via an AO procedure that uses Dinkelbach’s method for
power allocation, gradient ascent for RIS phase shift matrix optimization, and
beam alignment for analog precoder design. The results show that the proposed
method can achieve a better trade-off between EE and fairness compared to
methods that optimize only one of these metrics.

## Acknowledgement

The authors would like to thank Dr. Arman Farhang for helpful discussions.

## References

  * [1] W. Jiang, B. Han, M. A. Habibi, and H. D. Schotten, “The road towards 6G: A comprehensive survey,” _IEEE Open J. Commun. Soc._ , vol. 2, pp. 334–366, Feb. 2021.
  * [2] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M.-S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, Aug. 2019.
  * [3] X. Wang, L. Kong, F. Kong, F. Qiu, M. Xia, S. Arnon, and G. Chen, “Millimeter wave communication: A comprehensive survey,” _IEEE Commun. Surveys Tuts._ , vol. 20, no. 3, pp. 1616–1653, 2018.
  * [4] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4157–4170, Jun. 2019.
  * [5] Q. Wu, Y. Zhang, C. Huang, Y. Chau, Z. Yang, and M. Shikh-Bahaei, “Energy efficient intelligent reflecting surface assisted terahertz communications,” in _2021 IEEE ICC Workshops_ , Jun. 2021, pp. 1–6.
  * [6] Y. Liu, E. Liu, and R. Wang, “Energy efficiency analysis of intelligent reflecting surface system with hardware impairments,” in _2020 IEEE GLOBECOM_ , Dec. 2020, pp. 1–6.
  * [7] K. D. Gupta, R. Nigam, D. K. Sharma, and S. K. Dhurandher, “LSTM-based energy-efficient wireless communication with reconfigurable intelligent surfaces,” _IEEE Trans. Green Commun. Netw._ , vol. 6, no. 2, pp. 704–712, Dec. 2022.
  * [8] J. Woo, C. Song, and I. Lee, “Sum rate and fairness optimization for intelligent reflecting surface aided multiuser systems,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 12, pp. 13 436–13 440, Oct. 2021.
  * [9] D. Li, “Fairness-aware multiuser scheduling for finite-resolution intelligent reflecting surface-assisted communication,” _IEEE Commun. Lett._ , vol. 25, no. 7, pp. 2395–2399, Apr. 2021.
  * [10] K. Liu and Z. Zhang, “On the energy-efficiency fairness of reconfigurable intelligent surface-aided cell-free network,” in _IEEE VTC-Spring_ , Apr. 2021, pp. 1–6.
  * [11] C. Lin and G. Y. L. Li, “Terahertz communications: An array-of-subarrays solution,” _IEEE Commun. Mag._ , vol. 54, no. 12, pp. 124–131, Dec. 2016\.
  * [12] J. Branke, K. Deb, K. Miettinen, and R. Słowiński, “Multiobjective optimization: Interactive and evolutionary approaches,” _Springer_ , pp. 19–23, 2008.
  * [13] W. Dinkelbach, “On nonlinear fractional programming,” _Management Science_ , vol. 13, no. 7, pp. 492–498, 1967. [Online]. Available: http://www.jstor.org/stable/2627691
  * [14] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.
  * [15] Y. Nesterov, “Smooth minimization of non-smooth functions,” _Mathematical Programming_ , p. 127–152, 2005.
  * [16] M. R. Akdeniz, Y. Liu, M. K. Samimi, S. Sun, S. Rangan, T. S. Rappaport, and E. Erkip, “Millimeter wave channel modeling and cellular capacity evaluation,” _IEEE J. Sel. Areas Commun._ , vol. 32, no. 6, pp. 1164–1179, Jun. 2014.
  * [17] J. Ko, Y.-J. Cho, S. Hur, T. Kim, J. Park, A. F. Molisch, K. Haneda, M. Peter, D.-J. Park, and D.-H. Cho, “Millimeter-wave channel measurements and analysis for statistical spatial channel model in in-building and urban environments at 28 GHz,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 9, pp. 5853–5868, Jun. 2017.
  * [18] H. Shi, R. V. Prasad, E. Onur, and I. Niemegeers, “Fairness in wireless networks: Issues, measures and challenges,” _IEEE Commun. Surveys Tuts._ , vol. 16, no. 1, pp. 5–24, May 2014.


# quantum systems simulatability through classical networks

M. Caruso<EMAIL_ADDRESS>Fundación I+D del Software Libre, FIDESOL,
Granada (18100), España.  Campus de Fuentenueva, Universidad de Granada,
Granada (18071), España.

(August 27, 2024)

###### Abstract

We have shown that quantum systems on finite-dimensional Hilbert spaces are
equivalent under local transformations. Using these transformations give rise
to a gauge group that connects the hamiltonian operators associated with each
quantum system. Different quantum systems are connected in such way that
studying one of them allows to understand the other. This result can be
applied to the field of simulation of quantum systems, in order to mimic more
complicated quantum systems from another simulatable quantum system. Given
that there is a bridge that allows to simulate a particular quantum system on
this kind of Hilbert spaces using classical circuits we will provide a general
scenario to extend this bridge to simulate the time evolution, via Schrödinger
equation, of any of these quantum system using classical circuits. This
classical systems can be implemented and controlled more easily in the
laboratory than the quantum systems.

##  I Introduction

In this work, we have developed a procedure to connect a given pair of quantum
systems via a local transformation. We describe specifically a map among the
respective Hilbert spaces that connect its vector objects (which represent
quantum states) and its hamiltonian operators. We will studied the case in
which the corresponding Hilbert spaces are finite-dimensional, but this
results can be enunciated for infinite, but countable, dimensional Hilbert
spaces. This correspondence is a useful tool to map quantum systems in order
to study one of them through the other one.

A productive and promising field to apply these ideas may be the quantum
simulations. At the end of 20th century, R. Feynman asked the following
question: What kind of computer are we going to use to simulate physics? […] I
present that as another interesting problem: to work out the classes of
different kinds of quantum mechanical systems which are really
intersimulatable $-$which are equivalent$-$ […] The same way we should try to
find out what kinds of quantum mechanical systems are mutually
intersimulatable, and try to find a specific class, or a character of that
class which will simulate everything [1].

A quantum simulator is conceived as a controllable system whose aim is to
mimic the static or dynamical properties of another quantum system [2]. Today
this field can be roughly divided into digital quantum simulations, analog
quantum simulations, and a combination of both. The digital quantum
simulation, proposed by Lloyd [3], deals with the synthesis of a given
operator evolution in quantum gates. The advantage of this approach is the
flexibility, which introduce the quantum error correction and, then,
universality. One of its disadvantages is the large number of quantum gates
that may be involved in the synthesis, which implies huge technological effort
to maintain the coherence of the states. Whereas in the analog quantum
simulations there is a first quantum system that may be not experimentally
easy realizable or controllable and a second quantum system that mimics the
first one. The advantage of this approach is the possibility to describe
quantum systems in larger Hilbert spaces. One of the main issues concerns to
find a quantum system able to mimic certain aspects to simulate the first
quantum system. And finally, there is a win-win strategy that consists of a
hybrid digital-analog simulations to combine the best of both ideas. [4]

Focusing on the digital simulation of hamiltonian dynamics for quantum
systems, we see that it is inevitable to deal with a numerable and even finite
version of such quantum systems, i.e. the involved Hilbert spaces are finite-
dimensional. On the other hand in the analog quantum simulations, we will
focus on countable$-$dimensional Hilbert spaces. According to this, we must
study quantum systems on a countable dimensional (denumerable) Hilbert spaces
is essentially relevant.

We will show that any two quantum systems on respective Hilbert spaces which
are finite dimensional are connected via a gauge transformation. This includes
the case in which any of its corresponding hamiltonian may be time dependent.
We intend to deal with the topic of quantum simulation from an alternative
perspective, starting from a more well-known one, in the sense that the last
one can be analytically soluble and/or simulatable. We intend to open a way to
establish the equivalence class previously mentioned by Feynman [1].

On the other hand, from another formal equivalence between quantum and
classical systems proposed in [6] in order to simulate quantum systems through
specific circuits, it is possible to use such classical systems which their
controllability is simpler than for quantum systems in general, in order to
adequately describe its temporal evolution. We will use this equivalence
between quantum systems in order to show that a general simulation protocol is
possible to implement using classical circuits.

##  II Quantum systems on a denumerable hilbert space

Reviewing the basics aspects of quantum systems, let us consider a general
quantum system Q which can be described in a certain $n-$dimensional Hilbert
space $\mathscr{H}_{n}$. The deterministic temporal evolution of a quantum
system is driven by a hamiltonian operator $H$ (eventually time-dependent)
defined on $\mathscr{H}_{n}$. This operator modifies the vector state
$|\psi(t)\rangle$ at time $t{\in\mathcal{T}\subseteq}\mathbb{R}$, by the
equation

$i\partial_{t}|\psi(t)\rangle=H|\psi(t)\rangle,$ (1)

where $\partial_{t}$ represents the partial time derivative. Note that a
partial time derivative is used because $|\psi(t)\rangle$ can be dependent of
other quantities. The equation (1) is written in natural units, e.g.
$\hbar=1$. A solution of (1) is expressed as a $t-$parametrized curve on
$\mathscr{H}_{n}$ $|\psi\rangle{:}\mathcal{T}\longrightarrow\mathscr{H}_{n}$,
this vector curve can be represented using an orthonormal basis of $n-$states:
$\boldsymbol{\beta}=\\{|\beta_{k}\rangle\\}_{k\in I_{n}}$, where
$I_{n}{=}\\{1,\cdots,n\\}$ is the set of the first $n$ natural numbers. The
inner product defined in the Hilbert space $\mathscr{H}_{n}$ allow us to
express the state of the system at time $t$, $|\psi(t)\rangle$, in terms of
its coordinates in the basis $\boldsymbol{\beta}$ as
$\psi_{k}(t):=\langle\beta_{k}|\psi(t)\rangle$ where $k\in I_{n}$.

Note that the bra-ket notation is used to denote the inner product in
$\mathscr{H}_{n}$,
$\langle\star|\ast\rangle:\mathscr{H}_{n}\times\mathscr{H}_{n}\longrightarrow\mathbb{C}$.
Thus, we have a time-parametrized curve on $\mathbb{C}^{n}$,
$\boldsymbol{\psi}:\mathcal{T}\longrightarrow\mathbb{C}^{n}$, where
$\boldsymbol{\psi}$ is written in terms of the coordinates of $|\psi\rangle$
in base $\boldsymbol{\beta}$,
$\boldsymbol{\psi}{=}(\psi_{1},\cdots,\psi_{k},\cdots,\psi_{n})^{\mathfrak{t}}$,
where $\mathfrak{t}$ is the matrix transposition. Both time-parametrized
curves in the abstract Hilbert space and its matrix representation refers to
the same quantum system but belongs to different spaces $\mathscr{H}_{n}$ and
$\mathbb{C}^{n}$.

Also, the complex vector curve $\boldsymbol{\psi}$ satisfies another version
of the equation (1), given by

$id_{t}\boldsymbol{\psi}(t)=\boldsymbol{H\psi}(t),$ (2)

where $\boldsymbol{H}{\in}\mathbb{C}^{n\times n}$ is a complex matrix that
represents the hamiltonian operator $H$ in the basis $\boldsymbol{\beta}$ and
whose matrix elements are $H_{kl}=\langle\beta_{k}|H|\beta_{l}\rangle$. In
this manuscript, we refer to the hamiltonian operator, or hamiltonian matrix
simply as hamiltonian. Note that $d_{t}$ is a total time derivative is used,
because each coordinate $\psi_{k}$ depends on time only.

##  III Formal aspects of equivalent quantum systems

We considered a map $\boldsymbol{\Omega}_{\boldsymbol{\omega}}$, given a
nonsingular matrix $\boldsymbol{\omega}$, which transforms a matrix
$\boldsymbol{H}\in\mathbb{C}^{n\times n}$ as

$\boldsymbol{\Omega}_{\boldsymbol{\omega}}(\boldsymbol{H})=\boldsymbol{\omega}\boldsymbol{H}\boldsymbol{\omega}^{-1}+i(d_{t}\boldsymbol{\omega})\boldsymbol{\omega}^{-1},$
(3)

where $\boldsymbol{\omega}(t)$ is a $t-$differentiable non$-$singular matrix
of $n\times n$, i.e.
$\boldsymbol{\omega}:\mathcal{T}\longrightarrow\mathsf{GL}(n,\mathbb{C})$,
also $\boldsymbol{\Omega}_{\boldsymbol{\omega}}\in\mathcal{M}_{n\times
n}(\mathbb{C})$. The map $\boldsymbol{\Omega_{\omega}}$ (3) is composed by a
similarity transformation of $\boldsymbol{H}$, defined by
$\boldsymbol{\omega}$, plus another time-dependent term. The collection of
this transformations
$\\{\boldsymbol{\Omega}_{\boldsymbol{\omega}}:\forall\,\boldsymbol{\omega}\in\mathsf{GL}(n,\mathbb{C})\\}$
form a group of local (gauge) transformations, with the composition of maps as
a single associative binary operation. The locality of the transformation is
due to the $t-$dependence of $\boldsymbol{\omega}$. Note that
$\mathcal{M}_{n\times n}(\mathbb{K})$ are the $n\times n$ matrices over the
field $\mathbb{K}$.

This kind of mapping was studied in previous works from a pure mathematical
point of view for applications to differential equations in complex variables
with singular operators [7, 8, 9]. For a physical point of view the same kind
of mapping was presented in [10, 11, 12] in order to solve particular quantum
systems. Respect to that, in this section we studied the possibility to
connect any pair of hamiltonian $(H,H^{\prime})$ operators defined on their
respective $n-$dimensional Hilbert spaces
$(\mathscr{H}_{n},\mathscr{H}_{n}^{\prime})$; this hamiltonians are
represented by the matrices $(\boldsymbol{H},\boldsymbol{H}^{\prime})$
eventually time dependent. We proved that there is a non singular matrix
$\boldsymbol{\omega}$, $t-$dependent and differentiable, that connect
$\boldsymbol{H}$ and $\boldsymbol{H}^{\prime}$ in this way
$\boldsymbol{H}^{\prime}=\boldsymbol{\Omega}_{\boldsymbol{\omega}}(\boldsymbol{H})$.

If we composed two transformations
$\boldsymbol{\Omega}_{\boldsymbol{\omega}}\circ\,\boldsymbol{\Omega}_{\boldsymbol{\omega}^{\prime}}$
with $\boldsymbol{\omega}$ and $\boldsymbol{\omega}^{\prime}$ are nonsingular,
we see that
$\boldsymbol{\Omega}_{\boldsymbol{\omega}}\circ\,\boldsymbol{\Omega}_{\boldsymbol{\omega}^{\prime}}=\boldsymbol{\Omega}_{\boldsymbol{\omega.}\boldsymbol{\omega}^{\prime}}$,
thus
$[\,\boldsymbol{\omega,\omega}^{\prime}]=\boldsymbol{0}\Longrightarrow[\,\boldsymbol{\Omega_{\omega}},\boldsymbol{\Omega}_{\boldsymbol{\omega}^{\prime}}]$.
So that $\boldsymbol{\Omega_{I}}=\boldsymbol{I}$, where $\boldsymbol{I}$ is
the identity matrix. If we consider the composed property with
$\boldsymbol{\omega}^{\prime\prime}=\boldsymbol{\omega.\,\omega}^{\prime}$
such that
$\boldsymbol{\omega.\,\omega}^{\prime}=\boldsymbol{I}=\boldsymbol{\omega}^{\prime}\boldsymbol{.\,\omega}$
then
$\boldsymbol{\Omega_{\omega}}\circ\,\boldsymbol{\Omega}_{\boldsymbol{\omega}^{\prime}}=\boldsymbol{I}=\boldsymbol{\Omega}_{\boldsymbol{\omega}^{\prime}}\circ\,\boldsymbol{\Omega}_{\boldsymbol{\omega}}$,
we obtain a unique inverse of $\boldsymbol{\Omega}_{\boldsymbol{\omega}}$
given by
$(\boldsymbol{\Omega_{\omega}})^{{-}1}{=}\;\boldsymbol{\Omega}_{\boldsymbol{\omega}^{-1}}$.

We demonstrated that for any pair of $n\times n$, eventually $t-$dependent and
differentiable, matrices $\boldsymbol{H}$ and $\boldsymbol{H}^{\prime}$ there
exist a non-singular $n\times n$, $t-$dependent and differentiable matrix
$\boldsymbol{\omega}$ that connect them. For that we can define the following
equivalence relation:

$\displaystyle\boldsymbol{H}^{\prime}\sim\boldsymbol{H}\Longleftrightarrow\exists\;\boldsymbol{\omega}:\boldsymbol{H}^{\prime}=\boldsymbol{\Omega}_{\boldsymbol{\omega}}(\boldsymbol{H}).$
(4)

From the equivalence relation (4) then $\boldsymbol{\omega}$ satisfies the
differential equation:

$\displaystyle
id_{t}\,\boldsymbol{\omega}=\boldsymbol{H}^{\prime}\boldsymbol{\omega}-\boldsymbol{\omega
H}.$ (5)

First of all, the solution of (5) exists for the trivial cases
$\boldsymbol{H}=\boldsymbol{0}$ and $\boldsymbol{H}^{\prime}=\boldsymbol{0}$,
i.e. denoted by $\boldsymbol{\omega}_{1}$ and $\boldsymbol{\omega}_{2}$ such
that
$id_{t}\,\boldsymbol{\omega}_{1}{=}\boldsymbol{H}^{\prime}\boldsymbol{\omega}_{1}$
and
$id_{t}\,\boldsymbol{\omega}_{2}{=}-\boldsymbol{\omega}_{2}\boldsymbol{H}$.
Note that we can obtain $(\boldsymbol{\omega}_{1},\boldsymbol{\omega}_{2})$ as
iterative nonsingular solutions following [13]. The existence of these
solutions $\boldsymbol{\omega}_{1}$ and $\boldsymbol{\omega}_{2}$ in terms of
the defined equivalent relation (4) is written as
$\boldsymbol{H}^{\prime}\sim\boldsymbol{0}$ and
$\boldsymbol{0}\sim\boldsymbol{H}$, respectively. From transitivity of the
equivalence relation (4) we have $\boldsymbol{H}^{\prime}\sim\boldsymbol{H}$.
This means that there is a given $\boldsymbol{\omega}$ that
$\boldsymbol{H}^{\prime}=\boldsymbol{\Omega}_{\boldsymbol{\omega}}(\boldsymbol{H})$,
and using the composite property of
$\boldsymbol{\Omega}_{\boldsymbol{\omega}}$ finally express the solution
$\boldsymbol{\omega}$ as a function of the solutions
$(\boldsymbol{\omega}_{1},\boldsymbol{\omega}_{2})$. The equivalences
$\boldsymbol{H}^{\prime}\sim\boldsymbol{0}$ and
$\boldsymbol{0}\sim\boldsymbol{H}$ corresponds to
$\boldsymbol{H}^{\prime}=\boldsymbol{\Omega}_{\boldsymbol{\omega}_{1}}(\boldsymbol{0})$
and
$\boldsymbol{0}=\boldsymbol{\Omega}_{\boldsymbol{\omega}_{2}}(\boldsymbol{H})$
then
$\boldsymbol{H}^{\prime}=\boldsymbol{\Omega}_{\boldsymbol{\omega}_{1}}\boldsymbol{(}\boldsymbol{\Omega}_{\boldsymbol{\omega}_{2}}(\boldsymbol{H})\boldsymbol{)}$
and from the composition rule applied to
$\boldsymbol{H}^{\prime}=\boldsymbol{\Omega}_{\boldsymbol{\omega}_{1}\boldsymbol{\omega}_{2}}(\boldsymbol{H})$,
where the transitive solution is given by
$\boldsymbol{\omega}=\boldsymbol{\omega}_{1}\boldsymbol{\omega}_{2}$. The
solution $\boldsymbol{\omega}$ built in this way is called a transitive
solution, or composite solution.

We have demonstrated that for any pair of this kind of matrices
$(\boldsymbol{H},\boldsymbol{H}^{\prime})$, there is a nonsingular matrix
$\boldsymbol{\omega}$ that connects $\boldsymbol{H}$ and
$\boldsymbol{H}^{\prime}$ through the mapping $\boldsymbol{\Omega_{\omega}}$,
given by the expression (3), this is

$\boldsymbol{H}^{\prime}=\boldsymbol{\Omega_{\omega}}(\boldsymbol{H}).$ (6)

Suppose now that this pair of matrices $\boldsymbol{H}$ and
$\boldsymbol{H}^{\prime}$ are the hamiltonian operators of the following
differential equations

$\displaystyle id_{t}\boldsymbol{\psi}=\boldsymbol{H}\,\boldsymbol{\psi},\quad
id_{t}\boldsymbol{\psi}^{\prime}=\boldsymbol{H}^{\prime}\boldsymbol{\psi}^{\prime},$
(7)

finally, from (6) and (7) we have

$\displaystyle\boldsymbol{\psi}^{\prime}=\boldsymbol{\omega}\boldsymbol{\psi}.$
(8)

In summary, the connection between $\boldsymbol{H}$ and
$\boldsymbol{H}^{\prime}$ can be found at the level of the solutions of (7),
i.e. $\boldsymbol{\omega}$ connects both hamiltonian matrices via
$\boldsymbol{H}^{\prime}=\boldsymbol{\Omega_{\omega}}(\boldsymbol{H})$ (6) and
also both solutions via
$\boldsymbol{\psi}^{\prime}=\boldsymbol{\omega}\boldsymbol{\psi}$ (8). Of
course, if the natural units are not used, will be able to define a similar
map $\boldsymbol{\Omega_{\omega}}$, multiplying the second term on the right
side of expression (3) by $\hslash$.

An important comment about the map $\boldsymbol{\Omega_{\omega}}$: if
$\boldsymbol{\omega}$ is a unitary matrix (i.e.
$\boldsymbol{\omega}^{\dagger}=\boldsymbol{\omega}^{-1}$, where $\dagger$
corresponds to a complex transposition operation), then
$\boldsymbol{\Omega_{\omega}}$ is an endomorphism over the space of the self-
adjoint operators (each of those are defined over the respective Hilbert
spaces that have the same dimension $n$). This situation corresponds to a
mapping of self adjoint hamiltonian operators associated to closed quantum
systems.

In the general case these hamiltonians
$(\boldsymbol{H},\boldsymbol{H}^{\prime})$ are not necessarily hermitian and
even they can be time-dependent. We can express the evolution operator of each
quantum systems $(\boldsymbol{U},\boldsymbol{U}^{\prime})$ related via
$\boldsymbol{\omega}(t)$ as:
$\boldsymbol{U}^{\prime}(t,s)=\boldsymbol{\omega}(t)\,\boldsymbol{U}(t,s)\,\boldsymbol{\omega}^{-1}(s)$,
for all $t,s\in\mathbb{R}$. In Figure 1 a commutative diagram shows how is the
composition of this transformation. As we said for the case of
$(\boldsymbol{H},\boldsymbol{H^{\prime}})$ are hermitians we have

$\boldsymbol{U}^{\prime}(t,s)=\boldsymbol{\omega}(t)\,\boldsymbol{U}(t,s)\,\boldsymbol{\omega}^{\dagger}(s).$
(9)

where $\boldsymbol{\omega}$ is an unitary matrix.

${\boldsymbol{\psi}(s)}$${\boldsymbol{\psi}(t)}$${\boldsymbol{\psi}^{\prime}(s)}$${\boldsymbol{\psi}^{\prime}(t)}$$\scriptstyle{\boldsymbol{U}(t,s)}$$\scriptstyle{\boldsymbol{\omega}(s)}$$\scriptstyle{\boldsymbol{\omega}(t)}$$\scriptstyle{\boldsymbol{U}^{\prime}(t,s)}$
Figure 1: This commutativity diagram shows how is the composition of
transformations between equivalent quantum systems. The commutativity comes
from the existence of the inverse of $\boldsymbol{\omega}(x)$ or the inverse
of $\boldsymbol{U}(y,z)$ for all $x,y,z$.

A general strategy to simulate the arrival system Q′ consists the selection on
the departure system Q and the connection between them, given by
$\boldsymbol{\omega}$. The departure system Q is chosen in order to simulate
it analogically or digitally. The bridge between Q and Q′, namely
$\boldsymbol{\omega}$, is not necessarily associated with any particular
quantum system. For this reason the operator $\boldsymbol{\omega}$ will be
digitally simulated from a synthesis in quantum gates [14], which is an
universal method. Therefore the expression (9) could be implemented through a
hybrid kind of simulation, that involves a digital-analog protocol for its
implementation, as an indispensable requirement [4].

##  IV A connection with electrical networks

In previous works Rosner [15] has proposed an interesting analogy between
particular quantum system and a classical system of electrical oscillators.
This specifical proposal has been formalized and experimentally performed in
[5, 6]. In this paper we will do two things: one is to generalize this
formalization in order to include quantum systems over an arbitrary
dimensional Hilbert space. On the other hand, making use of the fact that all
these quantum systems are equivalent via $\boldsymbol{\Omega_{\omega}}$, the
route by which it is possible to implement the simulation of quantum systems
via classical circuits is well defined.

Given that the equivalence between quantum systems, we can consider without
loss of generality, a time-independent and self-adjoint hamiltonian operator,
represented by a constant hermitian matrix $\boldsymbol{H}$. Let’s start by
rewriting (2) using the decomplexification procedure a la Arnold [16, 17],
also called realification of the space $\mathbb{C}^{n}$ through the operator
$\mathfrak{D}:\mathbb{C}^{n}\longrightarrow\mathbb{R}^{2n}$, explicitly
defined as
$\mathfrak{D}(\boldsymbol{\psi}):=\boldsymbol{(}\mathrm{Re}(\boldsymbol{\psi}),\mathrm{Im}(\boldsymbol{\psi})\boldsymbol{)}^{\mathfrak{t}}=\boldsymbol{(}\boldsymbol{\varphi}_{1},\boldsymbol{\varphi}_{2}\boldsymbol{)}^{\mathfrak{t}}$;
which separates the real and the imaginary part of the complex
$\boldsymbol{\psi}$ and $\mathfrak{t}$ denotes the matrix transposition.

The equation (2) take the form of two separate equation for real and imaginary
part of $\boldsymbol{\psi}$.

$\displaystyle\begin{aligned}
\boldsymbol{\dot{\varphi}}_{1}&=\boldsymbol{H}_{2}\,\boldsymbol{\varphi}_{1}-\boldsymbol{H}_{1}\,\boldsymbol{\varphi}_{2}\\\
\boldsymbol{\dot{\varphi}}_{2}&=\boldsymbol{H}_{1}\,\boldsymbol{\varphi}_{1}+\boldsymbol{H}_{2}\,\boldsymbol{\varphi}_{2},\end{aligned}$
(10)

where $\boldsymbol{H}_{1}$ and $\boldsymbol{H}_{2}$ are the real and imaginary
part of the matrix $\boldsymbol{H}$, respectively. After a standard decoupling
procedure one gets the equations

$\ddot{\boldsymbol{\varphi}}_{l}(t)+\boldsymbol{A}_{\mathfrak{q}}\,\dot{\boldsymbol{\varphi}}_{l}(t)+\boldsymbol{B}_{\mathfrak{q}}\,\boldsymbol{\varphi}_{l}(t)=\boldsymbol{0},$
(11)

valid for $l=1,2$, it is clear that even if both the real and imaginary part
of $\boldsymbol{\psi}$ verify the same equation, one cannot leave out one of
them because the solution of (11) implies the knowledge of the initial
conditions. The $dots$ notations refers to the time derivative and the real
matrices $\boldsymbol{A}_{\mathfrak{q}}$ and $\boldsymbol{B}_{\mathfrak{q}}$
depends on $(\boldsymbol{H}_{1},\boldsymbol{H}_{2})$ [6], the subscript refers
to the quantum origin of this matrices.

Let us now go to a classical system. We start with a system of linear
differential equations of second order, entirely similar to (11)

$\ddot{\boldsymbol{q}}(t)+\mathbf{A}\>\dot{\boldsymbol{q}}(t)+\mathbf{B}\>\boldsymbol{q}(t)=\boldsymbol{0}$
(12)

with $\boldsymbol{q}{\in}\mathbb{R}^{n}$ are the generalized coordinates, and
$\mathbf{A},\mathbf{B}{\in}\mathcal{M}_{n\times n}(\mathbb{R})$. Following [5,
6] this kind of equations can be performed through a classical electric
network. We are particularly interested in lumped element model circuits where
voltage and current depend only upon time.

The corresponding dynamics of an electric network is defined by the
appropriate use of the Kirchhoff rules that take care of the topology of the
network. We restrict our analysis to passive networks, where the energy
provided by an external source is non negative. The network has exactly
$n-$ports: pairs of terminals that allow to exchange energy with the
surrounding and have a given voltage and current. One has the possibility of
choosing the voltage or the current as the representative state variable of
the excitation or the response of the network. We call $\boldsymbol{v}$ the
vector corresponding to the port voltage and $\boldsymbol{i}$ the port
current, in [5, 6] we choose the port voltage $\boldsymbol{v}$ to implement a
particular electric circuit.

The analysis provided in [5] of the time evolution of electric circuits
results in a system of linear differential equations as (12). Using the
synthesis methods from [18, 19] of all electric networks in a given family
also developed in [5] the proposed topology of this kind of electric network
is represented in the figure 2. There are $n-$dipole networks of kind $L\|C$
tandem circuit interconnected through an network $\mathcal{N}$.

$\boldsymbol{\mathcal{N}}$ $\boldsymbol{\mathcal{N}_{1}}$
$\boldsymbol{\mathcal{N}_{2}}$ $\boldsymbol{\mathcal{N}_{k}}$
$\boldsymbol{\mathcal{N}_{n}}$

Figure 2: There are $n$ dipole networks, denoted by
$\\{\mathcal{N}_{1},{\cdots},\mathcal{N}_{k},{\cdots},\mathcal{N}_{n}\\}$
interconnected through an interaction network $\mathcal{N}$.

We consider without loss of generality that the $n-$port electric network
$\mathcal{N}$ admit an admittance representation with admittance matrix
$\boldsymbol{Y}$, an $L\|C$ tandem circuit is connected to each port, the
elements of $\boldsymbol{Y}$ are obtained from the general method of $n-$port
network synthesis method presented in [19]. Using the Laplace transformation,
moving from the temporal domain to the $s-$domain variable it is possible to
obtain the elements of $\boldsymbol{Y}$ in order to obtain a differential
equation similar to (12), i.e.
$\boldsymbol{Y}(s)=\boldsymbol{\alpha}+\frac{1}{s}\boldsymbol{\beta}$ contains
resistors and inductors only, in the constant matrices $\boldsymbol{\alpha}$
and $\boldsymbol{\beta}$, respectively.

A concrete construction of such circuits with their corresponding experimental
measurements can be found in [6]. But in general from the configuration given
by the figure 2 each port voltage $v_{k}$ works as the real, or imaginary,
part of the $k-$component of the wave function $\boldsymbol{\psi}$. From this
generalization the classical matrices in (12) are given by
$\boldsymbol{A}=\boldsymbol{C}^{-1}\boldsymbol{\alpha}$ and
$\boldsymbol{B}=\boldsymbol{\omega}^{2}_{0}+\boldsymbol{C}^{-1}\boldsymbol{\beta}$,
where $\boldsymbol{C}=diag(C_{1},\cdots,C_{n})$ and
$\boldsymbol{L}=diag(L_{1},\cdots,L_{n})$ are the diagonal matrices which
contain the capacitors and inductors of each of the $n-$dipole networks
$\\{\mathcal{N}_{k}\\}_{k=1,\cdots n}$. Thus the matrix
$\boldsymbol{\omega}^{2}_{0}:=(\boldsymbol{LC})^{-1}$ contains the proper
frequencies of each $L\|C$ port tandem connected, see figure 2.

##  V Conclusion and final observations

The aim of present work it was to show that there is a way to modify the
behavior of a known quantum system, in order to get information of another
quantum system that, at least, has a difficulty to be resolved directly.

In summary, we have shown how for a given pair of quantum systems, finite-
dimensional Hilbert spaces and its respective hamiltonian:
$(\mathscr{H},\boldsymbol{H})$ and
$(\mathscr{H}^{\prime},\boldsymbol{H}^{\prime})$ they could be linked via
gauge (local) transformations $\boldsymbol{\omega}$, that allow us to obtain
$\boldsymbol{H}^{\prime}$ from $\boldsymbol{H}$, via
$\boldsymbol{\Omega_{\omega}}$.

Respect to the simulation of a quantum system Q′ we search for some other
system that imitates the behavior of Q as well as possible. In others words,
we must perform a casting call of quantum systems or actors which can be very
limited, because it is a hard task to find another Q one to simulate Q′. We
wanted to use this equivalence between quantum systems to simulate another
quantum system connected with Q′. But when we said another, we want to say any
other quantum system which is connected with Q′ through $\boldsymbol{\omega}$.
The map $\boldsymbol{\Omega_{\omega}}$ applied to a given hamiltonian
$\boldsymbol{H}$ in (3) works as makeup that allows any actor Q, to simulate
the first quantum system Q′, a priori, if Q is connected with Q′ through
$\boldsymbol{\omega}$. Following the metaphor, the equivalence between this
quantum systems expands that catalogue of actors that can make a good
performance in order to mimic another quantum system and becoming that casting
call, a priori more efficient.

A final comment in this regard could be the implementation of the formal
equivalence between quantum and classical systems proposed in [6] in order to
simulate quantum systems through specific circuits. The advantage of using
such classical systems is that their controllability is simpler than for
quantum systems in general, in order to adequately guide its temporal
evolution. This implementation open the possibility to expands this catalogue
of actors capable of simulating the quantum system $Q$ even more with
classical actors, who usually do not play that role.

We are convinced that this work enriches the field of quantum simulations and
paves the way for new protocols in hybrid (digital-analog) simulation. But
also opens the way to new approach to solve quantum systems. Studying the
solution of a simple one in order to obtain the solutions for another quantum
system, eventually a more complicated one. This approach allows simulating
quantum systems establish a strong connection between them and so build a
dictionary for interpreting concepts of a theory on the other.

## Acknowledgments

We thank to FIDESOL for the support and recall also the anonymous readers for
their constructive criticism to this work.

## Competing interest

The author declares that there are no competing interests.

## References

  * [1] Feynman, R.P.; Simulation Physics with Computers, Int. J. Theor. Phys. 21, $467$ (1982).
  * [2] Georgescu, I.M., Ashhab, S., Nori, F.; Quantum simulation, Rev. Mod. Phys. 86, $153-185$ (2014).
  * [3] Lloyd, S.; Universal Quantum Simulators, Science 273, (5278) $1073-1078$ (1996).
  * [4] Lamata, L., Parra-Rodriguez, A., Sanz M. and Solano, E. Digital-analog quantum simulations with superconducting circuits, Adv. in Phys. X, 3, 1457981 (2018).
  * [5] Caruso M., Fanchiotti H., García Canal C. A.;  Equivalence between classical and quantum dynamics. Neutral kaons and electric circuits, Ann. Phys. 326, 10, $2717-2736$ (2011).
  * [6] Caruso M., Fanchiotti H., García Canal C. A., Mayosky M. and Veiga A.; The quantum CP-violating kaon system reproduced in the electronic laboratory, Proc. R. Soc. A.472, 20160615 (2016).
  * [7] Varadarajan, V.S.; Formal Reduction Theory of Meromorphic Differential Equations: A Group Theoretic view, Pacific. J. Math. Vol. 109, 1, 1$-$80 (1986).
  * [8] Varadarajan, V.S.; Linear Meromorphic Differential Equations: A Modern Point of View, Am. Math. Soc. Vol. 33, 1$-$42 (1996).
  * [9] Varadarajan, V.S.; Vector Bundles and Connections in Physics and Mathematics: Some Historical Remarks, Trends Math, 502$-$541 (2003).
  * [10] Mostafazadeh, A.; Quantum canonical transformations and exact solution of the Schrödinger equation, J. Math. Phys., 38, $3489-3496$ (1997).
  * [11] Mostafazadeh, A.; Time dependent diffeomorphisms as quantum canonical transformations and the time dependent harmonic oscillator, J. Phys. A, 31, $6495-6503$, 1998).
  * [12] Mostafazadeh, A.; Geometric phases, symmetries of dynamical invariants, and exact solution of the Schrodinger equation, J. Phys. A, 34, $6325-6338$ (2001)
  * [13] Magnus, W.; On the exponential solution of differential equations for a linear operator. Comm. Pure and Appl. Math. VII 4, $649-673$ (1954).
  * [14] Giles, B.and Selinger, P.; Exact synthesis of multiqubit Clifford+$T$ circuits, Phys. Rev. A, 87, $032332$ (2013).
  * [15] J. L. Rosner, Am. J. Phys.64 (8), 982-985 (1996).
  * [16] Arnold, V.I., Ordinary differential equations, Springer-Verlag, (1992).
  * [17] V.I. Arnold, _Geometrical Methods in the Theory of Ordinary Differential Equations_ , Springer (1980)
  * [18] N. Balabanian, T.A. Bickart, _Linear network theory: analysis, properties, design and synthesis_ , Willey (1969)
  * [19] H. Carlin, A. Giordano, _Network Theory: An Introduction to Reciprocal and Nonreciprocal Circuits_ , Prentice Hall (1964)