# A hybrid volume-surface integral equation method for rapid electromagnetic
simulations in MRI

Ilias I. Giannakopoulos, Georgy D. Guryev, José E. C. Serrallés, Jan Paška,
Bei Zhang,
Luca Daniel, , Jacob K. White, , Christopher M. Collins,
Riccardo Lattanzi This work was supported in part by research grants from the
National Institutes of Health (NIH P41 EB017183, and R01 EB024536) and the
National Science Foundation (NSF CAREER 1453675) in the United States, and was
performed under the rubric of the Center for Advanced Imaging Innovation and
Research (CAI2R, www.cai2r.net). (Corresponding author: Ilias I.
Giannakopoulos.)Ilias I. Giannakopoulos, Jan Paška, Christopher M. Collins and
Riccardo Lattanzi are with the Center for Advanced Imaging Innovation and
Research (CAI2R), and with the Bernard and Irene Schwartz Center for
Biomedical Imaging, Department of Radiology, New York University Grossman
School of Medicine, NY 10016 USA (e-mail:
ilias.giannakopoulos@nyulangone.org).Bei Zhang is with the Advanced Imaging
Research Center, UT Southwestern Medical Center, 2201 Inwood Rd, Dallas, TX
75390 USA.Georgy D. Guryev, José E. C. Serrallés, Luca Daniel, and Jacob K.
White are with the Research Laboratory of Electronics, Department of
Electrical Engineering and Computer Science, Massachusetts Institute of
Technology, Cambridge, MA 02139 USA.

###### Abstract

Objective: We developed a hybrid volume surface integral equation (VSIE)
method based on domain decomposition to perform fast and accurate magnetic
resonance imaging (MRI) simulations that include both remote and local
conductive elements. Methods: We separated the conductive surfaces present in
MRI setups into two domains and optimized electromagnetic (EM) modeling for
each case. Specifically, interactions between the body and EM waves
originating from local radiofrequency (RF) coils were modeled with the
precorrected fast Fourier transform, whereas the interactions with remote
conductive surfaces (RF shield, scanner bore) were modeled with a novel cross
tensor train-based algorithm. We compared the hybrid-VSIE with other VSIE
methods for realistic MRI simulation setups. Results: The hybrid-VSIE was the
only practical method for simulation using 1 mm voxel isotropic resolution
(VIR). For 2 mm VIR, our method could be solved at least 23 times faster and
required 760 times lower memory than traditional VSIE methods. Conclusion: The
hybrid-VSIE demonstrated a marked improvement in terms of convergence times of
the numerical EM simulation compared to traditional approaches in multiple
realistic MRI scenarios. Significance: The efficiency of the novel hybrid-VSIE
method could enable rapid simulations of complex and comprehensive MRI setups.

###### Index Terms:

Cross approximation, integral equation, magnetic resonance imaging, matrix
compression, radiofrequency coil, tensor train.

## Nomenclature

Notation

Description

$a$

Scalar

$\bm{a}$

Vector in $\mathbb{C}^{3}$

$\mathbf{a}$

Vector in $\mathbb{C}^{n}$

$A$

Matrix in $\mathbb{C}^{n_{1}\times n_{2}}$

$A^{\rm T}$

Transpose of matrix

$A^{*}$

Conjugate transpose of matrix

$\bar{A}$

Conjugate of matrix

${\mathcal{A}}$

Tensor in $\mathbb{C}^{n_{1}\times n_{2}\times n_{3}\times\dots}$

${\cal{A}}$

Operator acting on vectors in $\mathbb{C}^{3}$

$\mathfrak{G}$

Dyad

${\mathrm{i}}$

Imaginary unit ${\mathrm{i}}^{2}=-1$

##  I Introduction

Accurate electromagnetic (EM) modeling of radiofrequency (RF) coils and their
interaction with biological tissue is critical in magnetic resonance (MR)
imaging (MRI) [1]. For example, due to the high prototyping cost of multi-
element coil arrays [2, 3, 4, 5], it is desirable to optimize coil design in
simulation before construction. Coil optimization is particularly important at
ultra-high-field (UHF) MRI, since poorly designed RF coils could worsen image
quality and compromise patient safety, due to the RF field inhomogeneities and
spatially non-uniform specific absorption rate amplifications [6, 7, 8, 9, 10,
11].

Although MRI simulation extensions have been incorporated into commercial EM
simulation software, these tools could have high memory requirements and long
solution times for problems involving multi-channel coils, because they are
based on finite difference and finite element approaches. The computational
burden further increases for UHF MRI because for accurate simulations one also
needs to include the RF shield, due to its effect on the EM field and signal-
to-noise ratio (SNR) distributions [12]. The MAgnetic Resonance Integral
Equation (MARIE) suite [13] was developed to address the limitations of
general-purpose commercial software. MARIE is an open-source software
specifically tailored to model EM interactions between biological tissue and
RF coils at MR frequencies. MARIE employs the volume-surface integral equation
(VSIE) method, which is especially useful for single-frequency problems like
MRI, since the arising linear system of equations can be easily adapted to
existing iterative numerical linear algebra routines, and solved rapidly.
Algorithms implemented on top of MARIE have demonstrated improved convergence
times [14, 15, 16] and superior accuracy [17, 18] compared with conventional
methods. However, the solution time in MARIE can still considerably increase
at UHF when using fine voxel resolutions and including an RF shield for
accurate coil simulations.

To address this, one could take advantage of the low-rank properties of the
VSIE coupling matrix [19]. For example, an algorithm based on the precorrected
fast Fourier transform (pFFT) was presented in [20, 21] to accelerate MARIE in
cases where the conductors and the object are close to each other. Techniques
for matrix compression have instead been proposed to improve simulation
efficiency in the case of conductors placed away from the sample (e.g., RF
shield) [22]. In this work, we developed a hybrid-VSIE method for MARIE to
rapidly simulate realistic experimental setups, with conductors positioned at
various distances from the sample. The proposed method assembles a low-
parametric representation of the VSIE by combining multiple approaches: pFFT
[20], adaptive cross approximation (ACA) [23], Tucker decomposition [16], and
tensor train (TT)-cross [24]. As a result, the memory footprint is
considerably reduced and simulations can be rapidly executed on a graphical
processing unit (GPU) even for fine voxel resolutions.

The remainder of this article is organized as follows. In Section II, we
review the technical background on surface and volume integral equations
methods, as well as TT and Tucker decomposition methods. In section III we
introduce our proposed hybrid-VSIE. Section IV and Section V describe the
numerical experiments and associated results, respectively. A discussion of
the results is provided in Section VI, whereas Section VII summarizes the main
conclusions of this work.

##  II Technical Background

###  II-A Integral Equation Method in MRI

The EM phenomena in MR can be thoroughly explained with Maxwell’s equations
[25]. The single operating frequency of the MR scanner allows for the
development of fast and customizable algorithms for the solution of the
integral equation form of Maxwell’s equations. In such a representation, the
coil and the sample are replaced with equivalent surface and volumetric
current sources [26, 27] that radiate the same electromagnetic waves as the
coil and the body, respectively:

$\displaystyle\bm{j}_{\rm c}\\!\left(\bm{r}\right)$
$\displaystyle=\hat{n}\\!\left(\bm{r}\right)\times\bm{h}_{\rm
c}\\!\left(\bm{r}\right)$ (1) $\displaystyle\bm{j}_{\rm
b}\\!\left(\bm{r}\right)$
$\displaystyle={\mathrm{i}}\omega\epsilon_{0}\chi_{e}\\!\left(\bm{r}\right)\bm{e}_{\rm
b}\\!\left(\bm{r}\right).$

Here, $\bm{j}_{\rm c}\\!\left(\bm{r}\right)$ and $\bm{j}_{\rm
b}\\!\left(\bm{r}\right)$ are the equivalent surface and polarization body
currents, respectively, $\bm{r}$ is the observation position vector,
$\hat{n}\\!\left(\bm{r}\right)\times\bm{h}_{\rm c}\\!\left(\bm{r}\right)$ is
the magnetic field impressed on the coil’s surface, $f=\omega/(2\pi)$ is the
Larmor frequency of interest, and $\epsilon_{0}$ is the vacuum permittivity.
$\chi_{e}\\!\left(\bm{r}\right)=\epsilon_{r}\\!\left(\bm{r}\right)-1$,
$\epsilon_{r}\\!\left(\bm{r}\right)$, and $\bm{e}_{\rm
b}\\!\left(\bm{r}\right)$ are the electric susceptibility, relative
permittivity, and the total electric field of the object, respectively.

The surface integral equation (SIE) used to compute $\bm{j}_{\rm
c}\\!\left(\bm{r}\right)$ is the electric field integral equation [28]:

$\frac{{\mathrm{i}}\hat{n}\\!\left(\bm{r}\right)}{\omega\mu_{0}}\times\bm{e}^{\rm
inc}_{\rm
c}\\!\left(\bm{r}\right)=\hat{n}\\!\left(\bm{r}\right)\times\iint\limits_{\rm
c}\mathfrak{G}\\!\left(\bm{r},\bm{r}^{\prime}\right)\bm{j}_{\rm
c}\\!\left(\bm{r}^{\prime}\right)d^{2}\bm{r}^{\prime},$ (2)

where $\rm c$ refers to the coil’s surface, $\hat{n}\\!\left(\bm{r}\right)$ is
the normal vector on $\rm c$, $\mu_{0}$ is the vacuum permeability, and
$\bm{r}^{\prime}$ is the source position vector. $\bm{e}^{\rm inc}_{\rm
c}\\!\left(\bm{r}\right)$ is the incident electric field on $\rm c$ and could
be interpreted as a voltage excitation at a one of the feeding ports of the
coil, as in [29, 30]. $\mathfrak{G}(\bm{r},\bm{r}^{\prime})$ is the dyadic
homogeneous Green’s function [31]:

$\displaystyle\mathfrak{G}(\bm{r},\bm{r}^{\prime})$
$\displaystyle=\frac{1}{4\pi}\\!\left(\mathfrak{I}+\frac{\nabla\nabla}{k_{0}^{2}}\right)g(\bm{r},\bm{r}^{\prime})$
(3) $\displaystyle g(\bm{r},\bm{r}^{\prime})$
$\displaystyle=\frac{e^{-{\mathrm{i}}k_{0}\lvert\bm{r}-\bm{r}^{\prime}\rvert}}{4\pi\lvert\bm{r}-\bm{r}^{\prime}\rvert},$

where $\mathfrak{I}$ is the unit dyad, $k_{0}$ is the vacuum wave number, and
$g(\bm{r},\bm{r}^{\prime})$ is the free-space Green’s function.

The volume integral equation (VIE) used to compute $\bm{j}_{\rm
b}\\!\left(\bm{r}\right)$ is the current-based volume integral equation [14]:

$\\!\left({\cal{M}}_{\epsilon_{r}\\!\left(\bm{r}\right)}-{\cal{M}}_{\chi_{e}\\!\left(\bm{r}\right)}{\cal{N}}\right)\bm{j}_{\rm
b}\\!\left(\bm{r}\right)={\mathrm{i}}\omega\epsilon_{0}{\cal{M}}_{\chi_{e}\\!\left(\bm{r}\right)}\bm{e}^{\rm
inc}_{\rm b}\\!\left(\bm{r}\right).$ (4)

Here $\rm b$ refers to the object (“body”), ${\cal{M}}$ is a multiplication
operator that multiplies quantities with the parameter indicated in subscript,
$\bm{e}^{\rm inc}_{\rm b}\\!\left(\bm{r}\right)$ is the incident electric
field on $\rm b$, and ${\cal{N}}$ is the following Green’s function integro-
differential operator.

${\cal{N}}\\!\left(\bm{j}_{\rm
b}\\!\left(\bm{r}\right)\right)\triangleq\nabla\times\nabla\times\iiint\limits_{\rm
b}g\\!\left(\bm{r},\bm{r}^{\prime}\right)\bm{j}_{\rm
b}\\!\left(\bm{r}^{\prime}\right)d^{3}\bm{r}^{\prime}.$ (5)

When the coil is loaded with a body, the incident electric field on the body
$\bm{e}^{\rm inc}_{\rm b}\\!\left(\bm{r}\right)$ is the electric field
scattered from the coil to the body, which is directly related to the current
distribution on the coil $\bm{j}_{\rm c}\\!\left(\bm{r}\right)$. Since such
current distribution is perturbed by the back-scattered field from the body,
it cannot be simply computed by solving equation (2), as in the unloaded case.
To account for the back-scattering effect due to the EM interactions between
coil and body, one needs to solve (2) and (4) in tandem, by constructing a
coupled VSIE.

###  II-B Discretized VSIE Linear System

SIEs and VIEs are usually solved with the Galerkin Method of Moments (MoM)
[32], where the coils and the body are discretized. In this work, we followed
the techniques presented in [15] and built our method on top of the open-
source software MARIE [13]. Next, we briefly review the construction of the
VSIE discretization system of linear equations.

Coils are modeled with a triangular tessellation and the unknown $\bm{j}_{c}$
is approximated with the Rao-Wilton-Glisson (RWG) basis functions [28]. After
applying the Galerkin method, (2) becomes the linear system

$Z_{\rm cc}\mathbf{j}_{\rm c}=\mathbf{v},$ (6)

where $Z_{\rm cc}\in\mathbb{C}^{m\times m}$ is the discretized counterpart of
the integral in the right-hand-side of (2), where the coil is discretized into
$m$ elements. $\mathbf{v}\in\mathbb{C}^{m\times 1}$ and $\mathbf{j}_{\rm
c}\in\mathbb{C}^{m\times 1}$ represent the voltage excitation vector and the
unknown coil currents in each of the $m$ discretization elements of the coil.
For coarse resolutions of the coil mesh, multiplications involving $Z_{\rm
cc}$ are fast, because the matrix is small and dense. When the coil’s mesh
resolution is fine, or when additional conductive structures, such as the
scanner’s bore, or an RF shield, are included in the simulation, the size of
the $Z_{\rm cc}$ matrix can become large. In such cases, since the off-
diagonal blocks are low-rank [33], one can compress them with a matrix
decomposition method, for example, a pseudo-skeleton approximation [34, 35,
36, 37] or ACA [23, 38], to reduce $Z_{\rm cc}$’s memory demands.

