present address: ]Physik-Institut, Universität Zürich, Winterthurerstrasse
190, CH-8057 Zürich, Switzerland

# Engineering Higher Order Van Hove singularities in two dimensions: the
example of the surface layer of Sr2RuO4

Anirudh Chandrasekaran  Department of Physics and Centre for the Science of
Materials, Loughborough University, Loughborough LE11 3TU, United Kingdom
Luke C. Rhodes  SUPA, School of Physics and Astronomy, University of St
Andrews, North Haugh, St Andrews, KY16 9SS, United Kingdom  Edgar Abarca
Morales  SUPA, School of Physics and Astronomy, University of St Andrews,
North Haugh, St Andrews, KY16 9SS, United Kingdom  Max Planck Institute for
Chemical Physics of Solids, Noethnitzer Strasse 40, 01187 Dresden, Germany
Carolina A. Marques  [  SUPA, School of Physics and Astronomy, University of
St Andrews, North Haugh, St Andrews, KY16 9SS, United Kingdom  Phil D.C. King
SUPA, School of Physics and Astronomy, University of St Andrews, North Haugh,
St Andrews, KY16 9SS, United Kingdom  Peter Wahl  SUPA, School of Physics and
Astronomy, University of St Andrews, North Haugh, St Andrews, KY16 9SS, United
Kingdom  Physikalisches Institut, Universität Bonn, Nussallee 12, 53115 Bonn,
Germany  Joseph J. Betouras  Department of Physics and Centre for the Science
of Materials, Loughborough University, Loughborough LE11 3TU, United Kingdom

###### Abstract

The properties of correlated electron materials are often intricately linked
to Van Hove singularities (VHs) in the vicinity of the Fermi energy. The class
of these VHs is of great importance, with higher order ones – with power-law
divergence in the density of states – leaving frequently distinct signatures
in physical properties. We use a new theoretical method to detect and analyse
higher order VHs (HOVHs) in two-dimensional materials and apply it to the
electronic structure of the surface layer of Sr2RuO4. We then constrain a low
energy model of the VHs of the surface layer of Sr2RuO4 against angle-resolved
photoemission spectroscopy and quasiparticle interference data to analyse the
VHs near the Fermi level. We show how these VHs can be engineered into HOVHs.

Introduction. The properties of strongly correlated electron materials can be
extremely sensitive to modest external stimuli, such as magnetic field,
pressure or uniaxial strain, often exhibiting complex phase diagrams as a
function of these tuning parameters Binz and Sigrist (2004); Hackl and Vojta
(2011); Shtyk et al. (2017); Steppke et al. (2017). In many cases, the
sensitivity of their properties is traced back to van Hove singularities (VHs)
in the vicinity of the Fermi energy ($E_{\mathrm{F}}$) Efremov et al. (2019);
Sherkunov and Betouras (2018); Yuan et al. (2019); Sunko et al. (2019). It has
been realised that such VHs can also be of a higher-order type (HOVHs), where
both the gradient and the Hessian determinant of the dispersion relation
$\varepsilon(k_{x},k_{y})$ vanish and unique thermodynamic signatures can be
observedEfremov et al. (2019); Chandrasekaran and Betouras (2022). A detailed
classification scheme of HOVHs is provided in Refs. Chandrasekaran et al.
(2020). Two of the authors have recently developed a new method to detect and
analyse a HOVHs from arbitrary two dimensional electronic structure models
Chandrasekaran and Betouras (2023); com (a), enabling us to uniquely
characterise the important band structure features in two dimensional systems.

Extending the method to real materials requires a detailed understanding of
the material’s low-energy electronic structure. Here, we focus on the surface
of Sr2RuO4, which provides an ideal test system to explore how the symmetry
and order of the VHs influences the electronic properties of multi-band and
spin-orbit coupled (SOC) systems. Indeed, strontium ruthenates provide an
ideal benchmark material class to explore the impact of VHs on macroscopic
properties Binz and Sigrist (2004). They host rich phase diagrams, with
metamagnetic phases, a putative quantum critical point in Sr3Ru2O7 Grigera et
al. (2001) and evidence for unconventional superconductivity in Sr2RuO4Maeno
et al. (1994). In both systems, there is strong evidence for the importance of
VHs in their electronic properties Allan et al. (2013); Shen et al. (2007);
Marques et al. (2021, 2022); Sunko et al. (2019); Marques et al. (2023). The
superconductivity in Sr2RuO4, for example, has been found to react sensitively
to uniaxial strain Steppke et al. (2017) which directly correlates with a VHs
crossing $E_{\mathrm{F}}$ Sunko et al. (2019), while in Sr3Ru2O7, an unusual
field dependence of the specific heat has been explained from Zeeman splitting
of HOVHs Efremov et al. (2019).

The surface layer of Sr2RuO4 provides a two-dimensional electronic system on
which high-resolution spectroscopic data is availableDamascelli et al. (2000);
Tamai et al. (2019); Wang et al. (2017); Marques et al. (2021), facilitating a
full assessment of the nature of its VHs. A $\sqrt{2}$x$\sqrt{2}$ structural
reconstruction due to octahedral rotations Damascelli et al. (2000); Matzdorf
et al. (2002); Marques et al. (2021); Morales et al. (2023) results in a unit
cell doubling that produces a four-fold saddle point in the electronic
structure at the Brillouin zone (BZ) corner. This type of saddle point is an
essential ingredient for the formation of a HOVHs (of classification type
$X_{9}$ Chandrasekaran et al. (2020)) in Sr3Ru2O7Efremov et al. (2019) with
power-law divergence in the density of states (DoS), distinct from the
logarithmic divergence of the $A_{1}$ saddle point at the BZ face in the bulk
of Sr2RuO4.

In this work, we apply analytical and numerical tools to study how the
octahedral rotation influences the formation of HOVHs in the surface layer of
Sr2RuO4. We first study the evolution of the VHs from density functional
theory (DFT) calculations and establish the impact of octahedral rotations on
the VHs, before constraining a minimal model directly against experimental
data from angle-resolved photoemission spectroscopy (ARPES) and scanning
tunnelling microscopy (STM). By combining these complementary data sets, we
are able to progressively refine the low-energy electronic structure and
obtain an accurate description of the dispersion relations in the vicinity of
the saddle point. From the experimentally-derived model, we classify multiple
VHs close to $E_{\mathrm{F}}$ and study how their order can be tuned,
providing new design principles for VHs engineering in two-dimensional
materials.

Figure 1: VHs in ruthenates. (a) top-down view of a single layer of Sr2RuO4
highlighting the surface-induced octahedral rotation $\theta$. Grey square
shows the unit cell without octahedral rotations, the blue shaded square
highlights the expanded unit cell with rotations. (b) Evolution of the band-
structure as a function of $\theta$. Darker lines correspond to larger Ru-O-Ru
rotation angles. (c) BZ of Sr2RuO4 with octahedral rotations (dashed) and
without (solid). (d) Fermi contours and 3D image of an $A_{3}$ HOVHs (and a
$\pi/2$ rotated copy required by symmetry) that can occur at the M-point due
to the reconstructed BZ.

Octahedral rotation-induced order of VHs. To analyse the impact of octahedral
rotation on the electronic structure of Sr2RuO4, we perform DFT calculations
of a monolayer of Sr2RuO4, shown in Fig. 1(a), with induced rotation of the
RuO6 octahedra, keeping the Ru-Ru distances constant. This produces a two-atom
unit cell (Fig. 1(a)). We then calculate the electronic structure for each
rotation angle and project onto a tight-binding model (TBM) Pizzi et al.
(2020). There are two dominant changes as a function of octahedral rotation.
First is the evolution of the saddle point VHs that lies close to
$E_{\mathrm{F}}$ at the $\bar{M}$ point. This band undergoes a Lifshitz
transition at an octahedral rotation angle of $\sim 5^{\circ}$ and evolves
from a concave dispersion to a convex dispersion as rotation angle increases.
The second feature is the evolution and Lifshitz transition of an electron
pocket at the $\bar{\Gamma}$ point. The energy of this band minimum is,
however, known to be located too low in energy in DFT, with ARPES measurements
placing it just above $E_{\mathrm{F}}$ Damascelli et al. (2000); Morales et
al. (2023). This only produces a linear offset to the DoS of the VHs around
$\bar{M}$ without affecting the order and symmetry of the relevant VHs around
$\bar{M}$ points. Here we focus only on the electronic structure of the low-
energy VHs around the $\bar{M}$ points com (b).

Figure 2: (a) Band-structure in the vicinity of the M point showing the
orbital character of VHs at a rotation angle $\theta=9^{\circ}$. Colours
encode orbital character (red: $d_{xz}$, green: $d_{yz}$, blue: $d_{xy}$). The
energy $E_{\mathrm{VHs}}$ of the $d_{xy}$ VHs is indicated by a horizontal
dashed line. (b) DoS for the band structure shown in (a), with the VHs on the
$d_{xy}$ band. (c) Type of the $d_{xy}$ VHs determined from the series
expansion Chandrasekaran et al. (2020) as function of angle and energy from
the VHs $E-E_{\mathrm{VHS}}$. At an angle of $\theta\sim 9^{\circ}$, the
behaviour exhibits a singular point at $E_{\mathrm{VHS}}$ where the VHs shows
power-law divergence with exponent $-\frac{1}{4}$. (d) Numerically determined
order of the tail of the VHs from $\frac{d\log\rho}{d\log E}$ for the
$d_{xy}$-derived VHs as function of angle and energy. The energy scale is
relative to the energy $E_{\mathrm{VHs}}$ at which the VHs occurs. The
numerically determined order is not reliable in the grey area, because it is
too close to the VHs. The dashed line shows the superimposed boundary between
logarithmic and polynomial order of the tail of the VHs determined
analytically.

Fig. 1(d)) and Figs 2(a) and (b) show the electronic structure close the
$\bar{M}$ point and throughout the surface BZ, respectively, for a $9^{\circ}$
octahedral rotation. The octahedral rotation significantly changes the
curvature of the electronic structure from convex to concave in the vicinity
of $\bar{M}$ point. This continues tuning of the band dispersion leads to
favorable conditions for the formation of HOVHs. By applying our method
Chandrasekaran and Betouras (2023) to extract the order and symmetry of the
VHs around the $\bar{M}$ point, we find that close to $\theta\sim 9^{\circ}$,
quadratic terms in the dispersion relation are suppressed in the
$\bar{\Gamma}-\bar{M}$ direction and the analytic expansion predicts a HOVHs
with an energy dependence in the tail of the DoS $\rho(E)\propto E^{-1/4}$,
classified as $A_{3}$-type Chandrasekaran et al. (2020); Yuan and Fu (2020).
For octahedral rotation angles below/above $9^{\circ}$, the dispersion is
dominated by quadratic terms resulting in a VHs with a logarithmic divergence.
Fig. 2(c) shows this phase diagram, as a function of energy and octahedral
rotation angle.

The result demonstrates that modification of the octahedral rotation provides
a method to engineer HOVHs in Sr2RuO4. However, in multi-band systems there is
a possibility that other electronic states could mask the clean divergent
signatures induced by HOVHs in thermodynamic properties. To check this
possibility and establish the energy range over which the HOVHs impacts
experimental observables, we numerically integrate the DoS $\rho(E)$
calculated using the TBM. The exponent of the VHs can then be determined from
the logarithmic derivative of $\rho(E)$,
$\frac{d\log\rho(E-E_{\mathrm{VHs}})}{d\log E}$, which provides the order of
the leading polynomial term in the tail of the VHs. The numerical analysis,
shown in Fig. 2(d), confirms that the DoS tail from the VHs changes
significantly as a function of rotation angle. While for small octahedral
rotation angles we observe a behaviour consistent with a logarithmic
divergence, as the angle $\theta$ approaches $9^{\circ}$, the VHs acquires an
exponent $-1/4$, consistent with the analytic theory.

Figure 3: Electronic structure from ARPES and fitting of TBM. (a) Band
dispersions of the bulk and surface electronic structure of Sr2RuO4 measured
by ARPES along the $\overline{\Gamma}-\overline{\mathrm{M}}$ direction of the
surface BZ, reproduced from Ref. Morales et al. (2023). (b) Zoom-in of the
measured electronic structure near $E_{\mathrm{F}}$, with the fitted TBM shown
atop (red lines). The best-fit parameters are given in com (c). (c) The TBM is
fitted to the bands extracted from the measurements shown in (b) and from
additional measurements (Ref. Morales et al. (2023)). (d) Corresponding Fermi
surface measured by ARPES and calculated from our TBM.

Fitting to ARPES. Having identified that HOVHs can be stabilised via
octahedral rotations in monolayer Sr2RuO4, we turn our attention to a system
that exhibits a series of VHs close to EF with the potential to tune them,
controlling their order. The surface of Sr2RuO4 is known to exhibit an
octahedral rotation, with low energy electron diffraction measurements placing
the angle at $\theta=8.5\pm 2.5^{\circ}$ Matzdorf et al. (2002), close to the
critical angle of 9∘ at which the $\bar{M}$-point VHs becomes higher order. To
establish the low energy model that describes the experimental electronic
structure in the vicinity of the VHs at the $\bar{M}$ point, we perform a
polynomial interpolation of the DFT-derived TBMs as a function of angle and a
least squares fitting to ARPES data of the surface electronic structure
Morales et al. (2023), presented in Fig. 3(a) com (c). We plot the obtained
TBM over the ARPES data in Fig. 3(b), and show this together with the
extracted data points from the ARPES measurements used for the fitting in Fig.
3(c). The model shows excellent agreement with the ARPES data (Fig. 3(d)), and
the octahedral rotation angle 8.03∘ is in agreement with measurements of the
surface structure Matzdorf et al. (2002).

Figure 4: Low-energy electronic structure from QPI. (a) Comparison of
experimental and calculated tunneling spectrum of the clean surface, for the
models fit to ARPES (grey, dotted) and to QPI (blue). (b) Left half: Quasi-
particle interference image at $E=2.8\mathrm{mV}$, right half: calculations
using model fitted to experimental data. (c-d) cuts through the atomic peaks
along $q_{y}$ (left) and $q_{x}$ (right) from experiment (c) and calculations
for the TBM fitted to the QPI (d). (e-f) Real-space image showing the
checkerboard pattern at the Sr positions, from experiment (e) and the model
(f). Blue and black spheres indicate the position of the Sr and Ru atoms,
respectively.

