# Answerability Fields: Answerable Location Estimation
via Diffusion Models

Daichi Azuma1, Taiki Miyanishi2, Shuhei Kurita3, Koya Sakamoto2,4 and Motoaki
Kawanabe2 1 Sony Semiconductor Solutions2 Advanced Telecommunications Research
Institute International (ATR), Kyoto, Japan3 RIKEN Center for Advanced
Intelligence Project, Tokyo, Japan4 Graduate School of Informatics, Kyoto
University, Japan

## References

  * [1] Ilge Akkaya, Marcin Andrychowicz, Maciek Chociej, Mateusz Litwin, Bob McGrew, Arthur Petron, Alex Paino, Matthias Plappert, Glenn Powell, Raphael Ribas, et al.  Solving rubik’s cube with a robot hand.  arXiv preprint arXiv:1910.07113, 2019.
  * [2] Peter Anderson, Qi Wu, Damien Teney, Jake Bruce, Mark Johnson, Niko Sünderhauf, Ian Reid, Stephen Gould, and Anton van den Hengel.  Vision-and-language navigation: Interpreting visually-grounded navigation instructions in real environments.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2018.
  * [3] Daichi Azuma, Taiki Miyanishi, Shuhei Kurita, and Motoaki Kawanabe.  Scanqa: 3d question answering for spatial scene understanding.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 19129–19139, June 2022.
  * [4] Dhruv Batra, Aaron Gokaslan, Aniruddha Kembhavi, Oleksandr Maksymets, Roozbeh Mottaghi, Manolis Savva, Alexander Toshev, and Erik Wijmans.  ObjectNav Revisited: On Evaluation of Embodied Agents Navigating to Objects.  In arXiv:2006.13171, 2020.
  * [5] Tim Brooks, Aleksander Holynski, and Alexei A. Efros.  Instructpix2pix: Learning to follow image editing instructions.  In CVPR, 2023.
  * [6] Devendra Singh Chaplot, Dhiraj Gandhi, Abhinav Gupta, and Ruslan Salakhutdinov.  Object goal navigation using goal-oriented semantic exploration.  In In Neural Information Processing Systems, 2020.
  * [7] Dave Zhenyu Chen, Angel X Chang, and Matthias Nießner.  Scanrefer: 3d object localization in rgb-d scans using natural language.  16th European Conference on Computer Vision (ECCV), 2020.
  * [8] Shoufa Chen, Peize Sun, Yibing Song, and Ping Luo.  Diffusiondet: Diffusion model for object detection.  arXiv preprint arXiv:2211.09788, 2022.
  * [9] Angela Dai, Angel X Chang, Manolis Savva, Maciej Halber, Thomas Funkhouser, and Matthias Nießner.  Scannet: Richly-annotated 3d reconstructions of indoor scenes.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2017.
  * [10] Abhishek Das, Samyak Datta, Georgia Gkioxari, Stefan Lee, Devi Parikh, and Dhruv Batra.  Embodied Question Answering.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2018.
  * [11] Abhishek Das, Georgia Gkioxari, Stefan Lee, Devi Parikh, and Dhruv Batra.  Neural Modular Control for Embodied Question Answering.  In Proceedings of the Conference on Robot Learning (CoRL), 2018.
  * [12] Samir Yitzhak Gadre, Mitchell Wortsman, Gabriel Ilharco, Ludwig Schmidt, and Shuran Song.  Clip on wheels: Zero-shot object navigation as object localization and exploration.  ArXiv, abs/2203.10421, 2022.
  * [13] Zhangxuan Gu, Haoxing Chen, Zhuoer Xu, Jun Lan, Changhua Meng, and Weiqiang Wang.  Diffusioninst: Diffusion model for instance segmentation.  arXiv preprint arXiv:2212.02773, 2022.
  * [14] Kaiming He, Georgia Gkioxari, Piotr Dollár, and Ross Girshick.  Mask r-cnn, 2018.
  * [15] Jonathan Ho and Tim Salimans.  Classifier-free diffusion guidance, 2022.
  * [16] Chenguang Huang, Oier Mees, Andy Zeng, and Wolfram Burgard.  Visual language maps for robot navigation.  In Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), London, UK, 2023.
  * [17] Krishna Murthy Jatavallabhula, Alihusein Kuwajerwala, Qiao Gu, Mohd Omama, Tao Chen, Shuang Li, Ganesh Iyer, Soroush Saryazdi, Nikhil Keetha, Ayush Tewari, Joshua B. Tenenbaum, Celso Miguel de Melo, Madhava Krishna, Liam Paull, Florian Shkurti, and Antonio Torralba.  Conceptfusion: Open-set multimodal 3d mapping.  Robotics: Science and Systems (RSS), 2023.
  * [18] Dimitrios Kanoulas Jianwei Liu, Maria Stamatopoulou.  Dipper: Diffusion-based 2d path planner applied on legged robots.  arXiv preprint arXiv:2309.14341, 2024.
  * [19] Dmitry Kalashnikov, Alex Irpan, Peter Pastor, Julian Ibarz, Alexander Herzog, Eric Jang, Deirdre Quillen, Ethan Holly, Mrinal Kalakrishnan, Vincent Vanhoucke, et al.  Scalable deep reinforcement learning for vision-based robotic manipulation.  In CoRL, pages 651–673, 2018.
  * [20] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes, 2022.
  * [21] Dongxu Li, Junnan Li, Hung Le, Guangsen Wang, Silvio Savarese, and Steven C.H. Hoi.  LAVIS: A one-stop library for language-vision intelligence.  In Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 3: System Demonstrations), pages 31–41, Toronto, Canada, July 2023. Association for Computational Linguistics.
  * [22] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization, 2019.
  * [23] Xiaojian Ma, Silong Yong, Zilong Zheng, Qing Li, Yitao Liang, Song-Chun Zhu, and Siyuan Huang.  Sqa3d: Situated question answering in 3d scenes.  arXiv preprint arXiv:2210.07474, 2022.
  * [24] Nelson Nauata, Sepidehsadat Hosseini, Kai-Hung Chang, Hang Chu, Chin-Yi Cheng, and Yasutaka Furukawa.  House-gan++: Generative adversarial layout refinement network towards intelligent computational agent for professional architects.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 13632–13641, 2021.
  * [25] Santhosh K. Ramakrishnan, Devendra Singh Chaplot, Ziad Al-Halah, Jitendra Malik, and Kristen Grauman.  Poni: Potential functions for objectgoal navigation with interaction-free learning.  In Computer Vision and Pattern Recognition (CVPR), 2022 IEEE Conference on. IEEE, 2022.
  * [26] Xuanchi Ren and Xiaolong Wang.  Look outside the room: Synthesizing a consistent long-term 3d scene video from a single image.  In CVPR, 2022.
  * [27] Marc Rigter, Jun Yamada, and Ingmar Posner.  World models via policy-guided trajectory diffusion.  arXiv preprint arXiv:2312.08533, 2023.
  * [28] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models, 2021.
  * [29] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation, 2015.
  * [30] Nur Muhammad Mahi Shafiullah, Chris Paxton, Lerrel Pinto, Soumith Chintala, and Arthur Szlam.  Clip-fields: Weakly supervised semantic fields for robotic memory.  arXiv preprint arXiv: Arxiv-2210.05663, 2022.
  * [31] Quan Sun, Yufeng Cui, Xiaosong Zhang, Fan Zhang, Qiying Yu, Zhengxiong Luo, Yueze Wang, Yongming Rao, Jingjing Liu, Tiejun Huang, and Xinlong Wang.  Generative multimodal models are in-context learners.  arXiv preprint arXiv:2312.13286, 2023.
  * [32] Zineng Tang, Ziyi Yang, Chenguang Zhu, Michael Zeng, and Mohit Bansal.  Any-to-any generation via composable diffusion.  In Thirty-seventh Conference on Neural Information Processing Systems, 2023.
  * [33] Peng Wang, An Yang, Rui Men, Junyang Lin, Shuai Bai, Zhikang Li, Jianxin Ma, Chang Zhou, Jingren Zhou, and Hongxia Yang.  Ofa: Unifying architectures, tasks, and modalities through a simple sequence-to-sequence learning framework.  CoRR, abs/2202.03052, 2022.
  * [34] Erik Wijmans, Samyak Datta, Oleksandr Maksymets, Abhishek Das, Georgia Gkioxari, Stefan Lee, Irfan Essa, Devi Parikh, and Dhruv Batra.  Embodied Question Answering in Photorealistic Environments with Point Cloud Perception.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [35] Andy Zeng, Pete Florence, Jonathan Tompson, Stefan Welker, Jonathan Chien, Maria Attarian, Travis Armstrong, Ivan Krasin, Dan Duong, Vikas Sindhwani, and Johnny Lee.  Transporter networks: Rearranging the visual world for robotic manipulation.  Conference on Robot Learning (CoRL), 2020.
  * [36] Lvmin Zhang, Anyi Rao, and Maneesh Agrawala.  Adding conditional control to text-to-image diffusion models, 2023.
  * [37] Kaizhi Zheng, Xuehai He, and Xin Eric Wang.  Minigpt-5: Interleaved vision-and-language generation via generative vokens, 2023.
  * [38] Xiyue Zhu, Vlas Zyrianov, Zhijian Liu, and Shenlong Wang.  Mapprior: Bird’s-eye view perception with generative models, 2023.


# The generalized Bernoulli numbers and its relation with the Riemann zeta
function at odd-integer arguments

Yayun Wu
School of Mathematics and Statistics, Hefei Normal University,
Hefei, 230601, People’s Republic of China E-mail address<EMAIL_ADDRESS>(Y.
Wu).

Abstract: By using the generalized Bernoulli numbers, we deduce new integral
representations for the Riemann zeta function at positive odd-integer
arguments. The explicit expressions enable us to obtain criteria for the
dimension of the vector space spanned over the rational by the
$\zeta(2n+1)/\pi^{2n}$, $n\geq 1$.

Key words: Riemann zeta function; Integral representation; Generalized
Bernoulli numbers; Hyperbolic functions

Mathematics Subject Classification: 11M06; 11C08; 11J72; 30A05

## 1\. Introduction

The Riemann zeta function $\zeta(s)$ is defined for complex values of $s$,
$\Re(s)>1$ as $\zeta(s):=\sum_{v=1}^{\infty}v^{-s}$, and its meromorphic
continuation over the whole complex $s-$plane, except for a simple pole at
$s=1$. The research of the arithmetical properties of the Riemann zeta
function at positive integers is one of the fascinating topics of the complex
analysis and number theory. Finding recurrence formulas and integral
representations of zeta function has become a significant issue. As is well
known, for Riemann’s zeta function $\zeta(s)$ at even points, we have the
formula (found by Euler)

$\displaystyle\zeta(2k)=\frac{(-1)^{k-1}(2k)^{2k}B_{2k}}{2(2k)!},~{}k\in\mathbb{N}^{+}.$

Here the $B_{k}=B_{k}(0)$ is the well-tabulated Bernoulli numbers defined by
means of the generating function

$\displaystyle\frac{ze^{zx}}{e^{z}-1}=\sum_{k=0}^{\infty}\frac{B_{k}(x)}{n!}z^{k},~{}|z|<2\pi.$

However, there is no analogous closed evaluation for the values of $\zeta(s)$
at positive odd integers or fractional points up to the present time. Using
the series of special functions and the multidimensional integral to
investigate $\zeta(s)$ are the main methods. Connecting Euler polynomials and
trigonometric functions, Srivastava and Choi[3] give the Riemann zeta function
at positive odd integers

$\displaystyle\zeta(2n+1)=\frac{(-1)^{n}\pi^{2n+1}}{4\delta(1-2^{-(2n+1)})(2n)!}\int_{0}^{\delta}E_{2n}(u)\csc(u\pi){\rm
d}u,~{}n\in\mathbb{N}^{+},~{}\delta=1~{}{\rm or}~{}1/2,$

where $E_{n}(u)$ is the Euler polynomials. Zurab[1], Cvijovi$\acute{c}$ and
Klinowski[2] establish similar results by using Euler polynomials and the
polylogarithm function[8]. Generalizing the trigonometric function to
hyperbolic version, Lima[6], D’Avanzo and Krylov[5] research the integral
result about $\zeta(n)$,

$\displaystyle\zeta(n)=\frac{2^{n}}{(2^{n}-1)(n-1)!}\int_{0}^{\infty}\ln^{n-1}(\coth(x)){\rm
d}x,~{}n\geq 2.$

These integral representations about $\zeta(s)$ hardly show the correlation
with the irrational problem of the zeta function values at odd integers. Since
the problem of the irrationality of the values of the zeta function at odd
integers is one of the most attractive topics in modern number theory. These
problems motivate us to modify the integral representations about $\zeta(s)$.

The article is organized as follows. In $\S$ 2, we reconsider the evaluation
of $\zeta(3)$ so that technical details of the general case do not obscure. In
$\S$ 3, we elaborate the general case by using the generalized Bernoulli
polynomials and numbers. In $\S$ 4, we evaluate the numerical results about
$\zeta(5)$ and $\zeta(7)$. In $\S$5, the formula for $\zeta(2n+1)$ enables us
to obtain linear forms with rational coefficients in the values of
$\zeta(2n+1)/\pi^{2n}$, $n\in\mathbb{N}^{+}$. We also ponder over the
arithmetical properties of the sequence
$\zeta(3)/\pi^{2},\zeta(5)/\pi^{4},\zeta(7)/\pi^{6}$,$\cdots$.

## 2\. Evaluation of $\zeta(3)$

In this part, we rediscover an integral representation for $\zeta(3)$. Let us
briefly describe the proof for a better exposition of the general case and its
consequences.

Lemma 2.1 For Riemann zeta function value $\zeta(3)$, we have

$\displaystyle\zeta(3)=\frac{2\pi^{2}}{7}\int_{-\infty}^{+\infty}\frac{e^{u}(e^{u}-1)}{(e^{u}+1)^{3}u}du.$
(1)

Proof. For the part sum of the series $\sum_{k=1}^{n}\frac{1}{(2k-1)^{3}}$,
the following expression is valid

$\displaystyle\sum_{k=1}^{n}\frac{1}{(2k-1)^{3}}=\frac{\pi^{2}}{4}\int_{\gamma_{n}}\frac{e^{z}(e^{z}-1)}{(e^{z}+1)^{3}z}{\rm
d}z,~{}n\in\mathbb{N}^{+},$ (2)