The discretization of the VIE in (4) also results in a linear system:

$Z_{\rm bb}\mathbf{j}_{\rm b}=\mathbf{e}.$ (7)

Here, $Z_{\rm bb}\in\mathbb{C}^{q\cdot n_{\rm v}\times q\cdot n_{\rm v}}$ is
the discretized version of the left-hand-side of (4), $\mathbf{j}_{\rm
b}\in\mathbb{C}^{q\cdot n_{\rm v}\times 1}$ are the unknown body currents in
each of the $n_{\rm v}$ discretization elements of the body, for each of the
$q$ components of the basis used to approximate the polarization currents, and
$\mathbf{e}\in\mathbb{C}^{q\cdot n_{\rm v}\times 1}$ is the discretized
external incident field that illuminates the body. To approximate $\bm{j}_{\rm
b}\\!\left(\bm{r}\right)$ in each voxel, one can employ polynomial basis
functions, either piecewise constant (PWC) [14], which use $q=3$ unknowns per
discretization element, or piecewise linear (PWL) [18], for which $q=12$. In
either case, since the body is a 3D object, the matrix $Z_{\rm bb}$ is usually
too large to be fully assembled. One can address this by exploiting the
translational invariance property of $g\\!\left(\bm{r},\bm{r}^{\prime}\right)$
to enforce a multilevel block-Toeplitz structure on $Z_{\rm bb}$. Precisely,
if we enclose the body in a uniform voxelized domain of size $n_{1}\times
n_{2}\times n_{3}=n_{\rm v}$, then only the first columns of each of the
Toeplitz blocks of $Z_{\rm bb}$ need to be assembled and the computations
required for the EM simulations can be rapidly executed using the fast Fourier
transform (FFT) [39, 40]. In addition, these defining columns of $Z_{\rm bb}$
can be reshaped into $q\times q$ 3D tensors, each of dimensions $n_{1}\times
n_{2}\times n_{3}$, which can be vastly compressed with the Tucker
decomposition [16].

Since conductive surfaces (e.g., RF coils) and body are both present in the
MRI simulation setup, equations (6), (7) must be solved together through a
coupled block system of equations:

$\begin{bmatrix}Z_{\rm cc}&\\!\left(Z_{\rm cb}\right)^{\rm T}\\\ Z_{\rm
cb}&Z_{\rm bb}\end{bmatrix}\begin{bmatrix}\mathbf{j}_{\rm c}\\\
\mathbf{j}_{\rm b}\end{bmatrix}=\begin{bmatrix}\mathbf{v}\\\ 0\end{bmatrix}.$
(8)

Note that $\mathbf{e}$ from (7) was discarded, since it can be expressed as
the product between the new matrix $Z_{\rm cb}\in\mathbb{C}^{q\cdot n_{\rm
v}\times m}$ and $\mathbf{j}_{\rm c}$. $Z_{\rm cb}$ is called the coupling
matrix and maps the equivalent surface electric currents on the conductors
onto the electric fields that illuminate the body, via the dyadic Green’s
function. In contrast with $Z_{\rm cc}$ and $Z_{\rm bb}$, the coupling matrix
is not the outcome of a Galerkin inner product, but of a Petrov-Galerkin one,
because the testing and basis functions used for its construction are
different [27].

###  II-C Compression of the Coupling Matrix

One can reshape the coupling matrix to a higher-order tensor based on the
dimensions of the simulation domain. Specifically, each column of the coupling
matrix is readily convertible to a 4D tensor of dimensions $n_{1}\times
n_{2}\times n_{3}\times q$, while the whole coupling matrix can be interpreted
as a 5D tensor with dimensions $n_{1}\times n_{2}\times n_{3}\times m\times
q$. In this way, the matrix can be compressed with a tensor decomposition-
based algorithm of choice and large-scale EM simulations can be practically
performed [22].

The $5$-th dimension of the reshaped coupling matrix corresponds to the number
of basis functions per voxel and is not compressible for the functions of
choice (PWC and PWL). However, $q$ is small comparing to the remaining
dimensions, therefore the 5D tensor can be treated as $q$ independent 4D
tensors with low-rank properties [19]. To exploit the low-rankness and
compress the coupling matrix, one can consider two cases: a) Conductors at a
distance from the body, such as in the case of the scanner bore, the scanner
body coil [41, 21] or other RF shields [12], and b) conductors placed close to
the body, such as in the case of local receive coils [42, 5].

For case a), a 4D tensor decomposition method can be employed for the
compression of the coupling matrix [22, 43]. For case b), the $4$-th dimension
$m$ cannot be compressed as much as the first three dimensions, so a 4D tensor
decomposition becomes suboptimal. One can instead use $m$ independent 3D
tensor decompositions [22], which, however, yields slow matrix-vector product
times for fine coil discretizations (high values of $m$). Alternatively, it is
possible to define an extended voxelized domain of dimensions $n^{\rm
ext}_{1}\times n^{\rm ext}_{2}\times n^{\rm ext}_{3}$, with $n^{\rm
ext}_{i}\geq n_{i},\>i=1,2,3$, which encloses both coil and body, and then
project each triangular pair of the coil mesh onto an expansion domain of
voxels [44]. Then, the VSIE problem could be transformed into an equivalent
VIE problem using the pFFT method as in [20]. The resulting extended VIE
domain is a uniform voxelized grid, like the original VIE one, so the Green’s
function operators can be compressed with the methods presented in [16].

###  II-D Tensor Decomposition Methods

Tensor decomposition methods are able to provide a low-parametric
representation of specific $d$-dimensional arrays that possess low-rank
properties. Available methods include the canonical polyadic model [45, 46],
the Tucker decomposition [47, 48], the tensor train [24, 49], the hierarchical
Tucker decomposition [50], and the tensor ring decomposition [51], with
various applications in electrodynamics [52, 53, 54] and MRI [55, 56]. Next we
briefly review TT and Tucker decomposition since they are used in this work.

####  II-D1 Tensor Train (TT)

A $d$ dimensional tensor ${\mathcal{A}}\in\mathbb{C}^{n_{1}\times
n_{2}\times\cdots\times n_{d}}$ can be approximated by a tensor
$\tilde{{\mathcal{A}}}$ of the same dimensions through TT as:

$\tilde{{\mathcal{A}}}_{i_{1}i_{2}\cdots
i_{d}}=\sum\limits_{j_{1},j_{2},\cdots,j_{d-1}}^{r_{1},r_{2},\cdots,r_{d-1}}G^{1}_{i_{1}j_{1}}{\mathcal{G}}^{2}_{j_{1}i_{2}j_{2}}\cdots{\mathcal{G}}^{d-1}_{j_{d-2}i_{d-1}j_{d-1}}G^{d}_{j_{d-1}i_{d}}.$
(9)

Here $r_{k}$, for $k=1,\dots,d-1$ are the TT-ranks and the matrices
$G^{1}\in\mathbb{C}^{n_{1}\times r_{1}}$, $G^{d}\in\mathbb{C}^{r_{d-1}\times
n_{d}}$ and tensors ${\mathcal{G}}^{k}\in\mathbb{C}^{r_{k-1}\times n_{k}\times
r_{k}}$, $k=2,\dots,d-1$ are the cores or carriages of the TT.

For any tensor ${\mathcal{A}}$ with low TT-rank, one can always find a low-
rank quasi-optimal TT approximation with a prescribed accuracy, using a
sequence of SVDs on the unfolding matrices of ${\mathcal{A}}$ [49]. One could
also use a TT-cross approximation method [24, 57] to construct the TT of
${\mathcal{A}}$. In this case, the tensor ${\mathcal{A}}$ has to be replaced
with a function that maps from indices in ${\mathcal{A}}$ to corresponding
values in ${\mathcal{A}}$. Then, the TT-cross method, iteratively calls this
function as a black-box procedure in order to construct the TT form of
${\mathcal{A}}$ [24, 50].

When the TT ranks are small, TT-cross can be considerably faster than TT-SVD,
because the former does not require assembly of all elements of
${\mathcal{A}}$. This can be particularly useful when compressing the $Z_{\rm
cb}^{\rm f}$ coupling matrix, because each of its elements is a 5D surface-
volume integral, which can be costly for a large number of quadrature
integration points.

####  II-D2 Tucker Decomposition

Tucker decomposition approximates the previously mentioned tensor
${\mathcal{A}}$ as:

$\tilde{{\mathcal{A}}}_{i_{1}i_{2}\cdots
i_{d}}\approx\sum\limits_{j_{1},j_{2},\cdots,j_{d}}^{r_{1},r_{2},\cdots,r_{d}}{\mathcal{G}}_{j_{1}j_{2}\cdots
j_{d}}U^{1}_{i_{1}j_{1}}U^{2}_{i_{2}j_{2}}\cdots U^{d}_{i_{d}j_{d}}.$ (10)

Here, ${\mathcal{G}}\in\mathbb{C}^{r_{1}\times r_{2}\times\cdots\times r_{d}}$
is the Tucker core, $U^{k}\in\mathbb{C}^{n_{k}\times r_{k}}$, $k=1,\dots,d$
are the Tucker factors, and $r_{k}$, $k=1,\dots,d$ are the Tucker ranks of
${\mathcal{A}}$. Tucker decomposition provides an excellent approximation of
${\mathcal{A}}$, with an upper error bound guarantee, if computed with the
higher-order SVD (HOSVD) [58] or the Cross3D algorithm [48, 59].

##  III Hybrid Volume-Surface Integral Equation

In this work, we propose a hybrid volume-surface integral equation approach
that combines pFFT [20], Tucker decomposition [16], TT-cross [24] and ACA [38]
to optimize the assembly and the solution of the VSIE system.

According to [60], the $2\times 2$ block system in (8) can be interpreted as
the outcome of a domain decomposition method. Building on that, we can expand
such system into a $3\times 3$ one, by separating the conductive surfaces into
a “near” ($\rm n$) and a “far” ($\rm f$) domain, based on their distance from
the sample. For a typical MRI simulation setup, the RF shield and the
scanner’s bore belong to the “far” domain, while the receive array is usually
in the “near” domain (Fig. 1). The domain where transmit arrays are placed can
vary, depending on their position relative to the body, and should be
judiciously chosen for each application.

Applying the proposed domain decomposition, the system of equations in (8)
becomes:

$\begin{bmatrix}Z^{\rm ff}_{\rm cc}&\\!\left(Z^{\rm fn}_{\rm cc}\right)^{\rm
T}&\\!\left(Z^{\rm f}_{\rm cb}\right)^{\rm T}\\\ Z^{\rm fn}_{\rm cc}&Z^{\rm
nn}_{\rm cc}&\\!\left(Z^{\rm n}_{\rm cb}\right)^{\rm T}\\\ Z^{\rm f}_{\rm
cb}&Z^{\rm n}_{\rm cb}&Z_{\rm bb}\end{bmatrix}\begin{bmatrix}\mathbf{j}^{\rm
f}_{\rm c}\\\ \mathbf{j}^{\rm n}_{\rm c}\\\ \mathbf{j}_{\rm
b}\end{bmatrix}=\begin{bmatrix}\mathbf{v}_{\rm f}\\\ \mathbf{v}_{\rm n}\\\
0\end{bmatrix}.$ (11)

Here, $\mathbf{j}_{\rm c}^{\rm n},\mathbf{v}_{\rm n}\in\mathbb{C}^{m_{n}\times
1}$ and $\mathbf{j}_{\rm c}^{\rm f},\mathbf{v}_{\rm
f}\in\mathbb{C}^{m_{f}\times 1}$ are the unknown currents and voltage
excitation vectors in the coils placed in the “near” and “far” domain,
respectively. Note that $\mathbf{v}_{\rm n}$ or $\mathbf{v}_{\rm f}$ can be
equal to zero, if no excitations exist in the corresponding domain. $Z^{\rm
ff}_{\rm cc}\in\mathbb{C}^{m_{f}\times m_{f}}$ and $Z^{\rm nn}_{\rm
cc}\in\mathbb{C}^{m_{n}\times m_{n}}$ are the Galerkin matrices that model the
interactions between the $m_{f}$ “far” and the $m_{n}$ “near” surface
discretization elements, respectively. $Z^{\rm fn}_{\rm
cc}\in\mathbb{C}^{m_{n}\times m_{f}}$ is the off-diagonal block of $Z_{\rm
cc}$ and models the coupling interactions between conductors in the “far” and
“near” domains. Since $Z^{\rm fn}_{\rm cc}$ has low-rank properties, it can be
compressed with the ACA method, as in [38]:

$UV^{*}=\text{ACA}\\!\left(Z^{\rm fn}_{\rm cc}\right).$ (12)

Returning to (11), $Z^{\rm n}_{\rm cb}\in\mathbb{C}^{m_{n}\times q\cdot
n_{v}}$ and $Z^{\rm f}_{\rm cb}\in\mathbb{C}^{m_{f}\times q\cdot n_{v}}$ are
the coupling matrices that model the interactions between the body and the
conductive surfaces placed on the “near” and “far” domain, respectively.

The $2\times 2$ sub-system in (11) that models the interactions between the
body and the conductors placed in the “near” domain can be solved efficiently
with the pFFT method [20]:

$Z^{\rm pFFT}_{\rm HOSVD}=\begin{bmatrix}Z^{\rm nn}_{\rm cc}&\\!\left(Z^{\rm
n}_{\rm cb}\right)^{\rm T}\\\ Z^{\rm n}_{\rm cb}&Z_{\rm bb}\end{bmatrix}.$
(13)