Comparison with QPI. The ARPES data reveal an important complication to the
low-energy electronic structure of the surface of Sr2RuO4. Whilst they allowed
us to locate the saddle point at around 18 meV below $E_{\mathrm{F}}$, the
bands in the vicinity of $E_{\mathrm{F}}$ are found to be dominated by a spin-
orbit gap in the electronic structure which can also give rise to VHs in the
DoS Morales et al. (2023). To gain further insight into the VHs mediated by
this spin-orbit gap, we refine our experimental model of the VHs by comparing
with QPI data obtained by STM at millikelvin temperatures and with sub-meV
energy resolution. QPI provides access to the electronic states above and
below $E_{\mathrm{F}}$. Low-energy low-temperature differential conductance
measurements of the surface of Sr2RuO4 were previously reported Marques et al.
(2021); Kreisel et al. (2021a) and established the existence of a partial gap
with four peaks in the vicinity of $E_{\mathrm{F}}$, indicating the presence
of four close-lying VHs. However, when we compare a theoretical continuum LDoS
calculation of the tunneling spectrum $g(V)$ from our ARPES-fit model to the
STM measurements (shown in grey and black, respectively, in Fig. 4(a)), we
find that although the model captures the partial gap, it only predicts two of
the observed four peaks, and the position of the partial gap is several meV
too low. The QPI data exhibits a small nematic term breaking the four-fold
rotation symmetry of order $\sim 1$ meV Marques et al. (2021); Kreisel et al.
(2021a). Incorporating this nematicity via an anisotropy of the nearest-
neighbour hopping into our model, and adjusting the chemical potential by
$\sim 3\mathrm{meV}$, we reproduce the experimental tunnelling conductance
with sub-meV resolution. To confirm the accuracy of our model, we perform
continuum QPI calculations, modelling scattering of electrons from a point-
like defect, and compare with experimental data Kreisel et al. (2021a) (see
supplementary for details). The calculation is in quantitative agreement with
the experimental data within a $\pm$5 meV window. Fig. 4(b) shows a comparison
between calculation and experiment of a representative constant energy cut,
while Figs.4(c,d) displays cuts as a function of energy.

We find that the formation of the gap-like structure close to the EF can be
understood completely as due to the near-$E_{\mathrm{F}}$ spin-orbit gap found
in ARPES. By calculating the atomic corrugation of the real space LDOS using
our new model, we find that SOC, coupled with a finite nematicity, produces
the experimentally observed checkerboard modulation (Fig. 4(e,f))Marques et
al. (2021), eliminating the need to introduce a checkerboard term into the TBM
as suggested previouslyMarques et al. (2021); Kreisel et al. (2021a). Thus our
model unifies ARPES and STM measurements of the low-energy electronic
structure of the surface of Sr2RuO4 with sub-meV energy resolution.

Figure 5: Analysis of spin-orbit induced VHs in the surface layer of Sr2RuO4.
(a) Band structure and (b) DoS around $E_{\mathrm{F}}$ for the TBM fitted to
QPI. Light grey lines in panel (a) show the band structure fitted to ARPES.
(c, d) phase diagrams of the VHs closest to the EF at the upper and lower edge
of the SOC-induced gap, marked by circles in (a). Solid lines indicate the
line where the VHs becomes an $A_{3}$ singularity. Similar results hold for
the set of singularities close to $\bar{M}_{1}$.

Order of the VHs. Our observation of a SOC-induced gap at EF raises the
question what the type of the VHs at the gap edge is. Fig. 5(a,b) presents the
band structure, and corresponding DoS, of our final optimised TBM. This model
is quantitatively accurate around the $\bar{M}$ point, now labelled by two
subscripts due to the small nematic splitting rendering the $M$ points
inequivalent along the Ru-Ru reciprocal lattice directions.

We proceed to characterise the VHs found within 20 meV of $E_{\mathrm{F}}$
around the $\bar{M}$ points taking into account the nematicity. The nematicity
changes the critical angle where the HOVHs with the divergent exponent $-1/4$
(and class still A3) is located. The critical line as a function of the
percentage change of the nematicity, for values experimentally achievable
through uniaxial strain, is computed (Fig. 5(c,d)). This line separates the
regions of maximum/minimum and saddle point VHs, therefore the system can be
tuned to an A3 HOVHs by applying strain as a tuning parameter.

Discussion. In this work, we analyse two sets of complementary experimental
data (STM and ARPES) of the benchmark material Sr2RuO4. We are able to
classify and understand the nature of HOVHs close to the Fermi surface. The
work refines a long-standing question on the appearance of HOVHs in Sr2RuO4,
as evidenced through ARPES measurements very early Lu et al. (1996); Yokoya et
al. (1996), before surface and bulk contributions were distinguished
Damascelli et al. (2000).

The combined theoretical and experimental characterization provides two key
results: (i) the model for the low-energy electronic structure, consistent
with STM and ARPES, captures related to the VHs. It provides an accurate
description of the gap-like structure seen in tunneling spectra that were
previously interpreted in terms of inelastic tunneling or correlation effects
Wang et al. (2017) or a hybridization gap Marques et al. (2021); Kreisel et
al. (2021b). This gap can be explained straightforwardly as the combined
result of SOC and octahedral rotation, leading to a suppression of DoS by
almost $50\%$. Our low-energy model promises an understanding of the, yet
unexplained, suppression of superconductivity in the surface layer of
Sr2RuO4Barker et al. (2003); Kambara et al. (2006); Marques et al. (2021).
(ii) From the series of TBMs used to fit the experimental data, we can deduce
‘design criteria’ for stabilizing HOVHs in Sr2RuO4.

To this end, the introduced set of analytical and numerical tools enables the
determination of the order of the HOVHs, leading to a phase diagram of their
order. The $d_{xy}$-VHs clearly changes it type as a function of rotation
angle, with a critical angle of $\sim 9.7^{\circ}$, surprisingly close to the
angle realized in the surface reconstruction of Sr2RuO4. Tunability of the
rotation angle can be achieved by biaxial strain, as, for example,
demonstrated in MBE-grown thin films Burganov et al. (2016). In combination
with doping or gating to control the chemical potential, thin films of
strontium ruthenates become a promising testbed to establish the influence and
importance of HOVHs for the physical properties. The theoretical and numerical
tools we introduce here provide new methods to study tunable VHs in a range of
existing quantum materials (e.g. Liu et al. (2023)) and in many more predicted
to exhibit almost flat bands Regnault et al. (2022). Even weak interactions
may lead to a plethora of quantum phases as a result of HOVHs Zervou et al.
(2023); Classen et al. (2020).

##  I Acknowledgements

We thank Andy Mackenzie and Andreas Rost for useful discussions. We are
grateful to the UK Engineering and Physical Sciences Resarch Council for
funding via Grant Nos EP/T034351/1 and EP/T02108X/1. This work used
computational resources of the Cirrus UK National Tier-2 HPC Service at EPCC
(http://www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC
(EP/P020267/1) and of the High-Performance Computing cluster Kennedy of the
University of St Andrews.

##  II Supplemental Material

###  II.1 Density Functional Theory Calculations

DFT calculations were performed using the augmented plane-wave code Quantum
Espresso Giannozzi et al. (2017) using the PBE for the exchange-correlation
functional. We used a monolayer of Sr2RuO4 with 15 $\AA$ of vacuum on either
side for each calculation with the RuO6 octahedra rigidly rotated by angles
between 0 and 12∘ to study the role of octahedral rotation on the electronic
structure. The energy cutoff for the wavefunction and charge density was
chosen at 90 Ry and 720 Ry respectively and a $\mathbf{k}$-grid of $8\times
8\times 1$ was used for each calculation.

###  II.2 Tight-binding models (TBMs)

These DFT band structures were then projected onto an orthonormal tight
binding basis consisting of the Ru $d_{xz}$ $d_{yz}$ and $d_{xy}$ orbitals in
a two-atom unit cell. This was achieved using a modified version of Wannier90
Pizzi et al. (2020) where we fix the phase of the Wannier functions against
the first orbital to preserve the sign of the wave functions between
neighbouring Ru atoms. This Wannierisation was performed on a $7\times 7\times
1$ $\mathbf{k}$-grid.

For fitting the ARPES band structure using the TBM, we perform a polynomial
interpolation on the hopping strengths of the thirteen individual models
obtained for angles $\theta=0^{\circ}$ through $12^{\circ}$ enabling fitting
of $\theta$ as a continuous variable. We add a local spin orbit term
$H_{\mathrm{SOC}}=\lambda\mathbf{L}\cdot\mathbf{S}$. Lastly, we allow for an
overall band renormalisation and shift of the chemical potential. Fitting the
TBM to the ARPES band structure therefore implies optimizing the following
free parameters: the SOC strength $\lambda$, RuO octahedral rotation angle
$\theta$, overall band renormalisation $Z$ and the chemical potential $\mu$.

For comparison with the QPI results, we introduce a nematicity of $0.5\%$
across all orbitals, multiplying the nearest neighbour hoppings in one
direction across all orbitals by $1.01$ and dividing it in the orthogonal
direction by $1.01$.

###  II.3 Continuum Local DoS calculations

For continuum local DoS (LDoS) calculations shown in fig. 4, we use the TBM
obtained from the fit to the ARPES data modified to be consistent with the STM
differential conductance. We use the localized wave functions obtained at an
angle close to the one obtained from the fit. Because the tunneling spectrum
is significantly better described with a slightly larger octahedral rotation
angle than the $8^{\circ}$ obtained from the fit, we use as a basis the model
and wave functions with $\theta=9^{\circ}$. The continuum LDoS calculations
are performed using the method introduced in Choubey et al. (2014); Kreisel et
al. (2015); Choubey et al. (2017) and recently applied to ruthenatesKreisel et
al. (2021a); Benedičič et al. (2022) using the St Andrews calcqpi
codeBenedičič et al. (2022); Rhodes et al. (2023); Marques et al. (2023);
Naritsuka et al. (2023). Calculations shown in fig. 4 were done on a
$\mathbf{k}$-space grid of $2048\times 2048$ points distributed uniformly
across the BZ of the two-atom unit cell. The real-space continuum LDoS is
calculated over an area corresponding to $128\times 128$ unit cells and with
an energy broadening of $0.2\mathrm{meV}$. We assume isotropic scattering with
a scattering potential of $V=1\mathrm{eV}$ across all orbitals and spin
directions.

###  II.4 Numerical Integration of DoS

We calculate the DoS $\rho(E)$ from the unperturbed Green’s function using
$\rho(E)=-\mathrm{Im}G_{0}(E)$ with an energy broadening of
$\Gamma=0.75\mathrm{meV}$, and using a $k$-space grid with $2088\times 2088$
over the first BZ. The DoS is calculated on a regular grid with an energy
spacing of $0.4\mathrm{meV}$. To obtain the logarithmic derivative, we use a
Savitzky-Golay algorithm with eleven points.

###  II.5 Tuning to a HOVHs

We briefly describe the procedure to tune the VHs in the TBMs into HOVHs. We
first choose a band and a critical point of interest to us. Let us denote them
by $n$ and $\mathbf{k}^{(0)}$ respectively. If $\mathbf{k}^{(0)}$ is a
critical point, then the gradient of the dispersion of band $n$ vanishes, i.e
$\gradient\varepsilon_{n}^{\,}(\mathbf{k})=0$ at
$\mathbf{k}=\mathbf{k}^{(0)}$. If the Hessian determinant
$\det[\partial_{i}^{\,}\partial_{j}^{\,}\varepsilon_{n}^{\,}(\mathbf{k})]$
also vanishes at $\mathbf{k}^{(0)}$, we then have a HOVHs. Now the tight-
binding Hamiltonian (and therefore the dispersion) also depend on the tuning
parameters which we lump together as a vector $\mathbf{t}$. Therefore, the
Hessian at any given $\mathbf{k}$ is a function of $\mathbf{t}$ as well. Thus,
if we set the Hessian determinant evaluated at $\mathbf{k}^{(0)}$ to zero, we
can solve for the critical value of $\mathbf{t}$ that gives us a HOVHs at
$\mathbf{k}^{(0)}$. To diagnose the type of the HOVHs, we will need the series
expansion of the dispersion to a higher degree than just quadratic order
(which is sufficient for computing the Hessian).

Now this procedure is fairly straightforward in theory: we simply have to
analytically diagonalise the band Hamiltonian $H(\mathbf{k})$ and obtain the
$n^{\text{th}}$ band-dispersion as a function of $\mathbf{k}$, denoted by
$\varepsilon_{n}^{\,}(\mathbf{k})$. We can then calculate the partial
derivatives of $\varepsilon_{n}^{\,}(\mathbf{k})$ or equivalently compute its
series expansion at the chosen $\mathbf{k}^{(0)}$ to any desired degree.
However, it is usually possible to implement this analytic procedure only for
two band Hamiltonians. Multi-band, detailed Hamiltonians necessitate the use
of an alternate, computationally implementable strategy for obtaining the
series expansion (and the partial derivatives).

The crux of this alternate method lies in the realisation that we only need to
compute the _numerical_ values of the partial derivatives
$\partial^{l+m}\varepsilon_{n}^{\,}(\mathbf{k})/\partial^{l}k_{x}^{\,}\partial^{m}k_{y}^{\,}$
evaluated _at_ $\mathbf{k}^{(0)}$ (for $1\leqslant l+m\leqslant N$), in order
to construct the Taylor expansion to degree $N$. In the one dimensional case,
the first derivative can be computed using the Feynman - Hellmann lemma which
reads

$\frac{d\varepsilon_{n}^{\,}(k)}{dk}\bigg{|}_{k^{(0)}}^{\,}=\left\langle
n,k^{(0)}\left|\frac{dH(k)}{dk}\bigg{|}_{k^{(0)}}\right|n,k^{(0)}\right\rangle.$
(1)

Note that this requires us to be able to 1) numerically obtain the eigenvector
$|n,k^{(0)}\rangle$ and 2) differentiate the Hamiltonian matrix
$H(\mathbf{k})$, both of which can be readily done, even in the case of large
multi-band Hamiltonians. Therefore, to compute the full Taylor expansion in
the two and three dimensional cases, we have to work with an appropriate
generalisation of the Feynman - Hellmann lemma which can also handle
complicated multi-band degeneracies that often prop up at high-symmetry
points. A comprehensive treatment of such an extension can be found in Ref
Chandrasekaran and Betouras (2023). The final formula for the Taylor expansion
up to second order reads