where
$\gamma_{n}=\Gamma_{1}^{(n)}+\Gamma_{2}^{(n)}+\Gamma_{3}^{(n)}+\Gamma_{4}^{(n)}+\Gamma_{\varepsilon}$
represent the contour with vertices $\pm R_{n}$, $R_{n}+{\rm i}R_{n}$,
$-R_{n}+{\rm i}R_{n}$ (see Fig. 1). $\Gamma_{\varepsilon}$ is a semicircle of
radius $\varepsilon\in(0,\pi/2)$ centered at the origin. $\Gamma_{1}^{(n)}$
represent $(-R_{n},0)\rightarrow(-\varepsilon,0)$ and
$(\varepsilon,0)\rightarrow(R_{n},0)$. $R_{n}=2n\pi$. In fact, since in the
neighborhood of the point $z=0$ and $z=z_{k}=(2k-1)\pi{\rm i}$,
$k=1,2,\cdots,n$, the following results are valid:

$\displaystyle\frac{e^{z}-1}{z}=O(1),\quad\frac{e^{z}}{(e^{z}+1)^{3}}=O((z-z_{k})^{-3}).$

Figure 1: The contour $\gamma_{n}$

Let $f(z)=\frac{e^{z}(e^{z}-1)}{(e^{z}+1)^{3}z}.$ The function $f(z)$ is
bounded on the boundary of the square with poles at the points $z_{k}$. It
follows that the integral in $(2)$ is equal to the sum of residues of the
integrand at the points $z_{1}=\pi{\rm i},z_{2}=3\pi{\rm
i},\cdots,z_{n}=(2n-1)\pi{\rm i}$,

$\displaystyle\int_{\gamma_{n}}f(z){\rm d}z=2\pi{\rm
i}\sum_{k=1}^{n}\underset{z=z_{k}}{\text{Res}}f(z)=\frac{4}{\pi^{2}}\sum_{k=1}^{n}\frac{1}{(2k-1)^{3}}.$
(3)

Furthermore, $|f(z)|=O(\frac{1}{R_{n}e^{R_{n}}})$,
$z\in\Gamma_{j}^{(n)},j=2,3,4$. In view of
$\sum\limits_{k=1}^{\infty}\frac{1}{(2k-1)^{3}}=7\zeta(3)/8$, and let
$R_{n}\rightarrow+\infty$, $\varepsilon\rightarrow 0$ in $\gamma_{n}$, we
obtain the required expression

$\displaystyle\int_{-\infty}^{+\infty}\frac{e^{u}(e^{u}-1)}{(e^{u}+1)^{3}u}{\rm
d}u=\frac{7}{2\pi^{2}}\zeta(3),$ (4)

which leads to the desired result. $\hfill\square$

Corollary 2.1 Remembering inverse hyperbolic secant function: ${\rm
asech}(u)=\log(1/u+\sqrt{1/u^{2}-1})$, $u\in(0,1]$, then we obtain

$\displaystyle\zeta(3)=-\frac{\pi^{2}}{14}\int_{0}^{\infty}\frac{1}{u}{\rm
d}{\rm sech^{2}(u)}=\frac{\pi^{2}}{7}\int_{0}^{1}\frac{u}{{\rm asech(u)}}{\rm
d}u.$

By using the integral transformation $\log u=x$ and remembering that
$\int_{0}^{\infty}2(1+x)^{-3}x^{t}{\rm d}x=\pi t(1-t)/\sin(\pi t)$,
$t\in(0,1)$, we find $(4)$ equals with $(22)$ in [1]:

$\displaystyle\zeta(3)=\frac{2\pi^{2}}{7}\int_{0}^{\infty}\frac{x-1}{(1+x)^{3}\log
x}{\rm
d}x=\frac{2\pi^{2}}{7}\int_{0}^{1}\int_{0}^{\infty}\frac{x^{t}}{(1+x)^{3}}{\rm
d}x{\rm d}t=\frac{1}{7}\int_{0}^{\pi}\frac{x(\pi-x)}{\sin x}{\rm d}x,$

which means the result $(4)$ is a special case of the more general result in
[2]. But for $\zeta(2n+1),n\geq 2$, we found no association between these
representations.

However, let us expand the integral representation of $\zeta(3)$ into the sum
of partial fractions with respect to $\exp(u)$:

$\displaystyle\zeta(3)=\frac{2\pi^{2}}{7}\sum_{l=1}^{3}\int_{-\infty}^{+\infty}\frac{w_{l}}{u(e^{u}+1)^{l}}{\rm
d}u,$ (5)

where $w_{1}=1,w_{2}=-3,w_{3}=2$. Since $\sum\limits_{l=1}^{3}w_{l}=0$, the
integrand $\sum\limits_{l=1}^{3}\frac{w_{l}}{u(e^{u}+1)^{l}}$ is an even
function in $\mathbb{R}$ except a removable singularity at origin. We wonder
whether the decomposition (5) is unique and can be generalized to the general
$\zeta(2n+1),n\geq 2$.

## 3\. Evaluation of $\zeta(m)$ at positive odd integers

### 3.1 A formula for $\zeta(m),m\geq 2$

The evaluation of $\zeta(3)$ can be straightforwardly generalized to
$\zeta(m)$, $m\geq 2$. Before considering the general case, we need the
following useful lemma.

Lemma 3.1 Let $m\geq 2$ be a positive integer, and $\gamma_{n}$ represent the
integral contour in Lemma 2.1. There exists unique $w_{l}\in\mathbb{Q}$, which
satisfying $\sum\limits_{l=1}^{m}w_{l}=0$, such that

$\displaystyle\sum_{l=1}^{m}\int_{\gamma_{n}}\frac{w_{l}}{z(e^{z}+1)^{l}}{\rm
d}z=\frac{2S_{m}}{({\rm
i}\pi)^{m-1}}\left(\zeta(m)\frac{2^{m}-1}{2^{m}}-\sum_{k=n+1}^{\infty}\frac{1}{(2k-1)^{m}}\right),~{}n=1,2,\cdots,$

wherein $S_{m}\in\mathbb{Q}$.

Proof. Step 1. Let
$f_{m}(z)=\sum\limits_{l=1}^{m}\frac{w_{l}}{z(e^{z}+1)^{l}}$. Note that the
integrand $f_{m}(z)$ over the vertical and top side $\Gamma^{(n)}_{j}$,
$j=2,3,4$ tend to $0$, when $R_{n}\rightarrow\infty$, since it can be bounded
by $Ce^{-R_{n}}/R_{n}$ and constant $C>0$. Let $z_{k}=(2k-1)\pi{\rm i}$,
$k=1,2,\cdots,n$. By means of the theorem of residues, we have

$\displaystyle\int_{\gamma_{n}}f_{m}(z){\rm d}z$ $\displaystyle=2\pi{\rm
i}\sum_{k=1}^{n}\lim_{z\rightarrow z_{k}}\frac{1}{(m-1)!}\left(\frac{{\rm
d}}{{\rm
d}z}\right)^{m-1}\sum_{l=1}^{m}\frac{w_{l}(z-z_{k})^{m}}{z(e^{z}+1)^{l}}$
$\displaystyle=2\pi{\rm i}\sum_{k=1}^{n}\lim_{z\rightarrow
z_{k}}\frac{1}{(m-1)!}\sum_{l=1}^{m}w_{l}\sum_{j=0}^{m-1}\binom{m-1}{j}\left[\frac{(z-z_{k})^{l}}{z(e^{z}+1)^{l}}\right]^{(j)}\left[(z-z_{k})^{m-l}\right]^{(m-j-1)}$
$\displaystyle=2\pi{\rm
i}\sum_{k=1}^{n}\sum_{l=1}^{m}\frac{w_{l}}{(l-1)!}\lim_{z\rightarrow
z_{k}}\left[\frac{(z-z_{k})^{l}}{z(e^{z}+1)^{l}}\right]^{(l-1)}.$ (6)

Step 2. Let us denote

$\displaystyle
g_{l}(z)=\frac{(z-z_{k})^{l}}{z(e^{z}+1)^{l}},~{}D_{l}=\lim_{z\rightarrow
z_{k}}\frac{g_{l}^{(l-1)}(z)}{(l-1)!},~{}l=1,2,\cdots,m.$ (7)

The conclusions

$g_{l}^{(l-1)}(z)=\sum_{j=0}^{l-1}\binom{l-1}{j}\left(\frac{1}{z}\right)^{(l-1-j)}\left[\left(\frac{z-z_{k}}{e^{z}+1}\right)^{l}\right]^{(j)},$

follow from the Leibniz rule for the differentiation of a product. For $j=0$
in $g_{l}^{(l-1)}(z)$, we have

$\displaystyle\lim_{z\rightarrow
z_{k}}\binom{l-1}{0}\left(\frac{1}{z}\right)^{(l-1)}\left(\frac{z-z_{k}}{e^{z}+1}\right)^{l}=-\frac{(l-1)!}{z_{k}^{l}}.$
(8)

According with
$\big{[}\frac{1}{(e^{z}+1)^{l}}\big{]}^{\prime}+\frac{l}{(e^{z}+1)^{l}}=\frac{l}{(e^{z}+1)^{l+1}}$,
which means for all $l$, the residues $c_{-1}^{(l)}$ of $(e^{z}+1)^{-l}$ at
$z=z_{k}$ are equal:

$c_{-1}^{(l)}=c_{-1}^{(1)}=\frac{1}{2\pi{\rm
i}}\int_{\Lambda_{k}}\frac{1}{(e^{z}+1)^{l}}{\rm d}z,~{}l=1,2,\cdots,m,$

where $\Lambda_{k}$ is the small circle centered at the poles $z_{k}$ and of
radius $\epsilon_{k}>0$. When $l=1$, we acquire
$c_{-1}^{(l)}=c_{-1}^{(1)}=\lim\limits_{z\rightarrow
z_{k}}(z-z_{k})/(e^{z}+1)=-1$. For $j=l-1$ in $g_{l}^{(l-1)}(z)$, by using
Laurent series of $(e^{z}+1)^{-l}$ at $z=z_{k}$, we obtain

$\displaystyle\lim_{z\rightarrow
z_{k}}\binom{l-1}{l-1}\frac{1}{z}\left[\left(\frac{z-z_{k}}{e^{z}+1}\right)^{l}\right]^{(l-1)}=\lim_{z\rightarrow
z_{k}}\frac{1}{z}[c_{-1}^{(l)}(z-z_{k})^{l-1}]^{(l-1)}=-\frac{(l-1)!}{z_{k}}.$
(9)

As a result of the relations $(8)$ and $(9)$, we can write out $D_{l}$:

$\displaystyle
D_{l}=-(\frac{b_{l}^{(l)}}{z_{k}^{l}}+\cdots+\frac{b_{j}^{(l)}}{z_{k}^{j}}+\cdots+\frac{b_{1}^{(l)}}{z_{k}}),~{}j=1,\cdots,l,$
(10)

in which $b_{l}^{(l)}=b_{1}^{(l)}=1$,
$b_{2}^{(l)},\cdots,b_{l-1}^{(l)}\in\mathbb{Q}$.

Step 3. In this part, we will prove that there exists unique $w_{l}$, which
satisfying

$\sum_{l=1}^{m}w_{l}D_{l}=\frac{S_{m}}{z_{k}^{m}},~{}S_{m}\in\mathbb{Q}.$

In fact, it should be noted that, by applying $(10)$ in Step 2, we obtain

$\displaystyle\sum_{l=1}^{m}w_{l}D_{l}=-\sum_{l=1}^{m}w_{l}(\frac{1}{z_{k}^{l}}+\cdots+\frac{b_{j}^{(l)}}{z_{k}^{j}}+\cdots+\frac{1}{z_{k}})=\frac{S_{m}}{z_{k}^{m}}.$
(11)

Comparing the coefficients on both sides of identity $(11)$ with respect to
$1/z^{l}_{k}$, we have the relation:

$\displaystyle\begin{pmatrix}1&1&\cdots&1\\\ 0&1&\cdots&b_{2}^{(m)}\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&1\end{pmatrix}_{m\times
m}\begin{pmatrix}w_{1}\\\ w_{2}\\\ \vdots\\\
w_{m}\end{pmatrix}=\begin{pmatrix}0\\\ 0\\\ \vdots\\\ -S_{m}\end{pmatrix}.$
(12)

Let $A_{m\times m}$ represents the coefficient matrix of $(12)$. Since
$A_{m\times m}$ is an upper-triangular matrix, and $\|A_{m\times m}\|=1,$ then
by Cramer’s Rule, the existence and unique of $w_{l}$ can be obtained. Set
$C^{(m)}=[0,0,\cdots,-S_{m}]^{T}$, $W^{(m)}=[w_{1},w_{2},\cdots,w_{m}]^{T}$,
then

$\displaystyle w_{l}=[W^{(m)}]_{l},~{}l=1,\cdots,m,$ (13)

in which $W^{(m)}=A_{m\times m}^{-1}C^{(m)}$ and $[\cdot]_{l}$ represents the
$l-$th element in the column vector.

Given that $z_{k}=(2k-1)\pi{\rm i}$, and substitute $w_{l}$ and $D_{l}$ into
$(6)$, we obtain

$\displaystyle\sum_{l=1}^{m}\int_{\gamma_{n}}\frac{w_{l}}{z(e^{z}+1)^{l}}{\rm
d}z=2\pi{\rm i}\sum_{k=1}^{n}\frac{S_{m}}{z_{k}^{m}}=\frac{2S_{m}}{({\rm
i}\pi)^{m-1}}\Big{[}\zeta(m)(1-\frac{1}{2^{m}})-\sum_{k=n+1}^{\infty}\frac{1}{(2k-1)^{m}}\Big{]},$

which gives the required equivalence. Lemma 3.1 is thus proved.
$\hfill\square$

Corollary 3.1 Let