Precisely, the pFFT is used to project the mesh of the coil located in the
“near” domain onto a uniform voxelized grid that encloses both the coil’s
triangular patches and the voxels that discretize the sample [20]. The
resulting $Z^{\rm pFFT}_{\rm HOSVD}$ matrix has a multilevel block-Toeplitz
structure and resembles the VIE integral equation matrix $Z_{\rm bb}$ of (7).
The subscript HOSVD indicates that the Tucker decomposition method can be
leveraged to compress the matrix and perform rapid multiplications in the case
of fine voxel resolutions, as in [16].

A previous study proposed performing a block-wise assembly of $Z_{\rm cb}^{\rm
f}$ and then compress each block with TT-SVD [43] or HOSVD [22]. In this work,
instead of assembling the coupling matrix, we used a function that maps from
indices to corresponding values in the matrix. Then, we calculated a low-
parametric approximation of the coupling matrix in $q$ 4D tensors, which
required only a fraction of its entries, using a tensor decomposition and a
cross approximation algorithm that iteratively calls the aforementioned
function as a black-box procedure. Precisely, we used the TT-cross method [24]
based on the adaptive density matrix renormalization group [57]:

$\displaystyle\tilde{{\mathcal{Z}}^{\chi}_{\rm cb}}$
$\displaystyle=\text{TT}\\!\left({\mathcal{Z}}^{\chi}_{\rm
cb}\right),\>\chi=1,\dots,q,$ (14)
$\displaystyle\tilde{{\mathcal{Z}}^{\chi}_{\rm cb}}_{,\>i_{1}i_{2}\cdots
i_{4}}$
$\displaystyle=\sum\limits_{j_{1},j_{2},j_{3}}^{r_{1},r_{2},r_{3}}G^{1}_{i_{1}j_{1}}{\mathcal{G}}^{2}_{j_{1}i_{2}j_{2}}{\mathcal{G}}^{3}_{j_{2}i_{3}j_{3}}G^{4}_{j_{3}i_{4}},$

where ${\mathcal{Z}}^{\chi}_{\rm cb}$ is the reshaped $\chi$-th 4D tensor
component of $Z^{\rm f}_{\rm cb}$, and $\tilde{{\mathcal{Z}}^{\chi}_{\rm cb}}$
is the TT compressed ${\mathcal{Z}}^{\chi}_{\rm cb}$ that is computed using
(9). Note that we never fully assemble $\tilde{{\mathcal{Z}}^{\chi}_{\rm
cb}}$, but we rather express it using the cores of its TT:
$G^{1},{\mathcal{G}}^{2},{\mathcal{G}}^{3}$, and $G^{4}$.

By combining all the aforementioned methods, we can transform (11) into a
hybrid-VSIE system with low memory footprint that can be solved rapidly:

$\scriptstyle\begin{bmatrix}Z^{\rm ff}_{\rm cc}&\bar{V}U^{\rm
T}&\left[\tilde{{\mathcal{Z}}^{1}_{\rm
cb}},\vdots,\tilde{{\mathcal{Z}}^{q}_{\rm cb}}\right]^{\rm T}\\\ UV^{*}&&\\\
\left[\tilde{{\mathcal{Z}}^{1}_{\rm cb}},\vdots,\tilde{{\mathcal{Z}}^{q}_{\rm
<EMAIL_ADDRESS>pFFT}_{\rm
HOSVD}$}\hfil\lx@intercol\end{bmatrix}\begin{bmatrix}\mathbf{j}^{\rm f}_{\rm
c}\\\ \mathbf{j}^{\rm n}_{\rm c}\\\ \mathbf{j}_{\rm
b}\end{bmatrix}=\begin{bmatrix}\mathbf{v}_{\rm f}\\\ \mathbf{v}_{\rm n}\\\
0\end{bmatrix}.$ (15)

###  III-A Solving the hybrid-VSIE

Since it is not practical to directly invert the hybrid-VSIE matrix in (15),
even for coarse resolutions, an iterative solver, like the generalized minimal
residual method (GMRES) [61], should be used. Multiplications involving
$Z^{\text{pFFT}}_{\text{HOSVD}}$ can be carried out efficiently using the
methods presented in [16]. Matrix-vector products involving $Z^{\rm ff}_{\rm
cc}$, $UV^{*}$, and $\bar{V}U^{\rm T}$ can be performed with standard matrix
multiplications. For the multiplications involving the TT-compressed $Z^{\rm
f}_{\rm cb}$ we propose the following approach.

Let us consider the product $[\tilde{{\mathcal{Z}}^{1}_{\rm
cb}},\smash{\vdots},\tilde{{\mathcal{Z}}^{q}_{\rm cb}}]\mathbf{j}_{\rm c}^{\rm
f}$. We can re-write it as:

$\left[\tilde{{\mathcal{Z}}^{1}_{\rm cb}},\vdots,\tilde{{\mathcal{Z}}^{q}_{\rm
cb}}\right]\mathbf{j}_{\rm c}^{\rm f}=\left[\tilde{{\mathcal{Z}}^{1}_{\rm
cb}}\mathbf{j}_{\rm c}^{\rm f},\vdots,\tilde{{\mathcal{Z}}^{q}_{\rm
cb}}\mathbf{j}_{\rm c}^{\rm f}\right].$ (16)

We can decompose $\tilde{{\mathcal{Z}}^{\chi}_{\rm cb}}$ using (14) and then
perform the product $\tilde{{\mathcal{Z}}^{\chi}_{\rm cb}}\mathbf{j}_{\rm
c}^{\rm f}$ in a multiple step process, for each $\chi\in{1,\dots,q}$. First
we can reshape ${\mathcal{G}}^{2}\in\mathbb{C}^{r_{1}\times n_{2}\times
r_{2}}\rightarrow G^{2}\in\mathbb{C}^{(r_{1}\cdot n_{2})\times r_{2}}$ and
${\mathcal{G}}^{3}\in\mathbb{C}^{r_{2}\times n_{3}\times r_{3}}\rightarrow
G^{3}\in\mathbb{C}^{(r_{2}\cdot n_{3})\times r_{3}}$, where $\rightarrow$
represents the reshape operation. Then, we can perform the following steps.

$\displaystyle\mathbf{y}$ $\displaystyle=G^{4}\mathbf{j}_{\rm c}^{\rm
f},\>\>\mathbf{y}\in\mathbb{C}^{r_{3}\times 1},$ (17)
$\displaystyle\mathbf{y}$
$\displaystyle=G^{3}\mathbf{y},\>\>\mathbf{y}\in\mathbb{C}^{r_{2}\cdot
n_{3}\times 1},\>\>\>\>\>\mathbf{y}\rightarrow Y\in\mathbb{C}^{r_{2}\times
n_{3}},$ $\displaystyle Y$
$\displaystyle=G^{2}Y,\>\>Y\in\mathbb{C}^{r_{1}\cdot n_{2}\times
n_{3}},\>\>Y\rightarrow Y\in\mathbb{C}^{r_{1}\times n_{2}\cdot n_{3}},$
$\displaystyle Y$ $\displaystyle=G^{1}Y,\>\>Y\in\mathbb{C}^{n_{1}\times
n_{2}\cdot n_{3}},\>\>Y\rightarrow\mathbf{y}\in\mathbb{C}^{n_{1}\cdot
n_{2}\cdot n_{3}\times 1}.$

Note that we implemented first the product $G^{4}\mathbf{x}$ to eliminate the
fourth and largest dimension of ${\mathcal{Z}}^{\chi}_{\rm cb}$, in order to
reduce the memory footprint and the operation cost of the remaining
multiplications.

To compute the product $[\tilde{{\mathcal{Z}}^{1}_{\rm
cb}},\smash{\vdots},\tilde{{\mathcal{Z}}^{q}_{\rm cb}}]^{\rm T}\mathbf{j}_{\rm
b}$, we applied the transpose operator:

$\mathbf{j}_{\rm b}^{\rm T}\left[\tilde{{\mathcal{Z}}^{1}_{\rm
cb}},\vdots,\tilde{{\mathcal{Z}}^{q}_{\rm
cb}}\right]=\sum\limits_{\chi=1}^{q}\\!\left(\\!\left(\mathbf{j}_{\rm
b}^{\chi}\right)^{\rm T}\tilde{{\mathcal{Z}}^{\chi}_{\rm cb}}\right)^{\rm T},$
(18)

where we separated $\mathbf{j}_{\rm b}$ into its vector components
$\mathbf{j}_{\rm b}=[\mathbf{j}_{\rm b}^{1},\smash{\vdots},\mathbf{j}_{\rm
b}^{q}]$, with each $\mathbf{j}_{\rm b}^{\chi}\in\mathbb{C}^{n_{1}\cdot
n_{2}\cdot n_{3}\times 1}$. The multiplications in (18) can be broken down
into concatenated steps as in (17):

$\displaystyle\mathbf{j}_{\rm b}^{\chi}$ $\displaystyle\rightarrow
Y\in\mathbb{C}^{n_{3}\cdot n_{2}\times n_{1}},$ (19) $\displaystyle Y$
$\displaystyle=YG^{1},\>\>Y\in\mathbb{C}^{n_{3}\cdot n_{2}\times
r_{1}},\>\>Y\rightarrow Y\in\mathbb{C}^{n_{3}\times r_{1}\cdot n_{2}},$
$\displaystyle Y$ $\displaystyle=YG^{2},\>\>Y\in\mathbb{C}^{n_{3}\times
r_{2}},\>\>\>\>\>\>\>Y\rightarrow\mathbf{y}\in\mathbb{C}^{1\times n_{3}\cdot
r_{2}},$ $\displaystyle\mathbf{y}$
$\displaystyle=\mathbf{y}G^{3},\>\>\mathbf{y}\in\mathbb{C}^{1\times r_{3}},$
$\displaystyle\mathbf{y}$
$\displaystyle=\mathbf{y}G^{4},\>\>\mathbf{y}\in\mathbb{C}^{1\times m_{f}}.$

##  IV Methods

###  IV-A Realistic Head Models

We evaluated the performance of our proposed hybrid-VSIE using the realistic
Billie and Duke head models from the virtual family population [62]. Both
models were discretized over a regular voxelized grid. Billie was discretized
with $5$, $2$, and $1$ mm voxel isotropic resolutions (VIR), corresponding to
a simulation domain of $34\times 38\times 45$, $84\times 96\times 116$, and
$168\times 188\times 222$ voxels, respectively. Duke was discretized with a
$5$ mm VIR ($38\times 47\times 46$ voxels).

###  IV-B Simulation Setups

####  IV-B1 Triangular head coil

We modeled the $7$ tesla transceiver triangular head array introduced in [63]
with a surrounding cylindrical perfectly electric conducting shield, loaded
with Billie head model. The shield had a $20$ cm radius, was $35$ cm long, and
was placed $1.7$ cm away in the $\hat{z}$ direction from the bottom end ring
of the coil array as shown in Fig. 1. Details on the coil’s geometry, tuning,
and decoupling can be found in [63]. For the case when Billie was discretized
with $5$ and $2$ mm VIR, the array and the RF shield were discretized with
$2140$ and $2662$ triangular elements, respectively. When Billie was
discretized with $1$ mm VIR, the RF shield’s discretization remained the same,
whereas the array’s mesh was refined to $6054$ triangular elements. The
discretization of the conductors was adjusted based on the body resolution
because for pFFT to be applicable, each pair of triangular elements must be
contained by the voxel expansion domain of pFFT, which was set to $5\times
5\times 5$ voxels.

Figure 1: Shielded 8-element triangular coil array loaded with the Billie head
model. The cyan volume defines the domain “near” the sample, which contains
the coil, whereas the shield belongs to the “far” domain, outside the cyan
region

.

####  IV-B2 Single loop

To qualitatively replicate previous results obtained with a commercial EM
simulation software [12], we modeled a single surface loop placed behind the
Duke head model (Fig. 2) at $3$, $7$, $9.4$, and $10.5$ tesla. The radius of
the loop was $4.15$ cm, the conductor’s width was $0.3$ cm and was discretized
with $1958$ elements. The loop was segmented with seven capacitors for tuning
and one capacitor connected in parallel to the feeding port for matching. The
loop and head were surrounded by a $50$ cm long RF shield, while the radius
was varied from $17.25$ to $32$ cm ($16$ positions, $1$ cm steps), for which
the number of discretization elements ranged from $3264$ to $5824$.

Figure 2: (left) RF shield encircling a loop loaded with the Duke head model.
(right) Details of the loop’s geometry. The figure shows a representative RF
shield radius.

####  IV-B3 Overlapping transceiver array

We modeled an eight-channel transceiver array with overlapping coil elements
(Fig. 3) at $9.4$ tesla. The array consisted of two sets of four identical
rectangular loops of $11.88$ cm width and $8$ cm length, arranged on two
cylindrical surfaces of length $8$ cm, and radii $12$ and $12.5$ cm. The four
loops of each set were interleaved by $45$ degrees. The neighboring coils
overlapped by $3$ cm to ensure geometric decoupling. The coil conductors were
$0.5$ cm wide and were discretized with $2792$ triangular elements. Each coil
was segmented with eight variable capacitors for tuning and one variable
capacitor connected in parallel at the feeding port for matching. The array
was surrounded by the same RF shields as in the single-loop setup and also
loaded with the Duke head model. For this setup, we modeled also the scanner’s
bore as an open cylindrical surface of $277.75$ cm length and $32.05$ cm
radius, which was discretized with $5596$ triangular patches.

Figure 3: (left) Transceiver array loaded with the Duke head, surrounded by an
RF shield and the scanner’s bore. Note that the bore’s length is cropped in
the figure for enhanced visualization and a representative radius is displayed
for the RF shield. (right) Details of the array’s geometry.

###  IV-C Numerical Experiments

All simulations were performed on a server running Ubuntu 18.04.2 LTS
operating system, with an Intel(R) Xeon(R) Gold 6248R CPU at 2.70GHz, 112
cores, 2 threads per core, and an NVIDIA A100 PCIe GPU with 40GB of memory.
All matrix-vector products were implemented in MATLAB 9.10 (The MathWorks
Inc., Natick, MA).

####  IV-C1 Comparison with Non-Hybrid Methods

We compared our proposed hybrid-VSIE (H) method with three non-hybrid-VSIE
approaches that solved the system of equations (8), for the setup in Fig. 1.
For the first non-hybrid method (I), the coupling matrix between the coil
array, the shield, and the head model was assembled in its full form. For the
second method (II), the coupling matrix was assembled in a compressed form
with the ACA method. Note that we used methods I and II for the comparison
because they are implemented in MARIE. For the third method (III), pFFT was
used to project both the coil and the shield to an extended voxelized domain
as in [20]. For the hybrid-VSIE method, the coil, tightly-fitting the head
model, was placed in the “near” domain, while the shield was placed in the
“far” domain. The polarization body currents were modeled with PWC and PWL
basis functions for comparison.

For all simulations, the system was solved with the GMRES, executed either on
GPU (G) or CPU (C) depending the memory requirements. The tolerance was set to
$10^{-5}$. GMRES’ restarts were set to $50$ for all cases, except for the $1$
mm VIR-PWL basis functions case, where they were reduced to $10$ to avoid GPU
memory overflows. The tolerances for TT-cross and ACA in the hybrid-VSIE
method were $10^{-3}$. The ACA tolerance in method II was set to $1$-$2$
orders of magnitude lower than the GMRES tolerance. The HOSVD tolerance was
$10^{-5}$ for all methods. For methods I,II, and III, we used $8$ volume and
$5$ surface quadrature integration points to approximate the coupling matrix
integrals. In the hybrid-VSIE, we used $8$ and $1$ volume, and $5$ and $1$
surface integration points for the computations of $Z_{\rm bc}^{\rm n}$ and
$Z_{\rm bc}^{\rm f}$, respectively.

####  IV-C2 Effect of the Shield

We used the hybrid-VSIE method to study the effect of the RF shield’s position
on the SNR distribution. For the loop in Fig. 2, we repeated $16$ simulations
with the hybrid-VSIE (one for each shield’s radius) for four main magnetic
field strengths and both PWC and PWL basis functions at $5$ mm VIR (128
simulations in total). For the array in Fig. 3, we used the hybrid-VSIE method
to perform simulations at $9.4$ tesla for all $16$ shield’s radii, both with
and without the scanner’s bore. To complement geometric decoupling between
first-order neighbors, the capacitor values were adjusted to further reduce
the relative scattering (S) parameters. We calculated the SNR [64] and
averaged its value over all voxels of the head. Both the loop and the multiple
channel array were re-tuned and matched for all shield’s positions and all
field strengths. Coil elements were placed in the “near” domain, while the RF
shield and the scanner’s bore were in the “far” domain. The polarization body
currents were modeled with PWC and PWL basis functions for the single-loop
case, and with PWC only for the array case. We used the same tolerances and
number of quadrature points as in the first experiment.

##  V Numerical Results

###  V-A Comparison with Non-Hybrid Methods

TABLE I compares the convergence time and iterations ($\left\lceil\rm
hh:mm:ss\right\rceil$ \- $\left\lceil\rm iter\right\rceil$), averaged over the
eight solutions of the system (one solution for each channel of the triangular
coil), for the hybrid-VSIE against methods I, II, and III. The bold font
indicates the faster solver in each case. Methods I, II, and III were not
applicable in all cases, either due to memory limitations or because the
resulting VSIE matrix from (8) was ill-conditioned.

TABLE I: Average GMRES Convergence Properties VIR | VSIE | PWC | PWL
---|---|---|---
hh:mm:ss - iter - PU | hh:mm:ss - iter - PU
$5$ mm | H |  00:00:03 \- 136 - G |  00:00:19 \- 165 - G
I | 00:00:06 - 136 - G | 00:04:25 - 165 - C
II | 00:00:04 - 136 - G | 00:00:22 - 165 - G
III | 00:00:04 - 170 - G | 00:00:24 - 199 - G
$2$ mm | H |  00:00:16 \- 177 - G |  00:01:37 \- 193 - G
I | 00:18:45 - 177 - C | N/A
II | 00:06:17 - 177 - C | N/A
$1$ mm | H |  00:02:57 \- 192 - G |  01:49:12 \- 327 - G

TABLE II shows the $\rm L_{2}$ norm relative differences between the solutions
of methods H, II, and III with respect to method I, which was assumed as the
gold-standard because it uses the fully assembled system. The comparison was
possible only for the $5$ mm VIR and both PWC and PWL basis functions, and $2$
mm VIR using only the PWC basis. In the remaining cases, method I could not be
implemented on our server due to memory overflow.

TABLE II: Relative differences with respect to Method I Method | PWC $5$ mm | PWL $5$ mm | PWC $2$ mm
---|---|---|---
H | $0.08\%$ | $2.76\%$ | $0.52\%$
II | $0.36\%$ | $0.50\%$ | $0.95\%$
III | $0.42\%$ | $0.38\%$ | N/A

A qualitatively comparison between PWC and PWL basis functions is shown in
Fig. 4 for $1$ mm VIR, for the circularly polarized (CP) mode of the array, in
which the eight transmit fields of the coil elements were combined to cancel
the phase at the center of the head.

Figure 4: Magnitude of transmit field of the triangular coil array in CP mode
combination. Results are shown for $1$ mm VIR, for PWC (top) and PWL (bottom)
basis functions, and three representative orthogonal slices of the head model.
The values are masked outside Billie for an enhanced visualization.

###  V-B Effect of the Shield

####  V-B1 Single Loop

The left axis in Fig. 5 shows the normalized average value of the SNR inside
the head as a function of the RF shield’s diameter. The right axis shows the
corresponding compression factor of $Z_{\rm bc}^{\rm f}$, computed as the
number of elements of the full matrix divided by the number of elements of the
TT-compressed one.

Figure 5: (left axis, blue) Normalized average SNR as a function of the radius
of the RF shield for different main magnetic field strengths and basis
functions (PWC and PWL). (right axis, red) Compression factor of $Z_{\rm
bc}^{\rm f}$. for each case. Note that the loop was re-tuned for each
simulation.

####  V-B2 8-element Array

The worst S parameter between first-order neighbors was $-17.3$ dB for $620$
mm shield diameter (w/o the bore), whereas coupling was slightly higher
between next-nearest and opposite neighbors, with the worst case being $-12.7$
dB between two opposite coils for $640$ mm shield diameter (w/o the bore).
Fig. 6 shows the normalized average SNR (left axis) and the compression factor
of $Z_{\rm bc}^{\rm f}$ (right axis), as a function of the shield’s diameter.

Figure 6: Normalized average SNR of the head voxels in respect to the radius
of the RF shield. The array was re-tuned and decoupled for all cases.

##  VI Discussion

The aim of this work was to introduce a hybrid-VSIE technique that exploits
the low-rank properties of the memory-intensive parts of the VSIE to enable
rapid EM simulations for realistic MRI scenarios. The separation of the
conductive surfaces into a “near” and a “far” domain yielded exceptional
numerical performances for the hybrid-VSIE case. In particular, the cross-TT
considerably reduced the memory footprint of the $Z^{\rm f}_{\rm cb}$ matrix.
This is useful when the “far” domain is represented by an RF shield, since the
compressed matrix can be stored and recycled for multiple EM simulations
involving different RF coils (for example, the same $Z^{\rm f}_{\rm cb}$ could
be used for the results in V.B.1 and V.B.2), or inside coil design
optimization pipelines [65].

While one could assemble $Z^{\rm f}_{\rm cb}$ using SVD or ACA instead of TT
for compression, the memory requirements would be larger [22]. TT-cross was
preferred over TT-SVD, because it does not require the full or a block-wise
assembly of $Z^{\rm f}_{\rm cb}$ prior to the decomposition - a task that can
be lengthy for fine VIR or PWL basis functions [43]. We chose to use ACA in
method II because such approach is implemented in MARIE. Alternatively, one
could use TT in method II to assemble the whole coupling matrix $Z_{\rm cb}$
as in [43]. However, the assembly time of the compressed matrix $Z_{\rm cb}$,
and the time to compute the relevant matrix-vector products (17), (18) could
increase, because the local interactions between the conductive surfaces and
the body could lead to higher matrix ranks, therefore the hybrid-VSIE was
deemed a superior choice. Further investigation is needed when the conductive
surfaces are neither near nor far from the body, for example, when simulating
body-transmit coils. In these cases, as in setups where all conductive
surfaces are near or far, the hybrid-VSIE would not be advantageous and one
could use either pFFT or TT-cross on (8) instead. The criteria to choose the
optimal method would depend on the distance of the coil from the body, and the
discretization chosen for the body and the coil. In future work, one could
also explore other approaches to achieve further acceleration, for example,
the use of the Calderón preconditioner [66] to reduce the number of iterations
of GMRES.

The hybrid-VSIE was the only practical method in the case of fine VIR. Method
III required a larger number of iterations for the convergence of GMRES even
in the case of $5$ mm VIR, because each triangular pair of the RF shield was
not small enough to be projected accurately onto the $5\times 5\times 5$ voxel
expansion domain of pFFT [44, 20], resulting in a ill-conditioned system. One
way to tackle this would be to enlarge the expansion domain of pFFT, at the
cost of additional assembly and solution times. At $1$ mm VIR, the hybrid-VSIE
became ill-conditioned like method III for the chosen pFFT expansion domain,
thus GMRES, required a larger number of iterations to converge. To address
this, instead of enlarging pFFT’s expansion domain, we refined the triangular
mesh of the array, leading to an $\sim\\!8$ times larger $Z_{\rm cc}^{\rm nn}$
matrix. Note that in method III, the same refinement would increase the size
of $Z_{\rm cc}$ by $\sim\\!1400$ times. Methods I and II required additional
memory for the coupling matrix assembly and multiplications, so they became
slow or impractical for fine VIR. In the case of $2$ mm VIR resolution and PWL
basis functions, one could assemble the coupling matrix $Z_{\rm cb}$ in
methods I and II using a point-matching integration (thus, neglecting the
contribution of the linear terms), instead of the Petrov-Galerkin
approximation, reducing its memory requirements by $4$, at the cost of lower
accuracy in the solution.

The hybrid-VSIE allowed us to use higher tolerances and fewer quadrature
integration points for the Galerkin matrices that involved remote interactions
of the body with conductors in the “far” domain. This was possible because
such interactions are less pronounced than the ones in the “near” domain. As a
result, the hybrid-VSIE method remained accurate ($2.76\%$ relative difference
of the solution with respect to method I) for PWL basis functions, even though
the contribution of the linear terms in $Z_{\rm bc}^{\rm f}$ was neglected. By
using more quadrature points to include linear term contributions, the
difference becomes smaller than $0.1\%$.

TABLE I shows that the PWL approximation required additional time to converge
to the desired tolerance. Nevertheless, as suggested by Fig. 3, PWL are needed
for accurate modeling of the magnetic field density [18]. Such level of
accuracy may not be needed for qualitative coil assessment, but is necessary
for other quantitative applications [67].

The marked improvement of the convergence time of hybrid-VSIE allowed us to
perform optimization of the RF shield’s position with respect to SNR. In
section V.B.1, we verified the accuracy of hybrid-VSIE by repeating the same
series of numerical simulations reported in [12]. The loop’s and shield’s
positions and geometries were exactly the same in both series of simulations,
while the head phantom was placed $1$ cm higher in the hybrid-VSIE simulation.
For the $3$T and $7$T cases, the results matched, with minor inconsistencies.
For $9.4$T and $10.5$T, the behavior of the average SNR as a function of the
shield diameter was similar, except for the location of the peak values, which
occurred for a $30$ mm smaller shield radius for the hybrid-VSIE. Such minor
discrepancies are expected, since in [12] the EM field was computed with a
different approach based on the finite integration technique.

In section V.B.2 we used the hybrid-VSIE to study the performance of a multi-
channel array with respect to the RF shield’s size at $9.4$ tesla. In
particular, we observed approximately $7\%$ difference between the highest
($345$ mm shield diameter) and lowest ($600$ mm shield diameter) average SNR
value. The inclusion of the scanner’s bore in the simulations had a negligible
effect on the average SNR ($\leq 3\%$). Note that, as in the single loop case,
the SNR value did not change monotonically with the shield diameter at $9.4$
tesla. Note that in the case of the array, the behavior of the average SNR as
a function of shield diameter was not smooth. We believe this is mainly due to
two reasons: 1) the RF interference patterns generated by the array elements
change based on the distance from the shield to the coils, as shown in [12] 2)
the tuning and coupling of each coil element were not identical for every
case. Specifically, after re-tuning for each shield position, the S parameters
of the array elements varied by up to a few dB, which affected the resulting
fields. This effect was instead negligible for the single loop, since the S
parameter was $\ll-50$ dB for all shield diameters.

The hybrid-VSIE is a domain-decomposition method and therefore carries the
limitations of the methods used to model each domain. Specifically, pFFT
requires that each triangular pair of the conductors in the “near” domain fits
in the pFFT voxel expansion domain in order to be performed. TT-cross
algorithm might overestimate the TT-ranks, which might lead to unnecessary
additional memory allocation. Note that this is not a practical limitation for
the simulation’s performance, since the coupling matrix could still be
compressed by a large factor with TT-cross. However, if one is interested in
assembling a large dataset consisting of multiple $Z^{\rm f}_{\rm cb}$
matrices for RF shield optimization pipelines, the overestimated TT-ranks
might lead to storage concerns. Finally, since the equivalent shield and coil
currents and the polarization body currents are solved in a coupled fashion in
our approach, the choice of the tolerance value for the different techniques
used in hybrid-VSIE is important to avoid false convergence of the iterative
solver. In fact, since the coil currents are normally larger than the body and
shield currents, by using a single tolerance, the accuracy of the latter could
be lower than the accuracy of the coil currents. While this is expected to
have a negligible effect on coil simulations, one could ensure the same
accuracy by solving first for the polarization currents, and then computing
the equivalent coil currents, as a post-processing step. Such approach would
allow further flexibility in the choice of the iterative’s solver tolerance.

##  VII Conclusion