$\displaystyle\varepsilon_{n}^{\,}$
$\displaystyle(\mathbf{k}^{(0)}+\mathbf{k})\approx\left\langle
n,\mathbf{k}^{(0)}\left|\frac{\partial
H(\mathbf{k}^{(0)}+\lambda\,\mathbf{k})}{\partial\lambda}\bigg{|}_{\lambda=0}\right|n,\mathbf{k}^{(0)}\right\rangle$
$\displaystyle+\left\langle
n,\mathbf{k}^{(0)}\left|\frac{\partial^{2}H(\mathbf{k}^{(0)}+\lambda\,\mathbf{k})}{\partial\lambda^{2}}\bigg{|}_{\lambda=0}\right|n,\mathbf{k}^{(0)}\right\rangle$
$\displaystyle+\frac{1}{2}\sum\limits_{m}\frac{\left|\left\langle
n,\mathbf{k}^{(0)}\left|\frac{\partial
H(\mathbf{k}^{(0)}+\lambda\,\mathbf{k})}{\partial\lambda}\bigg{|}_{\lambda=0}\right|m,\mathbf{k}^{(0)}\right\rangle\right|^{2}}{\varepsilon_{n}^{\,}-\varepsilon_{m}^{\,}}+\mathcal{O}\left(k^{3}\right),$

where $\mathbf{k}=(k_{x}^{\,},k_{y}^{\,})$ and the sum is over other the
eigenvectors $|m,\mathbf{k}^{(0)}\rangle$ that are not degenerate to
$|n,\mathbf{k}^{(0)}\rangle$. Note that only $\mathbf{k}$ and $\lambda$ are
symbolic variables while $\mathbf{k}^{(0)}$, $|n,\mathbf{k}^{(0)}\rangle$ and
$|m,\mathbf{k}^{(0)}\rangle$ are numerical vectors. This is what makes the
method viable for large and complicated models. By using the auxiliary
parameter $\lambda$ in the formula above, we will directly obtain the series
expansion as a multivariate polynomial in $k_{x}^{\,}$ and $k_{y}^{\,}$. The
partial derivatives can be extracted from the Taylor coefficients.

Furthermore, we can also readily track the effect of tuning parameters on the
low-energy expansion by using
$H(\mathbf{k}^{(0)}+\lambda\,\mathbf{k},\mathbf{t}^{(0)}+\lambda\,\delta\mathbf{t})$
in place of $H(\mathbf{k}^{(0)}+\lambda\,\mathbf{k})$ in the formula above,
where $\mathbf{t}^{(0)}$ is the set of tuning parameters (their chosen,
numerical values to be precise). The perturbations to $\mathbf{t}^{(0)}$ will
be incorporated through the symbolic vector $\delta\mathbf{t}=(\delta
t_{1}^{\,},\delta t_{2}^{\,},\cdots)$. In the present context, the tuning
parameters of interest are the RuO octahedral rotation angle $\theta$ and the
nematicity $\Delta_{\text{nem}}^{\,}$.

In this formulation, the numerical Taylor coefficients at the chosen
$(\mathbf{k}^{(0)},\theta^{(0)},\Delta_{\text{nem}}^{(0)})$ will be perturbed
by powers of $\delta\theta$, $\delta\Delta_{\text{nem}}^{\,}$. We can then
compute the Hessian determinant at
$(\mathbf{k}^{(0)},\theta^{(0)},\Delta_{\text{nem}}^{(0)})$. This will be a
polynomial in $\delta\theta$ and $\delta\Delta_{\text{nem}}^{\,}$. We set it
to zero to solve for $\delta\theta^{*}$ and $\delta\Delta_{\text{nem}}^{*}$.
The critical values of the tuning parameters that yield a HOVHs are then given
by $\theta^{(\text{c})}\approx\theta^{(0)}+\delta\theta^{*}$ and
$\Delta_{\text{nem}}^{(\text{c})}\approx\Delta_{\text{nem}}^{(0)}+\delta\Delta_{\text{nem}}^{*}$.
This process can be repeated a few times to obtain a more accurate value of
the tuning parameters that lead to a HOVHs at $\mathbf{k}^{(0)}$.

The strategy outlined above works quite well at high-symmetry points like the
$\overline{M}$ wherein the location of the critical point remains fixed even
as we change the tuning parameters. For the SOC gap VHs, we have to adopt a
slightly different strategy. In this case, the location of the critical point
$\mathbf{k}^{(0)}$ changes as we change $\theta$ and
$\Delta_{\text{nem}}^{\,}$. For $\Delta_{\text{nem}}^{\,}$ in the range of 0%
to 4%, we find that the Hessian determinant for the critical points in the SOC
gap changes sign in the range of $8^{\degree}$ to $12^{\degree}$. So for any
given $\Delta_{\text{nem}}^{\,}$, we perform binary search to find the value
of $\theta$ at which the Hessian determinant vanishes. At any intermediate
$\theta$, we have to compute the new location of $\mathbf{k}^{(0)}$. We start
at an approximate location and Taylor expand the band to cubic or quartic
order. This gives a polynomial approximation for
$\varepsilon_{n}^{\,}(\mathbf{k})$ near $\mathbf{k}^{(0)}$. We can compute the
gradient of this, set it to zero and solve the resultant system of
simultaneous polynomial equations numerically. Choosing the smallest real
solution, we can arrive at a more accurate estimate for $\mathbf{k}^{(0)}$. We
repeat this a few times to attain convergence and proceed with computing the
Hessian determinant at $\mathbf{k}^{(0)}$.

## References

  * Binz and Sigrist (2004) B. Binz and M. Sigrist, EPL (Europhysics Letters) 65, 816 (2004).
  * Hackl and Vojta (2011) A. Hackl and M. Vojta, Physical Review Letters 106, 137002 (2011).
  * Shtyk et al. (2017) A. Shtyk, G. Goldstein, and C. Chamon, Phys. Rev. B 95, 035137 (2017), ISSN 2469-9950, 2469-9969, URL https://link.aps.org/doi/10.1103/PhysRevB.95.035137.
  * Steppke et al. (2017) A. Steppke, L. Zhao, M. E. Barber, T. Scaffidi, F. Jerzembeck, H. Rosner, A. S. Gibbs, Y. Maeno, S. H. Simon, A. P. Mackenzie, et al., Science 355, eaaf9398 (2017).
  * Efremov et al. (2019) D. V. Efremov, A. Shtyk, A. W. Rost, C. Chamon, A. P. Mackenzie, and J. J. Betouras, Physical Review Letters 123, 207202 (2019).
  * Sherkunov and Betouras (2018) Y. Sherkunov and J. J. Betouras, Phys. Rev. B 98, 205151 (2018).
  * Yuan et al. (2019) N. F. Q. Yuan, H. Isobe, and L. Fu, Nat Commun 10, 5769 (2019), ISSN 2041-1723, URL http://www.nature.com/articles/s41467-019-13670-9.
  * Sunko et al. (2019) V. Sunko, E. A. Morales, I. Marković, M. E. Barber, D. Milosavljević, F. Mazzola, D. A. Sokolov, N. Kikugawa, C. Cacho, P. Dudin, et al., npj Quantum Materials 4, 46 (2019).
  * Chandrasekaran and Betouras (2022) A. Chandrasekaran and J. J. Betouras, Phys. Rev. B 105, 075144 (2022).
  * Chandrasekaran et al. (2020) A. Chandrasekaran, A. Shtyk, J. J. Betouras, and C. Chamon, Phys. Rev. Research 2, 013355 (2020), ISSN 2643-1564, URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013355.
  * Chandrasekaran and Betouras (2023) A. Chandrasekaran and J. J. Betouras, Advanced Physics Research 2, 2200061 (2023), ISSN 2751-1200, 2751-1200, URL https://onlinelibrary.wiley.com/doi/10.1002/apxr.202200061.
  * com (a) The method expands the dispersion relation directly in the vicinity of the VHs, and utilises the Feynmann-Hellmann theorem to evaluate higher-order derivatives.
  * Grigera et al. (2001) S. Grigera, R. Perry, A. Schofield, M. Chiao, S. Julian, G. Lonzarich, S. Ikeda, Y. Maeno, A. Millis, and A. Mackenzie, Science 294, 329 (2001).
  * Maeno et al. (1994) Y. Maeno, H. Hashimoto, K. Yoshida, S. Nishizaki, T. Fujita, J. G. Bednorz, and F. Lichtenberg, Nature 372, 532 (1994).
  * Allan et al. (2013) M. P. Allan, A. Tamai, E. Rozbicki, M. H. Fischer, J. Voss, P. D. C. King, W. Meevasana, S. Thirupathaiah, E. Rienks, J. Fink, et al., New Journal of Physics 15, 063029 (2013), URL https://dx.doi.org/10.1088/1367-2630/15/6/063029.
  * Shen et al. (2007) K. M. Shen, N. Kikugawa, C. Bergemann, L. Balicas, F. Baumberger, W. Meevasana, N. J. C. Ingle, Y. Maeno, Z.-X. Shen, and A. P. Mackenzie, Phys. Rev. Lett. 99, 187001 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.99.187001.
  * Marques et al. (2021) C. A. Marques, L. C. Rhodes, R. Fittipaldi, V. Granata, C.-M. Yim, R. Buzio, A. Gerbi, A. Vecchione, A. W. Rost, and P. Wahl, Advanced Materials 33, 2100593 (2021).
  * Marques et al. (2022) C. A. Marques, L. C. Rhodes, I. Benedičič, M. Naritsuka, A. B. Naden, Z. Li, A. C. Komarek, A. P. Mackenzie, and P. Wahl, Sci. Adv. 8, eabo7757 (2022), URL https://www.science.org/doi/10.1126/sciadv.abo7757.
  * Marques et al. (2023) C. A. Marques, P. A. E. Murgatroyd, R. Fittipaldi, W. Osmolska, B. Edwards, I. Benedičič, G.-R. Siemann, L. C. Rhodes, S. Buchberger, M. Naritsuka, et al., _Spin-orbit coupling induced Van Hove singularity in proximity to a Lifshitz transition in Sr4Ru3O_10_ (2023), arXiv:2303.05587 [cond-mat], URL http://arxiv.org/abs/2303.05587.
  * Damascelli et al. (2000) A. Damascelli, D. H. Lu, K. M. Shen, N. P. Armitage, F. Ronning, D. L. Feng, C. Kim, Z.-X. Shen, T. Kimura, Y. Tokura, et al., Phys. Rev. Lett. 85, 5194 (2000), URL https://link.aps.org/doi/10.1103/PhysRevLett.85.5194.
  * Tamai et al. (2019) A. Tamai, M. Zingl, E. Rozbicki, E. Cappelli, S. Riccò, A. de la Torre, S. McKeown Walker, F. Y. Bruno, P. D. C. King, W. Meevasana, et al., Phys. Rev. X 9, 021048 (2019), URL https://link.aps.org/doi/10.1103/PhysRevX.9.021048.
  * Wang et al. (2017) Z. Wang, D. Walkup, P. Derry, T. Scaffidi, M. Rak, S. Vig, A. Kogar, I. Zeljkovic, A. Husain, L. H. Santos, et al., Nature Physics 13, 799 (2017), ISSN 1745-2473, 1745-2481, URL http://www.nature.com/articles/nphys4107.
  * Matzdorf et al. (2002) R. Matzdorf, Ismail, T. Kimura, Y. Tokura, and E. W. Plummer, Phys. Rev. B 65, 085404 (2002), URL https://link.aps.org/doi/10.1103/PhysRevB.65.085404.
  * Morales et al. (2023) E. A. Morales, G.-R. Siemann, A. Zivanovic, P. A. E. Murgatroyd, I. Markovi, B. Edwards, C. A. Hooley, D. A. Sokolov, N. Kikugawa, C. Cacho, et al., Physical Review Letters 130, 096401 (2023), URL https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.130.096401.
  * Pizzi et al. (2020) G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, et al., Journal of Physics: Condensed Matter 32, 165902 (2020), URL https://doi.org/10.1088/1361-648x/ab51ff.
  * com (b) Known mechanisms to remove this pocket would require e.g DFT + U, or modification of the energies of the $e_{g}$ bands, neglected in these models.
  * Yuan and Fu (2020) N. F. Q. Yuan and L. Fu, Phys. Rev. B 101, 125120 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.101.125120.
  * com (c) Specifically, we perform a grid-based search in a parameter space that includes octahedral rotation, $\theta$, SOC $\lambda$, and a quasiparticle renormalisation factor, $Z$ sup . The best fit to the ARPES bands in the vicinity of the $\bar{M}$ point are $\theta=8.03^{\circ}$, $\lambda=0.17\text{ eV}$, and $Z=0.24$. For more details see Supplemental Material.
  * Kreisel et al. (2021a) A. Kreisel, C. A. Marques, L. C. Rhodes, X. Kong, T. Berlijn, R. Fittipaldi, V. Granata, A. Vecchione, P. Wahl, and P. J. Hirschfeld, npj Quantum Mater. 6, 100 (2021a), ISSN 2397-4648, URL https://www.nature.com/articles/s41535-021-00401-x.
  * Lu et al. (1996) D. H. Lu, M. Schmidt, T. R. Cummins, S. Schuppler, F. Lichtenberg, and J. G. Bednorz, Phys. Rev. Lett. 76, 4845 (1996).
  * Yokoya et al. (1996) T. Yokoya, A. Chainani, T. Takahashi, H. Katayama-Yoshida, M. Kasai, and Y. Tokura, Phys. Rev. Lett. 76, 3009 (1996).
  * Kreisel et al. (2021b) A. Kreisel, C. Marques, L. Rhodes, X. Kong, T. Berlijn, R. Fittipaldi, V. Granata, A. Vecchione, P. Wahl, and P. Hirschfeld, npj Quantum Materials 6, 100 (2021b).
  * Barker et al. (2003) B. Barker, S. Dutta, C. Lupien, P. McEuen, N. Kikugawa, Y. Maeno, and J. Davis, Physica B: Condensed Matter 329-333, 1334 (2003), ISSN 09214526, URL http://linkinghub.elsevier.com/retrieve/pii/S0921452602021580.
  * Kambara et al. (2006) H. Kambara, Y. Niimi, K. Takizawa, H. Yaguchi, Y. Maeno, and H. Fukuyama, in _AIP Conference Proceedings_ (AIP, Orlando, Florida (USA), 2006), vol. 850, pp. 539–540, URL http://aip.scitation.org/doi/abs/10.1063/1.2354823.
  * Burganov et al. (2016) B. Burganov, C. Adamo, A. Mulder, M. Uchida, P. D. C. King, J. W. Harter, D. E. Shai, A. S. Gibbs, A. P. Mackenzie, R. Uecker, et al., Phys. Rev. Lett. 116, 197003 (2016), URL http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.116.197003.
  * Liu et al. (2023) B. Liu, M.-Q. Kuang, Y. Luo, Y. Li, C. Hu, J. Liu, Q. Xiao, X. Zheng, L. Huai, S. Peng, et al., Phys. Rev. Lett. 131, 026701 (2023).
  * Regnault et al. (2022) N. Regnault, Y. Xu, M.-R. Li, D.-S. Ma, M. Jovanovic, A. Yazdani, S. S. P. Parkin, C. Felser, L. M. Schoop, N. P. Ong, et al., Nature 603, 824 (2022).
  * Zervou et al. (2023) A. Zervou, D. V. Efremov, and J. J. Betouras, Phys. Rev. Res. 5, L042006 (2023), URL https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.5.L042006.
  * Classen et al. (2020) L. Classen, A. V. Chubukov, C. Honerkamp, and M. M. Scherer, Phys. Rev. B 102, 125141 (2020).
  * Giannozzi et al. (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, et al., Journal of Physics: Condensed Matter 29, 465901 (2017), URL https://doi.org/10.1088/1361-648x/aa8f79.
  * Choubey et al. (2014) P. Choubey, T. Berlijn, A. Kreisel, C. Cao, and P. Hirschfeld, Phys. Rev. B 90, 134520 (2014), URL http://link.aps.org/doi/10.1103/PhysRevB.90.134520.
  * Kreisel et al. (2015) A. Kreisel, P. Choubey, T. Berlijn, W. Ku, B. Andersen, and P. Hirschfeld, Phys. Rev. Lett. 114, 217002 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.217002.
  * Choubey et al. (2017) P. Choubey, A. Kreisel, T. Berlijn, B. M. Andersen, and P. J. Hirschfeld, Phys. Rev. B 96, 174523 (2017), ISSN 2469-9950, 2469-9969, URL https://link.aps.org/doi/10.1103/PhysRevB.96.174523.
  * Benedičič et al. (2022) I. Benedičič, M. Naritsuka, L. C. Rhodes, C. Trainer, Y. Nanao, A. B. Naden, R. Fittipaldi, V. Granata, M. Lettieri, A. Vecchione, et al., Physical Review B 106, L241107 (2022), URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.106.L241107.
  * Rhodes et al. (2023) L. C. Rhodes, W. Osmolska, C. A. Marques, and P. Wahl, Phys. Rev. B 107, 045107 (2023), ISSN 2469-9950, 2469-9969, URL https://link.aps.org/doi/10.1103/PhysRevB.107.045107.
  * Naritsuka et al. (2023) M. Naritsuka, I. Benedičič, L. C. Rhodes, C. A. Marques, C. Trainer, Z. Li, A. C. Komarek, and P. Wahl, Proceedings of the National Academy of Sciences 120, e2308972120 (2023).
  * (47) Details can be found in the Supplemental material.


# On simultaneous similarity of families of commuting operators

Sherwin Kouchekian and Boris Shekhtman  Sherwin Kouchekian
Department of Mathematics & Statistics
University of South Florida
Tampa, FL 33620-5700<EMAIL_ADDRESS>Boris Shekhtman
Department of Mathematics
University of South Florida
Tampa, FL 33620-5700<EMAIL_ADDRESS>
###### Abstract.

Characterization of simultaneously similarity for commuting $m$-tuples of
operators is an open problem even in finite-dimensional spaces; known as “A
wild problem in linear algebra”. In this paper we offer a criteria for
simultaneously similarity of $m$-tuples of $k$-cyclic commuting operators on
an arbitrary Banach space. Moreover, we obtain an additional equivalence
condition in the case of finite dimensional Banach spaces, which extends the
result found in [7] for pairs of cyclic commuting matrices. We also present
two applications of our results, one in the case of general multiplication
operators on Banach spaces of analytic function, and one for $m$-tuples of
commuting square matrices.

2000 Mathematics Subject Classification. Primary: 15A21, 47A16, 47B99;
Secondary: 47L22, 15A30.

##  1\. Introduction

Characterization of simultaneously similarity for pairs of commuting square
matrices is a central problem in classifying algebras with wild type
representations, see [1], [2], and [5]. Gelfand and Ponomorev [4] showed that
characterization of simultaneously similarity for pairs of commuting square
matrices would provide a characterization of simultaneously similarity for
arbitrary pairs of square matrices. Since then, the problem of
characterization of simultaneously similarity for pairs of commuting square
matrices has become known as “A wild problem in linear algebra”. In [7], a
necessary and sufficient condition for cyclic pairs of commuting square
matrices was given in terms of vanishing ideals of polynomials. In this paper
we offer an extension of this result to $k$-cyclic commuting $m$-tuples of
operators on Banach spaces. We being by introducing some notations.

In what follows, $(\mathscr{X},\|\cdot\|)$ will always denote a complex Banach
space, and $\mathscr{B}(\mathscr{X})$ the space of all bounded linear
operators from $\mathscr{X}$ into $\mathscr{X}.$ An $m$-tuple
$(A_{1},\ldots,A_{m})$ is a vector in
$\left(\mathscr{B}(\mathscr{X})\right)^{m}$ and is called a commuting
$m$-tuple on $\mathscr{X}$ if $A_{i}A_{j}=A_{j}A_{i}$ for all $i,j=1,\dots,m.$
Two commuting $m$-tuples $(A_{1},\ldots,A_{m})$ and $(B_{1},\ldots,B_{m})$ on
$\mathscr{X},$ denoted by $\mathbf{A}$ and $\mathbf{B}$ respectively, are
called simultaneously similar if there exists an invertible operator $S$ in
$\mathscr{B}(\mathscr{X})$ such that $B_{j}=SA_{j}S^{-1}$ for all
$j=1,\ldots,m$. Moreover, we let $\mathbb{C}[{\bf
z}]:=\mathbb{C}[z_{1},...,z_{m}]$ stand for the algebra of polynomials in $m$
variables over the complex field $\mathbb{C}.$ Finally, if
$p({\mathbf{z}})=p(z_{1},...,z_{m})$ is a polynomial in $\mathbb{C}[{\bf z}]$
and if $\mathbf{A}=(A_{1},\ldots,A_{m})$ belongs to
$\left(\mathscr{B}(\mathscr{X})\right)^{m},$ then by $p(\mathbf{A})$ we mean
the operator $p(A_{1},...,A_{m})$ in $\mathscr{B}(\mathscr{X}).$

This paper is organized as follows. Section 2 contains our main result,
Theorem 2.3, which provides an equivalent condition for simultaneously
similarity between pairs of $k$-cyclic $m$-tuples of commuting operators on a
Banach space. In Section 3, we show that one can offer an additional useful
necessary and sufficient condition to Theorem 2.3 for finite dimensional case,
which is not true in general Banach spaces. The provided simple examples in
this section show the stark difference between finite and infinite dimensional
Banach spaces. Section 4 concludes our paper with two general applications of
Theorem 2.3 and Theorem 3.2.

##  2\. Main Result

To start with, we need the following definition.

######  Definition 2.1.

For a vector $\mathbf{u}=(u_{1},\ldots,u_{k})$ in $\mathscr{X}^{k}$ and an
$m$-tuple $\mathbf{A}=(A_{1},\ldots,A_{m})$ in
$\left(\mathscr{B}(\mathscr{X})\right)^{m}$ define

(2.2)
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})=\left\\{\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}:p_{j}\in\mathbb{C}[{\bf
z}]\right\\}.$