$S_{m}=\left\\{\begin{aligned}
&(-1)^{\frac{m-1}{2}}\Gamma(m),~{}~{}m=3,5,\cdots,\\\
&(-1)^{\frac{m+2}{2}}\Gamma(m),~{}~{}m=2,4,\cdots,\end{aligned}\right.$

where $\Gamma(\cdot)$ denotes the Gamma function and $\Gamma(m)=(m-1)!$. For
$n=1,2,\cdots$, by using Lemma 3.1, there exists unique $w_{l}\in\mathbb{Q}$,
$l=1,2,\cdots,m$, such that

$\displaystyle\sum_{l=1}^{m}\int_{\gamma_{n}}\frac{w_{l}}{z(e^{z}+1)^{l}}{\rm
d}z=\frac{2\Gamma(m)}{\pi^{m-1}{\rm
i}^{\kappa}}\Big{[}\zeta(m)(1-\frac{1}{2^{m}})-\sum_{k=n+1}^{\infty}\frac{1}{(2k-1)^{m}}\Big{]},$
(14)

in which $\kappa=0$ or $1$ with respect to $m$ is positive odd or even
integer.

### 3.2 A formula for $\zeta(m),m\geq 3$ at positive odd integers

The generalized Bernoulli polynomials[4] $B_{n}^{(l)}(x)$ of degree $n$ in $x$
are defined by the generating function:

$\displaystyle\left(\frac{z}{e^{z}-1}\right)^{l}e^{xz}=\sum_{n=0}^{\infty}B_{n}^{(l)}(x)\frac{z^{n}}{n!},~{}|z|<2\pi,1^{l}:=1,$
(15)

for arbitrary (real or complex) parameter $l$ and
$B_{n}^{(l)}(0)=B_{n}^{(l)}$. Let $m\geq 3$ be a positive odd integer. $w_{l}$
are solutions to $(12)$ and $z_{k}=(2k-1)\pi{\rm i}$, $k\geq 1$. Set
$x=z-z_{k}$ and $0<|x|<2\pi$. Using the property (22) of $B_{n}^{(l)}(x)$ in
[4], we obtain

$\displaystyle\frac{1}{(e^{z}+1)^{l}}$
$\displaystyle=\frac{(-1)^{l}}{(e^{x}-1)^{l}}=\sum_{n=0}^{\infty}(-1)^{l}B_{n}^{(l)}\frac{x^{n-l}}{n!},$
(16) $\displaystyle\frac{1}{(e^{-z}+1)^{l}}$
$\displaystyle=\frac{e^{xl}}{(e^{x}-1)^{l}}=\sum_{n=0}^{\infty}B_{n}^{(l)}(l)\frac{x^{n-l}}{n!}=\sum_{n=0}^{\infty}(-1)^{n}B_{n}^{(l)}\frac{x^{n-l}}{n!}.$
(17)

Since the Laurent series of $(e^{z}+1)^{-l}$ centered at $z=z_{k}$ have the
form

$\displaystyle\frac{1}{(e^{z}+1)^{l}}=\sum_{j=-l}^{\infty}c_{j}^{(l)}(z-z_{k})^{j},~{}l=1,2,\cdots,m,$
(18)

where

$\displaystyle c_{j}^{(l)}=\frac{1}{2\pi i}\int_{\Lambda_{k}}\frac{{\rm
d}z}{(e^{z}+1)^{l}(z-z_{k})^{j+1}},~{}j=-l,-(l-1),\cdots.$

Comparing the coefficients of $1/(z-z_{k})^{l}$ in $(16)$ and $(18)$, we
obtain

$\displaystyle
c_{j}^{(l)}=\frac{(-1)^{l}B_{l+j}^{(l)}}{\Gamma(l+j+1)},~{}j=-l,\cdots,-1.$
(19)

Taking $(7)$ and $(18)$ into account, we obtain

$\displaystyle
D_{l}=\sum_{j=0}^{l-1}(-1)^{l-j-1}\binom{l-1}{j}\frac{\Gamma(l-j)\Gamma(j+1)c_{-(l-j)}^{(l)}}{\Gamma(l)z_{k}^{l-j}}=\sum_{j=1}^{l}(-1)^{j+1}\frac{c_{-j}^{(l)}}{z_{k}^{j}}.$

Hence, it follows from $(10)$ that

$\displaystyle b_{j}^{(l)}=(-1)^{j}c_{-j}^{(l)}.$ (20)

Combining with $(19)$, we deduce that

$\displaystyle
b_{j}^{(l)}=\frac{(-1)^{l-j}}{\Gamma(l-j+1)}B_{l-j}^{(l)},~{}j=1,2,\cdots,l.$
(21)

With the help of $(17)$ and $(12)$, we obtain

$\displaystyle\begin{aligned} \frac{1}{2\pi{\rm
i}}\sum_{l=1}^{m}\int_{\gamma_{n}}\frac{w_{l}}{z(e^{-z}+1)^{l}}{\rm
d}z&=\frac{1}{2\pi{\rm
i}}\sum_{k=1}^{n}\sum_{j=1}^{m}\sum_{l=j}^{m}\int_{\Lambda_{k}}w_{l}(-1)^{l-j}\frac{B_{l-j}^{(l)}}{\Gamma(l-j+1)}\frac{1}{z(z-z_{k})^{j}}{\rm
d}z\\\
&=\sum_{k=1}^{n}\sum_{j=1}^{m}\sum_{l=j}^{m}w_{l}b_{j}^{(l)}\frac{(-1)^{j-1}}{z_{k}^{j}}=\sum_{k=1}^{n}\frac{(-1)^{m-1}w_{m}b_{m}^{(m)}}{z_{k}^{m}}\\\
&=-\sum_{k=1}^{n}\frac{S_{m}}{z_{k}^{m}}.\end{aligned}$

Hence, we have the relation

$\displaystyle\sum_{l=1}^{m}\int_{\gamma_{n}}\frac{w_{l}}{z(e^{z}+1)^{l}}{\rm
d}z$
$\displaystyle=\frac{1}{2}\sum_{l=1}^{m}\int_{\gamma_{n}}\frac{w_{l}}{z}\left[\frac{1}{(e^{z}+1)^{l}}-\frac{1}{(e^{-z}+1)^{l}}\right]{\rm
d}z$
$\displaystyle=\frac{1}{2}\sum_{l=1}^{m}\sum_{p=0}^{l-1}\int_{\gamma_{n}}\frac{(1-e^{z})w_{l}e^{zp}}{z(e^{z}+1)^{l}}{\rm
d}z.$ (22)

Since $z=0$ is a removable singularity pole for the integrand in $(22)$, let
$R_{n}\rightarrow+\infty$ and $\varepsilon\rightarrow 0$ in $\gamma_{n}$, then
we have

$\displaystyle
2\Gamma(m)(1-\frac{1}{2^{m}})\frac{\zeta(m)}{\pi^{m-1}}=\sum_{l=1}^{m}\int_{-\infty}^{+\infty}\frac{w_{l}}{u(e^{u}+1)^{l}}{\rm
d}u=\frac{1}{2}\sum_{l=1}^{m}\sum_{p=0}^{l-1}\int_{-\infty}^{+\infty}\frac{w_{l}(1-e^{u})e^{up}}{u(e^{u}+1)^{l}}{\rm
d}u.$

Note that the integrand
$\sum\limits_{l=1}^{m}\sum\limits_{p=0}^{l-1}\frac{w_{l}(1-e^{u})e^{up}}{u(e^{u}+1)^{l}}$
is an even function on $\mathbb{R}$ except a removable singularity at origin,
then we have

$\displaystyle\zeta(m)=\frac{(2\pi)^{m-1}}{(2^{m}-1)\Gamma(m)}\sum_{l=1}^{m}\sum_{p=0}^{l-1}\int_{0}^{+\infty}\frac{w_{l}(1-e^{u})e^{up}}{u(e^{u}+1)^{l}}{\rm
d}u,~{}~{}m=3,5,\cdots,$ (23)

in which $w_{l}$ are determined by $(13)$.

## 4\. Evaluation of $\zeta(5)$ and $\zeta(7)$ by using the generalized
Bernoulli numbers

In this part, we calculate the numerical results of $\zeta(5)$ and $\zeta(7)$.
In view of the Srivastava-Todorov formula [8, p.510, Eq.(3)], we know the
following representation for the generalized Bernoulli numbers:

$\displaystyle
B_{n}^{(l)}=\sum_{k=0}^{n}\binom{l+n}{n-k}\binom{l+k-1}{k}\frac{n!}{(n+k)!}\sum_{j=0}^{k}(-1)^{j}\binom{k}{j}j^{n+k}.$
(24)

Using the formula $(21)$ and $(24)$, it is verified by a straightforward
calculation that

$\displaystyle
D_{1}=-\frac{1}{z_{k}},~{}D_{2}=-\frac{1}{z_{k}^{2}}-\frac{1}{z_{k}},~{}D_{3}=-\frac{1}{2!}(\frac{2}{z_{k}^{3}}+\frac{3}{z_{k}^{2}}+\frac{2}{z_{k}}),$
$\displaystyle
D_{4}=-\frac{1}{3!}(\frac{3!}{z_{k}^{4}}+\frac{12}{z_{k}^{3}}+\frac{11}{z^{2}_{k}}+\frac{3!}{z_{k}}),~{}D_{5}=-\frac{1}{4!}(\frac{4!}{z^{5}_{k}}+\frac{60}{z_{k}^{4}}+\frac{70}{z_{k}^{3}}+\frac{50}{z_{k}^{2}}+\frac{4!}{z_{k}}),$
$\displaystyle
D_{6}=-\frac{1}{5!}(\frac{5!}{z_{k}^{6}}+\frac{360}{z_{k}^{5}}+\frac{510}{z_{k}^{4}}+\frac{450}{z_{k}^{3}}+\frac{274}{z_{k}^{2}}+\frac{5!}{z_{k}}),$
$\displaystyle
D_{7}=-\frac{1}{6!}(\frac{6!}{z_{k}^{7}}+\frac{2520}{z_{k}^{6}}+\frac{4200}{z_{k}^{5}}+\frac{4410}{z_{k}^{4}}+\frac{3248}{z_{k}^{3}}+\frac{1764}{z_{k}^{2}}+\frac{6!}{z_{k}}),$

where $z_{k}=(2k-1)\pi{\rm i},1\leq k\leq n$. According with $S_{m}$ in
Corollary $3.1$ and solving linear equation $(12)$, we obtain $w_{l}$ and the
following result.
Theorem 4.1 (1): For Riemann zeta function value $\zeta(5)$

$\displaystyle\zeta(5)=\frac{(2\pi)^{4}}{(2^{5}-1)\Gamma(5)}\sum_{l=1}^{5}\int_{-\infty}^{\infty}\frac{w_{l}}{(e^{u}+1)^{l}u}{\rm
d}z,$

where $w_{1}=-1,w_{2}=15,w_{3}=-50,w_{4}=60,w_{5}=-24$.

(2): For Riemann zeta function value $\zeta(7)$

$\displaystyle\zeta(7)=\frac{(2\pi)^{6}}{(2^{7}-1)\Gamma(7)}\sum_{l=1}^{7}\int_{-\infty}^{\infty}\frac{w_{l}}{(e^{u}+1)^{l}u}{\rm
d}z,$

where
$w_{1}=1,w_{2}=-63,w_{3}=602,w_{4}=-2100,w_{5}=3360,w_{6}=-2520,w_{7}=720$.

## 5\. The linear combination of $\zeta(2n+1)/\pi^{2n},n\geq 1$

In this part, we investigate the integral representation for the linear
combination of $\zeta(2n+1)/\pi^{2n}$. For all positive odd integer $m\geq 3$,
it follows from $(23)$ that

$\displaystyle
2\Gamma(m)(1-\frac{1}{2^{m}})\frac{\zeta(m)}{\pi^{m-1}}=-\sum_{l=1}^{m}\sum_{p=0}^{l-1}\int_{0}^{+\infty}\frac{w_{l}e^{(2p+1-l)u}}{u(e^{u}+e^{-u})^{l}}{\rm
d}(e^{u}+e^{-u}).$ (25)

Considering the symmetry of the power of $\exp(2p+1-l)$ with $p$ from $0$ to
$l-1$, we can draw the following conclusions:

$\displaystyle\sum_{p=0}^{l-1}\frac{e^{(2p+1-l)u}}{(e^{u}+e^{-u})^{l}}=\sum_{j=1}^{k}\frac{q_{j}^{(l)}}{(e^{u}+e^{-u})^{2j-1}},~{}~{}k=\lceil{l/2}\rceil,q_{j}^{(l)}\in\mathbb{Z},$
(26)

in which $\lceil{\cdot}\rceil$ means the ceiling function and
$q_{j}^{(l)}=1-\sum\limits_{\kappa=1}^{j-1}\binom{l+1-2\kappa}{j-\kappa}q_{\kappa}^{(l)},$
for $j\geq 2$, $q_{1}^{(l)}\equiv 1$. Remembering that
$\sum\limits_{l=1}^{m}w_{l}=0$ in Lemma 3.1, then we obtain

$\displaystyle 2\Gamma(m)(1-\frac{1}{2^{m}})\frac{\zeta(m)}{\pi^{m-1}}$
$\displaystyle=-\int_{0}^{\infty}\sum_{j=1}^{\lceil{m/2}\rceil}[\sum_{l=2j-1}^{m}w_{l}q_{j}^{(l)}]\frac{{\rm
d}(e^{u}+e^{-u})}{u(e^{u}+e^{-u})^{2j-1}}$
$\displaystyle=\int_{0}^{\infty}\sum_{j=2}^{\lceil{m/2}\rceil}[\sum_{l=2j-1}^{m}w_{l}q_{j}^{(l)}]\frac{2^{-2j+1}}{(j-1)u}{\rm
d}({\rm sech}{(u)})^{2j-2}$
$\displaystyle=-\int_{0}^{1}\sum_{j=2}^{\lceil{m/2}\rceil}\frac{1}{4^{j-1}}[\sum_{l=2j-1}^{m}w_{l}q_{j}^{(l)}]\frac{u^{2j-3}}{{\rm~{}asech}{(u)}}{\rm
d}u.$ (27)

Let us denote

$\displaystyle\tau_{j}^{(m)}=-\frac{2^{m-1}}{\Gamma(m)(2^{m}-1)}\sum\limits_{l=2j-1}^{m}\frac{w_{l}q_{j}^{(l)}}{4^{j-1}},$
(28)

then

$\displaystyle\frac{\zeta(m)}{\pi^{m-1}}=\sum_{j=2}^{\lceil{m/2}\rceil}\tau_{j}^{(m)}\int_{0}^{1}\frac{u^{2j-3}}{{\rm~{}asech}{(u)}}{\rm
d}u.$ (29)

Through above discussion, we have the following theorem.

Theorem 5.1 For any positive integer $n$, then there exists
$\theta_{k}\in\mathbb{Q}$, $k=1,2\cdots,n$, so that

$\sum_{k=1}^{n}\theta_{k}\frac{\zeta(2k+1)}{\pi^{2k}}=\theta_{n+1}I_{n},~{}\theta_{n+1}=0~{}or~{}1,$

where $I_{n}$ satisfy

$I_{n}=\int_{0}^{1}\frac{u^{2n-1}}{{\rm
asech}(u)}du>0,~{}~{}\lim_{n\rightarrow\infty}I_{n}=0.$

Proof. Since ${\rm asech(u)}=\log(1/u+\sqrt{1/u^{2}-1})$, then for $\vartheta$
near to $0$ from right side

$\displaystyle\frac{1}{{\rm asech(u)}}=-\frac{1}{\log u}-\frac{\log
2}{\log^{2}u}+o(1/\log^{2}u),~{}0<u<\vartheta<1,$ (30)

where $o(\cdot)$ means higher-order infinitesimal. Hence, we obtain $1/{\rm
asech}(u)\rightarrow 0$ as $u\rightarrow 0^{+}$. And for $\eta>0$ near to $1$
from left side, when $u\in(\eta,1)$,

$\displaystyle\frac{1}{{\rm
asech(u)}}=\frac{\sqrt{2}}{2\sqrt{1/u-1}}+\frac{\sqrt{2}}{24}\sqrt{1/u-1}+o(\sqrt{1/u-1}),$
(31)

then ${\rm asech(u)}/\sqrt{1/u-1}\rightarrow\sqrt{2}$ as $u\rightarrow 1^{-}$.
Since $\int_{0}^{1}(1/u-1)^{-1/2}du=\pi/2$, so we can deduce that
$\int_{0}^{1}1/{\rm asech(u)}{\rm du}$ converges. According with ${\rm
asech(u)}>0$, when $u\in(0,1)$, and $0<I_{n}\leq\int_{0}^{1}1/{\rm
asech(u)}{\rm du}$, then $\\{I_{n}\\}_{n\geq 1}$ is a monotonically decreasing
positive sequence with a lower bound, so it must converge.

From $(31)$, when $u\in(\eta,1)$, we can find proper $M>0$, so that ${\rm
1/asech(u)}\leq\frac{M}{\sqrt{1/u-1}}$. As a result, we have

$\displaystyle I_{n}$ $\displaystyle=\int_{0}^{1}\frac{u^{2n-1}}{{\rm
asech}(u)}du,$ $\displaystyle=\int_{0}^{\eta}\frac{u^{2n-1}}{{\rm
asech}(u)}du+\int_{\eta}^{1}\frac{u^{2n-1}}{{\rm asech}(u)}du$
$\displaystyle\leq\eta^{2n-1}\int_{0}^{1}\frac{1}{{\rm
asech}(u)}du+M\int_{0}^{1}u^{2n-1}(1/u-1)^{-1/2}du$
$\displaystyle\leq\eta^{2n-1}\int_{0}^{1}\frac{1}{{\rm
asech}(u)}du+M\sqrt{\pi}\frac{\Gamma(2n+1/2)}{\Gamma(2n+1)},$

since $\eta\in(0,1)$ and using Stirling’s approximation:
$\frac{\Gamma(2n+1/2)}{\Gamma(2n+1)}=O(1/\sqrt{2n+1})$ as
$n\rightarrow\infty$, so $\lim\limits_{n\rightarrow\infty}I_{n}=0$.

Using $(29)$, let $m=2n+1$, if for $\forall n\geq 1$, then
$\tau_{n+1}^{(2n+1)}\neq 0$, we can choose proper $\theta_{k}$, $1\leq k\leq
n$ so that $\theta_{n+1}=1$. If there exists $n\geq 2$ such that
$\tau_{n+1}^{(2n+1)}=0$, then we can choose proper $\theta_{k}$ (where
$\theta_{k}$ satisfying $\sum_{k=1}^{n}\theta_{k}^{2}\neq 0$) making
$\theta_{n+1}=0$. Therefore, there must exist proper $\theta_{k}\in\mathbb{Q}$
which make
$\sum\limits_{k=1}^{n}\theta_{k}\frac{\zeta(2k+1)}{\pi^{2k}}=\theta_{n+1}I_{n}$,
$\theta_{n+1}=0$ or $1$. Thus Theorem 5.1 is proved. $\hfill\square$

Since the problem of the irrationality of the values of the Riemann zeta
function at odd integers is one of the most attractive topics. Before Apery’s
results[10], the arithmetic properties of the zeta function at odd points
seemed inaccessible. According with Theorem $5.1$ in this paper and
Nesterenko’s theorem[9], we derive criteria for the dimension of the vector
space spanned over the rational by $\zeta(2n+1)/\pi^{2n}$, $n\geq 1$.

Theorem 5.2 Let $m=2n+1$ in $(29)$, if $\tau_{n+1}^{(2n+1)}\neq 0$ in $(28)$
for $n\geq 1$, then

$\displaystyle\dim_{\mathbb{Q}}(\mathbb{Q}\frac{\zeta(3)}{\pi^{2}}+\mathbb{Q}\frac{\zeta(5)}{\pi^{4}}+\cdots+\mathbb{Q}\frac{\zeta(2n+1)}{\pi^{2n}}+\cdots)=\infty;$
(32)

If there exits $n\geq 2$, so that $\tau_{n+1}^{(2n+1)}=0$, then

$\displaystyle\dim_{\mathbb{Q}}(\mathbb{Q}\frac{\zeta(3)}{\pi^{2}}+\mathbb{Q}\frac{\zeta(5)}{\pi^{4}}+\cdots+\mathbb{Q}\frac{\zeta(2n+1)}{\pi^{2n}})\leq
n-1.$ (33)

## 6\. Concluding remarks

Through numerical calculations of $\zeta(5)$, $\zeta(7)$ and $\zeta(9)$, we do
believe that only conclusion (32) is possible, which means the sequence
$\zeta(3)/\pi^{2},\zeta(5)/\pi^{4},\zeta(7)/\pi^{6}$,$\cdots$ contains
infinitely many irrational numbers. But it’s a pity that we can’t prove it
theoretically. This is also one of our future research contents.

## Conflict of interest

The author declares that there is no conflict of interest with respect to the
publication of this paper.

## Funding

Work funded by High-level Talent Research Project Fund of NHU (Grant Nos.
2022rcjj24).

## References

  * [1] Z. Silagadze, Sums of generalized harmonic series. Reson 20, 822-843 (2015).
  * [2] D. Cvijovi$\acute{c}$, J. Klinowski, Integral representations of the Riemann zeta function for odd-integer arguments. J. Comput. Appl. Math. 142, p.435-439 (2002).
  * [3] H. M. Srivastava, J. Choi, Zeta and Q-Zeta functions and associated series and integrals; Elsevier Inc., New York, NY,USA, 2012.
  * [4] J. ${\rm D}^{\prime}$Avanzo, N. A. Krylov, $\zeta(n)$ via hyperbolic functions, Involve, Vol.3, pp.289-296 (2010).
  * [5] F. M. S. Lima, New definite integrals and a two-term dilogarithm identity, Indag. Math., Vol.23, pp.1-9 (2012).
  * [6] F. Beukers, J.A.C. Kolk, E. Calabi, Sums of generalized harmonic series and volumes, Nieuw Archief voor Wiskunde, fourth series, Vol.11, pp.217-224 (1993).
  * [7] L. Lewin, Polylogarithms and associated functions, Elsevier (North-Holland), New York, London and Amsterdam, 1981.
  * [8] H. M. Srivastava, P. G. Todorov, An explicit formula for the generalized Bernoulli polynomials, J. Math. Anal. Appl. 130, pp.509-513 (1988).
  * [9] Yu. V. Nesterenko, On the linear independence of numbers, Vestnik Moskov. Univ. Ser. 1 Mat. Mekh, 1985, no.1, 46-54; English transl., Moscow Univ. Math. Bull. 40:1, 69-74 (1985).
  * [10] R. Ap$\acute{e}$ry, Irrationalit$\acute{e}$ de $\zeta(2)$ et $\zeta(3)$, Ast$\acute{e}$risque 61, 11-13 (1979).


# Heating a quantum dipolar fluid into a solid

J. Sánchez-Baena<EMAIL_ADDRESS>Center for Complex Quantum Systems,
Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C,
Denmark  Departament de Física, Universitat Politècnica de Catalunya, Campus
Nord B4-B5, 08034 Barcelona, Spain  Claudia Politi  Institut für Quantenoptik
und Quanteninformation, Österreichische Akademie der Wissenschaften,
Innsbruck, Austria  Institut für Experimentalphysik, Universität Innsbruck,
Austria  F. Maucher  Departament de Física, Universitat de les Illes Balears &
IAC-3, Campus UIB, E-07122 Palma de Mallorca, Spain  Center for Complex
Quantum Systems, Department of Physics and Astronomy, Aarhus University,
DK-8000 Aarhus C, Denmark  Francesca Ferlaino  Institut für Quantenoptik und
Quanteninformation, Österreichische Akademie der Wissenschaften, Innsbruck,
Austria  Institut für Experimentalphysik, Universität Innsbruck, Austria  T.
Pohl  Center for Complex Quantum Systems, Department of Physics and Astronomy,
Aarhus University, DK-8000 Aarhus C, Denmark<EMAIL_ADDRESS>
###### Abstract

Raising the temperature of a material enhances the thermal motion of
particles. Such an increase in thermal energy commonly leads to the melting of
a solid into a fluid and eventually vaporises the liquid into a gaseous phase
of matter. Here, we study the finite-temperature physics of dipolar quantum
fluids and find surprising deviations from this general phenomenology. In
particular, we describe how heating a dipolar superfluid from near-zero
temperatures can induce a phase transition to a supersolid state with a broken
translational symmetry. The predicted effect agrees with experimental
measurements on ultracold dysprosium atoms, which opens the door for exploring
the unusual thermodynamics of dipolar quantum fluids.

A supersolid is an exotic phase of matter in which particles develop regular
spatial order and simultaneously support the frictionless flow of a
superfluid. Having evaded experimental verification for several decades Chan
_et al._ (2013), supersolidity can now be observed in Bose-Einstein
condensates of ultracold atoms with finite-range interactions Léonard  _et
al._ (2017); Böttcher  _et al._ (2019a); Chomaz  _et al._ (2019); Tanzi  _et
al._ (2019a); Schuster  _et al._ (2020). Spontaneous symmetry breaking in
these systems occurs in the form of regular periodic patterns of the
condensate density as first predicted by Gross in 1957 Gross (1957). One would
thus expect the lowest possible temperatures to provide optimal conditions for
supersolidity by ensuring a high degree of phase coherence and maximal
population of the Bose-Einstein condensate. On the contrary, we demonstrate
here that thermal fluctuations in dipolar condensates do not merely diminish
global phase coherence but can instead facilitate the formation of periodic
modulations of the condensate density. This finding sheds light on recent
experimental observations and reveals an unusual fluid-solid phase transition,
whereby a supersolid state of matter emerges upon increasing the temperature.

Figure 1:  Heating a dipolar quantum fluid can lead to the emergence of a
supersolid phase of matter. (a) This is demonstrated in the thermodynamic
phase diagram for an infinitely elongated Bose-Einstein condensate in a radial
harmonic trap with no axial confinement, as illustrated in panel (c). In
between the superfluid (blue) and supersolid (red) region, both phases coexist
(purple region) as characteristic for a first-order phase transition. The
calculations were performed for a fixed chemical potential
$\mu/\epsilon_{d}=1$, where $\epsilon_{d}=\hbar^{2}/((12\pi)^{2}ma_{d}^{2})$
parametrizes the characteristic energy scale of the dipole-dipole
interactions. Our experiments in an ultracold gas of dysprosium atoms
demonstrate the temperature-driven emergence of supersolidity. The measured
contrast of axial density modulations is shown by the colored dots in panel
(b). The observations show that supersolidity at higher temperatures indeed
occurs for smaller atom numbers, $N_{\rm c}$ in the condensate, in agreement
with the theoretical transition line (purple line).

As we shall see below, this surprising behaviour arises from the anisotropic
nature of the dipole-dipole interaction

$V_{\rm dd}({\bf r})=\frac{C_{3}}{4\pi}\frac{1-3\cos^{2}\theta}{r^{3}},$ (1)

which has repulsive as well as attractive contributions, depending on the
angle $\theta$ between the atomic dipoles and the distance vector ${\bf r}$ of
two interacting atoms. The interaction strength $C_{3}$ and the atomic mass
$m$ define a length scale $a_{d}=mC_{3}/(12\pi\hbar^{2})$ that competes with
the scattering length $a$ of the short-range interaction between the atoms.
This competition between $a_{d}$ and $a>0$ can cause the condensate to
collapse when the stabilizing short-range repulsion is not sufficient to
overcome the attractive part of the dipole-dipole interaction between the
atoms Santos  _et al._ (2000); O’Dell  _et al._ (2004); Bortolotti  _et al._
(2006). Subsequent experiments Kadau  _et al._ (2016); Schmitt  _et al._
(2016); Chomaz  _et al._ (2016) have however found a higher level of
stability, which arises from quantum fluctuations Lima and Pelster (2011,
2012) that prevent the otherwise inevitable collapse of the condensate Petrov
(2015); Wächtler and Santos (2016); Cabrera  _et al._ (2018). In fact, the
balance of attraction and repulsion effectively enhances the role of quantum
fluctuations Wächtler and Santos (2016) beyond the semiclassical mean-field
physics of weakly interacting quantum gases. This yields a unique setting that
has revealed rich physics and a host of new quantum states, from self-bound
quantum droplets Schmitt  _et al._ (2016); Chomaz  _et al._ (2016); Böttcher
_et al._ (2019b) and supersolid phases Chomaz  _et al._ (2019); Böttcher  _et
al._ (2019a); Tanzi  _et al._ (2019a, b); Guo  _et al._ (2019) to complex
patterns in two-dimensional fluids Zhang  _et al._ (2021); Hertkorn  _et al._
(2021).

Given this striking role of quantum fluctuations in dipolar Bose-Einstein
condensates, one may also anticipate significant effects of thermal
fluctuations despite the ultralow temperatures that are required to reach
quantum degeneracy. To address this question we start from the grand canonical
potential $\Omega$ of the system at a finite temperature $T$. For a weakly
interacting gas with a high fraction of atoms in the condensate, one can use
Bogoliubov theory to determine $\Omega$. This yields simple expressions for
infinitely extended homogeneous systems Giorgini  _et al._ (1997) that can be
applied to describe trapped inhomogeneous gases within a local density
approximation. Hereby, one determines the Bogoliubov excitation spectrum and
all relevant observables for a homogeneous particle density $\rho$, which is
then identified as $\rho\equiv|\psi({\bf r})|^{2}$ with the local condensate
wave function $\psi({\bf r})$ at a given position ${\bf r}$. This permits to
express the grand canonical potential as

$\Omega=E_{0}+\frac{k_{\rm B}T}{(2\pi)^{3}}\int d{\bf r}\int d{\bf k}\text{
}\ln\\!\left(1-e^{-\frac{\varepsilon_{\bf k}\\!({\bf r})}{k_{\rm
B}T}}\right)\>,$ (2)

where $k_{\rm B}$ denotes the Boltzmann constant and $E_{0}$ is the zero-
temperature grand canonical energy that contains the mean-field interaction
energy and leading order corrections due to quantum fluctuations Lima and
Pelster (2011), i.e. small occupations of excited states above the formed
Bose-Einstein condensate. The dispersion $\varepsilon_{\bf k}=\sqrt{\tau_{\bf
k}(\tau_{\bf k}+2|\psi({\bf r})|^{2}\tilde{V}({\bf k}))}$ of these excitations
is determined by the kinetic energy $\tau_{\bf k}=\hbar^{2}k^{2}/(2m)$ of the
atoms and the Fourier transform $\tilde{V}({\bf
k})=\frac{4\pi\hbar^{2}a}{m}+\tilde{V}_{\rm dd}({\bf k})$ of their total
interaction potential.

Figure 2: Raising the temperature of a dipolar quantum fluid can induce a
pronounced roton-maxon spectrum of its collective excitations, as shown in
panel (a) for an infinitely elongated condensate along the $z$-axis [see
Fig.1(c)]. Heating the fluid tends to lower the energy of the roton minimum
and eventually softens the roton excitation as the temperature increases. This
effect can be traced back to the density dependence of the energy correction
caused by fluctuations, shown in panel (b). While quantum fluctuations yield
an energy $H_{\rm qu}$ (dashed line) that increases with a rising condensate
density $\rho=|\psi|^{2}$, the contribution $H_{\rm th}$ from thermal
fluctuations decreases (solid lines). The thermal energy correction $H_{\rm
th}({\bf r})$, therefore, acts as a focusing nonlinearity that supports the
formation of regular density modulations. This is illustrated in panel (c,d),
where we show the axial density $\rho_{||}(z)=\int{\rm d}x{\rm d}y\rho({\bf
r})$ along with the axial potential $\bar{H}_{\rm th}=\rho_{||}^{-1}\int{\rm
d}x{\rm d}y\rho({\bf r})\bar{H}_{\rm th}({\bf r})$, respectively. The
calculations are performed for $a/a_{\rm d}=0.7$ and $\mu=\varepsilon_{\rm
d}$.

Minimizing $\Omega$ with respect to $\psi({\bf r})$ then yields a nonlinear
wave equation that accounts for quantum as well as thermal fluctuations (see
Methods Section). At zero temperature, it describes the mean-field physics of
the condensate and captures leading-order effects of quantum fluctuations
through an effective density-dependent potential $H_{\rm qu}$ Wächtler and
Santos (2016) that increases the energy of the system. The second term in
Eq.(2) yields an additional potential

$\displaystyle H_{\rm th}({\bf r})$ $\displaystyle=\\!\int\\!\frac{{\rm d}{\bf
k}}{(2\pi)^{3}}\tilde{V}({\bf k})f_{\bf k}({\bf r})\frac{\tau_{\bf
k}}{\varepsilon_{\bf k}({\bf r})}\ ,$ (3)

that accounts for finite-temperature effects. It describes the interaction
between the condensate and thermally created excitations that populate
Bogoliubov modes according to the Bose distribution $f_{\bf
k}=1/(e^{\varepsilon_{\bf k}/k_{\rm B}T}-1)$. The resulting form of the
finite-temperature extended Gross-Pitaevskii equation (TeGPE) agrees with the
result of Hartree-Fock Bogoliubov theory Giorgini  _et al._ (1997); Aybar and
Oktel (2019), and includes relevant fluctuation terms that are commonly
neglected within the Popov approximation Griffin (1996) (see Methods Section).

Let us first use this framework to study an elongated atomic gas that is
confined harmonically in the $x-y$ plane and extends infinitely in the
$z$-direction without confinement along the $z$-axis. Figure 1(a) shows the
thermodynamic phase diagram obtained by simulating the imaginary time
evolution of the TeGPE at a fixed chemical potential $\mu$ (see Methods
Section). At zero temperature, we find a superfluid-supersolid quantum phase
transition, with a co-existence region that is expected for a first-order
phase transition Blakie  _et al._ (2020). While increasing the temperature may
generally be expected to melt the supersolid phase Cinti  _et al._ (2014), we
find instead that it shifts the transition towards weaker dipole-dipole
interactions. As a result, heating the system effectively drives a phase
transition from a fluid into a solid phase.

We can understand this effect from the excitation spectrum of the condensate
in the superfluid phase. To this end, we solve the time-dependent TeGPE within
linear response theory to find the excitation spectrum $\omega_{k_{z}}$ for
periodic plane-wave excitations along the $z$-direction. As shown in Fig.2,
the obtained dispersion exhibits the expected roton-maxon form Santos  _et
al._ (2003); Wilson  _et al._ (2008); Petter  _et al._ (2019); Schmidt  _et
al._ (2021), known from low-temperature helium Landau (1949) and Bose-Einstein
condensates with finite-range interactions O’Dell  _et al._ (2003); Henkel
_et al._ (2010); Mottl  _et al._ (2012); Zhang  _et al._ (2018). The local
minimum at finite momenta supports the formation of roton quasiparticles,
which were introduced by Landau as elementary vortices to describe
superfluidity in 4He Landau (1949). Experiments show that the roton minimum in
helium decreases with increasing temperature Dietrich  _et al._ (1972) due to
roton-roton scattering Kebukawa (1973). Yet, the roton energy remains sizable
at the transition to a normal-fluid phase Dietrich  _et al._ (1972), beyond
which it only varies weakly with temperature. The presence of a Bose-Einstein
condensate in dilute dipolar superfluids, however, enhances the effect of
thermal fluctuations due to the larger energy scale of the interaction between
Bogoliubov excitations and the condensate. A similar effect is found for atoms
with light-induced interactions and predicted to lower the roton minimum and
cause enhanced condensate depletion Mazets  _et al._ (2004). In the present
case, we find a thermal softening of the roton mode that can drive an
instability of the superfluid and thereby cause the formation of a supersolid
phase with increasing temperature.

Figure 3: Our measured axial density $|\psi_{||}(z)|^{2}$ (filled area),
observed for $T=76.5$ nK and $N_{\rm c}=13400$ condensate atoms, demonstrates
the formation of a supersolid state and agrees well with the numerical
simulation of the TeGPE (red line). On the other hand, equivalent simulations
at zero temperature (blue line) disagree qualitatively and instead yield an
unmodulated superfluid phase.

We can gain further intuition about the underlying mechanism by closer
inspection of the two fluctuation energies $H_{\rm qu}$ and $H_{\rm th}$ that
both contribute a local nonlinearity to the wave equation for $\psi({\bf r})$.
$H_{\rm qu}>0$ is the Lee-Huang-Yang correction to the equation of state Lima
and Pelster (2011, 2012), and raises the ground state energy due to the small
condensate depletion caused by the atomic interactions. It therefore increases
for higher particle densities and stronger interactions, as shown in Fig.
2(b). Consequently, $H_{\rm qu}$ generates an effective repulsion that
stabilizes the condensate against collapse Wächtler and Santos (2016), and
shifts the roton instability towards higher densities and stronger dipole-
dipole interactions. On the contrary, $H_{\rm th}$ increases as we lower the
density of the condensate [see Fig.2(b)]. This behaviour is readily understood
as follows. Decreasing the condensate density increases the fraction of
thermally excited, non-condensed atoms Aybar and Oktel (2019). In the limit
where this fraction remains small, such an increase implies a larger potential
energy due to interactions with the thermal atoms. It therefore contributes a
positive energy correction that decreases upon increasing the density
$\rho=|\psi|^{2}$ of the condensate. As a result, thermal fluctuations
energetically favour higher condensate densities, such that $H_{\rm th}$ acts
as a focusing nonlinearity which lowers the roton energy and facilitates the
formation of a density-modulated phase, as illustrated in Fig.2c.

We recently observed experimental signatures of this effect by studying the
cooling-heating lifecycle of bosonic dysprosium atoms at ultralow temperatures
Sohmen  _et al._ (2021). The experiment starts from a thermal cloud of
$~{}10^{5}$ atoms in a cigar-shaped optical dipole trap and traces the time
evolution of the gas as it is cooled evaporatively to quantum degeneracy by
lowering the depth of the trap. During the continual cooling and
thermalization we observed the expected emergence of supersolidity, and
studied the equilibrium states of the quantum fluid across the supersolid
phase transition. The measured density profiles indicate a higher degree of
modulation at higher temperatures. While this has cast mystery on the origin
of the observations, they can now be used to corroborate and benchmark our
theoretical understanding . Figure 1b shows our measured contrast of the axial
density modulations for different temperatures and condensed-atom numbers (see
Methods Section). The results confirm the formation of a supersolid phase with
increasing temperature in good agreement with the theoretical transition line
obtained numerically by the TeGPE. Moreover, Fig. 3 compares our measured
axial density, $\rho_{||}(z)$, to the theoretical prediction. The predicted
zero-temperature ground state corresponds to an unstructured superfluid and
deviates qualitatively from the observed supersolid state. The result of our
finite-temperature TeGPE simulation, however, agrees with the experiment and
reproduces quantitatively the period and amplitude of the measured density
modulations. This remarkable level of agreement offers strong indication that
the observed supersolid has indeed been generated by the finite temperature of
the atoms.

The possibility to make detailed comparisons between theory and experiments
opens up several directions for exploring the surprising thermodynamic
behaviour of quantum ferrofluids. Already, the ground state phase diagram
exhibits a rich structure, including first-order as well as second-order
quantum phase transitions in one and two-dimensional systems Blakie  _et al._
(2020); Zhang  _et al._ (2019); Biagioni  _et al._ (2022). This offers a
promising starting point for investigating how thermal fluctuations influence
the nature of the fluid-solid phase transition and may affect the physics of
higher dimensional supersolids Norcia  _et al._ (2021); Tanzi  _et al._
(2021), which can come in a diverse range of complex patterns Zhang  _et al._
(2021); Hertkorn  _et al._ (2021); Zhang  _et al._ (2019). Our present
findings motivate future experiments and first-principle simulations Macia
_et al._ (2016); Cinti  _et al._ (2017) to expand the phase diagram of Fig.1
into the high-temperature domain and draw a direct connection to the more
familiar physics of liquid-solid phase transitions in the absence of
superfluidity. Such numerical approaches may also reveal how the present
phenomenology extends into the regime of strong interactions, which is
becoming accessible in experiments with ultracold polar molecules Anderegg
_et al._ (2021); Schindewolf  _et al._ (2022); Schmidt  _et al._ (2022).
Equally important, an improved understanding of finite-temperature effects in
dipolar quantum fluids could help resolving current questions about
quantitative discrepancies between measurements and theory Böttcher  _et al._
(2019a); Petter  _et al._ (2019).

We thank Yongchang Zhang, Georg Bruun, and Jordi Boronat for valuable
discussions and the Er-Dy team in Innsbruck for experimental support. This
work was supported by the DNRF through the Center of Excellence ”CCQ” (Grant
agreement no.: DNRF156) and the Carlsberg Foundation through the ’Semper
Ardens’ Research Project QCooL. JB acknowledges funding by the European Union,
the Spanish Ministry of Universities and the Recovery, Transformation and
Resilience Plan through a grant from Universitat Politècnica de Catalunya. FF
and CP acknowledge support through an ERC Consolidator Grant (RARE, No.
681432), the QuantERA grant MAQS (No I4391-N) and the FOR grant (2247/PI2790)
by the Austrian Science Fund FWF.

##  Appendix A The nonlinear wave equation

The grand canonical potential is minimal in equilibrium such that we can
minimize Eq.(2) with respect to the condensate wave function $\psi({\bf r})$.
This yields the nonlinear wave equation

$\displaystyle\mu\psi({\bf r})=\left(-\frac{\hbar^{2}\nabla^{2}}{2m}+U({\bf
r})+\frac{4\pi\hbar^{2}a}{m}\absolutevalue{\psi({\bf r})}^{2}\right.$ (4)
$\displaystyle\left.+\\!\\!\int\\!{\rm d}{\bf r}^{\prime}V_{\rm dd}({\bf
r}-{\bf r}^{\prime})\absolutevalue{\psi({\bf r}^{\prime})}^{2}+H_{\rm qu}({\bf
r})+H_{\rm th}({\bf r})\right)\psi({\bf r}),$

which determines the equilibrium state of the condensate for a given chemical
potential $\mu$. We consider a harmonic trapping potential $U({\bf
r})=\frac{m}{2}(\omega_{x}x^{2}+\omega_{y}y^{2}+\omega_{z}z^{2})$, with
trapping frequencies $\omega_{x,y,z}$ along the three cartesian axes. The
first four terms correspond to the Gross-Pitaevskii equation that describes
the mean-field physics of the condensate at zero temperature. The next term is
given by $H_{\rm qu}({\bf r})=\gamma_{\rm qu}|\psi({\bf r})|^{3}$ and accounts
for leading order effects of quantum fluctuations with a strength $\gamma_{\rm
qu}$ that increases with $a$ and $a_{d}$ Lima and Pelster (2011, 2012);
Wächtler and Santos (2016). Finite-temperature effects are captured by the
last term as given in Eq.(3). We note here that the applied local-density
approximation can cause an infrared divergence of the momentum integral in
Eq.(3). However, the finite system size of trapped systems yields a natural
momentum cutoff that ensures converged results. Indeed, we find that our
calculated condensate wave functions are not sensitive to the precise choice
of the momentum cutoff for relevant trap geometries.

##  Appendix B Finite-temperature simulations

We have calculated the condensate wave function at finite temperatures by
simulating the imaginary time evolution of the wave equation (4). More
concretely, we replace $\mu\psi$ by $-\partial_{t}\psi$ in Eq.((4)) and
simulate the time evolution until $\psi({\bf r},t)$ reaches a steady state for
a given norm $N_{\rm c}=\int{\rm d}{\bf r}|\psi({\bf r},t)|^{2}$. $N_{\rm c}$
corresponds to the number of condensate atoms under 3D confinement as
considered in Figs.1(b) and 3, and yields the axial density $N_{\rm
c}/L=L^{-1}\int{\rm d}x{\rm d}y\int_{0}^{L}{\rm d}z|\psi({\bf r},t)|^{2}$ for
a given length $L$ of the periodic simulation box as considered in Figs.1(a)
and 2. Finally, we determine the chemical potential from Eq.(4) in order to
construct the thermodynamic phase diagram shown in Fig.1(a). The results shown
in Figs.1(b) and 3 are obtained for the experimental trap parameters
$\omega_{x}/2\pi=88$Hz, $\omega_{y}/2\pi=141$Hz, and $\omega_{z}/2\pi=36$Hz,
and a scattering length of $89a_{0}$, where $a_{0}$ is the Bohr radius. The
simulations of Figs.1(a) and 2 have been performed for
$w_{x}=0.0717\varepsilon_{\rm d}$, $w_{y}=0.142\varepsilon_{\rm d}$, and
$w_{z}=0$. In all cases, the dipoles are considered to be polarized along the
$y$-axis .

##  Appendix C Experimental determination of the average axial density and the
density contrast

We probe the emergence of a supersolid state via time-of-flight measurements
and in-situ Faraday phase contrast imaging Sohmen  _et al._ (2021). The former
provides information on the global phase coherence and the latter allows us to
extract the modulation contrast from the in-situ atomic density. In
particular, to obtain the modulation contrast of the experimental points in
Fig. 1c, we proceed as follows. For each $\\{N_{c},\text{ }T\\}$, we record
the in-trap density distribution and integrate along the direction orthogonal
to the droplet chain to get a 1D density profile $\rho_{||}(z)$. We repeat the
measurement under the same experimental conditions for 10-20 times. For each
profile, we then calculate its Fourier Transform, $\tilde{\rho}_{||}(k)=\int
e^{-ik_{z}z}\rho_{||}(z){\rm d}z$. Finally, we determine the average of
$|\tilde{\rho}_{||}(k)|$ and obtain the modulation contrast as the ratio
between the Fourier component at the modulation wavelength and the density
$|\tilde{\rho}_{||}(0)|$ at $k_{z}=0$.

The experimental profile shown in Fig. 3 is obtained by averaging the density
distributions from the repeated measurements. The central maximum in each
profile is shifted to the origin, $z=0$, to correct for unavoidable center-of-
mass fluctuations in the experiment. Moreover, we remove negative density
contributions in each density image, which are caused by small but inevitable
misalignments of the imaging objective.

## References

  * Chan  _et al._ (2013) M. H. W. Chan, R. B. Hallock, and Reatto L., “Overview on solid 4He and the issue of supersolidity,” Journal of Low Temperature Physics 172, 317 (2013).
  * Léonard  _et al._ (2017) J. Léonard, A. Morales, P. Zupancic, T. Esslinger, and T. Donner, “Supersolid formation in a quantum gas breaking a continuous translational symmetry,” Nature 543, 87 (2017).
  * Böttcher  _et al._ (2019a) F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, “Transient supersolid properties in an array of dipolar quantum droplets,” Phys. Rev. X 9, 011051 (2019a).
  * Chomaz  _et al._ (2019) L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, “Long-lived and transient supersolid behaviors in dipolar quantum gases,” Phys. Rev. X 9, 021012 (2019).
  * Tanzi  _et al._ (2019a) L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, “Observation of a dipolar quantum gas with metastable supersolid properties,” Phys. Rev. Lett. 122, 130405 (2019a).
  * Schuster  _et al._ (2020) S. C. Schuster, P. Wolf, S. Ostermann, S. Slama, and C. Zimmermann, “Supersolid properties of a Bose-Einstein condensate in a ring resonator,” Phys. Rev. Lett. 124, 143602 (2020).
  * Gross (1957) Eugene P. Gross, “Unified theory of interacting bosons,” Phys. Rev. 106, 161–162 (1957).
  * Santos  _et al._ (2000) L. Santos, G. V. Shlyapnikov, P. Zoller, and M. Lewenstein, “Bose-Einstein condensation in trapped dipolar gases,” Phys. Rev. Lett. 85, 1791 (2000).
  * O’Dell  _et al._ (2004) D. H. J. O’Dell, S. Giovanazzi, and C. Eberlein, “Exact hydrodynamics of a trapped dipolar Bose-Einstein condensate,” Phys. Rev. Lett. 92, 250401 (2004).
  * Bortolotti  _et al._ (2006) D. C. E. Bortolotti, S. Ronen, J. L. Bohn, and D. Blume, “Scattering length instability in dipolar Bose-Einstein condensates,” Phys. Rev. Lett. 97, 160402 (2006).
  * Kadau  _et al._ (2016) H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier, I. Ferrier-Barbut, and T. Pfau, “Observing the Rosensweig instability of a quantum ferrofluid,” Nature 530, 194 (2016).
  * Schmitt  _et al._ (2016) M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut, and T. Pfau, “Self-bound droplets of a dilute magnetic quantum liquid,” Nature 539, 259 (2016).
  * Chomaz  _et al._ (2016) L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, “Quantum-fluctuation-driven crossover from a dilute Bose-Einstein condensate to a macrodroplet in a dipolar quantum fluid,” Phys. Rev. X 6, 041039 (2016).
  * Lima and Pelster (2011) Aristeu R. P. Lima and Axel Pelster, “Quantum fluctuations in dipolar bose gases,” Phys. Rev. A 84, 041604 (2011).
  * Lima and Pelster (2012) A. R. P. Lima and A. Pelster, “Beyond mean-field low-lying excitations of dipolar Bose gases,” Phys. Rev. A 86, 063609 (2012).
  * Petrov (2015) D. S. Petrov, “Quantum mechanical stabilization of a collapsing Bose-Bose mixture,” Phys. Rev. Lett. 115, 155302 (2015).
  * Wächtler and Santos (2016) F. Wächtler and L. Santos, “Quantum filaments in dipolar Bose-Einstein condensates,” Phys. Rev. A 93, 061603 (2016).
  * Cabrera  _et al._ (2018) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, “Quantum liquid droplets in a mixture of Bose-Einstein condensates,” Science 359, 301 (2018).
  * Böttcher  _et al._ (2019b) F. Böttcher, M. Wenzel, J.-N. Schmidt, M. Guo, T. Langen, I. Ferrier-Barbut, T. Pfau, R. Bombín, J. Sánchez-Baena, J. Boronat, and F. Mazzanti, “Dilute dipolar quantum droplets beyond the extended Gross-Pitaevskii equation,” Phys. Rev. Research 1, 033088 (2019b).
  * Tanzi  _et al._ (2019b) L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Famà, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, “Supersolid symmetry breaking from compressional oscillations in a dipolar quantum gas,” Nature 574, 382 (2019b).
  * Guo  _et al._ (2019) M. Guo, F. Böttcher, J. Hertkorn, J.-N. Schmidt, M. Wenzel, H. P. Büchler, T. Langen, and T. Pfau, “The low-energy goldstone mode in a trapped dipolar supersolid,” Nature 574, 386 (2019).
  * Zhang  _et al._ (2021) Y.-C. Zhang, T. Pohl, and F. Maucher, “Phases of supersolids in confined dipolar Bose-Einstein condensates,” Phys. Rev. A 104, 013310 (2021).
  * Hertkorn  _et al._ (2021) J. Hertkorn, J.-N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, T. Langen, M. Zwierlein, and T. Pfau, “Pattern formation in quantum ferrofluids: From supersolids to superglasses,” Phys. Rev. Research 3, 033125 (2021).
  * Giorgini  _et al._ (1997) S. Giorgini, L. P. Pitaevskii, and S. Stringari, “Thermodynamics of a trapped Bose-condensed gas,” Journal of Low Temperature Physics 109, 309 (1997).
  * Aybar and Oktel (2019) E. Aybar and M. Ö. Oktel, “Temperature-dependent density profiles of dipolar droplets,” Phys. Rev. A 99, 013620 (2019).
  * Griffin (1996) A. Griffin, “Conserving and gapless approximations for an inhomogeneous Bose gas at finite temperatures,” Phys. Rev. B 53, 9341 (1996).
  * Blakie  _et al._ (2020) P. B. Blakie, D. Baillie, L. Chomaz, and F. Ferlaino, “Supersolidity in an elongated dipolar condensate,” Phys. Rev. Research 2, 043318 (2020).
  * Cinti  _et al._ (2014) F. Cinti, M. Boninsegni, and T. Pohl, “Exchange-induced crystallization of soft-core Bosons,” New Journal of Physics 16, 033038 (2014).
  * Santos  _et al._ (2003) L. Santos, G. V. Shlyapnikov, and M. Lewenstein, “Roton-maxon spectrum and stability of trapped dipolar Bose-Einstein condensates,” Phys. Rev. Lett. 90, 250403 (2003).
  * Wilson  _et al._ (2008) R. M. Wilson, S. Ronen, J. L. Bohn, and H. Pu, “Manifestations of the roton mode in dipolar Bose-Einstein condensates,” Phys. Rev. Lett. 100, 245302 (2008).
  * Petter  _et al._ (2019) D. Petter, G. Natale, R. M. W. van Bijnen, A. Patscheider, M. J. Mark, L. Chomaz, and F. Ferlaino, “Probing the roton excitation spectrum of a stable dipolar bose gas,” Phys. Rev. Lett. 122, 183401 (2019).
  * Schmidt  _et al._ (2021) J.-N. Schmidt, J. Hertkorn, M. Guo, F. Böttcher, M. Schmidt, K. S. H. Ng, S. D. Graham, T. Langen, M. Zwierlein, and T. Pfau, “Roton excitations in an oblate dipolar quantum gas,” Phys. Rev. Lett. 126, 193002 (2021).
  * Landau (1949) L. Landau, “On the theory of superfluidity,” Phys. Rev. 75, 884 (1949).
  * O’Dell  _et al._ (2003) D. H. J. O’Dell, S. Giovanazzi, and G. Kurizki, “Rotons in gaseous Bose-Einstein condensates irradiated by a laser,” Phys. Rev. Lett. 90, 110402 (2003).
  * Henkel  _et al._ (2010) N. Henkel, R. Nath, and T. Pohl, “Three-dimensional roton excitations and supersolid formation in Rydberg-excited Bose-Einstein condensates,” Phys. Rev. Lett. 104, 195302 (2010).
  * Mottl  _et al._ (2012) R. Mottl, F. Brennecke, K. Baumann, R. Landig, T. Donner, and T. Esslinger, “Roton-type mode softening in a quantum gas with cavity-mediated long-range interactions,” Science 336, 1570 (2012).
  * Zhang  _et al._ (2018) Y.-C. Zhang, V. Walther, and T. Pohl, “Long-range interactions and symmetry breaking in quantum gases through optical feedback,” Phys. Rev. Lett. 121, 073604 (2018).
  * Dietrich  _et al._ (1972) O. W. Dietrich, E. H. Graf, C. H. Huang, and L. Passell, “Neutron scattering by rotons in liquid helium,” Phys. Rev. A 5, 1377 (1972).
  * Kebukawa (1973) T. Kebukawa, “The Temperature Dependence of Phonon Velocity and Roton Minimum in Liquid He II,” Progress of Theoretical Physics 49, 388. (1973).
  * Mazets  _et al._ (2004) I. E. Mazets, D. H. J. O'Dell, G. Kurizki, Davidson N., and W. P. Schleich, “Depletion of a Bose-Einstein condensate by laser-induced dipole-dipole interactions,” Journal of Physics B: Atomic, Molecular and Optical Physics 37, S155 (2004).
  * Sohmen  _et al._ (2021) M. Sohmen, C. Politi, L. Klaus, L. Chomaz, M. J. Mark, M. A. Norcia, and F. Ferlaino, “Birth, life, and death of a dipolar supersolid,” Phys. Rev. Lett. 126, 233401 (2021).
  * Zhang  _et al._ (2019) Y.-C. Zhang, F. Maucher, and T. Pohl, “Supersolidity around a critical point in dipolar bose-einstein condensates,” Phys. Rev. Lett. 123, 015301 (2019).
  * Biagioni  _et al._ (2022) G. Biagioni, N. Antolini, A. Alaña, M. Modugno, A. Fioretti, C. Gabbanini, L. Tanzi, and G. Modugno, “Dimensional crossover in the superfluid-supersolid quantum phase transition,” Phys. Rev. X 12, 021019 (2022).
  * Norcia  _et al._ (2021) M. A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M. J. Mark, R. N. Bisset, L. Santos, and . Ferlaino, “Two-dimensional supersolidity in a dipolar quantum gas,” Nature 596, 357 (2021).
  * Tanzi  _et al._ (2021) L. Tanzi, J. G. Maloberti, G. Biagioni, A. Fioretti, C. Gabbanini, and G. Modugno, “Evidence of superfluidity in a dipolar supersolid from nonclassical rotational inertia,” Science 371, 1162 (2021).
  * Macia  _et al._ (2016) A. Macia, J. Sánchez-Baena, J. Boronat, and F. Mazzanti, “Droplets of trapped quantum dipolar bosons,” Phys. Rev. Lett. 117, 205301 (2016).
  * Cinti  _et al._ (2017) F. Cinti, A. Cappellaro, L. Salasnich, and T. Macrì, “Superfluid filaments of dipolar bosons in free space,” Phys. Rev. Lett. 119, 215302 (2017).
  * Anderegg  _et al._ (2021) L Anderegg, S. Burchesky, Y. Bao, S. S. Yu, T. Karman, E. Chae, K.-K. Ni, W. Ketterle, and J. M. Doyle, “Observation of microwave shielding of ultracold molecules,” Science 373, 779 (2021).
  * Schindewolf  _et al._ (2022) A. Schindewolf, R. Bause, X.-Y. Chen, M. Duda, T. Karman, I. Bloch, and X.-Y. Luo, “Evaporation of microwave-shielded polar molecules to quantum degeneracy,” (2022), arXiv:2201.05143 [cond-mat.quant-gas] .
  * Schmidt  _et al._ (2022) M. Schmidt, L. Lassablière, G. Quéméner, and T. Langen, “Self-bound dipolar droplets and supersolids in molecular Bose-Einstein condensates,” Phys. Rev. Research 4, 013235 (2022).


# NLOS Error Mitigation Using Weighted Least Squares and Kalman Filter in UWB
Positioning

Ruixin Fan and Xin Du Ruixin Fan and Xin Du are with the College of
Information Science & Electronic Engineering, Zhejiang University, China.
(E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In wireless positioning systems, non-line-of-sight (NLOS) is a challenging
problem. NLOS causes great ranging bias and location error, so NLOS mitigation
is essential for high accuracy positioning. In this letter, we propose the
Weighted-Least-Squares Robust Kalman Filter (WLS-RKF) for NLOS identification
and mitigation. WLS-RKF employs a hypothesis test based on Mahalanobis
distance for NLOS identification, and updates the corresponding Kalman filter
using the WLS solution. It requires no prior knowledge about NLOS distribution
or signal features. We perform simulations and experiments for ultra-wideband
(UWB) positioning in various scenarios. The results confirm that WLS-RKF
effectively mitigates NLOS error and achieves 5cm positioning accuracy.

###### Index Terms:

Wireless Positioning, Non-line-of-sight (NLOS), Kalman Filter, Weighted-Least-
Squares (WLS), Ultra-wideband (UWB).

##  I Introduction

Wireless positioning has been widely used for robot navigation and other
scenarios. For indoor positioning where the Global Navigation Satellite System
(GNSS) is unavilable, ultra-wideband (UWB) is often adopted for its high
accuracy and low cost. In a typical UWB positioning system, the distances
between the tag and the anchors can be measured based on time of flight (TOF),
and then the tag’s position can be estimated. However, if the signal is
obstructed, it propagates slower in the obstacle than in the air or arrives at
the receiver from reflections, which makes the measured distance larger than
the true distance [1]. This is referred to as the non-line-of-sight (NLOS)
effect, which leads to a significant location error.

Currently, there are two kinds of approaches for NLOS mitigation. The first
kind is to identify NLOS measurements based on the features of the received
signal, such as confidence metric, delay spread, and change of SNR [2]. [3]
proposes a decision tree-based method using the UWB channel quality
indicators. After NLOS identification, an accurate position is calculated by
the weighted-least-squares (WLS) technique [4]. Because the signal features
highly depend on the propagation environment, such algorithms are difficult to
be widely used [5]. The second kind is based on the measured distances. For
example, the residual-weighting (Rwgh) algorithm adopts distance residual as a
measure of the presence of NLOS[6]. Algorithms such as RCC-Rwgh [7] are
proposed to reduce the computational complexity of Rwgh. However, [8] shows
that these methods work only for location estimates within the convex-hull
region enclosed by the anchors. Based on the assumption that the measured
distance is positively biased when NLOS occurs, [9] uses the standard
deviation of measured distances for NLOS identification. Recently, [10]
proposes the robust Kalman Filter (RKF). RKF uses the Mahalanobis distance to
detect outliers of the measurements, which is probably affected by NLOS. Then
the effect on the posterior estimation can be mitigated by increasing the
covariance of the filter. However, we find that the positioning error
increases when the tag moves in the NLOS condition.

In this letter, we focus mainly on the second kind of approaches. We propose
the Weighted-Least-Squares Robust Kalman Filter (WLS-RKF) method for NLOS
mitigation. In WLS-RKF, the NLOS is identified by a hypothesis test based on
Mahalanobis distance. Then WLS technique is adopted where the NLOS measurement
gets a smaller weight, and the corresponding filter is updated using the WLS
solution. It uses only the measured distances and requires no prior knowledge
about NLOS distribution or signal features. We evaluate the method by
simulations and experiments in UWB positioning. The results show that the
location error is reduced by over 90% and positioning accuracy is up to 5cm.

The rest of this letter is organized as follows. In Section II, we give a
brief introduction of wireless positioning model, the Kalman filter and the
NLOS error model. In Section III, the WLS-RKF is presented and discussed in
detail. The simulation and experiment results for the scenario of UWB
positioning are shown in Section IV and V respectively. Finally, Section IV
gives a conclusion.

##  II Preliminaries

###  II-A Wireless Positioning Model

We consider a wireless positioning system with $N$ anchors $A_{1}$, $A_{2}$,
…, $A_{N}$ and a tag $T$. The position vector of $A_{i}$ is
$\bm{Z}_{i}=[x_{i}\ y_{i}]^{T}$, and the position vector of $T$ is $\bm{Z}=[x\
y]^{T}$. $r_{i}$ is the measured distance between $T$ and $A_{i}$. Once all
the distances are available, we have the following equations:

$||\bm{Z}-\bm{Z}_{i}||=\sqrt{(x-x_{i})^{2}+(y-y_{i})^{2}}=r_{i},\
i=1,2,...,N.$ (1)

Let $w_{i}$ denote the weight of the $i$th equation in (1). Then the WLS
estimate of the tag’s position can be expressed as

$\hat{\bm{Z}}=\arg\min_{\bm{Z}}\sum_{i=1}^{N}w_{i}^{2}(r_{i}-||\bm{Z}-\bm{Z}_{i}||)^{2},$
(2)

which can be solved by Gauss-Newton techniques [11] or other numerical search
methods.

###  II-B Kalman Filter

The Kalman filter is used to smooth the measurement data. From [9], the state
data vector satisfies:

$\bm{X}_{k}=\bm{A}\bm{X}_{k-1}+\bm{\Gamma}\bm{W}_{k-1},$ (3)

where $\bm{X}_{k}$ is the state vector at the time $t_{k}$. $\bm{W}_{k-1}$ is
the driving noise vector with covariance matrix $\bm{Q}=\sigma_{u}^{2}\bm{I}$
and

$\bm{A}=\left[\begin{array}[]{cc}\bm{I}&\Delta t\bm{I}\\\
\bm{O}&\bm{I}\end{array}\right],\quad\bm{\Gamma}=\left[\begin{array}[]{c}\bm{O}\\\
\Delta t\bm{I}\end{array}\right].$ (4)

The measurement process is

$\bm{Y}_{k}=\bm{H}\bm{X}_{k}+\bm{U}_{k},$ (5)

where $\bm{Y}_{k}$ is the measured data vector, $\bm{H}=[\bm{I}\ \bm{O}]$, and
$\bm{U}_{k}$ is the measurement noise vector with covariance matrix
$\bm{R}=\sigma_{x}^{2}\bm{I}$.

In each iteration, the Kalman filter includes the following steps:

$\displaystyle\hat{\bm{X}}_{k,k-1}$
$\displaystyle=\bm{A}\hat{\bm{X}}_{k-1,k-1}$ (6) $\displaystyle\bm{P}_{k,k-1}$
$\displaystyle=\bm{A}\bm{P}_{k-1,k-1}\bm{A}^{T}+\bm{\Gamma}\bm{Q}\bm{\Gamma}^{T}$
(7) $\displaystyle\bm{K}_{k}$
$\displaystyle=\bm{P}_{k,k-1}\bm{H}^{T}(\bm{H}\bm{P}_{k,k-1}\bm{H}^{T}+\bm{R})$
(8) $\displaystyle\hat{\bm{X}}_{k,k}$
$\displaystyle=\hat{\bm{X}}_{k,k-1}+\bm{K}_{k}(\bm{Y}_{k}-\bm{H}\hat{\bm{X}}_{k,k-1})$
(9) $\displaystyle\bm{P}_{k,k}$
$\displaystyle=(\bm{I}-\bm{K}_{k}\bm{H})\bm{P}_{k,k-1}$ (10)

where $\bm{K}_{k}$ is the Kalman gain and $\bm{P}_{k,k}$ is the covariance
matrix of $\hat{\bm{X}}_{k,k}$.

In wireless positioning, Kalman filter can be used to smooth distance
measurement data. The state vector is $\bm{X}_{k}=[r_{k}\ \dot{r}_{k}]^{T}$,
where $\dot{r}_{k}$ denotes the first derivative of $r_{k}$.

###  II-C NLOS Ranging Model

Due to the NLOS effect, distance measurement is positively biased. The
measured distance $r$ can be expressed as [1]:

$r=d+b+\epsilon,$ (11)

where $d$ is the true distance, $b$ is the ranging bias caused by NLOS
propagation, and $\epsilon$ is the measurement noise.

Measurement noise can be modeled as Gaussian distributed $\epsilon\sim
N(0,\sigma_{m})$. NLOS bias can be modeled in different ways such as
exponentially distributed, uniformly distributed and Gaussian distributed [4].
In [12], a simplified through-the-wall (TTW) NLOS ranging model is proposed.
When an RF wave traverses a wall, the ranging bias is

$b\approx w(\sqrt{\epsilon_{r}}-1)+0.31w\theta_{i}^{2},$ (12)

where $\epsilon_{r}$ is the wall’s real relative permittivity, $w$ is the
wall’s thickness, and $\theta_{i}$ is the incidence angle.

In general, the NLOS bias is much larger than the measurement noise, that is,
$b>>\epsilon$. Data from UWB range measurements shows that $\sigma_{m}$ is
only a few centimeters[13], while $b$ can be up to 1.5m [14].

##  III Weighted-Least-Squares
Robust Kalman Filter

In this section we detail the proposed algorithm, Weighted-Least-Squares
Robust Kalman Filter (WLS-RKF). The main idea is to identify the positive bias
of the NLOS measurements, and mitigate the position error using the other
line-of-sight (LOS) measurements.

###  III-A NLOS Identification

Kalman filter can be used to identify outliers in the stream of data. Similar
to [10], we use a Kalman filter KFi to handle the distance measurements to the
$i$th anchor respectively. In each iteration, the Kalman filter predicts the
distance ${d}_{ki}=\bm{H}\hat{\bm{X}}_{k,k-1}$. In LOS condition, the measured
distance $\bm{Y}_{k}=r_{ki}$ obeys the Gaussian distribution, with the mean of
${d}_{ki}$ and the variance of ($\bm{H}\bm{P}_{k,k-1}\bm{H}+\bm{R}$).
Therefore, the square of the Mahalanobis distance from $r_{ki}$ to ${d}_{ki}$
obeys the $\chi^{2}$ distribution[15] with the degree of freedom of 1, that is

$\gamma_{ki}=(r_{ki}-{d}_{ki})^{T}(\bm{H}\bm{P}_{k,k-1}\bm{H}+\bm{R})^{-1}(r_{ki}-{d}_{ki})\sim\chi_{1}^{2},$
(13)

For the significance level $\alpha$, we have

$\text{Pr}(\gamma_{ki}\leq\chi_{1,\alpha}^{2})=1-\alpha,$ (14)

where Pr() is the probability of a random event, and $\chi_{1,\alpha}^{2}$ is
the $\alpha$-quantile of the $\chi_{1}^{2}$ distribution. In this letter, we
choose $\alpha=0.001$, and accordingly, $\chi_{1,\alpha}^{2}=6.2$.

If the measured distance does not satisfy the condition
$\gamma_{ki}\leq\chi_{1,\alpha}^{2}$, it is considered as an outlier, which is
probably caused by NLOS. Since the NLOS bias is always positive, we add
another constraint $r_{ki}>{d}_{ki}$, which can reduce the probability of
false alarm. Therefore, the following hypothesis test is employed to identify
NLOS measurements:

$\displaystyle H_{0}:\gamma_{ki}>\chi_{1,\alpha}^{2}\ \text{and}\
r_{ki}>{d}_{ki}\quad\text{NLOS case}$ (15) $\displaystyle
H_{1}:\gamma_{ki}\leq\chi_{1,\alpha}^{2}\ \text{ or }\
r_{ki}\leq{d}_{ki}\quad\text{LOS case}$

###  III-B NLOS Mitigation

After NLOS measurements are identified, we use the WLS technique to mitigate
the position error. According to (2), the weights and the distances from the
tag to the anchors should be determined.

The distance $\hat{r}_{ki}$ is the output of the $i$th Kalman filter:

$\hat{r}_{ki}=\left\\{\begin{array}[]{cc}{d}_{ki},&\text{NLOS case}\\\
\bm{H}\hat{\bm{X}}_{k,k},&\text{ LOS case }\\\ \end{array}\right.$ (16)

where ${d}_{ki}$ denotes the predicted distance, and $\hat{\bm{X}}_{k,k}$ is
the estimated state vector by (9). $\bm{H}\hat{\bm{X}}_{k,k}$ denotes the
first element of $\hat{\bm{X}}_{k,k}$, which is the estimated distance. In LOS
case, the Kalman filter gives an accurate estimation of the distance. In NLOS
case, however, the measured distance ${r}_{ki}$ is biased and not reliable. So
we use the predicted distance.

The weight $w_{i}$ is selected as:

$w_{i}=\left\\{\begin{array}[]{cc}\dfrac{1}{\sqrt{{\gamma_{ki}}/{\chi_{1,\alpha}^{2}}}},&\text{NLOS
case}\\\ 1,&\text{ LOS case }\\\ \end{array}\right.$ (17)

where $\sqrt{{\gamma_{ki}}/{\chi_{1,\alpha}^{2}}}$ is the Mahalanobis distance
normalized with the threshold $\chi_{1,\alpha}^{2}$.

In NLOS case, according to (15), $w_{i}$ is less than 1 and inversely
proportional to the Mahalanobis distance, which is determined by the NLOS
bias. Therefore, distance measurement with a large NLOS bias will get a small
weight. The WLS solution is mainly determined by LOS measurements, and the
location error caused by NLOS bias is mitigated.

Then the estimated position of the tag $\hat{\bm{Z}}_{k}$ is the WLS solution:

$\hat{\bm{Z}}_{k}=\mathop{\arg\min}\limits_{\bm{Z}}\sum_{i=1}^{N}w_{i}^{2}(\hat{r}_{ki}-||\bm{Z}-\bm{Z}_{i}||)^{2}.$
(18)

In LOS case, we update the state of Kalman filter with the measured distance
$\bm{Y}_{k}=r_{ki}$. In NLOS case, however, this will make the filter will
enter the wrong state, since the measured distance $r_{ki}$ is biased.
Therefore, after we get $\hat{\bm{Z}}_{k}$, we update the state of Kalman
filter with the distance between $\hat{\bm{Z}}_{k}$ and the anchor’s position
$\bm{Z}_{i}$:

$\bm{Y}_{k}=||\hat{\bm{Z}}_{k}-\bm{Z}_{i}||.$ (19)

As discussed before, $\hat{\bm{Z}}_{k}$ is mainly determined by the LOS
measurements. So the state of the Kalman filter can be better estimated with
the information of the LOS measurements. This ensures that the filter gives an
accurate distance prediction for NLOS identification in the next iteration. By
comparison, in RKF [10], distance measurements to each anchor are handled with
an independent Kalman filter. Although the efect of NLOS can be mitigated by
increasing the covariance, the state of the filter still diverges slowly.

Algorithm 1 WLS-RKF

Input: number of anchors $N$, the threshold $\chi_{1,\alpha}^{2}$, measured
distances $r_{ki}$ at time $t_{k}$, $i=1,2,...,N,k=1,2,...$
Output: estimated positions $\hat{\bm{Z}}_{k}$ at time $t_{k}$, $k=1,2,...$

1: Initialization for Kalman filters KFi, $i=1,2,...,N$

2: for k = 1, 2, … do

3: for $i=1,2,...,N$ do

4: Predict distance ${d}_{ki}$ by KFi

5: Calculate $\gamma_{ki}$ by (13)

6: if $\gamma_{ki}>\chi_{1,\alpha}^{2}$ and $r_{ki}>{d}_{ki}$ then

7: $NLOS_{i}$ = true

8: $\hat{r}_{ki}={d}_{ki}$

9: $w_{i}=\sqrt{{\chi_{1,\alpha}^{2}}/{\gamma_{ki}}}$

10: else

11: $NLOS_{i}$ = false

12: Update state for KFi with $\bm{Y}_{k}=r_{ki}$

13: $\hat{r}_{ki}=\bm{H}\hat{\bm{X}}_{k,k}$

14: $w_{i}=1$

15: end if

16: end for

17:
$\hat{\bm{Z}}_{k}=\mathop{\arg\min}\limits_{\bm{Z}}\sum_{i=1}^{N}w_{i}^{2}(\hat{r}_{ki}-||\bm{Z}-\bm{Z}_{i}||)^{2}$

18: for $i=1,2,...,N$ do

19: if $NLOS_{i}$ == true then

20: Update state for KFi with $\bm{Y}_{k}=||\hat{\bm{Z}}_{k}-\bm{Z}_{i}||$

21: end if

22: end for

23: end for

The operation of the proposed WLS-RKF is summarized in Algorithm 1. For each
distance measurement $r_{ki}$, LOS/NLOS case is identified in Line 6. The NLOS
case is handled in Line 8-9 and 20. The LOS case is handled in Line 12-14.

It should be noted that enough LOS measurements is required to ensure the
position accuracy. For 2-D positioning, the number of LOS measurements should
be no less than 2. We suggest that the positioning system includes at least 4
anchors to achieve a better performance.

##  IV Simulations and Results

###  IV-A Simulation Setup

We consider 4 representative cases, as shown in Fig.1.

  * •

Case 1: There are 4 anchors with coordinates $A_{1}$(0, 0), $A_{2}$(10, 0),
$A_{3}$(10, 10) and $A_{4}$(0, 10). The tag $T$ moves from $(0,3)$ to $(10,3)$
in the straight line at a constant velocity of 0.5m/s.

  * •

Case 2: The settings are the same as case 1, and another anchor with
coordinate $A_{5}$(5, 15) is added.

  * •

Case 3: There are 4 anchors with the same coordinates as in case 1. The tag
$T$ moves around the rectangle at a constant velocity of 0.5m/s for two laps.
The rectangle is 8m long and 6m wide, and its corner is round with a radius of
0.5m.

  * •

Case 4: The settings are the same as case 3, and another anchor with
coordinate $A_{5}$(5, 15) is added.

(a)  (b)

Figure 1: Simulation cases. (a) Case 1 & 2\. (b) Case 3 & 4\. The anchor
$A_{5}$ is excluded in Case 1 and 3. The obstacle will introduce NLOS bias to
the distance measurements.

In case 1 and 2, there is a wall with the length $L$ and width $W$. The wall
makes the path $TA_{3}$, $TA_{5}$ and $TA_{4}$ change from LOS to NLOS
successively as $T$ moves along the line. In case 2 and 4, two walls with the
length $L_{1}$, $L_{2}$ and width $W$ is placed orthogonally. The walls
introduce a more severe obstruction: at least one path $TA_{i}$ is NLOS at
almost any point of the trajectory. $L$, $L_{1}$, $L_{2}$ and $W$ are randomly
chosen from 3-8m, 4-7m, 2-5m and 0.3-0.7m. The settings ensure that at least 2
LOS measurements exist at any time.

To create the measurement data for simulation, we calculate the true distances
from $T$ to $A_{i}$ with sampling interval 0.05s and add the measurement noise
and NLOS bias to the true distances. The measurement noise is Gaussian
distribution with standard deviation $\sigma_{m}=0.02m$ and the NLOS bias is
calculated with the TTW model 12, where $\epsilon_{r}=6$ and $\theta_{i}$ is
calculated from the geometry of $TA_{i}$ and the walls.

For the proposed WLS-RKF, we select $\chi_{1,\alpha}^{2}$=6.2. The parameters
of the Kalman filters are $\sigma_{x}=\sigma_{m}=0.02$ and $\sigma_{u}=0.5$.
And the initial value of $\bm{P}$ is $\bm{P}_{0}=diag(\sigma_{m}^{2},0)$. We
also simulate with the RKF [10] for comparison, and with unweighted LS as the
benchmark.

###  IV-B Simulation Results

(a)  (b)

Figure 2: Simulation results for LS, RKF and the proposed WLS-RKF. (a) Case 2.
(b) Case 4. The results of Case 1 and 3 are similar to (a) and (b).

(a)  (b)

Figure 3: The experiment setting. (a) Photo. (b) Schematic. For Line 1, a person obstructs the signal path between the anchor and the tag. For Line 2, the statue in the lobby is used as an obstacle. TABLE I: Location Error(cm) for Simulations Algorithm | Case 1 | Case 2 | Case 3 | Case 4
---|---|---|---|---
|  RMS |  90% |  RMS |  90% |  RMS |  90% |  RMS |  90%
LS |  41.2 |  61.6 |  44.7 |  65.9 |  78.5 |  137.6 |  68.3 |  115.2
RKF |  40.4 |  59.7 |  43.8 |  64.4 |  79.6 |  140.2 |  69.1 |  116.4
WLS-RKF |  1.7 |  2.1 |  1.9 |  2.0 |  1.9 |  3.3 |  1.8 |  3.0

We run 20 simulations for each case. Fig. 2 shows the trajectories of Case 2
and 4. From Fig. 2(a), we can see that the trajectories are obviously biased
due to NLOS when the tag passes by the obstacle, and the position error
increases as more paths to the anchors are obstructed. RKF slows down the
error growth but the trajectory still deviates gradually from the ground
truth. The proposed WLS-RKF effectively mitigates the NLOS error and achieves
accurate location results. From Fig. 2(b), we can see that even the initial
position A is biased due to NLOS. As the tag moves aside the obstacle, the
paths to the anchors changes from LOS to NLOS or otherwise. WLS-RKF would
adjust the NLOS judgement as the change occurs. The filtered trajectory
”converges” to the true position at B, indicating that the NLOS identification
is right and the error is migigated successfully. From then on, the position
results are always accurate.

We calculate the root-mean-square (RMS) and 90% location error, as shown in
Table I. For Case 3 and 4, we ignore results of the first lap since the filter
does not converge at the beginning. With the WLS-RKF, the position error is
reduced by over 95% in all cases compared to unweighted LS.

##  V Experiments and Results

###  V-A Experiment setup

We carry out experiments in the lobby of the laboratory building. As shown in
Fig. 3(a), there are 4 UWB anchors located at $A_{1}$(0.63, 0), $A_{2}$(7.57,
0), $A_{3}$(7.57, 5.40) and $A_{4}$(0.63, 5.40). An Unmanned Ground Vehicle
(UGV) with an UWB tag on the gimbal is controlled to move along a line at the
speed of about 0.5m/s. Each anchor or tag includes a DW1000 and an STM32F103
microcontroller. DW1000 is a low-power low-cost UWB transceiver. Through
double-sided two-way ranging (DS-TWR) [16], distances between the tag and the
anchors can be measured from the TOF of UWB signals, and calibrated with the
method in [13]. Preliminary experiments shows that the ranging interval is
about 0.05s and the standard deviation measurement noise is about 2cm.

(a)  (b)  (c)  (d)

Figure 4: Experiment results of LS, RKF and the proposed WLS-RKF. (a) Line
1-1. (b) Line 1-2. (c) Line 2-1. (d) Line 2-2.

We consider 2 cases, where the UGV moves along Line 1 and Line 2, as shown in
Fig. 3(b). Line 1 and 2 are parallel to the x-axis and y-axis respectively. In
the first case, a person stands inside the area as an obstacle. In the second
case, the UGV passes by a statue, which introduces a great NLOS bias.

At every sample time, we collect the measured distances and apply the proposed
WLS-RKF to get the position. All the parameters in the algorithm are the same
as in the simulations. For comparison, the positions with unweighted LS and
RKF are also calculated. The ground truth of the trajectory is measured by a
Sokkia SRX1 total station, which has an accuracy of 1mm.

###  V-B Experiment Results

In each case, we perform experiments where the UGV passes on both sides of the
obstacle, noted as Line 1-1, 1-2, 2-1 and 2-2. Fig. 4 shows the trajectories
of each algorithm. The proposed WLS-RKF achieves a good performance in all the
cases. In comparison, the error grows more slowly with RKF but still
accumulates as the UGV moves in NLOS condition. In Fig. 4 (a) and (d), the
filtered position is still biased after NLOS effect no longer exists.

Table II shows the RMS and 90% location error, and Fig. 5 shows the cumulative
distribution function (CDF). As the trajectories are roughly parallel to the
coordinate axes, we calculate the error in the y-axis for Line 1-1, 1-2, and
in the x-axis for Line 2-1, 2-2. Note that we only consider the section where
the NLOS effect occurs, that is, $x\in(5.5,7)$ for Line 1-1, $x\in(2,4)$ for
Line 1-2, and $x\in(0,2)$ for Line 2-1 and 2-2. We can see that WLS-RKF
mitigates the positioning error by over 90% for all the cases. The RMS error
is within 5cm.

TABLE II: Location Error(cm) for Experiments Algorithm | Line 1-1 | Line 1-2 | Line 2-1 | Line 2-2
---|---|---|---|---
|  RMS |  90% |  RMS |  90% |  RMS |  90% |  RMS |  90%
LS |  31.8 |  47.6 |  60.5 |  44.8 |  31.9 |  32.8 |  78.3 |  155.8
RKF |  20.6 |  37.5 |  26.1 |  50.1 |  14.7 |  26.2 |  94.6 |  177.8
WLS-RKF |  2.5 |  4.4 |  3.1 |  5.1 |  2.7 |  4.1 |  4.1 |  6.3

(a)  (b)  (c)  (d)

Figure 5: CDF of position error for LS, RKF and the proposed WLS-RKF (a) Line
1-1. (b) Line 1-2. (c) Line 2-1. (d) Line 2-2.

##  VI Conclusion

In this letter, we propose WLS-RKF, an NLOS mitigation method in UWB
positioning. It leverages the measured distances and mitigates NLOS error
using WLS technique and Kalman Filter without any prior knowledge about NLOS
distribution or signal features. The simulation and experiment results show
that WLS-RKF achieves good performances in various scenarios.

For future work, it may be possible for multiple sensor fusion. Besides UWB
positioning, other sensors such as inertial measurement unit (IMU) and
odometer can also be used to increase the robustness of the positioning
system.

## References

  * [1] D. B. Jourdan, D. Dardari, and M. Z. Win, “Position error bound for UWB localization in dense cluttered environments,” _IEEE transactions on aerospace and electronic systems_ , vol. 44, no. 2, pp. 613–628, 2008.
  * [2] J. Schroeder, S. Galler, K. Kyamakya, and K. Jobmann, “Nlos detection algorithms for ultra-wideband localization,” in _2007 4th Workshop on Positioning, Navigation and Communication_. IEEE, 2007, pp. 159–166.
  * [3] A. Musa, G. D. Nugraha, H. Han, D. Choi, S. Seo, and J. Kim, “A decision tree‐based NLOS detection method for the UWB indoor location tracking accuracy improvement,” _International Journal of Communication Systems_ , vol. 32, no. 13, p. e3997, 2019.
  * [4] İ. Güvenç, C.-C. Chong, F. Watanabe, and H. Inamura, “Nlos identification and weighted least-squares localization for uwb systems using multipath channel statistics,” _EURASIP Journal on Advances in Signal Processing_ , vol. 2008, pp. 1–14, 2007.
  * [5] L. Wang, R. Chen, L. Chen, L. Shen, P. Zhang, Y. Pan, and M. Li, “A robust filter for toa based indoor localization in mixed los/nlos environment,” in _2018 Ubiquitous Positioning, Indoor Navigation and Location-Based Services (UPINLBS)_. IEEE, 2018, pp. 1–9.
  * [6] P.-C. Chen, “A non-line-of-sight error mitigation algorithm in location estimation,” in _WCNC. 1999 IEEE Wireless Communications and Networking Conference (Cat. No. 99TH8466)_ , vol. 1. IEEE, 1999, pp. 316–320.
  * [7] L. Jiao, J. Xing, and F. Y. Li, “Performance comparison of residual related algorithms for toa positioning in wireless terrestrial and sensor networks,” in _2009 1st International Conference on Wireless Communication, Vehicular Technology, Information Theory and Aerospace & Electronic Systems Technology_. IEEE, 2009, pp. 278–283.
  * [8] D. Liu, M.-C. Lee, C.-M. Pun, and H. Liu, “Analysis of wireless localization in nonline-of-sight conditions,” _IEEE transactions on vehicular technology_ , vol. 62, no. 4, pp. 1484–1492, 2013.
  * [9] B. L. Le, K. Ahmed, and H. Tsuji, “Mobile location estimator with nlos mitigation using kalman filtering,” in _2003 IEEE Wireless Communications and Networking, 2003. WCNC 2003._ , vol. 3. IEEE, 2003, pp. 1969–1973.
  * [10] X. Li, Y. Wang, and K. Khoshelham, “Comparative analysis of robust extended kalman filter and incremental smoothing for uwb/pdr fusion positioning in nlos environments,” _Acta Geodaetica et Geophysica_ , vol. 54, no. 2, pp. 157–179, 2019.
  * [11] F. Gustafsson and F. Gunnarsson, “Mobile positioning using wireless networks: possibilities and fundamental limitations based on available wireless network measurements,” _IEEE Signal processing magazine_ , vol. 22, no. 4, pp. 41–53, 2005.
  * [12] B. Silva and G. P. Hancke, “Ranging error mitigation for through-the-wall non-line-of-sight conditions,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 11, pp. 6903–6911, 2020.
  * [13] J. González, J.-L. Blanco, C. Galindo, A. Ortiz-de Galisteo, J.-A. Fernández-Madrigal, F. A. Moreno, and J. L. Martinez, “Mobile robot localization based on ultra-wide-band ranging: A particle filter approach,” _Robotics and autonomous systems_ , vol. 57, no. 5, pp. 496–507, 2009.
  * [14] D. B. Jourdan, J. J. Deyst, M. Z. Win, and N. Roy, “Monte Carlo localization in dense multipath environments using UWB ranging,” in _2005 IEEE International Conference on Ultra-Wideband_. IEEE, 2005, pp. 314–319.
  * [15] G. Chang, “Kalman filter with both adaptivity and robustness,” _Journal of Process Control_ , vol. 24, no. 3, pp. 81–87, 2014.
  * [16] Decawave, “DW1000 user mannual how to use, configure and program the DW1000 UWB transceiver,” 2017.