We presented a new hybrid-VSIE method suitable for rapid and comprehensive EM
simulations in MRI. Our proposed method combines TT, pFFT, Tucker
decomposition, and the ACA to heavily compress the coupling interactions
between coils, shields, bore, and body. The achievable compression can enable
GPU programming to accelerate the integral equation’s solution in the case of
clinically relevant voxel resolutions, for simulation setups that involve both
local and remote interactions between conductive surfaces and tissue. Since
the compressed coupling matrices can be stored and re-used without repeating
the assembly process, the proposed method could facilitate iterative coil
optimization pipelines that include the scanner bore and an RF shield.
Finally, the idea of domain decomposition in the hybrid-VSIE method could be
easily extended to wire integral equations methods [68] to facilitate
simulations in the case of coils made with wires rather than flat conductors
[5].

## References

  * [1] T. S. Ibrahim, A. Kangarlu, and D. W. Chakeress, “Design and performance issues of RF coils utilized in ultra high field MRI: experimental and numerical evaluations,” _IEEE Transactions on Biomedical Engineering_ , vol. 52, no. 7, pp. 1278–1284, 2005.
  * [2] J. H. Duyn, “The future of ultra-high field MRI and fMRI for study of the human brain,” _Neuroimage_ , vol. 62, no. 2, pp. 1241–1248, 2012.
  * [3] H. Fujita _et al._ , “RF surface receive array coils: the art of an LC circuit,” _Journal of magnetic resonance imaging_ , vol. 38, no. 1, pp. 12–25, 2013.
  * [4] D. Hernandez and K.-N. Kim, “A review on the RF coil designs and trends for ultra high field magnetic resonance imaging,” _Investigative Magnetic Resonance Imaging_ , vol. 24, no. 3, pp. 95–122, 2020.
  * [5] N. Tavaf _et al._ , “A self-decoupled 32-channel receive array for human-brain MRI at 10.5 T,” _Magnetic resonance in medicine_ , 2021.
  * [6] J. Jin and J. Chen, “On the SAR and field inhomogeneity of birdcage coils loaded with the human head,” _Magnetic resonance in medicine_ , vol. 38, no. 6, pp. 953–963, 1997.
  * [7] R. Lattanzi _et al._ , “Electrodynamic constraints on homogeneity and radiofrequency power deposition in multiple coil excitations,” _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ , vol. 61, no. 2, pp. 315–334, 2009.
  * [8] X. Zhang _et al._ , “From complex B1 mapping to local SAR estimation for human brain MR imaging using multi-channel transceiver coil at 7T,” _IEEE transactions on medical imaging_ , vol. 32, no. 6, pp. 1058–1067, 2013\.
  * [9] O. Kraff _et al._ , “MRI at 7 Tesla and above: demonstrated and potential capabilities,” _Journal of Magnetic Resonance Imaging_ , vol. 41, no. 1, pp. 13–33, 2015.
  * [10] M. A. Ertürk _et al._ , “Toward imaging the body at 10.5 tesla,” _Magnetic resonance in medicine_ , vol. 77, no. 1, pp. 434–443, 2017.
  * [11] M. E. Ladd _et al._ , “Pros and cons of ultra-high-field MRI/MRS for human application,” _Progress in nuclear magnetic resonance spectroscopy_ , vol. 109, pp. 1–50, 2018.
  * [12] B. Zhang _et al._ , “Effect of radiofrequency shield diameter on signal-to-noise ratio at ultra-high field MRI,” _Magnetic resonance in medicine_ , vol. 85, no. 6, pp. 3522–3530, 2021.
  * [13] J. F. Villena _et al._ , “MARIE: A MATLAB-based open dource MRI electromagnetic analysis software,” _Biological, Medical Devices, and Systems_ , p. 3, 2016.
  * [14] A. G. Polimeridis _et al._ , “Stable FFT-JVIE solvers for fast analysis of highly inhomogeneous dielectric objects,” _Journal of Computational Physics_ , vol. 269, pp. 280–296, 2014.
  * [15] J. F. Villena _et al._ , “Fast electromagnetic analysis of MRI transmit RF coils based on accelerated integral equation methods,” _IEEE Transactions on Biomedical Engineering_ , vol. 63, no. 11, pp. 2250–2261, 2016.
  * [16] I. I. Giannakopoulos, M. S. Litsarev, and A. G. Polimeridis, “Memory footprint reduction for the FFT-based volume integral equation method via tensor decompositions,” _IEEE Transactions on Antennas and Propagation_ , vol. 67, no. 12, pp. 7476–7486, 2019.
  * [17] A. A. Tambova _et al._ , “On the generalization of DIRECTFN for singular integrals over quadrilateral patches,” _IEEE Transactions on Antennas and Propagation_ , vol. 66, no. 1, pp. 304–314, 2017.
  * [18] I. P. Georgakis _et al._ , “A fast volume integral equation solver with linear basis functions for the accurate computation of EM fields in MRI,” _IEEE Transactions on Antennas and Propagation_ , vol. 69, no. 7, pp. 4020–4032, 2020.
  * [19] W. Chai and D. Jiao, “Theoretical study on the rank of integral operators for broadband electromagnetic modeling from static to electrodynamic frequencies,” _IEEE Transactions on Components, Packaging and Manufacturing Technology_ , vol. 3, no. 12, pp. 2113–2126, 2013.
  * [20] G. D. Guryev _et al._ , “Fast field analysis for complex coils and metal implants in MARIE 2.0.” in _Proc. ISMRM_ , 2019, p. 1035.
  * [21] E. Milshteyn _et al._ , “Individualized SAR calculations using computer vision-based MR segmentation and a fast electromagnetic solver,” _Magnetic Resonance in Medicine_ , vol. 85, no. 1, pp. 429–443, 2021.
  * [22] I. I. Giannakopoulos _et al._ , “Compression of Volume-Surface Integral Equation Matrices via Tucker Decomposition for Magnetic Resonance Applications,” _IEEE Transactions on Antennas and Propagation_ , vol. 70, no. 1, pp. 459–471, 2022.
  * [23] M. Bebendorf and S. Rjasanow, “Adaptive low-rank approximation of collocation matrices,” _Computing_ , vol. 70, no. 1, pp. 1–24, 2003.
  * [24] I. Oseledets and E. Tyrtyshnikov, “TT-cross approximation for multidimensional arrays,” _Linear Algebra and its Applications_ , vol. 432, no. 1, pp. 70–88, 2010.
  * [25] J. C. Maxwell, “VIII. A dynamical theory of the electromagnetic field,” _Philosophical transactions of the Royal Society of London_ , no. 155, pp. 459–512, 1865.
  * [26] W. C. Chew, M. S. Tong, and B. Hu, “Integral equation methods for electromagnetic and elastic waves,” _Synthesis Lectures on Computational Electromagnetics_ , vol. 3, no. 1, pp. 1–241, 2008.
  * [27] P. Yla-Oijala _et al._ , “Surface and volume integral equation methods for time-harmonic solutions of Maxwell’s equations,” _Progress In Electromagnetics Research_ , vol. 149, pp. 15–44, 2014.
  * [28] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” _IEEE Transactions on antennas and propagation_ , vol. 30, no. 3, pp. 409–418, 1982.
  * [29] D. Jiao and J.-M. Jin, “Fast frequency-sweep analysis of RF coils for MRI,” _IEEE transactions on biomedical engineering_ , vol. 46, no. 11, pp. 1387–1390, 1999.
  * [30] Y. H. Lo _et al._ , “Finite-width gap excitation and impedance models,” in _2011 IEEE International Symposium on Antennas and Propagation (APSURSI)_. IEEE, 2011, pp. 1297–1300.
  * [31] C.-T. Tai, _Dyadic Green functions in electromagnetic theory_. Institute of Electrical & Electronics Engineers (IEEE), 1994.
  * [32] R. F. Harrington, _Field computation by moment methods_. Wiley-IEEE Press, 1993.
  * [33] J.-G. Wei, Z. Peng, and J.-F. Lee, “A fast direct matrix solver for surface integral equation methods for electromagnetic wave scattering from non-penetrable targets,” _Radio Science_ , vol. 47, no. 05, pp. 1–9, 2012\.
  * [34] S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, “A theory of pseudoskeleton approximations,” _Linear algebra and its applications_ , vol. 261, no. 1-3, pp. 1–21, 1997.
  * [35] S. A. Goreinov, N. L. Zamarashkin, and E. E. Tyrtyshnikov, “Pseudo-skeleton approximations of matrices,” in _Doklady Akademii Nauk_ , vol. 343, no. 2. Russian Academy of Sciences, 1995, pp. 151–152.
  * [36] S. A. Goreinov and E. E. Tyrtyshnikov, “The maximal-volume concept in approximation by low-rank matrices,” _Contemporary Mathematics_ , vol. 280, pp. 47–52, 2001.
  * [37] S. Kurz, O. Rain, and S. Rjasanow, “The adaptive cross-approximation technique for the 3D boundary-element method,” _IEEE transactions on Magnetics_ , vol. 38, no. 2, pp. 421–424, 2002.
  * [38] K. Zhao, M. N. Vouvakis, and J.-F. Lee, “The adaptive cross approximation algorithm for accelerated method of moments computations of EMC problems,” _IEEE transactions on electromagnetic compatibility_ , vol. 47, no. 4, pp. 763–773, 2005.
  * [39] M. F. Catedra, E. Gago, and L. Nuno, “A numerical scheme to obtain the RCS of three-dimensional bodies of resonant size using the conjugate gradient method and the fast Fourier transform,” _IEEE transactions on antennas and propagation_ , vol. 37, no. 5, pp. 528–537, 1989.
  * [40] P. Zwamborn and P. M. Van Den Berg, “The three dimensional weak form of the conjugate gradient FFT method for solving scattering problems,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 40, no. 9, pp. 1757–1766, 1992.
  * [41] J. Vaughan _et al._ , “Efficient high-frequency body coil for high-field MRI,” _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ , vol. 52, no. 4, pp. 851–859, 2004.
  * [42] J. R. Corea _et al._ , “Screen-printed flexible MRI receive coils,” _Nature communications_ , vol. 7, no. 1, pp. 1–7, 2016.
  * [43] I. I. Giannakopoulos _et al._ , “A tensor train compression scheme for remote volume-surface integral equation interactions,” _ACES_ , 2021.
  * [44] J. R. Phillips and J. K. White, “A precorrected-FFT method for electrostatic analysis of complicated 3-D structures,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 16, no. 10, pp. 1059–1072, 1997.
  * [45] J. D. Carroll and J.-J. Chang, “Analysis of individual differences in multidimensional scaling via an N-way generalization of “Eckart-Young” decomposition,” _Psychometrika_ , vol. 35, no. 3, pp. 283–319, 1970.
  * [46] R. A. Harshman _et al._ , “Foundations of the PARAFAC procedure: Models and conditions for an" explanatory" multimodal factor analysis,” _UCLA Working Papers in Phonetics_ , 1970.
  * [47] L. R. Tucker, “Some mathematical notes on three-mode factor analysis,” _Psychometrika_ , vol. 31, no. 3, pp. 279–311, 1966.
  * [48] I. V. Oseledets, D. Savostianov, and E. E. Tyrtyshnikov, “Tucker dimensionality reduction of three-dimensional arrays in linear time,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 30, no. 3, pp. 939–956, 2008.
  * [49] I. V. Oseledets, “Tensor-train decomposition,” _SIAM Journal on Scientific Computing_ , vol. 33, no. 5, pp. 2295–2317, 2011.
  * [50] J. Ballani, L. Grasedyck, and M. Kluge, “Black box approximation of tensors in hierarchical Tucker format,” _Linear algebra and its applications_ , vol. 438, no. 2, pp. 639–657, 2013.
  * [51] Q. Zhao _et al._ , “Tensor ring decomposition,” _arXiv preprint arXiv:1606.05535_ , 2016.
  * [52] Z. Chen _et al._ , “Sparsity-aware precorrected tensor train algorithm for fast solution of 2-D scattering problems and current flow modeling on unstructured meshes,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 67, no. 12, pp. 4833–4847, 2019.
  * [53] M. Wang _et al._ , “VoxCap: FFT-accelerated and Tucker-enhanced capacitance extraction simulator for voxelized structures,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 68, no. 12, pp. 5154–5168, 2020.
  * [54] C. Qian, M. Wang, and A. C. Yucel, “Compression of far-fields in the fast multipole method via tucker decomposition,” _IEEE Transactions on Antennas and Propagation_ , 2021.
  * [55] B. Yaman _et al._ , “Low-rank tensor models for improved multidimensional MRI: application to dynamic cardiac $t_{1}$ mapping,” _IEEE transactions on computational imaging_ , vol. 6, pp. 194–207, 2019.
  * [56] Y. Yu _et al._ , “Multidimensional compressed sensing MRI using tensor decomposition-based sparsifying transform,” _PloS one_ , vol. 9, no. 6, p. e98441, 2014.
  * [57] D. Savostyanov and I. Oseledets, “Fast adaptive interpolation of multi-dimensional arrays in tensor train format,” in _The 2011 International Workshop on Multidimensional (nD) Systems_. IEEE, 2011, pp. 1–8.
  * [58] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” _SIAM journal on Matrix Analysis and Applications_ , vol. 21, no. 4, pp. 1253–1278, 2000.
  * [59] I. I. Giannakopoulos, M. S. Litsarev, and A. G. Polimeridis, “3D cross-Tucker approximation in FFT-based volume integral equation methods,” in _2018 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting_. IEEE, 2018, pp. 2507–2508.
  * [60] K. Zhao, V. Rawat, and J.-F. Lee, “A domain decomposition method for electromagnetic radiation and scattering analysis of multi-target problems,” _IEEE Transactions on Antennas and Propagation_ , vol. 56, no. 8, pp. 2211–2221, 2008.
  * [61] Y. Saad and M. H. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” _SIAM Journal on scientific and statistical computing_ , vol. 7, no. 3, pp. 856–869, 1986.
  * [62] A. Christ _et al._ , “The Virtual Family-development of surface-based anatomical models of two adults and two children for dosimetric simulations,” _Physics in Medicine & Biology_, vol. 55, no. 2, p. N23, 2009\.
  * [63] I. Giannakopoulos _et al._ , “Magnetic-resonance-based electrical property mapping using Global Maxwell Tomography with an 8-channel head coil at 7 Tesla: a simulation study,” _IEEE Transactions on Biomedical Engineering_ , vol. 68, no. 1, pp. 236–246, 2021.
  * [64] R. Lattanzi _et al._ , “Performance evaluation of a 32-element head array with respect to the ultimate intrinsic SNR,” _NMR in Biomedicine: An International Journal Devoted to the Development and Application of Magnetic Resonance In vivo_ , vol. 23, no. 2, pp. 142–151, 2010.
  * [65] J. E. C. Serrallés _et al._ , “Parametric Coil Optimization via Global Optimization,” in _Proc. ISMRM_ , 2021, p. 1399.
  * [66] F. P. Andriulli _et al._ , “A multiplicative Calderon preconditioner for the electric field integral equation,” _IEEE Transactions on Antennas and Propagation_ , vol. 56, no. 8, pp. 2398–2412, 2008.
  * [67] J. E. Serrallés _et al._ , “Noninvasive estimation of electrical properties from magnetic resonance measurements via Global Maxwell Tomography and match regularization,” _IEEE Transactions on Biomedical Engineering_ , vol. 67, no. 1, pp. 3–15, 2019.
  * [68] M. Y. Monin _et al._ , “A hybrid volume-surface-wire integral equation for the anisotropic forward problem in electroencephalography,” _IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology_ , vol. 4, no. 4, pp. 286–293, 2020.