If $\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ is dense in $\mathscr{X},$ we call
$\mathbf{u}$ a cyclic $k$-tuple for $\mathbf{A}$, or $\mathbf{A}$ is
$k$-cyclic with respect to the $k$-tuple $\mathbf{u}.$ Moreover if
$\mathbf{u}$ is understood, $\mathbf{A}$ is simply called a $k$-cyclic
$m$-tuple. It is clear that $\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ is a
(linear) subspace of $\mathscr{X},$ and our definition extends the standard
definition of cyclicity to an $m$-tuple of operators.

We are now in the position to state our main result.

######  Theorem 2.3.

Let $\mathbf{A}=(A_{1},\ldots,A_{m})$ be a commuting $m$-tuple on
$\mathscr{X},$ and suppose $\mathbf{u}=(u_{1},...,u_{k})$ is a cyclic
$k$-tuple for $\mathbf{A}$. A commuting $m$-tuple
$\mathbf{B}=(B_{1},\ldots,B_{m})$ on $\mathscr{X}$ is simultaneously similar
to $\mathbf{A}$ if and only if there exists a cyclic $k$-tuple
$\mathbf{v}=(v_{1},...,v_{k})$ for $\mathbf{B}$ and a positive constant $c>0$
such that

(2.4)
$c^{-1}\left\|\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\right\|\leq\left\|\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}\right\|\leq
c\left\|\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\right\|$

for all polynomials $p_{j}$ in $\mathbb{C}[{\bf z}]$.

###### Proof.

First, suppose that $\mathbf{A}$ and $\mathbf{B}$ are simultaneously similar.
Thus, there exists an invertible operator $S$ in $\mathscr{B}(\mathscr{X})$
such that $B_{j}=SA_{j}S^{-1}$ for all $j=1,\ldots,m$. Define $v_{j}=Su_{j}$
for $j=1,\dots,k$, and let $\mathbf{v}=(v_{1},\dots,v_{k}).$ For any $p$ in
$\mathbb{C}[{\bf z}],$ since $\mathbf{B}$ is a commuting $m$-tuple, we have

$p(\mathbf{A})=p(A_{1},\dots,A_{m})=p(S^{-1}B_{1}S,\dots,S^{-1}B_{m}S)=S^{-1}p(B_{1},\dots,B_{m})S=S^{-1}p(\mathbf{B})S.$

Therefore,

(2.5)
$\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}=\sum_{j=1}^{k}S^{-1}p_{j}(\mathbf{B})Su_{j}=S^{-1}\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\quad\text{
for all }\quad p_{j}\in\mathbb{C}[{\bf z}].$

Recalling the definition (2.2), it follows from (2.5) that
$\mathcal{L}_{\mathbf{B}}(\mathbf{v})=S\mathcal{L}_{\mathbf{A}}(\mathbf{u}).$
Since $S$ is onto and $\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ is dense in
$\mathscr{X}$ by the assumption, it follows that
$\mathcal{L}_{\mathbf{B}}(\mathbf{v})$ is dense in $\mathscr{X}$ as well.
Thus, $\mathbf{v}=(v_{1},...,v_{k})$ is a cyclic $k$-tuple for $B$. Moreover,
(2.5) also implies that

(2.6)
$\left\|\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}\right\|\leq\left\|S^{-1}\right\|\left\|\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\right\|\quad\text{
for all }\quad p_{j}\in\mathbb{C}[{\bf z}].$

Applying $S$ from the left to (2.5) and taking the norm, we obtain

(2.7)
$\frac{1}{\|S\|}\left\|\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\right\|\leq\left\|\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}\right\|\quad\text{
for all }\quad p_{j}\in\mathbb{C}[{\bf z}].$

Now inequality (2.4) follows from (2.6) and (2.7) with
$c=\max\\{\|S\|,\|S^{-1}\|\\}.$ This establishes the proof of the necessity
part.

Conversely, suppose that $\mathbf{B}$ is a commuting $m$-tuple, and let
$\mathbf{v}=(v_{1},...,v_{k})$ be a cyclic $k$-tuple for $\mathbf{B}$ which
satisfies (2.4). Define $\mathfrak{L}$ from
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})\subseteq\mathscr{X}$ onto
$\mathcal{L}_{\mathbf{B}}(\mathbf{v})\subseteq\mathscr{X}$ as

(2.8)
$\mathfrak{L}\Bigl{(}\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}\Bigr{)}=\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\quad\text{
for all }\quad p_{j}\in\mathbb{C}[{\bf z}].$

First we show that $\mathfrak{L}$ is well-defined. To see this, suppose that
$x\in\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ has two representations
$x=\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}$ and
$x=\sum_{j=1}^{k}q_{j}(\mathbf{A})u_{j}$ for some
$p_{j},q_{j}\in\mathbb{C}[{\bf z}]$ and $1\leq j\leq k.$ In other words,
$\sum_{j=1}^{k}(p_{j}(\mathbf{A})-q_{j}(\mathbf{A}))u_{j}=0.$ Now the first
inequality of (2.4) implies that
$\sum_{j=1}^{k}(p_{j}(\mathbf{B})-q_{j}(\mathbf{B}))v_{j}=0;$ or equivalently,
$\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}=\sum_{j=1}^{k}q_{j}(\mathbf{B})v_{j}.$
In view of the definition for $\mathfrak{L},$ the last equality is equivalent
to
$\mathfrak{L}\,x=\mathfrak{L}(\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j})=\mathfrak{L}(\sum_{j=1}^{k}q_{j}(\mathbf{A})u_{j}).$
Thus, $\mathfrak{L}$ is well-defined.

Furthermore, $\mathfrak{L}$ is clearly linear. It follows trivially from the
first inequality of (2.4) that $\mathfrak{L}$ is also bounded on
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ with $\|\mathfrak{L}\|\leq c.$ Since
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ is dense in $\mathscr{X},$ the bounded
linear transformation (BLT) theorem implies that $\mathfrak{L}$ has a unique
norm-preserving extension to $\mathscr{X};$ that is, there exists a unique $S$
in $\mathscr{B}(\mathscr{X})$ such that $Sx=\mathfrak{L}x$ for all $x$ in
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ and $\|S|=\|\mathfrak{L}\|.$ We prove
that $S$ is the desired operator which provides the simultaneous similarity
between $\mathbf{A}$ and $\mathbf{B}.$ To show $S$ is invertible, in view of
the open mapping theorem, it suffices to prove that $S$ is a bijection.

To do this, first let $x\in\mathcal{L}_{\mathbf{A}}(\mathbf{u}).$ Then $x$ has
the form $x=\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}$ for some $p_{j}$ in
$\mathbb{C}[{\bf z}].$ Now, suppose $x\in\mathcal{L}_{\mathbf{A}}(\mathbf{u})$
such that $Sx=0.$ Using the fact that $Sx=\mathfrak{L}x,$ the definition of
$\mathfrak{L}$ together with the second inequality of (2.4) imply that $x=0.$
Since $\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ is dense in $\mathscr{X},$ it
follows from the continuity of $S$ that $\text{Ker\,}(S)=\\{0\\},$ where
$\text{Ker\,}(S)$ stands for the kernel of $S.$ This proves that $S$ is one-
to-one.