# Zeros of Ramanujan-type Polynomials

Bibekananda Maji  Bibekananda Maji
Department of Mathematics
Indian Institute of Technology Indore
Indore, Simrol, Madhya Pradesh 453552, India<EMAIL_ADDRESS><EMAIL_ADDRESS>and  Tithi Sarkar  Tithi Sarkar
Department of Mathematics
Indian Institute of Technology Indore
Indore, Simrol, Madhya Pradesh 453552, India<EMAIL_ADDRESS>
###### Abstract.

Ramanujan’s notebooks contain many elegant identities and one of the
celebrated identities is a formula for $\zeta(2k+1)$. In 1972, Grosswald gave
an extension of the Ramanujan’s formula for $\zeta(2k+1)$, which contains a
polynomial of degree $2k+2$. This polynomial is now well-known as the
Ramanujan polynomial $R_{2k+1}(z)$, first studied by Gun, Murty, and Rath.
Around the same time, Murty, Smith and Wang proved that all the non-real zeros
of $R_{2k+1}(z)$ lie on the unit circle. Recently, Chourasiya, Jamal, and the
first author found a new polynomial while obtaining a Ramanujan-type formula
for Dirichlet $L$-functions and named it as Ramanujan-type polynomial
$R_{2k+1,p}(z)$. In the same paper, they conjectured that all the non-real
zeros of $R_{2k+1,p}(z)$ lie on the circle $|z|=1/p$. The main goal of this
paper is to present a proof of this conjecture.

2010 Mathematics Subject Classification. Primary 26C10; Secondary 11B68.
Keywords and phrases. Bernoulli numbers, odd zeta values, Ramanujan
polynomials, roots on the unit circle, reciprocal polynomials

##  1\. Introduction

In the Second Notebook [33, p. 173, Ch. 14, Entry 21(i)] as well as in the
Lost Notebook [34, p. 319, Entry (28)], Ramanujan mentioned the following
remarkable identity for odd zeta values. Considering $\alpha>0$ and
$\beta=\frac{\pi^{2}}{\alpha}$ and non-zero integer $k$, he proved that

$\displaystyle
G_{k}(\alpha)=(-1)^{k}G_{k}(\beta)-2^{2k}\sum_{j=0}^{k+1}(-1)^{j-1}\frac{B_{2j}}{(2j)!}\frac{B_{2k+2-2j}}{(2k+2-2j)!}\alpha^{k+1-j}\beta^{j},$
(1.1)

where

$G_{k}(x)=x^{-k}\left\\{\frac{1}{2}\zeta(2k+1)+\sum_{n=0}^{\infty}\frac{n^{-2k-1}}{e^{2xn}-1}\right\\}.$

In 1972, Emil Grosswald [17] gave a notable extension of Ramanujan’s formula
(1.1). Let us recall the divisor function $\sigma_{z}(m)=\sum_{d|m}d^{z}$ for
any $z\in\mathbb{C}$. For $\xi\in\mathbb{H}$, $k\in\mathbb{Z}-\\{0\\}$, we
define

$\mathfrak{F}_{k}(\xi):=\sum_{n=1}^{\infty}\sigma_{-k}(n)e^{2{\pi}in\xi}.$

Then we have

$\displaystyle\mathfrak{F}_{2k+1}(\xi)-\xi^{2k}\mathfrak{F}_{2k+1}\left(-\frac{1}{\xi}\right)$
$\displaystyle=\frac{1}{2}\zeta(2k+1)(\xi^{2k}-1)$
$\displaystyle+\frac{(2{\pi}i)^{2k+1}}{2\xi}\sum_{j=0}^{k+1}\xi^{2k+2-2j}\frac{B_{2j}}{(2j)!}\frac{B_{2k+2-2j}}{(2k+2-2j)!}.$
(1.2)

Substituting $\xi=i\beta/\pi,$ $\alpha\beta={\pi}^{2},$ with $\alpha,\beta>0,$
the identity (1) becomes same as the Ramanujan’s identity (1.1) for
$\zeta(2k+1).$ Later, Gun, Murty, and Rath [18] studied the finite sum
involving Bernoulli numbers present on the right hand side of (1) and named it
as the Ramanujan polynomial. Thus, for $k\in\mathbb{N}$ and $z\in\mathbb{C}$,
Ramanujan polynomials are defined as follows:

$R_{2k+1}(z):=\sum_{j=0}^{k+1}z^{2k+2-2j}\frac{B_{2j}}{(2j)!}\frac{B_{2k+2-2j}}{(2k+2-2j)!}.$
(1.3)

One can easily verify that Ramanujan polynomials are reciprocal or self-
inverse. In 2011, Murty, Smith, and Wang [20] proved that all the non-real
roots of $R_{2k+1}(z)$ lie on the unit circle $|z|=1$. In 2013, Lalin and
Rogers [22] observed similar behaviour for a few variant of Ramanujan
polynomials involving Bernoulli numbers and Euler numbers. One of the
polynomials they [22, Theorem 2.2] have studied is the following:

$\displaystyle
R_{2k+1,2}(z)=\sum_{j=1}^{k}(2^{2j}-1)(2^{2k+2-2j}-1)\frac{B_{2j}B_{2k+2-2j}}{{(2j)!(2k+2-2j)!}}z^{2k+2-2j}.$
(1.4)

Interestingly, an existence of the above polynomial can be seen from a
Ramanujan-type formula studied by Berndt [4, Theorem 4.7] and Malurkar [19].
Recently, Chourasiya, Jamal, and the first author [11] found the following
Grosswald-type formula for the Dirichlet $L$-function $L(s,\chi)$. Let $p$ be
a prime number and $\chi_{1}$ be the principal Dirichlet character modulo $p$.
For any $k\in\mathbb{N}$ and $z\in\mathbb{H}$, we define

$\displaystyle\mathfrak{F}_{k,\chi_{1}}(z):=\sum_{n=1}^{\infty}a_{n}\sigma_{-k,\chi_{1}}(n)\exp(2\pi
inz),\quad\text{with}\quad a_{n}=\begin{cases}1,&\text{if}\,\,\gcd(n,p)=1,\\\
1-p,&\text{if}\,\,\gcd(n,p)=p.\end{cases}$

Then we have

$\displaystyle\left(pz\right)^{2k}$
$\displaystyle\mathfrak{F}_{2k+1,\chi_{1}}\left(-\frac{1}{p^{2}z}\right)-\mathfrak{F}_{2k+1,\chi_{1}}(z)=\frac{p-1}{2}L(2k+1,\chi_{1})\left\\{\left(pz\right)^{2k}-1\right\\}$
$\displaystyle+$ $\displaystyle\frac{(2\pi
i)^{2k+1}}{2z\,p^{2k+2}}\sum_{j=1}^{k}{(p^{2j}-1)(p^{2k+2-2j}-1)\frac{B_{2j}}{(2j)!}\frac{B_{2k+2-2j}}{(2k+2-2j)!}}(pz)^{2k+2-2j}.$
(1.5)

The above identity motivated them to define Ramanujan-type polynomials as

$R_{2k+1,p}(z):=\sum_{j=1}^{k}(p^{2j}-1)(p^{2k+2-2j}-1)\frac{B_{2j}B_{2k+2-2j}}{{(2j)!(2k+2-2j)!}}(pz)^{2k+2-2j},$
(1.6)

where $k\in\mathbb{N}$ and $p$ is a prime number. In the same paper, they [11,
Conjecture 2.9] have stated the below conjecture about the location of non-
real roots of $R_{2k+1,p}(z)$.

######  Conjecture 1.1.

Let $k\in\mathbb{N}$ and $p$ be any prime number, the only real root of
Ramanujan-type polynomial $R_{2k+1,p}(z)$ is $z=0$ of multiplicity $2$.
Moreover, the non-real roots are simple and lie on the circle
$|z|=\frac{1}{p}$.

At the end of their paper [11, p. 748], they also mentioned that the above
conjecture might be true for any nature number $n\geq 2$ instead of any prime
number $p$. The main goal of this article is to prove this conjecture. Our
main result is the following.

######  Theorem 1.2.

Let $k$ and $n\geq 2$ be two fixed natural numbers. Then the following
Ramanujan-type polynomial

$R_{2k+1,n}(z):=\sum_{j=1}^{k}(n^{2j}-1)(n^{2k+2-2j}-1)\frac{B_{2j}B_{2k+2-2j}}{{(2j)!(2k+2-2j)!}}(nz)^{2k+2-2j},$
(1.7)

has only one real root at $z=0$ of multiplicity $2$. Again, the non-real roots
are simple and lie on the circle $|z|=\frac{1}{n}$.

In the next section, we collect a few well-known results about the location of
zeros of reciprocal polynomials.

##  2\. Known results on reciprocal polynomials

Let $p(x)\in\mathbb{R}[x]$ be a non-zero polynomial of degree $n$ with

$p(x)=a_{0}+a_{1}x+a_{2}x^{2}+\ldots+a_{n}x^{n}.$

We say that $p(x)$ is a reciprocal or self-inverse polynomial if and only if

$x^{n}p\left(\frac{1}{x}\right)=a_{n}+a_{n-1}x+a_{n-2}x^{2}+\ldots+a_{0}x^{n}=p(x).$

That implies that the coefficients of real reciprocal polynomials satisfy

$a_{k}=a_{n-k},~{}\forall~{}k=0,1,\cdots,n.$

######  Theorem 2.1 (Lakatos [21]).

If a reciprocal polynomial

$p(x)=\sum_{i=0}^{k}a_{i}x^{i}$

of degree $k$ with real coefficients satisfy the condition

$|a_{k}|\geq\sum_{i=0}^{k}|a_{i}-a_{k}|,$ (2.1)

then the polynomial $p(x)$ certainly bears all its zeros on the unit circle.
If the inequality (2.1) is strict, then the multiplicity of all the zeros of
$p(x)$ is one.

Later, the above theorem has been modified by many mathematicians, thus the
condition (2.1) for reciprocal polynomials has many improvements. Schinzel
[23] gave one generalization of Theorem 2.1 for any reciprocal polynomial over
$\mathbb{C}$.

######  Theorem 2.2 (Schinzel [23]).

Let $p(z)$ be a reciprocal polynomial, say,

$p(z)=\sum_{j=0}^{k}a_{j}z^{j}$

of degree $k$ with $a_{j}\in\mathbb{C},\forall j$. If the coefficients of
$p(z)$ satisfy the condition

$|a_{k}|\geq\inf_{\begin{subarray}{c}c,d\in\mathbb{C}\\\
{|d|=1}\end{subarray}}\sum_{j=0}^{k}|ca_{j}-d^{k-j}a_{k}|,$ (2.2)

then all the zeros of $p(z)$ lie on $|z|=1$. If the inequality (2.2) is
strict, then all the zeros of $p(x)$ must be simple.

Next, we collect a few results related to Bernoulli numbers. The following
bounds hold for Bernoulli numbers.

######  Lemma 2.3.

[1, p. 805] For $n\in\mathbb{N}$, we have

$\displaystyle\frac{2(2n)!}{(2\pi)^{2n}}<|B_{2n}|<\frac{2(2n)!}{(2\pi)^{2n}(1-2^{1-2n})}.$
(2.3)

We are now going to mention one useful identity about the convolution of
Bernoulli polynomials.

######  Lemma 2.4.

[12, p. 31, eq. (3.2)]. For $m\in\mathbb{N}$, we have

$\displaystyle\sum_{j=0}^{m}{{m}\choose{j}}B_{j}(a)B_{m-j}(b)=m(a+b+1)B_{m-1}(a+b)-(m-1)B_{m}(a+b).$

In the next section, we present a proof of our main result.

##  3\. Proof of Theorem 1.2

Before proving our main theorem, we need following lemmas.

######  Lemma 3.1.

Let $k\geq 1$ and $n\geq 2$ be two fixed natural numbers. Let $f_{n}(x)$ be a
function $f_{n}:[0,~{}2k-2]\rightarrow\mathbb{R}$ defined as

$\displaystyle f_{n}(x):=(n^{x+2}-1)(n^{2k-x}-1).$ (3.1)

The function $f_{n}(x)$ has global maxima at $x=k-1$ and its minima at
$x=0,2k-2$, with the same minimum value $(n^{2}-1)(n^{2k}-1)$ at both of the
end points.

  * Proof.

For fixed two natural numbers $n\geq 2$ and $k\geq 1$, it is given that

$f_{n}(x):=(n^{x+2}-1)(n^{2k-x}-1).$