Next, let $y\in\mathscr{X}.$ By density of
$\mathcal{L}_{\mathbf{B}}(\mathbf{v}),$ there exists a sequence
$\\{y_{n}\\}\subseteq\mathcal{L}_{\mathbf{B}}(\mathbf{v})$ such that $y_{n}\to
y$ in $\mathscr{X};$ that is, $\lim_{n\to\infty}\|y_{n}-y\|=0,$ where
$y_{n}=\sum_{j=1}^{k}p_{j}^{(n)}(\mathbf{B})v_{j}$ for some sequence of
polynomials $\\{p_{1}^{(n)},\dots,p_{k}^{(n)}\\}.$ Letting
$x_{n}=\sum_{j=1}^{k}p_{j}^{(n)}(\mathbf{A})u_{j},$ we have from the
definition of $\mathfrak{L}$ that $\mathfrak{L}x_{n}=y_{n}$ for all $n\geq 1.$
Since $\\{y_{n}\\}$ converges in $\mathscr{X},$ and hence a Cauchy sequence,
the second inequality of (2.4) implies that $\\{x_{n}\\}$ is also a Cauchy
sequence in $\mathscr{X}.$ Therefore, there exists an $x\in\mathscr{X}$ such
that $x_{n}\to x$ in $\mathscr{X}.$ Now, it follows from the continuity of $S$
that

(2.9)
$Sx=\lim_{n\to\infty}Sx_{n}=\lim_{n\to\infty}\mathfrak{L}\,x_{n}=\lim_{n\to\infty}y_{n}=y.$

Thus, $S$ is also surjective.

It remains to show $\mathbf{A}$ and $\mathbf{B}$ are simultaneously similar.
We start by fixing $x\in\mathscr{X}$ and use the density of
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ to get a sequence
$\\{x_{n}\\}\subseteq\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ such that $x_{n}\to
x$ in $\mathscr{X},$ where $x_{n}=\sum_{j=1}^{k}p_{j}^{(n)}(\mathbf{A})u_{j}$
for some sequence of polynomials $\\{p_{1}^{(n)},\dots,p_{k}^{(n)}\\}$ in
$\mathbb{C}[{\bf z}].$ For each fixed $n\geq 1$ and $1\leq j\leq k,$ let
$q_{\ell,j}^{(n)}=z_{\ell}p_{j}^{(n)},$ where $1\leq\ell\leq m.$ Clearly for
each fixed $\ell,$

$A_{\ell}x_{n}=A_{\ell}\sum_{j=1}^{k}p_{j}^{(n)}(\mathbf{A})u_{j}=\sum_{j=1}^{k}A_{\ell}p_{j}^{(n)}(\mathbf{A})u_{j}=\sum_{j=1}^{k}q_{\ell,j}^{(n)}(\mathbf{A})u_{j}.$

Thus, it follows from the definitions of $\mathfrak{L}$ and $q_{\ell,j}^{(n)}$
that

(2.10)
$\mathfrak{L}A_{\ell}x_{n}=\mathfrak{L}\sum_{j=1}^{k}q_{\ell,j}^{(n)}(\mathbf{A})u_{j}=\sum_{j=1}^{k}q_{\ell,j}^{(n)}(\mathbf{B})v_{j}=B_{\ell}\sum_{j=1}^{k}p_{j}^{(n)}(\mathbf{B})v_{j}=B_{\ell}\mathfrak{L}x_{n}.$

Finally, the equation (2.10), together with a similar argument as in (2.9),
implies

(2.11)
$SA_{\ell}x=\lim_{n\to\infty}SA_{\ell}x_{n}=\lim_{n\to\infty}\mathfrak{L}A_{\ell}x_{n}=\lim_{n\to\infty}B_{\ell}\mathfrak{L}x_{n}=\lim_{n\to\infty}B_{\ell}Sx_{n}=B_{\ell}Sx.$

Since $x$ is arbitrary, it follows from (2.11) that $SA_{\ell}=B_{\ell}S$ for
all $\ell=1,\dots,m.$ This finishes the proof of the theorem. ∎

We make a couple of observations here. First of all, we note that the positive
constants $c^{-1}$ and $c$ in (2.4) could have been equivalently substituted
by two arbitrary positive constants $c_{1}$ and $c_{2},$ respectively. The
given format, however, makes (2.4) invariant under the substitution of
$(\mathbf{A},\mathbf{u})$ with $(\mathbf{B},\mathbf{v}),$ or vice versa. This
also agrees with the statement of Theorem 2.3, which is symmetric with respect
to $\mathbf{A}$ and $\mathbf{B}.$ Next, the proof of necessity part of Theorem
2.3 clearly shows that if $\mathbf{A}$ and $\mathbf{B}$ are two simultaneously
similar $m$-tuples, then $A$ is $k$-cyclic if and only if $\mathbf{B}$ is
$k$-cyclic. For instance, if $\mathbf{u}=(u_{1},\dots,u_{k})$ is a cyclic
$k$-tuple for $\mathbf{A},$ then $\mathbf{v}=(Su_{1},\dots,Su_{k})$ is a
cyclic $k$-tuple for $\mathbf{B},$ where $\mathbf{u}$ and $\mathbf{v}$ both
satisfy (2.4). In general, however, the converse is not true; that is, there
are simultaneously similar $k$-cyclic commuting $m$-tuples $\mathbf{A}$ and
$\mathbf{B},$ where the corresponding cyclic $k$-tuples $\mathbf{u}$ and
$\mathbf{v}$ do not satisfy (2.4). At the end of the next section, we provide
a simple example by utilizing Theorem 3.2.

##  3\. The Finite Dimensional Case

In this section, we assume $\mathscr{X}$ has a finite dimension. Therefore, it
may be assumed without loss of generality that $\mathscr{X}=\mathbb{C}^{N}$
for some natural number $N.$ It is also clear that an $m$-tuple
$\mathbf{A}=(A_{1},\ldots,A_{m})$ in $(\mathscr{B}(\mathscr{X}))^{m}$ is now
simply an $m$-tuple of $N\times N$ matrices on $\mathbb{C}^{N}.$ Our main goal
here is Theorem 3.2, where we show that one can add another very useful
equivalent condition to Theorem 2.3. This equivalent condition, however, does
not hold in the infinite dimensional case. A simple counterexample will be
provided at the end of this section. We start with a general lemma.

######  Lemma 3.1.

Suppose $\mathbf{A}=(A_{1},\ldots,A_{m})$ and
$\mathbf{B}=(B_{1},\ldots,B_{m})$ are two $m$-tuples of commuting $N\times N$
matrices on $\mathbb{C}^{N}.$ If $\mathbf{u}=(u_{1},...,u_{k})$ and
$\mathbf{v}=(v_{1},...,v_{k})$ are two arbitrary vectors in
$\left(\mathbb{C}^{N}\right)^{k},$ then the following conditions are
equivalent.

  1. (I)

There exists a positive constant $c>0$ such that

$c^{-1}\left\|\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\right\|\leq\left\|\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}\right\|\leq
c\left\|\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\right\|,$

for all polynomials $p_{j}$ in $\mathbb{C}[{\bf z}]$.

  2. (II)

$\left\\{p_{1},\dots,p_{k}\in\mathbb{C}[{\bf
z}]:\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}=0\right\\}=\left\\{p_{1},\dots,p_{k}\in\mathbb{C}[{\bf
z}]:\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}=0\right\\}.$ In other words if
$p_{1},\dots,p_{k}\in\mathbb{C}[{\bf z}],$ then

$\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}=0\quad\text{if and only
if}\quad\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}=0.$

###### Proof.

Since all norms on a finite a dimensional space are equivalent, we may assume
the norm in (I) is any arbitrary norm on $\mathbb{C}^{N}.$ Now if
$\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}=0$ for some $p_{1},\dots,p_{k}$ in
$\mathbb{C}[{\bf z}],$ then it follows from the first inequality of (I) that
$\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}=0,$ and vice versa. Thus (I) implies
(II).

To prove (II) implies (I), define $\mathfrak{L}$ from
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ onto
$\mathcal{L}_{\mathbf{B}}(\mathbf{v})$ as in (2.8). Using an exact similar
argument as in the proof of Theorem 2.3, replacing inequality (2.4) with the
condition (II), one concludes that $\mathfrak{L}$ is well defined. Since
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ and
$\mathcal{L}_{\mathbf{B}}(\mathbf{v})$ are now subspaces of $\mathbb{C}^{N},$
they are both closed, and thus complete. This implies that $\mathfrak{L}$ is a
bounded operator from $\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ onto
$\mathcal{L}_{\mathbf{B}}(\mathbf{v})$. Moreover if
$x=\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}$ belongs to
$\mathcal{L}_{\mathbf{A}}(\mathbf{u})$ such that $\mathfrak{L}\,x=0;$ or
equivalently $\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}=0$, then it follows from
(II) that $x=0;$ that is $\text{Ker\,}(S)=\\{0\\}.$ This shows that
$\mathfrak{L}$ is also one to one, and thus it is invertible. Now using a
similar argument given in the first part of the proof of Theorem 2.3, we
obtain the inequalities in (II) with
$c=\max\\{\|\mathfrak{L}\|,\|\mathfrak{L}^{-1}\|\\}.$ ∎

An immediate consequence of Lemma 3.1 when combined with Theorem 2.3 is the
following result for the finite dimensional case.

######  Theorem 3.2.

Let $\mathbf{A}=(A_{1},\ldots,A_{m})$ and $\mathbf{B}=(B_{1},\ldots,B_{m})$ be
two $m$-tuples of commuting $N\times N$ matrices on $\mathbb{C}^{N}.$ If
$\mathbf{A}$ is $k$-cyclic with respect to the $k$-tuple
$\mathbf{u}=(u_{1},...,u_{k}),$ then the following statements are equivalent.

  1. (a)

$\mathbf{B}$ is simultaneously similar to $\mathbf{A}.$

  2. (b)

There exists a cyclic $k$-tuple $\mathbf{v}=(v_{1},...,v_{k})$ for
$\mathbf{B}$ and a positive constant $c>0$ such that

$c^{-1}\left\|\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\right\|\leq\left\|\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}\right\|\leq
c\left\|\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}\right\|,$

for all polynomials $p_{j}$ in $\mathbb{C}[{\bf z}]$.

  3. (c)

There exists a cyclic $k$-tuple $\mathbf{v}=(v_{1},...,v_{k})$ for
$\mathbf{B}$ such that

$\Bigl{\\{}p_{1},\dots,p_{k}\in\mathbb{C}[{\bf
z}]:\sum_{j=1}^{k}p_{j}(\mathbf{A})u_{j}=0\Bigr{\\}}=\Bigl{\\{}p_{1},\dots,p_{k}\in\mathbb{C}[{\bf
z}]:\sum_{j=1}^{k}p_{j}(\mathbf{B})v_{j}=0\Bigr{\\}}.$

As mentioned in the beginning of this section, condition (c) of Theorem 3.2 is
not equivalent to simultaneous similarity in infinite dimensional case. Here
is a simple example.

######  Example 3.3.

Let $(\mathscr{X},\|\cdot\|)=(H^{2},\|\cdot\|_{2}),$ where $H^{2}$ denotes the
Hardy space of analytic functions on the open unit disk
$\mathbb{D}=\\{z\in\mathbb{C}:|z|<1\\}$ defined as (see [3])

$H^{2}=\Biggl{\\{}f:f(z)=\sum_{n=0}^{\infty}a_{n}z^{n},\ \ z\in\mathbb{D},\ \
a_{n}\in\mathbb{C},\ \ \text{ and }\ \
\|f\|_{2}^{2}=\sum_{n=0}^{\infty}|a_{n}|^{2}<\infty\Biggr{\\}}.$

Let $A=M_{z}$ be the operator of multiplication by the independent variable
$z$ on $H^{2}$ defined by $M_{z}:f(z)\mapsto zf(z).$ Clearly
$A\in\mathscr{B}(\mathscr{X}).$ Since polynomials are dense in $H^{2},$ $A$ is
cyclic with the cyclic vector, say, $u=1.$ Next, we set $B=2A=M_{2z}.$ It
follows that $B\in\mathscr{B}(\mathscr{X})$ and $B$ is also cyclic with
respect to the same cyclic vector $u=1.$

Now if $p\in\mathbb{C}[z],$ then it follows from the definition that
$p(A)u=p(z).$ Thus, $p(A)u=0$ if and only if $p(z)\equiv 0;$ or equivalently,
$\\{p\in\mathbb{C}[z]:p(A)u=0\\}=\\{0\\}.$ A similar argument also implies
that $\\{p\in\mathbb{C}[z]:p(B)u=0\\}=\\{0\\}.$ Therefore, condition (c) of
Theorem 3.2 holds. However, $A$ and $B=2A$ are not (simultaneously) similar to
each other. For the sake of completeness, we provide a proof here. So suppose
that there exists an invertible $S\in\mathscr{B}(\mathscr{X})$ such that
$2A=SAS^{-1.}$ Itteration of the last equality for $n\geq 1$ gives
$2^{n}A^{n}=SA^{n}S^{-1}.$ Consequently,

$2^{n}\|A^{n}\|=\|2^{n}A^{n}\|=\|SA^{n}S^{-1}\|\leq\|S\|\|A^{n}\|\|S^{-1}\|\quad\text{
for all }n\geq 1.$

Since $\|A^{n}\|>0$ ($n\geq 1$), the above inequality implies that
$\|S\|\|S^{-1}\|\geq 2^{n}$ for all $n\geq 1,$ which is absurd as both $\|S\|$
and $\|S^{-1}\|$ are finite. Therefore, no such $S$ can exist.

We should mention that in the infinite dimensional case, such as in Example
3.3, it may very well happen that the condition in part (c) of Theorem 3.2 is
trivially verified because both sets described in this part turn out to be the
zero ideal $\\{0\\}.$ In the finite dimensional case, however, this phenomena
can never occur. To see this, recall that by the Cayley-Hamilton theorem every
$N\times N$ matrix $A$ satisfies its own characteristic polynomial; that is,
there always exists a polynomial $p,$ and hence infinitely many, such that
$p(A)$ vanishes identically. This is the underlying reason why the ideal
approach in the infinite dimensional case cannot serve as a fruitful strategy.

We conclude this section with an example of two cyclic (simultaneously)
similar operators $A$ and $B$ with a common cyclic $k$-tuple for which the
condition (2.4) of Theorem 2.3 does not hold, see the remarks at the end of
Section 2.

######  Example 3.4.

Let $\mathscr{X}=\mathbb{C}^{2},$ considered as the usual 2-dimensional vector
space over $\mathbb{C}$ with
$\\{e_{1},e_{2}\\}=\bigl{\\{}(1,0)^{T},(0,1)^{T}\bigr{\\}}$ as its standard
basis. Note also that $(e_{1},e_{2})$ is trivially a cyclic 2-tuple for any
$2\times 2$-matrix. Consider $A=\begin{pmatrix}0&1\\\ 0&0\end{pmatrix}$ and
$B=\begin{pmatrix}-1&1\\\ -1&1\end{pmatrix}.$ A routine check shows that
$B=SAS^{-1},$ where $S=\begin{pmatrix}-1&2\\\ -1&1\end{pmatrix}.$ Thus, $A$
and $B$ are (simultaneously) similar with a common cyclic $2$-tuple
$\mathbf{u}=\mathbf{v}=(e_{1},e_{2}).$ Now if we let $p_{1}(z)=1$ and
$p_{2}(z)=-z,$ then clearly $\sum_{j=1}^{2}p_{j}(A)u_{j}=e_{1}-Ae_{2}=0,$
whereas $\sum_{j=1}^{2}p_{j}(B)v_{j}=e_{1}-Be_{2}\neq 0.$ This shows that
condition (c) of Theorem 3.2 does not hold. But in light of Lemma 3.1,
conditions (b) and (c) of Theorem 3.2 are always equivalent in the finite
dimensional case. Thus, (2.4) does not hold either. This completes our
argument. We make a final remark that the provided trivial example was only
possible due to the added equivalence condition in the finite dimensional
case, which can be easily checked.

##  4\. Examples and Illustrations

We conclude this paper with some applications of the obtained result. As a
first application, let $\Omega$ be a non-empty open subset of
$\mathbb{C}^{m},$ and denote by
$(\mathscr{X}(\Omega),\|\cdot\|_{\mathscr{X}})$ the Banach space of
holomorphic functions on $\Omega$ for which $\mathbb{C}[{\bf z}]$ is a dense
subset. For example, one could consider Dirichlet-type spaces on the
$m$-dimensional unit polydisc

$\mathbb{D}^{m}=\\{(z_{1},\dots,z_{m})\in\mathbb{C}^{m}:|z_{j}|<1,\ \
j=1\dots,m\\},$

which includes the well-known Hardy and Bergman spaces over the polydiscs and
the polydisc algebra $A(\mathbb{D}^{m})$, see [6]. Furthermore let $M_{z_{j}}$
denote the operator of multiplication by the $j^{th}$-coordinate $z_{j}$ on
$\mathscr{X}(\Omega)$ defined as $M_{z_{j}}:f({\mathbf{z}})\mapsto
z_{j}f({\mathbf{z}}).$ It follows that
$\mathbf{M}=(M_{z_{j}},\dots,M_{z_{m}})$ is a cyclic $m$-tuple of commuting
operators on $\mathscr{X}(\Omega)$ with the cyclic vector $1.$ Since
$p(\mathbf{M})1=p({\mathbf{z}})$ for all $p\in\mathbb{C}[{\bf z}],$ the
following result is an immediate consequence of Theorem 2.3.

######  Corollary 4.1.

A commuting $m$-tuple $A=(A_{1},\dots,A_{m})$ on $\mathscr{X}(\Omega)$ is
simultaneously similar to $\mathbf{M}=(M_{z_{1}},\dots,M_{z_{m}})$ if and only
if $\mathbf{A}$ is cyclic and there exists a cyclic vector $u$ in
$\mathscr{X}(\Omega)$ for $\mathbf{A}$ such that

$c^{-1}\left\|p\right\|_{\mathscr{X}}\leq\left\|p(\mathbf{A})u\right\|_{\mathscr{X}}\leq
c\left\|p\right\|_{\mathscr{X}},$

for some constant $c>0$ and all polynomials $p$ in $\mathbb{C}[{\bf z}].$

Our next and final result is concerned with the finite dimensional case
generalizing the idea already presented in Example 3.4. So let
$\mathscr{X}=\mathbb{C}^{N}$ with the standard basis

$e_{1}=(1,0,\dots,0)^{T},\quad e_{2}=(0,1,\dots,0)^{T},\quad
e_{N}=(0,0,\dots,1)^{T}.$

As noted in Example 3.4, $\mathbf{u}=(e_{1},\dots,e_{N})$ is clearly a cyclic
$m$-tuple for any $m$-tuple of commuting $N\times N$ matrices
$\mathbf{A}=(A_{1},\dots,A_{m}).$ Therefore, as a consequence of Theorem 3.2,
we have the following corollary which provides an answer to a ”A wild problem
in Linear Algebra” for the cyclic case.

######  Corollary 4.2.

Two $m$-tuples of commuting $N\times N$ matrices
$\mathbf{A}=(A_{1},\dots,A_{m})$ and $\mathbf{B}=(B_{1},\dots,B_{m})$ are
simultaneously similar if and only if there exists a basis
$\mathbf{v}=(v_{1},...,v_{N})$ for $\mathbb{C}^{N}$ such that

$\Biggl{\\{}p_{1},\dots,p_{N}\in\mathbb{C}[{\bf
z}]:\sum_{j=1}^{N}p_{j}(\mathbf{A})e_{j}=0\Biggr{\\}}=\Biggl{\\{}p_{1},\dots,p_{N}\in\mathbb{C}[{\bf
z}]:\sum_{j=1}^{N}p_{j}(\mathbf{B})v_{j}=0\Biggr{\\}}.$

## References

  * [1] P. Donovan and M.R. Freislich, Some evidence for an extension of the Brauer-Thrall conjecture, Sonderforschungsbereich Theor. Math., 40 (1972), 24–26.
  * [2] Y. A. Drozd, Tame and wild matrix problems, Representations and Quadratic Forms, (1979), 39–74.
  * [3] P. L. Duren, Theory of Hp Spaces, Academic Press, New York, 1970.
  * [4] I.M. Gelfand and V.A. Ponomarev, Remarks on the classification of a pair of commuting linear transformations in a finite-dimensional space, Funktsional. Anal. i Prilozhen., 3 (1969) no. 4, 81–82. [English transl.: Functionul Anal. Appl. 3 (1969), 325–326]
  * [5] Michiel Hazewinkel, Nadiya Gubareni, and V. V. Kirichenko, Algebras, rings and modules. Vol. 2, Springer (2007).
  * [6] W. Rudin, Function Theory in Polydisks, W. A. Benjamin, Inc., New York-Amsterdam, 1969.
  * [7] B. Shekhtman, On Cyclicity and Regularity of Commuting Matrices, Mathematics and Statistics 1 (2013) 157–161.


In 5G/6G wireless systems, reconfigurable intelligent surfaces (RIS) can play a role as a passive anchor to enable and enhance localization in various scenarios.
However, most existing RIS-aided localization works assume that the geometry of the RIS is perfectly known, which is not realistic in practice due to calibration errors.
In this work, we derive the misspecified Cramér-Rao bound (MCRB) for a single-input-single-output RIS-aided localization system with RIS geometry mismatch.

Specifically, unlike most existing works that use numerical methods, we propose a closed-form solution to the pseudo-true parameter determination problem for MCRB analysis.

Simulation results demonstrate the validity of the derived pseudo-true parameters and MCRB,
and show that the RIS geometry mismatch causes performance saturation in the high signal-to-noise ratio regions.

Localization, RIS, 5G/6G, geometry mismatch, calibration error, MCRB

§ INTRODUCTION

ris-aided localization has been extensively studied recently [1, 2, 3]. One of the merits is that ris enables localization in extreme scenarios. For example, a siso system consisting of a bs and a ue can perform communication, but achieving localization is impossible. With the introduction of the RIS channel and the corresponding delay and aod measurements, however, joint ue localization and synchronization can be completed [4]. RIS is also beneficial in other scenarios such as wireless fingerprinting localization [5], signal strength-based localization [6], localization under mobility [7], user tracking [8], terahertz band localization [9], etc.

Unfortunately, almost all the mentioned systems assume that the ris position and orientation are perfectly known. However, the geometry mismatch of the ris is likely to be introduced in reality due to calibration errors. These calibration errors propagate to the estimated channel parameters and cause the model mismatch in the position estimation process.

For the mismatched model, the mcrb is a tool to quantify the impact of the ris calibration error [10].

When using mcrb, the assumed channel model is different from the true model, and a misspecified performance bound can be derived with the model mismatch considered.
The mcrb analysis for radio localization under hardware impairment [11] and channel model mismatch [12] have been reported in previous works. In the mcrb derivation, one of the most essential steps is determining the pseudo-true parameters [10], which is to find a solution that minimizes the kld between the true and mismatched statistical models and is usually accomplished using numerical methods [11, 12].

In this work, we aim to derive the mcrb to evaluate the impact of ris geometry mismatch on RIS-assisted localization.
The main contribution of this work is that we derive a closed-form solution to the pseudo-true parameter determination problem. The corresponding geometrical interpretation is also given.

The simulation code of this paper is available at  <https://github.com/ZPinjun/RISgeoMCRB2023ICASSP>.

§ SYSTEM MODEL

§.§ Geometrical Relations

We consider a downlink siso wireless system with a bs, a ue, and a ris, as shown in Fig. <ref>.
We indicate the position of the bs and the ue by $\mathbf{p}_\text{b}\in\mathbb{R}^3$ and $\mathbf{p}\in\mathbb{R}^3$.
The position and orientation of the ris are denoted as $\mathbf{p}_\text{r}\in\mathbb{R}^3$ and $\mathbf{R}_\text{r}\in \text{SO(3)}$ [13], respectively.
The entries of $\mathbf{p}_\text{b}$, $\mathbf{p}_\text{r}$, and $\mathbf{R}_\text{r}$ are assumed to be known.

Besides, we assume an unknown clock bias $\Delta\in\mathbb{R}$ exists between the ue and the bs [14, 15].

In the ris's lcs, the aoa from the bs consists of an azimuth angle $\theta_{\text{az}}$ and an elevation angle $\theta_{\text{el}}$,
while the aod towards the ue consists of an azimuth angle $\phi_{\text{az}}$  and an elevation angle $\phi_{\text{el}}$.
For compactness, we define $\bm{\theta}=[\theta_{\text{az}},\theta_{\text{el}}]^\mathsf{T} $ and $\bm{\phi}=[\phi_{\text{az}},\phi_{\text{el}}]^\mathsf{T} $.
Since we assume the geometry of the bs and the ris to be known in the ideal model, we focus on the aod towards the ue with an unknown position, which is given by

\begin{align}
\label{eq_phi}
\phi_\text{az} &= \arctan 2\left([\mathbf{R}_\text{r}^\mathsf{T} (\mathbf{p}-\mathbf{p}_\text{r})]_2, [\mathbf{R}_\text{r}^\mathsf{T} (\mathbf{p}-\mathbf{p}_\text{r})]_1\right),
\\
\phi_\text{el} &= \arcsin\left({[\mathbf{R}_\text{r}^\mathsf{T} (\mathbf{p}-\mathbf{p}_\text{r})]_3}/{\|\mathbf{p}-\mathbf{p}_\text{r}\|}\right),
\end{align}

where $[\cdot]_i$ represents the $i$-th element of a vector.
The delays over the bs-ue los and the bs-ris-ue nlos paths are given by

\begin{align}
\tau_\text{b} &= {\|\mathbf{p}_\text{b}-\mathbf{p}\|}/{c} + \Delta, \\
\tau_\text{r} &= {\|\mathbf{p}_\text{b}-\mathbf{p}_\text{r}\|/{c} + \|\mathbf{p}_\text{r}-\mathbf{p}\|}/{c} + \Delta,
\label{eq_tau}
\end{align}

where $c$ is the speed of light.

Illustration of a ris-aided localization system.

§.§ Channel Model

We consider the transmission of $L$ ofdm pilot symbols with $K$ subcarriers.
The frequency of the $k$-th subcarrier is denoted as $f_k= f_c + \frac{(2k-1-K)}{2}\Delta_f, k=1,\dots,K$, where $\Delta_f=B/K$ is the subcarrier spacing, $f_c$ is the carrier frequency, and $B$ is the bandwidth.
The received baseband signal in the ofdm block with index $\ell$ is given by [14]

\begin{equation}\label{eq_channel}
\mathbf{y}_\ell = \underbrace{ \overbrace{g_\text{b}\mathbf{d}(\tau_\text{b})\odot\mathbf{x}}^{\text{\ac{los} channel}} + \overbrace{g_\text{r}\mathbf{b}(\bm{\phi})^\mathsf{T}\bm{\gamma}_\ell(\mathbf{d}(\tau_\text{r})\odot\mathbf{x})}^{\text{\ac{ris} reflected channel}}}_{\bm{\mu}_\ell} + \mathbf{n},
\vspace{-1em}
\end{equation}

where $g_\text{b}$ and $g_\text{r}$ indicate the complex channel gain for the los path and the ris reflected path, $\bm{\gamma}_\ell$ is the ris phase profile, $\mathbf{x}$ is the transmitted signal,
$\mathbf{n}\sim\mathcal{CN}(\mathbf{0},\sigma^2\mathbf{I}_K)$ is a complex awgn vector, and $\bm{\mu}_\ell$ is the noise-free version of the received signal.
The vector $\mathbf{d}(\tau)$ is the delay steering vector defined as

\begin{equation}
\mathbf{d}(\tau) = [1, e^{-j2\pi\Delta_f\tau}, \dots, e^{-j2\pi(K-1)\Delta_f\tau}]^\mathsf{T}.
\end{equation}

Moreover, $\mathbf{b}(\bm{\phi})\triangleq\mathbf{a}(\bm{\theta})\odot\mathbf{a}(\bm{\phi})$,
where $\mathbf{a}(\bm{\theta})$ and $\mathbf{a}(\bm{\phi})$ are respectively the array response vectors for the aoa and aod of the ris.
Under the narrowband far-field model, the array response vector of the ris can be described as
\left[\mathbf{a}\left(\bm{\alpha}\right)\right]_i = e^{j\frac{2\pi f_c}{c}\mathbf{t}(\bm{\alpha})^\mathsf{T}\mathbf{p}_i},
where $\mathbf{p}_i$ is the position of the $i$-th ris element given in its lcs, and $\mathbf{t}(\bm{\alpha})$ is the direction vector defined as

\begin{equation}\notag
\mathbf{t}(\bm{\alpha}) \triangleq
\begin{bmatrix}
\cos(\alpha_{\text{az}})\cos(\alpha_{\text{el}}),\
\sin(\alpha_{\text{az}})\cos(\alpha_{\text{el}}),\
\sin(\alpha_{\text{el}})
\end{bmatrix}^\mathsf{T} .
\end{equation}

§.§ Geometry Mismatch