Differentiating with respect to $x$, one can see that

$f_{n}^{\prime}(x)=-n^{x+2}\log(n)+n^{2k-x}\log(n)\Rightarrow
f_{n}^{\prime\prime}(x)=-\left(n^{x+2}+n^{2k-x}\right)\log^{2}(n).$

Further, we can easily check that $f_{n}^{\prime\prime}(x)$ is always negative
in the given domain, so all the extrema points obtained by equating
$f^{\prime}_{n}(x)$ to zero will correspond to the local maximum values. Now
we find extrema points of the function $f_{n}(x)$ by equating
$f^{\prime}_{n}(x)$ to zero, that is,

$f_{n}^{\prime}(x)=0\Rightarrow n^{x+2}=n^{2k-x}\Rightarrow x=k-1.$

Thus, $x=k-1$ is the only local maxima point in the interval $[0,2k-2]$. Hence
it must global maxima as well. Since the function is continuously
differentiable, so we conclude that $x=2k-2~{}\text{and}~{}x=0$ are the only
points of minima. Moreover, we have

$f_{n}(0)=f_{n}(2k-2)=(n^{2}-1)(n^{2k}-1).$

Hence the result follows. ∎

######  Lemma 3.2.

Let $k\geq 1$ and $n\geq 2$ be two fixed natural numbers. The function
$g_{n}:[0,2k-2]\rightarrow\mathbb{R}$ defined as

$g_{n}(x):=\frac{(n^{x+2}-1)}{(2^{x+2}-1)}\frac{(n^{2k-x}-1)}{(2^{2k-x}-1)},$
(3.2)

has an upper bound $\frac{(n^{k+1}-1)^{2}}{3(2^{2k}-1)}.$

  * Proof.

We can write $g_{n}(x)$ in terms of $f_{n}(x)$ as follows:

$g_{n}(x)=\frac{f_{n}(x)}{f_{2}(x)}.$

Therefore, the ratio of maximum value of $f_{n}(x)$ and the minimum value of
$f_{2}(x)$ will provide us an upper bound for $g_{n}(x).$ From Lemma 3.1, we
know that the maxima of $f_{n}(x)$ is attained at $x=k-1$ and the
corresponding maximum value is,

$\displaystyle f_{n}(k-1)=(n^{k+1}-1)^{2}.$ (3.3)

Also the minimum value of $f_{n}(x)$ is attained at $x=0$ and in particular
for $n=2$ the minimum value will be,

$\displaystyle f_{2}(0)=(2^{2}-1)(2^{2k}-1)=3(2^{2k}-1).$ (3.4)

Hence, using (3.3) and (3.4), it follows that

$g_{n}(x)\leq\frac{(n^{k+1}-1)^{2}}{3(2^{2k}-1)}.$

This completes the proof. ∎

Next, we state a result which will be crucial for the proof of our main
theorem.

######  Proposition 3.3.

Let $R_{2k+1,n}(z)$ be the Ramanujan-type polynomial defined as in (1.7) and
$A_{j}$ be the $(2k-2-2j)^{th}$ coefficient of
$z^{-2}R_{2k+1,n}\left(\frac{z}{n}\right).$ Then, we have

$|A_{k-1}|\geq\inf_{\begin{subarray}{c}c,d\in\mathbb{C}\\\
{|d|=1}\end{subarray}}\sum_{j=0}^{k-1}|cA_{j}-dA_{k-1}|.$

  * Proof.

By appropriate shifting, from the definition (1.7) of $R_{2k+1,n}(z)$, we can
write

$R_{2k+1,n}(z)=\sum_{j=0}^{k-1}(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{{(2j+2)!(2k-2j)!}}(nz)^{2k-2j}.$
(3.5)

Replacing $z$ by $\frac{z}{n}$, we have

$\displaystyle R_{2k+1,n}\left(\frac{z}{n}\right)$
$\displaystyle=z^{2}\sum_{j=0}^{k-1}(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{{(2j+2)!(2k-2j)!}}z^{2k-2-2j}$
$\displaystyle=z^{2}H(z),$ (3.6)

where $H(z)$ is defined as

$\displaystyle
H(z):=\sum_{j=0}^{k-1}(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{{(2j+2)!(2k-2j)!}}z^{2k-2-2j}.$
(3.7)

It is given that the coefficient of $z^{2k-2-2j}$ in (3.7) as

$\displaystyle
A_{j}=(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{{(2j+2)!(2k-2j)!}}.$
(3.8)

Now we show that the polynomial $H(z)$ satisfies Theorem 2.2, that is, the
coefficients of $H(z)$ satisfy the following inequality:

$\displaystyle|A_{k-1}|\geq\inf_{\begin{subarray}{c}c,d\in\mathbb{C}\\\
{|d|=1}\end{subarray}}\sum_{j=0}^{k-1}|cA_{j}-d^{k-1-j}A_{k-1}|.$ (3.9)

In particular, considering $d=1$, we only need to show

$|A_{k-1}|\geq\sum_{j=0}^{k-1}|cA_{j}-A_{k-1}|.$ (3.10)

Using the lower bound for Bernoulli numbers i.e., Lemma 2.3 in (3.7), we have

$\displaystyle|A_{j}|$
$\displaystyle=(n^{2j+2}-1)(n^{2k-2j}-1)\frac{|B_{2j+2}||B_{2k-2j}|}{{(2j+2)!(2k-2j)!}}$
$\displaystyle>\frac{4(n^{2j+2}-1)(n^{2k-2j}-1)}{(2\pi)^{2k+2}}.$ (3.11)

Similarly, using the upper bound for $B_{2k}$, we get

$\displaystyle|A_{k-1}|$
$\displaystyle=(n^{2k}-1)(n^{2}-1)\frac{|B_{2k}||B_{2}|}{(2k)!2!}$
$\displaystyle<(n^{2k}-1)(n^{2}-1)\frac{2.(2k)!}{(2\pi)^{2k}(1-2^{1-2k})(2k)!}\frac{2.2!}{(2\pi)^{2}(1-2^{1-2})2!}$
$\displaystyle=\frac{8(n^{2k}-1)(n^{2}-1)}{(2\pi)^{2k+2}(1-2^{1-2k})}.$ (3.12)

Note that the sign of $B_{2n}$ is $(-1)^{n+1}$, which implies that the sign of
$A_{j}$ will be $(-1)^{k+1},\forall~{}j$. Therefore,

$|cA_{j}-A_{k-1}|=(-1)^{k+1}(cA_{j}-A_{k-1}).$

Now we shall try to find a value of positive $c$ for which

$\displaystyle|cA_{j}-A_{k-1}|=(-1)^{k+1}(cA_{j}-A_{k-1})>0$
$\displaystyle\Leftrightarrow\quad c|A_{j}|-|A_{k-1}|>0.$ (3.13)

For any $n\in\mathbb{N}$, it is well-known that

$\displaystyle B_{n}\left(\frac{1}{2}\right)=(2^{1-n}-1)B_{n},\quad
B_{n}(0)=B_{n},~{}\forall n\neq 1.$

We use these properties of Bernoulli numbers to see that

$\displaystyle B_{2j+2}\left(\frac{1}{2}\right)-B_{2j+2}(0)$
$\displaystyle=B_{2j+2}(2^{1-2j-2}-1)-B_{2j+2}=B_{2j+2}(2^{1-2j-2}-2)$
$\displaystyle\Rightarrow B_{2j+2}$
$\displaystyle=\frac{B_{2j+2}(\frac{1}{2})-B_{2j+2}(0)}{2(2^{-2j-2}-1)}.$
(3.14)

Similarly, one can show that

$B_{2k-2j}=\frac{B_{2k-2j}(\frac{1}{2})-B_{2k-2j}(0)}{2(2^{-2k+2j}-1)}.$
(3.15)

Now employing the above identities (3) and (3.15) in (3.8), we obtain

$\displaystyle A_{j}=$
$\displaystyle(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{{(2j+2)!(2k-2j)!}}$
$\displaystyle=$
$\displaystyle{{2k+2}\choose{2j+2}}(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{(2k+2)!}$
$\displaystyle=$
$\displaystyle{{2k+2}\choose{2j+2}}\frac{(n^{2j+2}-1)(n^{2k-2j}-1)}{(2k+2)!}\left(\frac{B_{2j+2}(\frac{1}{2})-B_{2j+2}(0)}{2(2^{-2j-2}-1)}\right)\left(\frac{B_{2k-2j}(\frac{1}{2})-B_{2k-2j}(0)}{2(2^{-2k+2j}-1)}\right).$
(3.16)

Recall that the sum on the right hand side of (3.10) is

$\sum_{j=0}^{k-1}|cA_{j}-A_{k-1}|=(-1)^{k+1}\sum_{j=0}^{k-1}\left(cA_{j}-A_{k-1}\right).$

Using (3), we can write

$\displaystyle(-1)^{k+1}\sum_{j=0}^{k-1}cA_{j}=\sum_{j=0}^{k-1}c|A_{j}|$
$\displaystyle=\sum_{j=0}^{k-1}c{{2k+2}\choose{2j+2}}\frac{(n^{2j+2}-1)(n^{2k-2j}-1)}{(2k+2)!}\Bigg{|}\left(\frac{B_{2j+2}(\frac{1}{2})-B_{2j+2}(0)}{2(2^{-2j-2}-1)}\right)\left(\frac{B_{2k-2j}(\frac{1}{2})-B_{2k-2j}(0)}{2(2^{-2k+2j}-1)}\right)\Bigg{|}$
$\displaystyle=\sum_{j=0}^{k-1}D(j){{2k+2}\choose{2j+2}}\Bigg{|}\left(B_{2j+2}\left(\frac{1}{2}\right)-B_{2j+2}(0)\right)\left(B_{2k-2j}\left(\frac{1}{2}\right)-B_{2k-2j}(0)\right)\Bigg{|},$
(3.17)

where

$\displaystyle D(j)$
$\displaystyle=\frac{c}{4(2k+2)!}\frac{(n^{2j+2}-1)(n^{2k-2j}-1)}{(2^{-2j-2}-1)(2^{-2k+2j}-1)}=\frac{2^{2k}c}{(2k+2)!}\frac{(n^{2j+2}-1)}{(2^{2j+2}-1)}\frac{(n^{2k-2j}-1)}{(2^{2k-2j}-1)}.$
(3.18)

Using Lemma 3.2, one can find the following upper bound for $D(j)$:

$\displaystyle
D(j)=\frac{2^{2k}c}{(2k+2)!}g_{n}(2j)<\frac{2^{2k}c}{(2k+2)!}\frac{(n^{k+1}-1)^{2}}{3(2^{2k}-1)}.$
(3.19)

Utilize (3.19) in (3) to see that

$\displaystyle\sum_{j=0}^{k-1}c|A_{j}|$
$\displaystyle<\frac{2^{2k}c}{(2k+2)!}\frac{(n^{k+1}-1)^{2}}{3(2^{2k}-1)}\sum_{j=0}^{k-1}{{2k+2}\choose{2j+2}}$
$\displaystyle\Bigg{|}\left(B_{2j+2}\left(\frac{1}{2}\right)-B_{2j+2}(0)\right)\left(B_{2k-2j}\left(\frac{1}{2}\right)-B_{2k-2j}(0)\right)\Bigg{|}.$
(3.20)

Now we shall use Lemma 2.4 to calculate the value of the sum given in (3).
That is,

$\displaystyle\sum_{j=0}^{k-1}{{2k+2}\choose{2j+2}}\Bigg{|}\left(B_{2j+2}\left(\frac{1}{2}\right)-B_{2j+2}(0)\right)\left(B_{2k-2j}\left(\frac{1}{2}\right)-B_{2k-2j}(0)\right)\Bigg{|}$
$\displaystyle=(-1)^{k-1}\sum_{j=0}^{k-1}{{2k+2}\choose{2j+2}}\Bigg{[}B_{2j+2}\left(\frac{1}{2}\right)B_{2k-2j}\left(\frac{1}{2}\right)-B_{2j+2}\left(\frac{1}{2}\right)B_{2k-2j}(0)$
$\displaystyle-
B_{2j+2}(0)B_{2k-2j}\left(\frac{1}{2}\right)+B_{2j+2}(0)B_{2k-2j}(0)\Bigg{]}$
$\displaystyle=(-1)^{k}4(2k+1)(1-2^{-2k-2})B_{2k+2}.$ (3.21)

Employing (3) in (3), we get

$\sum_{j=0}^{k-1}c|A_{j}|<\frac{2^{2k+2}c}{(2k+2)!}\frac{(n^{k+1}-1)^{2}}{3(2^{2k}-1)}(2k+1)(1-2^{-2k-2})|B_{2k+2}|.$
(3.22)

Our main aim is to find a positive $c$ such that (3.10) holds. From (3.10), we
have

$\displaystyle|A_{k-1}|>\sum_{j=0}^{k-1}|cA_{j}-A_{k-1}|$
$\displaystyle\Leftrightarrow(-1)^{k+1}A_{k-1}>(-1)^{k+1}\sum_{j=0}^{k-1}(cA_{j}-A_{k-1})$
$\displaystyle\Leftrightarrow(1+k)(-1)^{k+1}A_{k-1}>(-1)^{k+1}\sum_{j=0}^{k-1}cA_{j}$
$\displaystyle\Leftrightarrow(1+k)|A_{k-1}|>\sum_{j=0}^{k-1}c|A_{j}|.$ (3.23)

To prove (3), we shall show that the left hand side expression of (3) is in
fact bigger than an upper bound of the sum $\sum_{j=0}^{k-1}c|A_{j}|$, which
we have already found in (3.22). So we need to find a positive $c$ such that
the following inequality holds:

$\displaystyle(1+k)|A_{k-1}|>\frac{c}{(2k+2)!}\frac{(n^{k+1}-1)^{2}}{3(2^{2k}-1)}(2k+1)(2^{2k+2}-1)|B_{2k+2}|.$
(3.24)

This implies that $c$ should satisfy the below inequality:

$\displaystyle\frac{3(1+k)(2^{2k}-1)(2k+2)!|A_{k-1}|}{(2k+1)(n^{k+1}-1)^{2}(2^{2k+2}-1)|B_{2k+2}|}>c.$
(3.25)

Moreover, using the lower bound (3.11) of $|A_{k-1}|$ and upper bound (2.3)
for $|B_{2k+2}|$, one can show that

$\displaystyle\frac{3(1+k)(2^{2k}-1)(2k+2)!|A_{k-1}|}{(2k+1)(n^{k+1}-1)^{2}(2^{2k+2}-1)|B_{2k+2}|}>c_{n,k}$

where

$\displaystyle
c_{n,k}:=\frac{(1+k)(n^{2k}-1)(n^{2}-1)(2^{2k}-1)(2^{2}-1)(2^{2k+1}-1)}{2^{2k}(2k+1)(n^{k+1}-1)^{2}(2^{2k+2}-1)}.$

Note the constant $c_{n,k}$ is a positive constant. Thus, for any positive $c$
such that $c<c_{n,k}$ will satisfy the inequality (3.25) and hence the
equation (3). This proves the inequality (3.10) for the polynomial $H(z)$ and
consequently the result follows. ∎

Now we are ready for proving the main theorem.

  * Proof of Theorem 1.2.

We are interested to show that all non-real roots of $R_{2k+1,n}(z)$ lie on
the circle $z=1/n$, which is equivalent to show that all non-real roots of
$R_{2k+1,n}\left(\frac{z}{n}\right)$ lie on the unit circle $|z|=1$. From the
equation (3), we know that

$\displaystyle R_{2k+1,n}\left(\frac{z}{n}\right)$
$\displaystyle=z^{2}\sum_{j=0}^{k-1}(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{{(2j+2)!(2k-2j)!}}z^{2k-2-2j}$
$\displaystyle=z^{2}H(z),$ (3.26)

where $H(z)$ is defined as

$\displaystyle H(z)=\sum_{j=0}^{k-1}A_{j}z^{2k-2-2j},$

and

$\displaystyle
A_{j}=(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{{(2j+2)!(2k-2j)!}}.$

We can easily check that $A_{j}$ and $A_{k-1-j}$ are same, which indicates
that the polynomial $R_{2k+1,n}\left(\frac{z}{n}\right)$ is reciprocal and so
$R_{2k+1,n}(z)$. Now from ( Proof of Theorem 1.2.), we can easily see that
$R_{2k+1,n}\left(\frac{z}{n}\right)$ has degree $2k$ containing only even
powers of $z$ with $z^{2}$ being the least power of $z$. Therefore,
$R_{2k+1,n}\left(\frac{z}{n}\right)$ has a real zero at $z=0$ of multiplicity
$2$. Utilizing Proposition 3.3, we can see that $H(z)$ satisfies Theorem 2.2,
so we can conclude that $H(z)$ has all its roots on the unit circle. Moreover,
as we have proved that the inequality (3) is strict so all the zeros of $H(z)$
are simple. Thus, the polynomial $R_{2k+1,n}\left(\frac{z}{n}\right)$ has all
the zeros on the unit circle except $z=0$. Next, we show that there is no real
root except zero. Let us suppose that the polynomial
$R_{2k+1,n}\left(\frac{z}{n}\right)$ has a real root except $z=0$, then it
must be on the unit circle. Consequently, the possibility of the non-trivial
real root will be either $z=1~{}\text{or}~{}z=-1$ or both. So, putting $z=1$
in (3.7), we get

$\displaystyle H(1)=0,$ $\displaystyle\Rightarrow$
$\displaystyle\sum_{j=0}^{k-1}(n^{2j+2}-1)(n^{2k-2j}-1)\frac{B_{2j+2}B_{2k-2j}}{{(2j+2)!(2k-2j)!}}=0,$
$\displaystyle\Rightarrow$ $\displaystyle\sum_{j=0}^{k-1}A_{j}=0.$ (3.27)

We have already seen that the sign of $A_{j}$ is $(-1)^{k+1}$ for all $j$.
Note that all even Bernoulli numbers are never zero. Therefore, depending on
the parity of $k$ the sum $\sum_{j=0}^{k-1}A_{j}$ is either positive or
negative. This is a contradiction to ( Proof of Theorem 1.2.). This completes
the proof that $z=1$ is not a root of $H(z)$. Again, as $H(z)$ is an even
function, so by the same argument we can conclude that $z=-1$ is not a root of
$H(z)$ as well. Hence, all roots of $H(z)$ are complex and simple, and lie on
the unit circle. This completes the proof. ∎

##  4\. Concluding Remarks

Motivated from the Ramanujan’s formula (1.1) as well as Grosswald’s identity
(1), Gun, Murty, and Rath [18] defined Ramanujan polynomial as

$R_{2k+1}(z)=\sum_{j=0}^{k+1}z^{2k+2-2j}\frac{B_{2j}}{(2j)!}\frac{B_{2k+2-2j}}{(2k+2-2j)!}.$
(4.1)

Around the same time, Murty, Smith, and Wang [20] showed that all the non-real
roots of $R_{2k+1}(z)$ lie on the unit circle. In this paper, we studied the
following Ramanujan-type polynomial:

$R_{2k+1,n}(z):=\sum_{j=1}^{k}(n^{2j}-1)(n^{2k+2-2j}-1)\frac{B_{2j}B_{2k+2-2j}}{{(2j)!(2k+2-2j)!}}(nz)^{2k+2-2j},$
(4.2)

We showed that it has only one real root is at $z=0$ of multiplicity $2$.
Moreover, the non-real roots of $R_{2k+1,n}(z)$ are simple and all of them lie
on the circle $|z|=\frac{1}{n}$, which says that all the non-real roots of
$R_{2k+1,n}\left(\frac{z}{n}\right)$ lie on the unit circle. Ramanujan’s
formula (1.1) has been generalized by many mathematicians, over the years.
Recently, Dixit and Gupta [13, Theorem 2.1] found an interesting analogue of
Ramanujan’s formula (1.1) for squares of odd zeta values and in the process
they encountered the following finite sum involving Bernoulli numbers. For any
$k\in\mathbb{N}$,

$\displaystyle\sum_{j=0}^{k+1}(-1)^{j}\frac{B_{2j}^{2}}{{(2j)!}^{2}}\frac{B_{2k+2-2j}^{2}}{{(2k+2-2j)!}^{2}}\alpha^{2j}\beta^{2k+2-2j}.$
(4.3)

Very recently, Banerjee and Sahani [3, Theorem 1.1] generalized the work of
Dixit and Gupta and found a formula for higher power of odd zeta values. While
obtaining this generalization, they found the following finite sum involving
higher powers of Bernoulli numbers. For any two fixed natural numbers $k$ and
$\ell$,

$\displaystyle\sum_{j=0}^{k+1}(-1)^{j}\frac{B_{2j}^{\ell}}{{(2j)!}^{\ell}}\frac{B_{2k+2-2j}^{\ell}}{{(2k+2-2j)!}^{\ell}}\alpha^{\ell
j}\beta^{\ell(k+1-j)}.$ (4.4)

It is interesting to note that, in both of the equations (4.3) and (4.4),
$\alpha$ and $\beta$ are positive and satisfy $\alpha\beta=\pi^{2}$, which is
exactly the same condition also required in Ramanujan’s formula (1.1). This
observation motivated us to define a one-variable generalization of Ramanujan
polynomial. Let $k$ and $\ell$ be two fixed natural numbers, then we define

$\displaystyle
R_{2k+1}^{(\ell)}(z):=\sum_{j=0}^{k+1}(-1)^{(\ell+1)j}\frac{B_{2j}^{\ell}}{{(2j)!}^{\ell}}\frac{B_{2k+2-2j}^{\ell}}{{(2k+2-2j)!}^{\ell}}z^{2\ell
j}.$ (4.5)

For $\ell=1$, this polynomial coincides with the Ramanujan polynomial (4.1).
For simplicity, replace $z^{2\ell}$ by $Z$ to write

$\displaystyle
R_{2k+1}^{(\ell)}(Z):=\sum_{j=0}^{k+1}(-1)^{(\ell+1)j}\frac{B_{2j}^{\ell}}{{(2j)!}^{\ell}}\frac{B_{2k+2-2j}^{\ell}}{{(2k+2-2j)!}^{\ell}}Z^{j}.$
(4.6)

We can easily verify that this is a reciprocal polynomial. Based on numerical
evidences checked in Mathematica software, we conjecture the following
statement.

######  Conjecture 4.1.

All the non-real roots of $R_{2k+1}^{(\ell)}(Z)$ are simple and lie on the
unit circle $|Z|=1$.

Ramanujan’s formula for $\zeta(2k+1)$ has been generalized by many
mathematicians for different $L$-functions. For examples, readers are
encouraged to see [2, 5, 6, 7, 9, 10, 14, 15, 16] and references therein.
Interestingly, all of these generalizations consist of a finite sum involving
Bernolli numbers or Euler numbers or some special values of $L$-functions.
This suggests that we can always find a generalized Ramanujan polynomial from
any analogue of Ramanujan’s formula for $\zeta(2k+1)$ and that might be
interest of study. In this spirit, Berndt and Straub studied a character
analogue of Ramanujan polynomial and gave a more general conjecture [8,
Conjecture 7.3], which is different from our Conjecture 4.1.

##  5\. Acknowledgement

The first author wants to thank SERB for the MATRICS Grant MTR/2022/000545.
This article is based on the M.Sc. thesis work of the second author under the
guidance of the first author. Both authors sincerely thank Bhaskaracharya
Mathematics Laboratory and Brahmagupta Mathematics Library of the Department
of Mathematics at IIT Indore, supported by DST FIST Project (File No.:
SR/FST/MS-I/2018/26).

## References

  * [1] M. Abramowitz and I. A. Stegun, Handbook of mathematical functions with formulas, graphs, and mathematical tables. National Bureau of Standards Applied Mathematics Series, Washington, D.C. 1964.
  * [2] S. Banerjee, R. Gupta, and R. Kumar, _A note on odd zeta values over number field and extended Eisenstien series_ , arXiv:2304.07870.
  * [3] S. Banerjee and V. Sahani, _Transformation formulas for the higher power of odd zeta values and generalized Eisenstein series_ , arXiv:2206.13331.
  * [4] B. C. Berndt, _Analytic Eisenstein series, theta-functions, and series relations in the spirit of Ramanujan_ , J. Reine Angew. Math. 303 (1978), 332–365.
  * [5] B. C. Berndt, _On Eisenstein series with characters and the values of Dirichlet $L$-functions_, Acta Arith. 28 (1975), 299–320.
  * [6] B. C. Berndt, _Modular transformations and generalizations of several formulae of Ramanujan_ , Rocky Mountain J. Math. 7 (1977), 147–189.
  * [7] B. C. Berndt and A. Straub, _Ramanujan’s formula for $\zeta(2n+1)$_, Exploring the Riemann zeta function, Eds. H. Montgomery, A. Nikeghbali, and M. Rassias, pp. 13–34, Springer, 2017.
  * [8] B. C. Berndt and A. Straub, _On a secant Dirichlet series and Eichler integrals of Eisenstein series_ , Math. Z. 284 (2016), 827–852.
  * [9] D. M. Bradley, _Series acceleration formulas for Dirichlet series with periodic coefficients_ , Ramanujan J. 6 (2002), 331–346.
  * [10] P. Chavan, S. Chavan, C. Vignat, and T. Wakhare, _Dirichlet series under standard convolutions: Variations on Ramanujan’s identity for odd zeta values_ , Ramanujan J. 59 (2022), 1245–1285.
  * [11] S. Chourasiya, Md. K. Jamal, and B. Maji, _A new Ramanujan-type identity for $L(2k+1,\chi_{1}$)_, Ramanujan J. 60 (2023), 729–750.
  * [12] K. Dilcher, Sums of products of Bernoulli numbers J. Number Theory 60(1996), 23–41.
  * [13] A. Dixit and R. Gupta, _On squares of odd zeta values and analogues of Eisenstein series_ , Adv. Appl. Math. 110 (2019), 86–119.
  * [14] A. Dixit, R. Gupta, R. Kumar, and B. Maji, _Generalized Lambert series, Raabe’s cosine transform and a generalization of Ramanujan’s formula for $\zeta(2m+1)$_, Nagoya Math. J. 239 (2020), 232–293.
  * [15] A. Dixit and B. Maji, _Generalized Lambert series and arithmetic nature of odd zeta values_ , Proc. Roy. Soc. Edinburg Sect. A: Mathematics 150 (2020), 741–769.
  * [16] A. Gupta and B. Maji, _On Ramanujan’s formula for $\zeta(1/2)$ and $\zeta(2m+1)$_, J. Math. Anal. Appl. 507 (2022), 125738.
  * [17] E. Grosswald, _Remarks concerning the values of the Riemann zeta function at integral, odd arguments_ , J. Number Theory 4 (1972), 225–235.
  * [18] S. Gun, M. R. Murty and P. Rath, _Transcendental values of certain Eichler integrals,_ Bull. Lond. Math. Soc., 43 (2011), 939–952.
  * [19] S. L. Malurkar, _On the application of Herr Mellin’s integrals to some series,_ J. Indian Math. Soc. 16 (1925/26), 130–138.
  * [20] R. Murty, C. Smyth, and R. Wang, _Zeros of Ramanujan’s polynomials,_ J. Ramanujan Math. Soc. 26 (2011), 107–125.
  * [21] P. Lakatos,  _On zeros of reciprocal polynomials,_ C. R. Math. Rep. Acad. Sci. Canada 24 (2002), 91–96
  * [22] M. Lalin and M. Rogers _Variations of the Ramanujan polynomials and remarks on $\zeta(2j+1)/\pi^{2j+1}$_, Funct. Approx. Comment. Math. 48 (2013), 91–111.
  * [23] A. Schinzel, Self-inversive polynomials with all zeros on the unit circle, Ramanujan J. 9 (2005), 19-23.