Taking the bs's position as a reference, we can focus on the calibration errors in $\mathbf{p}_\text{r}$ and $\mathbf{R}_\text{r}$. Thus, the accessible prior geometry information of the ris can be modeled as

\begin{align}
\tilde{\mathbf{p}}_\text{r} = {\mathbf{p}}_\text{r} + \mathbf{u},\qquad
\tilde{\mathbf{R}}_\text{r} = {\mathbf{R}}_\text{p}(\mathbf{v}){\mathbf{R}}_\text{r},
\end{align}

where $\mathbf{u}$
and $\mathbf{v}$ are the calibration errors that cause the mismatch in the model.
Here, ${\mathbf{R}}_\text{p}(\mathbf{v})= \mathbf{R}_z(v_3)\mathbf{R}_y(v_2)\mathbf{R}_x(v_1)$,
where $\mathbf{R}_x(v_1)$ denotes a rotation of $v_1$ degree around the X-axis, and likewise for $\mathbf{R}_y(v_2)$ and $\mathbf{R}_z(v_3)$.

§ LOCALIZATION LOWER BOUND

§.§ Misspecified Cramér-Rao Bound

Suppose a two-stage localization framework is employed [16, 17], that consists of a channel parameter estimation process followed by a UE location estimation from the obtained channel parameters.
For the channel parameter estimation stage, we estimate the following channel parameters from the received signal $\mathbf{y}_\ell, \ell=1,\dots,L$ given in (<ref>):

\begin{equation}\label{eq_eta_ch}
\bm{\eta}_{\text{ch}} \triangleq [ {\phi}_{\text{az}}, {\phi}_{\text{el}}, \tau_\text{b}, {\tau}_\text{r},\mathfrak{R}(g_\text{b}),\mathfrak{I}(g_\text{b}),\mathfrak{R}(g_\text{r}),\mathfrak{I}(g_\text{r}) ]^\mathsf{T},
\end{equation}

where $\mathfrak{R}(\cdot) $ / $\mathfrak{I}(\cdot)$ denote the operations of taking the real / imaginary part.
The fim of $\bm{\eta}_\text{ch}$ can be calculated by the Slepian-Bangs formula <cit.>

\begin{equation}\label{eq_FIMeta}
\mathbf{J}({\bm{\eta}_\text{ch}}) = \frac{2}{\sigma^2}\sum_{\ell = 1}^{L} \mathfrak{R}\left( \left(\frac{\partial \bm{\mu}_\ell}{\partial \bm{\eta}_\text{ch}}\right)^\mathsf{H}  \frac{\partial \bm{\mu}_\ell}{\partial \bm{\eta}_\text{ch}}\right).
\vspace{-0.2em}
\end{equation}

By removing the nuisance parameters related to $g_\text{b}$ and $g_\text{r}$, we further define $\bm{\eta}\triangleq[{\phi}_{\text{az}}, {\phi}_{\text{el}}, \tau_\text{b}, {\tau}_\text{r}]^\mathsf{T} $ used for position and orientation estimation.

We can compute the fim for $\bm{\eta}$ using Schur's complement: we partition $\mathbf{J}(\bm{\eta}_\text{ch})=[\mathbf{X},\mathbf{Y};\mathbf{Y}^\mathsf{T},\mathbf{Z}]$,
where $\mathbf{X}\in\mathbb{R}^{4\times 4}$ so that $\mathbf{J}(\bm{\eta})=\mathbf{X}-\mathbf{Y}\mathbf{Z}^{-1}\mathbf{Y}^\mathsf{T}$.

For the ue location estimation,
we define a parameter vector in the location domain as
\mathbf{r} \triangleq [ \mathbf{p}^\mathsf{T}, \Delta ]^\mathsf{T}\in\mathbb{R}^4.
The objective here is to estimate $\mathbf{r}$ from the estimated $\hat{\bm{\eta}}$, where the ris geometry mismatch occurs.
We express $\bm{\eta}$ as a function of $\mathbf{r}$, i.e., $\bm{\eta} = \mathbf{g}(\mathbf{r}|\mathbf{p}_\text{b},\mathbf{p}_\text{r},\mathbf{R}_\text{r})$, a relationship that is defined by (<ref>)-(<ref>).
Thus, the distribution of the estimated $\hat{\bm{\eta}}$ can be obtained as $\hat{\bm{\eta}}\sim\mathcal{N}\left(\mathbf{g}(\mathbf{r}|\mathbf{p}_\text{b},\mathbf{p}_\text{r},\mathbf{R}_\text{r}),\bm{\Sigma}\right)$, where $\bm{\Sigma} = \mathbf{J}^{-1}(\bm{\eta})$ if an efficient channel estimator is applied in the first stage.
To estimate $\mathbf{r}$, however, we adopt the mismatched model by using $\hat{\bm{\eta}}\sim\mathcal{N}\left(\mathbf{g}({\mathbf{r}}|\mathbf{p}_\text{b},\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{R}}_\text{r}),\bm{\Sigma}\right)$ to obtain an estimate $\hat{\mathbf{r}}$ from $\hat{\bm{\eta}}$.
Then, we have the true likelihood function ($f_\text{T}$) and the mismatched likelihood function ($f_\text{M}$) as

\begin{align}\notag
\ln f_\text{T} &= -\frac{1}{2}(\hat{\bm{\eta}}-\mathbf{g}(\mathbf{r}|\mathbf{p}_\text{r},\mathbf{R}_\text{r}))^\mathsf{T} \bm{\Sigma}^{-1} (\hat{\bm{\eta}}-\mathbf{g}(\mathbf{r}|\mathbf{p}_\text{r},\mathbf{R}_\text{r})),\\
\ln f_\text{M} &= -\frac{1}{2}(\hat{\bm{\eta}}-\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{R}}_\text{r}))^\mathsf{T} \bm{\Sigma}^{-1} (\hat{\bm{\eta}}-\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{R}}_\text{r})).\notag
\end{align}

Note that we omit the constant term and the parameter $\mathbf{p}_\text{b}$ as they do not affect the estimation.

Now, the lower bound matrix of the estimation mse  based on $f_\text{M}$ can be obtained as [12]

\begin{equation}\label{eq_MCRB_matrix}
\text{LBM}(\hat{\mathbf{r}},\bar{\mathbf{r}}) = \underbrace{\mathbf{A}_{\mathbf{r}_0}^{-1} \mathbf{B}_{\mathbf{r}_0} \mathbf{A}_{\mathbf{r}_0}^{-1}}_{\text{MCRB}(\mathbf{r}_0)} + \underbrace{(\bar{\mathbf{r}}-\mathbf{r}_0)(\bar{\mathbf{r}}-\mathbf{r}_0)^\mathsf{T}}_{\text{Bias}(\mathbf{r}_0)},
\end{equation}

where $\bar{\mathbf{r}}$ is the true parameter vector, $\mathbf{r}_0$ is the pseudo-true parameter vector that minimizes the kld between $f_\text{T}$ and $f_\text{M}$. $\mathbf{A}_{\mathbf{r}_0}$ and $\mathbf{B}_{\mathbf{r}_0}$ are two generalizations of the FIMs as [10]

\begin{equation}\label{eq_KLD}
\mathbf{r}_0 = \argmin_\mathbf{r} D(f_\text{T}(\hat{\bm{\eta}}|\bar{\mathbf{r}})\|f_\text{M}(\hat{\bm{\eta}}|\mathbf{r})),
\end{equation}

\begin{align}
\left[\mathbf{A}_{\mathbf{r}_0}\right]_{i,j} &= \mathbb{E}_{f_\text{T}}\left\{\frac{\partial^2}{\partial r_i\partial r_j}\ln f_\text{M}(\hat{\bm{\eta}}|\mathbf{r})\Big|_{\mathbf{r}=\mathbf{r}_0} \right\},\\
&= \left(\frac{\partial^2\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{o}}_\text{r})}{\partial r_i\partial r_j}\right)^\mathsf{T}
\bm{\Sigma}^{-1}(\bm{\eta}-\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{o}}_\text{r}))\Big|_{\mathbf{r}=\mathbf{r}_0} \notag\\
&\qquad\quad -\left(\frac{\partial\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{o}}_\text{r})}{\partial r_i}\right)^\mathsf{T}  \bm{\Sigma}^{-1}\frac{\partial\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{o}}_\text{r})}{\partial r_j}\bigg|_{\mathbf{r}=\mathbf{r}_0},\notag
\end{align}

\begin{align}
\mathbf{B}_{\mathbf{r}_0} &= \mathbb{E}_{f_\text{T}}\left\{
\frac{\partial\ln f_\text{M}(\hat{\bm{\eta}}|\mathbf{r})}{\partial \mathbf{r}}\Big|_{\mathbf{r}=\mathbf{r}_0}
\cdot\left(\frac{\partial\ln f_\text{M}(\hat{\bm{\eta}}|\mathbf{r})}{\partial \mathbf{r}}\Big|_{\mathbf{r}=\mathbf{r}_0}\right)^\mathsf{T}
\right\},\notag\\
&=\left(\frac{\partial\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{o}}_\text{r})}{\partial \mathbf{r}}\right)^\mathsf{T} \bm{\Sigma}^{-1}
\tilde{\bm{\Sigma}}(\mathbf{r})
\bm{\Sigma}^{-1}  \frac{\partial\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{o}}_\text{r})}{\partial \mathbf{r}}  \bigg|_{\mathbf{r}=\mathbf{r}_0}\notag,
\end{align}

where $D(\cdot\|\cdot)$ denotes the kld and $\tilde{\bm{\Sigma}}(\mathbf{r})=\bm{\Sigma}+(\bm{\eta}-\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{o}}_\text{r}))  (\bm{\eta}-\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{o}}_\text{r}))^\mathsf{T}$.
Therefore, the remaining problem is solving (<ref>) to obtain the pseudo-true parameter $\mathbf{r}_0$.
As mentioned, most previous works on localization mismatch analysis solve (<ref>) using numerical methods, e.g., gradient descent, which is time-consuming and offers no guarantees to obtain the global minimum.
The next subsection presents a closed-form solution to the pseudo-true parameter estimation problem with the global minimum guarantee.

§.§ The Closed-form Pseudo-true Parameter Vector

In this subsection, we derive a closed-form solution for (<ref>).
According to the definition of kld, we have

\begin{align}
&D(f_\text{T}(\hat{\bm{\eta}}|\bar{\mathbf{r}})\|f_\text{M}(\hat{\bm{\eta}}|\mathbf{r})) = \mathbb{E}_{f_\text{T}}\left\{\ln f_\text{T}(\hat{\bm{\eta}}|\bar{\mathbf{r}}) - \ln f_\text{M}(\hat{\bm{\eta}}|\mathbf{r})\right\},\notag\\
&= -\frac{1}{2}\mathbb{E}_{f_\text{T}}\left\{ (\hat{\bm{\eta}}-\mathbf{g}(\bar{\mathbf{r}}|\mathbf{p}_\text{r},\mathbf{R}_\text{r}))^\mathsf{T} \bm{\Sigma}^{-1} (\hat{\bm{\eta}}-\mathbf{g}(\bar{\mathbf{r}}|\mathbf{p}_\text{r},\mathbf{R}_\text{r}))\right\}\notag \\
&\qquad + \frac{1}{2}\mathbb{E}_{f_\text{T}}\left\{(\hat{\bm{\eta}}-\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{R}}_\text{r}))^\mathsf{T} \bm{\Sigma}^{-1} (\hat{\bm{\eta}}-\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{R}}_\text{r}))\right\},\notag \\
&\stackrel{(a)}{=}\frac{1}{2}\mathbf{h}(\mathbf{r})^\mathsf{T} \bm{\Sigma}^{-1} \mathbf{h}(\mathbf{r}),\label{eq_quad}
\end{align}

where $\mathbf{h}(\mathbf{r})\triangleq\mathbf{g}(\bar{\mathbf{r}}|\mathbf{p}_\text{r},\mathbf{R}_\text{r})-\mathbf{g}(\mathbf{r}|\tilde{\mathbf{p}}_\text{r},\tilde{\mathbf{R}}_\text{r})$
and step $(a)$ can be obtained by using <cit.>.
Since $\bm{\Sigma}$ is a covariance matrix, then $\bm{\Sigma}^{-1}$ is positive definite.
Hence, the quadratic form (<ref>) always satisfies

\begin{equation}
D(f_\text{T}(\hat{\bm{\eta}}|\bar{\mathbf{r}})\|f_\text{M}(\hat{\bm{\eta}}|\mathbf{r})) = \frac{1}{2}\mathbf{h}(\mathbf{r})^\mathsf{T} \bm{\Sigma}^{-1} \mathbf{h}(\mathbf{r}) \geqslant 0, \label{eq_minimum}
\end{equation}

where the equality holds if and only if $\mathbf{h}(\mathbf{r})=\mathbf{0}$.
The inequality in (<ref>) implies that if there exists a vector $\mathbf{r}_0$ that satisfies $\mathbf{h}(\mathbf{r}_0)=\mathbf{0}$, then it must be the global minimum of (<ref>).
The following steps aim to find such a vector $\mathbf{r}_0$, which is possible in the considered problem since the number of unknowns ($\mathbf{r}\in\mathbb{R}^4$) equals the number of observations ($\mathbf{h}(\mathbf{r})\in\mathbb{R}^4$).

Let $\mathbf{r}_0 = [{\mathbf{p}}^\mathsf{T}_0, {\Delta}_0]^\mathsf{T}$ be the pseudo-true parameters and $[\bar{\mathbf{p}}^\mathsf{T}, \bar{\Delta}]^\mathsf{T} =\bar{\mathbf{r}}$, then $\mathbf{h}(\mathbf{r}_0)=\mathbf{0}$ can be rewritten based on (<ref>)–(<ref>) as

\begin{align}\label{eq_1}
&{\|\mathbf{p}_\text{b}-\bar{\mathbf{p}}\|} + {c}\bar{\Delta} = {\|\mathbf{p}_\text{b}-\mathbf{p}_0\|} + {c}\Delta_0 , \\
&{\|\mathbf{p}_\text{b}-\mathbf{p}_\text{r}\| + \|\mathbf{p}_\text{r}-\bar{\mathbf{p}}\|} + {c}\bar{\Delta} = \|\mathbf{p}_\text{b}-\tilde{\mathbf{p}}_\text{r}\|\notag \\
&\qquad\qquad\qquad\qquad\qquad\qquad + \|\tilde{\mathbf{p}}_\text{r}-\mathbf{p}_0\| + {c}\Delta_0,\\
& {\mathbf{R}_\text{r}^\mathsf{T} (\bar{\mathbf{p}}-\mathbf{p}_\text{r})}/{\|\bar{\mathbf{p}}-\mathbf{p}_\text{r}\|} = {\tilde{\mathbf{R}}_\text{r}^\mathsf{T} (\mathbf{p}_0-\tilde{\mathbf{p}}_\text{r})}/{\|\mathbf{p}_0-\tilde{\mathbf{p}}_\text{r}\|}.\label{eq_3}
\end{align}

Combining (<ref>)–(<ref>), we obtain

\begin{equation}\label{eq_intersec}
(\alpha + \|\mathbf{p}_0-\tilde{\mathbf{p}}_\text{r}\| - \|\mathbf{p}_0-{\mathbf{p}}_\text{b}\|)\frac{\mathbf{p}_0-\tilde{\mathbf{p}}_\text{r}}{\|\mathbf{p}_0-\tilde{\mathbf{p}}_\text{r}\|} = \tilde{\mathbf{R}}_\text{r}\mathbf{R}_\text{r}^\mathsf{T}(\bar{\mathbf{p}}-\mathbf{p}_\text{r}),
\end{equation}

where $\alpha=\|{\mathbf{p}}_\text{b}-\tilde{\mathbf{p}}_\text{r}\| + \|{\mathbf{p}}_\text{b}-\bar{\mathbf{p}}\| - \|{\mathbf{p}}_\text{b}-{\mathbf{p}}_\text{r}\|$.
From (<ref>), it can be inferred that a pseudo-true ue position can be obtained as the intersection of a line $s_l$ and a hyperboloid $s_h$ given by

\begin{align}
s_l:\quad & \mathbf{p} = x\tilde{\mathbf{R}}_\text{r}\mathbf{R}_\text{r}^\mathsf{T}(\bar{\mathbf{p}}-\mathbf{p}_\text{r}) + \tilde{\mathbf{p}}_\text{r},\label{eq_line}\\
s_h:\quad & \|{\mathbf{p}}-\tilde{\mathbf{p}}_\text{r}\| - \|{\mathbf{p}}-{\mathbf{p}}_\text{b}\| = \beta,\label{eq_hyperboloid}
\end{align}

where $x$ is a positive scalar representing the length of the line segment, and $\beta=\|{\bar{\mathbf{p}}}-{\mathbf{p}}_\text{r}\|-\alpha$.
As a result, the pseudo-true ue position can be determined by solving $x$, which is

\begin{equation}
x_0 = \frac{\beta^2-\|\tilde{\mathbf{p}}_\text{r}-\mathbf{p}_\text{b}\|^2}{2[\mathbf{a}^\mathsf{T}(\tilde{\mathbf{p}}_\text{r}-\mathbf{p}_\text{b}) + \beta\|\mathbf{a}\| ]},
\end{equation}

with $\mathbf{a}=\tilde{\mathbf{R}}_\text{r}\mathbf{R}_\text{r}^\mathsf{T}(\bar{\mathbf{p}}-\mathbf{p}_\text{r})$. Then, the pseudo-true ue position ($\mathbf{p}_0$) is obtained by substituting $x_0$ into (<ref>).
An estimate of pseudo-true clock bias $\Delta_0$ can be directly determined from (<ref>) given $\mathbf{p}_0$.

§.§ Mismatched Bound and Estimator

Based on (<ref>), the expected rmse of the ue position estimation (under ris geometry mismatch) can be lower bounded as

\begin{equation}\label{eq_LB}
\sqrt{\mathbb{E}\{\|\mathbf{p}-\hat{\mathbf{p}}\|^2\}} \geq \sqrt{\text{tr}\big(\left[\text{LBM}(\hat{\mathbf{r}},\bar{\mathbf{r}})\right]_{1:3,1:3}\big)} \triangleq \text{LB},
\end{equation}

where $\text{tr}(\cdot)$ returns the trace of a matrix.

To verify the derived mcrb, we formulate a ml-based localization estimator.
To simplify the analysis, we assume that an efficient estimator is applied to obtain the channel parameters $\hat{\bm{\eta}}$. As such, we focus on analyzing the ue position estimation stage where the mismatch impact appears.

Hence we have the following misspecified maximum likelihood (MML) estimation:

\begin{equation}\label{eq_ML}
\hat{\mathbf{r}}_\text{MML} = \argmax_{\mathbf{r}}\  \ln f_\text{M}(\hat{\bm{\eta}}_\text{ch}|\mathbf{r}),
\end{equation}

which can be solved using the gradient descent method.

§ SIMULATION RESULTS

§.§ Simulation Setup

We consider a localization system that consists of a bs at $[5,0,3]^\mathsf{T}$,
a ris at $[0,-5,2.5]^\mathsf{T}$ with orientation $\mathbf{o}_\text{r}=[0^\circ,0^\circ,90^\circ]$,
and a ue with the default location $[-2.5,2.5,0]^\mathsf{T} $.
The ris is of size $64\times 64$ elements with half-wavelength spacing. The other parameters are as follows: average transmission power $P = \unit[10]{dBm}$, carrier frequency $f_c = \unit[28]{GHz}$, bandwidth $W = \unit[400]{MHz}$, number of transmissions $L = 32$, number of subcarriers $K = 3000$, noise PSD $N_0 = \unit[-173.855]{dBm/Hz}$ and noise figure $N_f = \unit[10]{dB}$.

§.§ Results Analysis

The geometric relationship between the pseudo-true ue position $\mathbf{p}_0$ and true ue position $\bar{\mathbf{p}}$.

We first present the geometry relationship between the pseudo-true ue position $\mathbf{p}_0$ and the true ue position $\bar{\mathbf{p}}$.
We set a fixed ris position and orientation mismatch error as [$\mathbf{u}=0.2\times\mathbf{1}$]m and [$\mathbf{v}=3\times\mathbf{1}$]deg, respectively.
Fig. <ref> shows the projection of the corresponding geometry onto the XOY plane.
We can see that the ris geometry mismatch causes a deviation between the true and pseudo-true ue position,
and the pseudo-true ue position $\mathbf{p}_0$ (obtained by solving (<ref>) with a gradient descent method) coincides with the intersection of the line $s_l$ and hyperboloid $s_h$,
which affirms the solution stated in Subsection <ref>.
We can also observe that depending on the position of the ue, either $\beta>0$ or $\beta<0$ can occur. An example of each case is shown in Fig. <ref>.

ML-RMSE, LB, and bias term versus transmitted power (SNR) for RIS position mismatch ([$\mathbf{u}=0.01\times\mathbf{1}$]m) and orientation mismatch ([$\mathbf{v}=0.5\times\mathbf{1}$]deg).

Next, we evaluate the position estimation performance of the estimator (<ref>) and compare it with the theoretical bound  LB (<ref>) to validate the derived pseudo-true parameter vector $\mathbf{r}_0$ and the mcrb.
Fig. <ref> shows the rmse for the ml estimator [20], LB, bias term, and the mismatch-free peb versus different transmitted power for a ris position mismatch ([$\mathbf{u}=0.01\times\mathbf{1}$]m) and orientation mismatch ([$\mathbf{v}=0.5\times\mathbf{1}$]deg) separately.
We can observe that at low transmit power (i.e., low snr given a fixed noise) levels, the LB and the mismatch-free peb coincide, implying that the ris geometry mismatch is not the main source of error.
At higher signal snr, however, LB deviates from the mismatch-free peb and saturates, which reveals the positioning performance is thus more severely affected by the ris geometry mismatch.
The rmse of the ml estimator closely follows the LB, which demonstrates the validity of our derivation.

LB and PEB versus different levels of ris position mismatch ($\sigma_\text{p}=\unit[\{1,2,3,4,5,6\}\times 10^{-2}]{m}$) and orientation mismatch ($\sigma_\text{o}=\unit[\{1,2,3,4,5,6\}\times 10^{-1}]{deg}$).

Finally, in Fig.<ref>, we evaluate the localization performance as a function of the standard deviation of the mismatch errors.
We assume $\mathbf{u}\sim\mathcal{N}(\mathbf{0},\sigma_{\text{p}}^2\mathbf{I}_3)$, $\mathbf{v}\sim\mathcal{N}(\mathbf{0},\sigma_{\text{o}}^2\mathbf{I}_3)$, and generate 100 statistical realizations of the mismatch for each standard deviation value and collect the minimum, maximum, and mean values of the LBs.
It can be seen that, for both position and orientation mismatches, a larger mismatch standard deviation results in an average higher performance bound and produces a larger perturbation in the LB.

§ CONCLUSION

This paper considered a siso ris-aided localization system with ris geometry mismatch.

We proposed a closed-form solution to determine the pseudo-true parameters, which is used for mcrb derivation.
The derived mcrb is validated using the empirical rmse of a ml estimator. A phenomenon is observed whereby the bound saturates as the snr increases.

[1]
Henk Wymeersch, Jiguang He, Benoit Denis, Antonio Clemente, and Markku Juntti,
“Radio localization and mapping with reconfigurable intelligent
surfaces: Challenges, opportunities, and research directions,”
IEEE Veh. Technol. Mag., vol. 15, no. 4, pp. 52–61, 2020.

[2]
Hongliang Zhang, Boya Di, Kaigui Bian, Zhu Han, H. Vincent Poor, and Lingyang
“Toward ubiquitous sensing and localization with reconfigurable
intelligent surfaces,”
Proc. IEEE, vol. 110, no. 9, pp. 1401–1422, 2022.

[3]
Rawan Alghamdi, Reem Alhadrami, Dalia Alhothali, Heba Almorad, Alice Faisal,
Sara Helal, Rahaf Shalabi, Rawan Asfour, Noofa Hammad, Asmaa Shams, Nasir
Saeed, Hayssam Dahrouj, Tareq Y. Al-Naffouri, and Mohamed-Slim Alouini,
“Intelligent surfaces for 6G wireless networks: A survey of
optimization and performance analysis techniques,”
IEEE Access, vol. 8, pp. 202795–202818, 2020.

[4]
Kamran Keykhosravi, Musa Furkan Keskin, Gonzalo Seco-Granados, and Henk
“SISO RIS-enabled joint 3D downlink localization and
in Proc. IEEE Int. Conf. Commun. (ICC), 2021.

[5]
Cam Ly Nguyen, Orestis Georgiou, Gabriele Gradoni, and Marco Di Renzo,
“Wireless fingerprinting localization in smart environments using
reconfigurable intelligent surfaces,”
IEEE Access, vol. 9, pp. 135526–135541, 2021.

[6]
Haobo Zhang, Hongliang Zhang, Boya Di, Kaigui Bian, Zhu Han, and Lingyang Song,
“Metalocalization: Reconfigurable intelligent surface aided
multi-user wireless indoor localization,”
IEEE Trans. Wireless Commun., vol. 20, no. 12, pp. 7743–7757,

[7]
Kamran Keykhosravi, Musa Furkan Keskin, Gonzalo Seco-Granados, Petar Popovski,
and Henk Wymeersch,
“RIS-enabled SISO localization under user mobility and
spatial-wideband effects,”
IEEE J. Sel. Top. Signal Process., vol. 16, no. 5, pp.
1125–1140, 2022.

[8]
Boyu Teng, Xiaojun Yuan, Rui Wang, and Shi Jin,
“Bayesian user localization and tracking for reconfigurable
intelligent surface aided MIMO systems,”
IEEE J. Sel. Top. Signal Process., vol. 16, no. 5, pp.
1040–1054, 2022.

[9]
Hadi Sarieddeen, Mohamed-Slim Alouini, and Tareq Y. Al-Naffouri,
“An overview of signal processing techniques for terahertz
Proc. IEEE, vol. 109, no. 10, pp. 1628–1665, 2021.

[10]
Stefano Fortunati, Fulvio Gini, Maria S Greco, and Christ D Richmond,
“Performance bounds for parameter estimation under misspecified
models: Fundamental findings and applications,”
IEEE Signal Process Mag., vol. 34, no. 6, pp. 142–157, 2017.

[11]
Cüneyd Öztürk, Musa Furkan Keskin, Henk Wymeersch, and Sinan
“On the impact of hardware impairments on RIS-aided
in Proc. IEEE Int. Conf. Commun. (ICC), 2022.

[12]
Hui Chen, Ahmed Elzanaty, Reza Ghazalian, Musa Furkan Keskin, Riku Jäntti,
and Henk Wymeersch,
“Channel model mismatch analysis for XL-MIMO systems from a
localization perspective,”
in Proc. IEEE GLOBECOM, 2022.

[13]
P-A Absil, Robert Mahony, and Rodolphe Sepulchre,
“Optimization algorithms on matrix manifolds,”
in Optimization Algorithms on Matrix Manifolds. Princeton
University Press, 2009.

[14]
Emil Björnson, Henk Wymeersch, Bho Matthiesen, Petar Popovski, Luca
Sanguinetti, and Elisabeth de Carvalho,
“Reconfigurable intelligent surfaces: A signal processing
perspective with wireless applications,”
IEEE Signal Process Mag., vol. 39, no. 2, pp. 135–158, 2022.

[15]
Pinjun Zheng, Tarig Ballal, Hui Chen, Henk Wymeersch, and Tareq Y Al-Naffouri,
“Coverage analysis of joint localization and communication in THz
systems with 3D arrays,”
TechRxiv preprint, 2022, doi: 10.36227/techrxiv.21385080.v1.

[16]
Arash Shahmansoori, Gabriel E. Garcia, Giuseppe Destino, Gonzalo Seco-Granados,
and Henk Wymeersch,
“Position and orientation estimation through millimeter-wave MIMO
in 5G systems,”
IEEE Trans. Wireless Commun., vol. 17, no. 3, pp. 1822–1835,

[17]
Hui Chen, Hadi Sarieddeen, Tarig Ballal, Henk Wymeersch, Mohamed-Slim Alouini,
and Tareq Y. Al-Naffouri,
“A tutorial on terahertz-band localization for 6g communication
IEEE Commun. Surv. Tutorials, vol. 24, no. 3, pp. 1780–1815,

[18]
Steven M Kay,
Fundamentals of statistical signal processing: estimation
Prentice-Hall, Inc., 1993.

[19]
Kaare Brandt Petersen, Michael Syskind Pedersen, et al.,
“The matrix cookbook,”
Technical University of Denmark, vol. 7, no. 15, 2008.

[20]
N. Boumal, B. Mishra, P.-A. Absil, and R. Sepulchre,
“Manopt, a Matlab toolbox for optimization on manifolds,”
Journal of Machine Learning Research, vol. 15, no. 42, pp.
1455–1459, 2014.