4pt 4pt

# No More Potentially Dynamic Objects: Static Point Cloud Map Generation based
on 3D Object Detection and Ground Projection

Soojin Woo*, Donghwi Jung*, and Seong-Woo Kim * indicates equal contribution.
All authors are with Seoul National University, Seoul, South Korea.
<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose an algorithm to generate a static point cloud map
based on LiDAR point cloud data. Our proposed pipeline detects dynamic objects
using 3D object detectors and projects points of dynamic objects onto the
ground. Typically, point cloud data acquired in real-time serves as a snapshot
of the surrounding areas containing both static objects and dynamic objects.
The static objects include buildings and trees, otherwise, the dynamic objects
contain objects such as parked cars that change their position over time.
Removing dynamic objects from the point cloud map is crucial as they can
degrade the quality and localization accuracy of the map. To address this
issue, in this paper, we propose an algorithm that creates a map only
consisting of static objects. We apply a 3D object detection algorithm to the
point cloud data which are obtained from LiDAR to implement our pipeline. We
then stack the points to create the map after performing ground segmentation
and projection. As a result, not only we can eliminate currently dynamic
objects at the time of map generation but also potentially dynamic objects
such as parked vehicles. We validate the performance of our method using two
kinds of datasets collected on real roads: KITTI and our dataset. The result
demonstrates the capability of our proposal to create an accurate static map
excluding dynamic objects from input point clouds. Also, we verified the
improved performance of localization using a generated map based on our
method. Our code is available at https://github.com/woo-
soojin/no_more_potentially_dynamic_objects.

##  I Introduction

As the use of 3D-related technologies spreads, research related to detecting
3D space with sensors and representing it as 3D data is actively underway.
There are various methods for representing space in 3D, however, among them,
point cloud is a representation method that depicts space with countless
points consisting of three values: $x$, $y$, and $z$. Therefore, as the number
of points increases, the information about space becomes richer, and it
becomes possible to grasp more advanced spatial information, such as
calculating normal vectors through clustering with surrounding points, in
addition to basic coordinate information. Furthermore, as mentioned in
references [1, 2, 3, 4, 5], the point cloud is essential in various fields
including autonomous driving and robotics applications. Specifically, A point
cloud can be utilized for the representation of 3D objects, the creation of
meshes, a 3D scene reconstruction[1], and the digital twins that replicate
real-world spaces in virtual environments.[2] In autonomous driving and
robotics fields, point cloud data is also applied for various applications
such as map construction based on stacked point clouds [3], localization that
supports vehicle and robot to find their current position [4], and the path
planning [5]. For the path planning task, the point cloud can be adopted to
decide safe routes avoiding obstacles by applying an object detection
algorithm that identifies surrounding objects.

Figure 1: Changed distribution of potentially dynamic objects over time at the
same location.

Figure 2: The process of our static map creation method.

In addition, a point cloud map is also one of the common utilization of point
clouds. A point cloud map is generated by the accumulation of point clouds
from sensors such as LiDAR using the estimated pose based on the algorithm
[3]. This point cloud map represents the 3D space based on the distributions
of points containing information such as position and intensity. Among many
kinds of point cloud maps, the point cloud maps used for autonomous driving
are mostly generated from the data collected on the roads. Therefore, the data
is filled with dynamic objects such as vehicles and pedestrians. And the
distribution of these objects varies over time as shown in Fig. 1 by changing
their locations. The movement of dynamic objects not only degrades the
accuracy of scan-to-map matching in localization but also can lead to
inefficient path planning by mistakenly perceiving the areas as impassable
even though the place is actually empty. For these reasons, the representation
of surroundings only with static objects having fixed positions is important
for the accuracy of algorithms including localization and path planning. In
other words, dynamic objects in the 3D point cloud map must be identified and
eliminated.
Removing dynamic objects, as mentioned above, is a crucial issue in autonomous
driving and robotics. However, at the moment of map creation, it is difficult
to distinguish static objects and dynamic objects. In previous research that
removes dynamic objects from the map includes studies such as [6, 7, 8, 9, 10,
11, 12, 13, 14]. These studies differentiate dynamic objects from static
objects, and remove them predicting the dynamic objects that moved at the time
of data acquisition. Therefore, this approach can only identify and eliminate
objects that changed their positions at the time of data collection. However,
there is a limitation in that potentially dynamic objects, such as parked
cars, can not be recognized as dynamic objects. Consequently, they are not
removed from the map because these objects are identified as static objects
during the mapping process. As shown in Fig. 1, the potentially dynamic
objects that can change their locations over time are indistinguishable from
static objects such as buildings and trees during the data collecting process.
However, if the dynamic objects move, they might not be visible when the data
is acquired next time. Therefore, it is necessary to remove potentially
dynamic objects from the map, however this is a challenging task to resolve
with existing methods. This limitation not only decreases the quality of the
static map but also the performance of algorithms including localization and
path planning which rely on the static map.
In this paper, we propose a pipeline that detects objects based on their
features in a single frame rather than referring to its time-dependent
features and then removes dynamic objects by projecting them to the calculated
ground. Methods for detecting objects using their traits have been extensively
studied in various fields, including camera-based approaches [15, 16] and
LiDAR-based approaches [17, 18]. When using data collected by the camera,
objects are found in images, classified, and then 2D bounding boxes are
obtained. Considering the sparsity of LiDAR data, the methods based on LiDAR
data convert point cloud into various formats to improve the detection
performance. The [17] converts point cloud data into voxels for 3D object
detection. In [18], the input point cloud data is converted to a 2D Birds-Eye-
View (BEV) RGB-map to get the coordinates of the object’s bounding box.
Among these methods, we employed a voxel-based method for 3D object detection
which adopts LiDAR point cloud as input data. Using this approach, we detect
surrounding 3D objects and distinguish dynamic and static objects based on
returned classes on each object. The objects classified as dynamic are
projected to the ground plane calculated by the ground segmentation algorithm.
The datasets used in this paper were collected on the road which means, the
majority of dynamic objects such as cars exist over the road. Considering this
trait, we conducted ground projection on dynamic objects to replace the car
with the road as if there was no car on the road in the first place.
Subsequently, we create a static point cloud map with point cloud data without
dynamic objects. As a consequence, the static map constructed by our method
only includes static objects by removing both dynamic and potentially dynamic
objects from the map which makes our approach independent from data
acquisition timing.
To the best of our knowledge, our research is the very first to propose a
process that utilizes 3D object detection and ground projection algorithms on
a point cloud to generate a static point cloud map by removing potentially
dynamic objects. We confirmed the superior localization performance when we
used static maps created by our methods compared to maps generated by existing
approaches.
The contributions of our paper are as follows:

  * •

_A Static Map Creation Process_ : We proposed a process that connects 3D
object detection and ground projection to remove potentially dynamic objects
and generate a static point cloud map.

  * •

_Public Release of Code_ : We made the code used in our research publicly
available, enabling other researchers to utilize it.

Figure 3: The process of voxel-based 3D object detection.

##  II Related Works

There are methods for detecting and removing dynamic objects that exploit
occupancy maps. This approach uses ray tracing to evaluate the occupancy of a
specific area and determines the presence of dynamic objects by inspecting the
continuity of occupancy between consecutive frames. While this approach
demonstrates high accuracy in dynamic object removal, it imposes a significant
computational burden, particularly in occupancy calculations using ray
tracing. To address this computational burden, previous research, such as
Hornung _et al._ ’s work [6], attempted to enhance optimization. However, as
seen in the paper [7], computational challenges remained. Additionally, Pagad
_et al._ [8] improved performance by applying object detection to Octomap [6],
however computational burdens still persisted.
To mitigate the computational challenges posed by occupancy map-based methods,
visibility-based methods were proposed. Unlike occupancy map-based methods,
this approach focuses on the visibility of objects, reducing the computational
load. For instance, when a specific object is occluded by another object, it
is considered a dynamic point, reducing the computational requirements
compared to occupancy map-based methods. However, this simplification leads to
decreased accuracy, resulting in frequent false positives. To address this
issue, Kim _et al._ [9] introduced an iterative refinement process to remove
false positives. Nevertheless, challenges, such as less accuracy, still exist
in visibility-based methods. Therefore, Fan _et al._ [10] proposed a method
that integrates both visibility-based and map-based approaches to achieve high
accuracy and less computational burden in dynamic object removal. However,
their method is limited to detecting currently dynamic objects. In other
words, Dynamicfilter [10] could not detect the potentially dynamic objects
such as parked cars.
There are also model-based methods for removing dynamic objects using a model.
Among these, there is research that utilizes ground models including the paper
[11]. In this paper, the authors removed dynamic objects by assuming that all
dynamic objects are in contact with the ground. However, this research also
relies on the motion characteristics of dynamic objects to identify them,
potentially dynamic objects that may not change their positions over a short
period. Additionally, previous methods [12, 13] that adopted deep learning-
based approaches for dynamic object removal identified objects based on
features such as shape rather than their motion. Therefore, these method could
detect dynamic objects including potentially dynamic objects. However, these
methods projected the point cloud into 2D range images and identified dynamic
objects in these images. This projection selected one point per pixel,
resulting in information loss and reduced object detection performance. To
address this issue, in our paper, instead of projecting the point cloud into
2D range images, we use the shape of 3D voxels to ensure that more points are
kept. Moreover, Chen _et al._ [14] proposed a LiDAR SLAM algorithm using deep
learning-based dynamic object removal. In their study, they utilized the
PointPillars [19] to find dynamic objects in 3D point clouds without 2D
transformations. However, they removed points corresponding to dynamic
objects, which reduced accuracy of LiDAR SLAM due to sparsity of point clouds.
Therefore, in our paper, we improve mapping performance by utilizing ground
projection through ground segmentation, preserving points instead of removing
them.
In this paper, we develop a model-based approach capable of detecting both
currently dynamic objects and potentially dynamic objects. To achieve higher
accuracy in obtaining the 3D bounding box, we opt for a voxel-based object
detection method in our process. Subsequently, after projecting the points
associated with dynamic objects onto the ground, we generate a highly accurate
static map through mapping.

##  III Methods

In this section, we will describe the process, as depicted in Fig. 2, which
removes dynamic components from a given 3D point cloud map to create a static
point cloud map. First, we transform the point cloud data from the LiDAR
sensor into voxel features. We perform dynamic object detection by providing
transformed data as input to a CNN-based neural network. Additionally, we
apply a ground segmentation algorithm to project data to the ground after
calculating the ground plane based on raw point cloud data. Specifically, we
use 3D bounding boxes corresponding to dynamic objects returned by a 3D object
detector to find out the points inside the boxes to project them on the
ground. Finally, we stack these projected point clouds through mapping to
create the 3D static point cloud map. The following steps will provide
sequential explanations of our process.

###  III-A 3D Object Detection

There are several methods for 3D object detection, including the use of
mathematical equations. Among them, in this paper, we utilize the neural
network-based method that has been widely studied recently. The reason for
using the neural networks is that we can leverage GPUs to process incoming
point clouds in real-time with high accuracy. However, unlike data such as
camera image data, the point clouds coming from LiDAR do not have a consistent
size for each data sequence. Therefore, applying data embedding using neural
networks, such as convolution layers, is challenging. To overcome this
limitation, in this paper, we adopt a 3D object detection method that utilizes
spaced voxels based on subdivided 3D space. Voxels are data with fixed
dimensions in the shape of 3D cubes, allowing for embedding using 3D
convolution. As a result of this embedding, 2D features are obtained, and
applying 2D convolution layers ultimately allows us to predict 3D bounding
boxes. At this point, we set the prediction values through the network as the
center position, size, and orientation angle of the bounding box, ultimately
obtaining the 3D bounding box. Then, we use this obtained bounding box to
filter the point cloud. As a result, we can identify the points corresponding
to dynamic objects. The process of 3D object detection is described in Fig. 3
and formulated as follows:

$\displaystyle p=(x,y,z),$ (1) $\displaystyle P=\\{p\\},$ (2)
$\displaystyle\bar{P}=\\{\bar{p}|\theta^{min}\leq\bar{p}\leq\theta^{max}\\},$
(3) $\displaystyle\bar{p}\in
V_{ijk}\;(\text{if}\;\tau^{min}_{ijk}\leq\bar{p}\leq\tau^{max}_{ijk}),$ (4)
$\displaystyle(x_{v},y_{v},z_{v})=\left(x_{c}\pm\frac{w}{2},y_{c}\pm\frac{l}{2},z_{c}\pm\frac{h}{2}\right),$
(5) $\displaystyle R=\begin{pmatrix}\cos{\theta}&-\sin{\theta}&0\\\
\sin{\theta}&\cos{\theta}&0\\\ 0&0&1\end{pmatrix},$ (6)
$\displaystyle(\bar{x_{e}},\bar{y_{e}},\bar{z_{e}})^{\top}=R\cdot(x_{v},y_{v},z_{v})^{\top},$
(7)

where $P$ represents the incoming point clouds, and $p$ represents the
individual points composing the point clouds, with coordinates of $x$, $y$,
and $z$. In addition, $\theta^{min}$ and $\theta^{max}$ indicate the minimum
and maximum thresholds for filtering the point clouds in $x$, $y$, and $z$ for
voxelization. The points that are included in the voxel being passed through
the filtering represent $\bar{p}$, while $\bar{P}$ signifies the point clouds
composed of these $\bar{p}$ points. Additionally, $V_{ijk}$ corresponds to the
$i$th, $j$th, and $k$th voxel in the $x$, $y$, and $z$ directions,
respectively. $\tau^{min}_{ijk}$ and $\tau^{max}_{ijk}$ represent the
thresholds for the $x$, $y$, and $z$ values that determine the points included
in that specific voxel. Furthermore, $x_{c}$, $y_{c}$, and $z_{c}$ denote the
$x$, $y$, and $z$ values which are the center position of the bounding box.
Moreover, $w$, $l$, and $h$ indicate the width, length, and height of the 3D
bounding box. On the $x-y$ plane, $\theta$ is the orientation of the bounding
box. $x_{c}$, $y_{c}$, $z_{c}$, $w$, $l$ and $h$ are obtained as an output of
3D object detection model. On top of that, $x_{v}$, $y_{v}$, and $z_{v}$ mean
the $x$, $y$, and $z$ values expressing the positions of the bounding box’s
vertices. Additionally, $R$ indicates the rotation matrix based on $\theta$.
Lastly, $\bar{x}_{e}$, $\bar{y}_{e}$ and $\bar{z}_{e}$ are the transformed
positions of bounding box’s vertices.

###  III-B Ground Segmentation

Simply removing the points of dynamic objects consequently reduces the number
of points that are matched to the previously accumulated scans during the
scan-to-map matching stage, and it leads to the decreased accuracy of LiDAR
mapping by diminishing the number of points that can be useful for LiDAR
mapping. Therefore, we maintained the total number of points by projecting the
points corresponding to dynamic objects onto the ground plane in the middle of
our dynamic object removal process to sustain the performance of the LiDAR
mapping. To execute the projection, it is crucial to determine the equation on
the target plane, which in this case is the ground plane. Assuming that the
vehicle is on the ground and the slope is minimal during motion, the ground
plane will be parallel to the $x$-$y$ plane in the world coordinate system. As
a result, the normal vector of the ground plane will have values only in the
$z$ direction. To project dynamic object points onto the ground plane, the
average $z$ value of the ground plane denoted as $\bar{z}$ needs to be
calculated. In our paper, we achieve the ground plane calculation by
performing ground segmentation targeting the entire input point clouds,
isolating the points identified as ground, and then calculating the average
$z$ value of these points. This can be mathematically expressed as follows:

$\displaystyle P_{gr}=\\{p_{i}|i\in I_{gr}\\},$ (8) $\displaystyle
P=P_{ng}\cup P_{gr},$ (9)
$\displaystyle\bar{z}=\frac{1}{n}\cdot\Sigma^{n}_{k=1}P_{gr},$ (10)

where $P_{gr}$ represents the set of points in the point clouds that belong to
the ground, and $I_{gr}$ is the set of indices of points in $P$ which are part
of $P_{gr}$. On the other hand, $P_{ng}$ denotes the set of points in the
point clouds which are not included in the ground. Additionally, $n$ expresses
the number of ground points.

###  III-C Projection and Mapping

We utilize the bounding boxes obtained from 3D object detection process in
Sec. III-A to find out the point cloud and extract points corresponding to
dynamic objects. Subsequently, we use the mean $z$ value of the ground,
$\bar{z}$, calculated in Sec. III-B, to project the labeled points as dynamic
objects onto the ground. The appearance before and after projection is shown
in Fig. 4. The points corresponding to the dynamic objects are projected onto
the ground through ground segmentation, and then we create the final 3D point
cloud map using 3D LiDAR Simultaneous Localization and Mapping (SLAM).
Specifically, the data which is composed of projected points, and the map
which is accumulated until the previous step are used as input data of SLAM.
During this process, loop closure is used to correct prediction errors
resulting from long-distance driving. This procedure can be expressed
mathematically as follows:

$\displaystyle\hat{p_{i}}=\begin{cases}(x_{i},y_{i},z_{i}),&\text{if }i\in
I_{gr},\\\ (x_{i},y_{i},\bar{z}),&\text{otherwise},\end{cases},$ (11)
$\displaystyle\hat{P}=\\{\hat{p}\\},$ (12) $\displaystyle\hat{P}_{k+1}\oplus
M_{k}=M_{k+1},$ (13)

where $\hat{p}$ represents points that have been projected onto the ground
plane, and $\hat{P}$ denotes the point cloud composed of $\hat{p}$.
Additionally, $M$ represents the point cloud map. Furthermore, $\oplus$
signifies the mapping process utilizing LiDAR SLAM.

Figure 4: The projection result of the point cloud onto the ground plane
within the 3D bounding box obtained by the 3D object detection model. Dynamic
objects (pink), Static objects (gray).

##  IV Experiments

###  IV-A Implementation

In this paper, we utilized VoxelNeXt [20] for 3D object detection. For ground
segmentation, we applied the Patchwork++ algorithm [21]. Point cloud mapping
was performed using LOAM [3], while loop closure utilized Scan Context [22].
We assessed the performance of localization based on the accuracy of scan-to-
map matching using Iterative Closest Point (ICP) registration.

###  IV-B Data

We evaluated the performance of the proposed algorithm that generates the
static map both in a qualitative and quantitative way using datasets, such as
KITTI [23] and ours. Additionally, our dataset was acquired following the same
path on different dates to validate the performance of mapping and
localization according to the existence of dynamic objects on the same path.
For the evaluation process, we selected two types of datasets. One contained a
number of dynamic objects, while the other one had comparatively fewer dynamic
objects by applying a dynamic object removal algorithm. The reason for varying
the distribution of dynamic objects is that our research aims to verify
whether autonomous driving algorithms, such as mapping and localization, are
affected by changing distributions of dynamic objects over time. In the
mapping process, we utilized the dataset with abundant dynamic objects to
prove the ability of our algorithm that remove dynamic objects from the map.
For the localization task, we used the data that dynamic objects have been
removed to verify whether the algorithm performs well in different
environments, apart from the mapping stage.

###  IV-C Evaluation

We proved the capability of our algorithm in quantitative and qualitative
ways. In Sec. IV-C1, we qualitatively compared the results of mapping using a
map created by our algorithm and by other methods. Ground truth denotes the
map created using GPS information, and without removal expresses the map
constructed without applying a dynamic object removal algorithm.

For quantitative analysis, in Sec. IV-C2, we measured the accuracies of
localization on each static map generated by dynamic object removal algorithms
including our methods. As reference methods, we applied Removert[9] and
Erasor[11]. Additionally, we established the baseline by not applying any
dynamic object removal method for performance comparisons. For the
localization performance evaluation, we measured the accuracy of localization
at time $k+1$ by performing ICP registration between the 3D point cloud map
and the point cloud at time $k+1$. We used the ground truth pose at time $k$
as the initial transformation for global registration. Then, we calculated the
Root Mean Squared Error (RMSE) of the relative transformation between the
ground truth poses at time $k$ and time $k+1$ for the entire path. The results
are shown in Table II.

TABLE I: Comparison of SLAM odometry accuracy.
(Average error; meter). |  seq |  w/o removal |  Our method
---|---|---|---
KITTI[23] |  00 |  1.60 |  2.68
|  05 |  2.40 |  2.48
|  07 |  1.55 |  1.52
|  08 |  6.42 |  5.97
Our data |  |  3.56 |  4.33
TABLE II: Comparison of localization performance.
(Average error; meter, radian). |  |  w/o removal |  Our method
---|---|---|---
Our data |  x-y |  0.950 |  0.900
|  yaw |  0.070 |  0.066

####  IV-C1 Mapping

Fig. 5 and 6 show the results that applied our algorithm on KITTI and the
dataset of our campus. Firstly, we detected dynamic objects and then projected
points composing dynamic objects to the ground. Finally, we adapted LiDAR SLAM
for the mapping process to validate the capability of our proposal. As shown
in Table I, we compared the mapping results for each dataset with the results
of mapping using the raw point cloud before ground projection. By comparing
the mapping results by our method with the raw point cloud-based result, we
demonstrated that the performance of mapping is maintained even after ground
projection. As demonstrated in Fig. 5, 6 and Table I, we showed high accuracy
of mapping results for all two datasets including KITTI and our campus
dataset. Therefore, we confirmed that ground projection for dynamic objects
maintains mapping performance.

Figure 5: Comparison of SLAM odometry prediction performance on the 00
sequence of KITTI dataset.

Figure 6: Comparison of SLAM odometry prediction performance on our dataset.

Figure 7: Scan to map matching. Scan (yellow), Map (Cyan).

####  IV-C2 Localization

For localization, we calculated the performance using the accuracy of ICP-
based registration between the scan and map, as shown in Fig. 7. Table II
shows the localization error with the unit of meter. We demonstrated the
superior performance of our algorithm by showing lower localization errors
compared to the baseline (without removal). We achieved higher localization
accuracy by using our proposal to remove potentially dynamic objects from the
data used as input for the localization.

##  V Conclusion

In our paper, we proposed a novel and robust method for removing dynamic
objects from a given point cloud map through 3D object detection and ground
projection. The previous research relied on temporal information, therefore
they could only eliminate currently moving objects from their map. However,
our method not only removed moving objects at the time of data acquisition but
also potentially dynamic objects that change their location across time.
Rather than simply removing the points from the map that compose dynamic
objects, we projected corresponding points to the ground plane during the
mapping stage to maintain the performance of mapping and localization. In
particular, we found the dynamic objects based on the results of the 3D object
detection algorithm. To calculate the average height value of the ground
plane, we applied ground segmentation to the raw point clouds. With an
estimated value by the ground segmentation process, we projected the points
corresponding to the dynamic objects based on the results of the detection
algorithm. Subsequently, we used previously gained point cloud data by
projecting dynamic points as input of the LiDAR mapping process to create a
static map, dynamic objects are completely erased. The static map constructed
by our method does not possess both currently moving and potentially dynamic
objects. We evaluated the mapping performance of our proposal with two types
of datasets. One is publicly available and the other one was acquired by
ourselves. We also verified the localization performance using the generated
map with our method. Through evaluation, we confirmed the robustness of our
pipeline.

## Acknowledgement

This work is supported by the National Research Foundation of Korea (NRF)
through the Ministry of Science and ICT under Grant 2021R1A2C1093957, by the
Korean Ministry of Land, Infrastructure and Transport(MOLIT) as ”Innovative
Talent Education Program for Smart City”, and by Korea Institute for
Advancement of Technology(KIAT) grant funded by the Korea Government(MOTIE)
(P0020536, HRD Program for Industrial Innovation). The Institute of
Engineering Research at Seoul National University provided research facilities
for this work.

## References

  * [1] R. Wang, J. Peethambaran, and D. Chen, “Lidar point clouds to 3-d urban models $:$ a review,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 11, no. 2, pp. 606–627, 2018.
  * [2] L. Ladicky, O. Saurer, S. Jeong, F. Maninchedda, and M. Pollefeys, “From point clouds to mesh using regression,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 3893–3902.
  * [3] J. Zhang and S. Singh, “Loam: Lidar odometry and mapping in real-time.” in _Robotics: Science and systems_ , vol. 2, no. 9. Berkeley, CA, 2014, pp. 1–9.
  * [4] D. Rozenberszki and A. L. Majdik, “Lol: Lidar-only odometry and localization in 3d point cloud maps,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 4379–4385.
  * [5] Y. Peng, D. Qu, Y. Zhong, S. Xie, J. Luo, and J. Gu, “The obstacle detection and obstacle avoidance algorithm based on 2-d lidar,” in _2015 IEEE international conference on information and automation_. IEEE, 2015, pp. 1648–1653.
  * [6] A. Hornung, K. M. Wurm, M. Bennewitz, C. Stachniss, and W. Burgard, “Octomap: An efficient probabilistic 3d mapping framework based on octrees,” _Autonomous robots_ , vol. 34, pp. 189–206, 2013.
  * [7] J. Schauer and A. Nüchter, “The peopleremover—removing dynamic objects from 3-d point cloud data by traversing a voxel occupancy grid,” _IEEE robotics and automation letters_ , vol. 3, no. 3, pp. 1679–1686, 2018.
  * [8] S. Pagad, D. Agarwal, S. Narayanan, K. Rangan, H. Kim, and G. Yalla, “Robust method for removing dynamic objects from point clouds,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 10 765–10 771.
  * [9] G. Kim and A. Kim, “Remove, then revert: Static point cloud map construction using multiresolution range images,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 10 758–10 765.
  * [10] T. Fan, B. Shen, H. Chen, W. Zhang, and J. Pan, “Dynamicfilter: an online dynamic objects removal framework for highly dynamic environments,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 7988–7994.
  * [11] H. Lim, S. Hwang, and H. Myung, “Erasor: Egocentric ratio of pseudo occupancy-based dynamic object removal for static 3d point cloud map building,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 2272–2279, 2021.
  * [12] A. Milioto, I. Vizzo, J. Behley, and C. Stachniss, “Rangenet++: Fast and accurate lidar semantic segmentation,” in _2019 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2019, pp. 4213–4220.
  * [13] X. Chen, S. Li, B. Mersch, L. Wiesmann, J. Gall, J. Behley, and C. Stachniss, “Moving object segmentation in 3d lidar data: A learning-based approach exploiting sequential data,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 6529–6536, 2021.
  * [14] Y. Chen, S. Sun, H. Yin, and M. H. Ang, “Exploring the effect of 3d object removal using deep learning for lidar-based mapping and long-term vehicular localization,” in _2022 IEEE 25th International Conference on Intelligent Transportation Systems (ITSC)_. IEEE, 2022, pp. 1730–1735.
  * [15] R. Girshick, “Fast r-cnn,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 1440–1448.
  * [16] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, “You only look once: Unified, real-time object detection,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 779–788.
  * [17] Y. Zhou and O. Tuzel, “Voxelnet: End-to-end learning for point cloud based 3d object detection,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 4490–4499.
  * [18] M. Simony, S. Milzy, K. Amendey, and H.-M. Gross, “Complex-yolo: An euler-region-proposal for real-time 3d object detection on point clouds,” in _Proceedings of the European conference on computer vision (ECCV) workshops_ , 2018, pp. 0–0.
  * [19] A. H. Lang, S. Vora, H. Caesar, L. Zhou, J. Yang, and O. Beijbom, “Pointpillars: Fast encoders for object detection from point clouds,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2019, pp. 12 697–12 705.
  * [20] Y. Chen, J. Liu, X. Zhang, X. Qi, and J. Jia, “Voxelnext: Fully sparse voxelnet for 3d object detection and tracking,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2023, pp. 21 674–21 683.
  * [21] S. Lee, H. Lim, and H. Myung, “Patchwork++: Fast and robust ground segmentation solving partial under-segmentation using 3d point cloud,” 2022.
  * [22] G. Kim and A. Kim, “Scan context: Egocentric spatial descriptor for place recognition within 3d point cloud map,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 4802–4809.
  * [23] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for autonomous driving? the kitti vision benchmark suite,” in _2012 IEEE conference on computer vision and pattern recognition_. IEEE, 2012, pp. 3354–3361.


††footnotetext:  File: main.tex, printed: 2024-8-28, 1.35

# Landau-Bloch type theorem for elliptic and $K$-quasiregular harmonic
mappings

Vasudevarao Allu  Vasudevarao Allu, School of Basic Science, Indian Institute
of Technology Bhubaneswar, Bhubaneswar-752050, Odisha, India.
<EMAIL_ADDRESS>and  Rohit Kumar  Rohit Kumar, School of Basic Science,
Indian Institute of Technology Bhubaneswar, Bhubaneswar-752050, Odisha, India.
<EMAIL_ADDRESS>
###### Abstract.

In this paper, we establish an improved coefficient bounds for quasiregular
and elliptic harmonic mappings and using these bounds we establish Landau-
Bloch type theorem for $(K,K^{\prime})$-elliptic and K-quasiregular harmonic
mappings in plane. Furthermore, we prove the coefficient estimates for
$K$-quasiconformal harmonic self maps defined on the unit disk $\mathbb{D}$.

###### Key words and phrases:

Harmonic Mapping, Quasiregular mapping, $(K,K^{\prime})$-elliptic mapping,
Quasiconformal mapping, Landau-Bloch type theorem.

###### AMS Subject Classification: Mathematics Subject Classification:

Primary 30C50, 31A05; Secondary 32A18, 30C62, 30C45, 33E05.

##  1\. Introduction

A continuous complex-valued function $f=u+iv$ is a complex-valued harmonic
function in a domain $\Omega\subset\mathbb{C}$, if both $u$ and $v$ are real-
valued harmonic functions in $\Omega$ (see [7]). The inverse function theorem
and a result of Lewy [13] shows that a harmonic function $f$ is locally
univalent (i.e., one-to-one) in $\Omega$ if, and only if, the Jacobian $J_{f}$
is non-zero in $\Omega$, where

$J_{f}(z)=|f_{z}(z)|^{2}-|f_{\overline{z}}(z)|^{2}.$

A harmonic function $f$ is said to be sense-preserving if $J_{f}(z)>0$ for $z$
in $\Omega$. Let $\mathcal{H}$ denote the class of complex-valued harmonic
functions $f$ in the unit disk $\mathbb{D}=\\{z\in\mathbb{C}:|z|<1\\}$,
normalized by $f(0)=0=f_{z}(0)-1$. Each function $f$ in $\mathcal{H}$ can be
expressed as $f=h+\overline{g}$, where $h$ and $g$ are analytic functions in
$\mathbb{D}$ (see [7]). Here $h$ and $g$ are called the analytic and co-
analytic parts of $f$ respectively, and have power series representations

$h(z)=z+\sum_{n=2}^{\infty}a_{n}z^{n}\quad\mbox{and}\quad
g(z)=\sum_{n=1}^{\infty}b_{n}z^{n}.$

Let $\mathcal{S_{H}}$ be the subclass of $\mathcal{H}$ consisting of
univalent, and sense-preserving harmonic mappings on $\mathbb{D}$. Let
$\mathcal{H}^{0}=\\{f\in\mathcal{H}:f_{\overline{z}}(0)=0\\}$, and
$\mathcal{S_{H}}^{0}=\\{f\in\mathcal{S_{H}}:f_{\overline{z}}(0)=0\\}$. Hence
for any function $f=h+\overline{g}$ in $\mathcal{H}^{0}$, its analytic and co-
analytic parts can be represented by

$h(z)=z+\sum_{n=2}^{\infty}a_{n}z^{n}\quad\mbox{and}\quad
g(z)=\sum_{n=2}^{\infty}b_{n}z^{n}$

respectively. It is interesting to note that $\mathcal{S_{H}}$ reduces to
$\mathcal{S}$, the class of normalized analytic and univalent functions in
$\mathbb{D}$, if the co-analytic part of functions in the class
$\mathcal{S_{H}}$ is zero. In 1984, Clunie and Sheil-Small [6] investigated
the class $\mathcal{S_{H}}$, together with some geometric subclasses.
Subsequently, the class $\mathcal{S_{H}}$ and its subclasses have been
extensively studied by several authors (see [1, 2, 6, 11, 16]).

For a continuously differentiable function $f$, denote

$\Lambda_{f}=\max_{0\leq\theta\leq
2\pi}\left|f_{z}+e^{-2i\theta}f_{\bar{z}}\right|=\left|f_{z}\right|+\left|f_{\bar{z}}\right|,$
$\lambda_{f}=\min_{0\leq\theta\leq
2\pi}\left|f_{z}+e^{-2i\theta}f_{\bar{z}}\right|=||f_{z}|-|f_{\bar{z}}||.$

It is easy to see that $\left|J_{f}\right|=\Lambda_{f}\lambda_{f}$. A harmonic
mapping $f=h+\overline{g}$ defined in the unit disc $\mathbb{D}$ can be
expressed as

$f\left(re^{i\theta}\right)=\sum_{n=0}^{\infty}a_{n}r^{n}e^{in\theta}+\sum_{n=}^{\infty}\bar{b}_{n}r^{n}e^{-in\theta},\quad
0\leq r<1,$

where

$h(z)=\sum_{n=0}^{\infty}{a}_{n}z^{n},\ g(z)=\sum_{n=1}^{\infty}{b}_{n}z^{n}.$

The classical Landau theorem (see [12]) asserts that if $f$ is a holomorphic
mapping with $f(0)=0$, $f^{\prime}(0)=1$ and $|f(z)|<M$ for $z\in\mathbb{D}$,
then $f$ is univalent in $\mathbb{D}_{r_{0}}$, and $f(\mathbb{D}_{r_{0}})$
contains a disc $\mathbb{D}_{\sigma_{0}}$, where

$r_{0}=\frac{1}{M+\sqrt{M^{2}-1}}\ \ \text{and}\ \ \sigma_{0}=Mr_{0}^{2}.$

The quantities $r_{0}$ and $\sigma_{0}$ can not be improved. The extremal
function is $f_{0}(z)=Mz\left(\frac{1-Mz}{M-z}\right)$.

In 2000, Chen et al. [3] established the following two versions of Landau-type
theorems for bounded harmonic mapping on the unit disc under a suitable
restriction.

Theorem A. [3] Let $f$ be a harmonic mapping of the unit disc $\mathbb{D}$
such that $f(0)=0$, $f_{\bar{z}}(0)=0,f_{z}(0)=1$, and $|f(z)|<M$ for
$z\in\mathbb{D}$. Then, $f$ is univalent on a disc $\mathbb{D}_{\rho_{0}}$
with

$\rho_{0}=\frac{\pi^{2}}{16mM},$

and $f\left(\mathbb{D}_{\rho_{0}}\right)$ contains a schlicht disc
$\mathbb{D}_{R_{0}}$ with

$R_{0}=\rho_{0}/2=\frac{\pi^{2}}{32mM},$

where $m\approx 6.85$ is the minimum of the function $(3-r^{2})/(r(1-r^{2}))$
for $0<r<1$.

Theorem B. [3] Let $f$ be a harmonic mapping of the unit disc $\mathbb{D}$
such that $f(0)=0$, $\lambda_{f}(0)=1$ and $\Lambda_{f}(z)\leq\Lambda$ for
$z\in\mathbb{D}$. Then, $f$ is univalent on a disc $\mathbb{D}_{\rho_{0}}$
with

$\rho_{0}=\frac{\pi}{4(1+\Lambda)},$

and $f\left(\mathbb{D}_{\rho_{0}}\right)$ contains a schlicht disc
$\mathbb{D}_{R_{0}}$ with

$R_{0}=\frac{1}{2}\rho_{0}=\frac{\pi}{8(1+\Lambda)}.$

Theorems A and B are not sharp. In 2006, better estimates for theorem A and B
were given by Grigoyan [9]. In fact, Grigoyan has proved the following lemma.

Lemma C. [9] Assume that $f=h+\overline{g}$ with $h$ and $g$ analytic in the
unit disc $\mathbb{D}$ and let $h(z)=\sum_{n=1}^{\infty}a_{n}z^{n}$ and
$g(z)=\sum_{n=1}^{\infty}b_{n}z^{n}$ for $z\in\mathbb{D}$.

(a) If $|f(z)|<M$ for $z\in\mathbb{D}$, then

$\left|a_{n}\right|,\left|b_{n}\right|\leq M,\quad n=1,2,\ldots.$

(b) If $\Lambda_{f}(z)\leq\Lambda$ for $z\in\mathbb{D}$, then

$\left|a_{n}\right|+\left|b_{n}\right|\leq\frac{\Lambda}{n},\quad
n=1,2,\ldots.$

The following result has been proved by Grigoyan using Lemma C, which improves
the estimates in Theorems A and B.

Theorem D. [9] Let $f$ be harmonic mapping of the unit disc $\mathbb{D}$ such
that $f(0)=0,J_{f}(0)=1$ and $|f(z)|<M$ for $z\in\mathbb{D}$. Then, $f$ is
univalent on a disc $\mathbb{D}_{\rho_{1}}$ with

$\rho_{1}=1-\frac{2\sqrt{2}M}{\sqrt{\pi+8M^{2}}}$

and $f\left(\mathbb{D}_{\rho_{1}}\right)$ contains a schlicht disc with

$R_{1}=\frac{\pi}{4M}+4M-4M\sqrt{1+\frac{\pi}{8M^{2}}}.$

Theorem E. [9] Let $f$ be harmonic mapping of the unit disc $\mathbb{D}$ such
that $f(0)=0,\lambda_{f}(0)=1$ and $\Lambda_{f}(z)\leq\Lambda$ for all
$z\in\mathbb{D}$. Then, $f$ is univalent on a disc $\mathbb{D}_{\rho_{2}}$
with

$\rho_{2}=\rho_{2}(\Lambda)=\frac{1}{1+\Lambda}$

and $f\left(\mathbb{D}_{\rho_{2}}\right)$ contains a schlicht disc with

$R_{2}(\Lambda)=1-\Lambda\ln\left(1+\frac{1}{\Lambda}\right).$

######  Definition 1.1.

Let $\Omega\subset\mathbb{C}$ be a doamain. A mapping
$f:\Omega\rightarrow\mathbb{C}$ is said to be absolutely continuous on lines,
$ACL$ in brief, in a domain $\Omega$ if for every closed rectangle
$R\subset\Omega$ with sides parallel to the axes $x$ and $y,f$ is absolutely
continuous on almost every horizontal line and almost every vertical line in
$R$. Such a mapping $f$ has partial derivatives $f_{x}$ and $f_{y}$ a.e. in
$\Omega$. Moreover, we say $f\in ACL^{2}$ if $f\in ACL$ and its partial
derivatives are locally $L^{2}$ integrable in $\Omega$.

######  Definition 1.2.

[8] A sense-preserving and continuous mapping $f$ of $\mathbb{D}$ onto
$\mathbb{C}$ is said to be $(K,K^{\prime})$-elliptic mapping if

  1. (1)

$f$ is $ACL^{2}$ in $\mathbb{D}$, $J_{f}\neq 0$ a.e. in $\mathbb{D}$,

  2. (2)

$\exists\ K\geq 1$ and $K^{\prime}\geq 0$ such that

$\left\lVert D_{f}\right\rVert^{2}\leq KJ_{f}+K^{\prime}\ {\it{a.e.}}\
\text{in}\ \mathbb{D}.$

In particular, if $K^{\prime}=0$, then a $(K,K^{\prime})$-elliptic mapping
becomes $K$-quasiregular mapping. For more details on elliptic mappings, we
refer to ([5], [8], [14]).

It is easy to see that $\left\lVert D_{f}(z)\right\rVert=\Lambda(z)$.
Therefore, a sense-preserving harmonic mapping $f$ is said to be
$K$-quasiregular harmonic $(K\geq 1)$ on $\mathbb{D}$ if $\Lambda_{f}^{2}\leq
KJ_{f}$, or equivalently $\Lambda_{f}\leq K\lambda_{f}$ holds everywhere on
$\mathbb{D}$.

If $f:G\to G^{\prime}$ is $K$\- quasiregular harmonic homeomorphism then $f$
is called $K$-quasiconformal harmonic mapping.

For $\theta\in[0,2\pi]$ and $z\in\mathbb{D}$. Let

$P\left(z,e^{i\theta}\right)=\frac{1}{2\pi}\frac{1-|z|^{2}}{\left|1-ze^{-i\theta}\right|^{2}}$

be the Poisson kernel. For a mapping $F\in L^{1}(\partial\mathbb{D})$, the
Poisson integral of $F$ is defined by

$f(z)=P[F](z)=\int_{0}^{2\pi}P\left(z,e^{i\theta}\right)F\left(e^{i\theta}\right)d\theta.$

In 2007, Partyka and Sakan [15] proved following result.

Theorem F. [15] Given $K\geq 1$, and let $f(z)=P[F](z)$ be a harmonic
$K$-quasiconformal mapping of $\mathbb{D}$ onto itself, with the boundary
function $F(t)$. If $f(0)=0$, then for a.e. $z=e^{it}\in\partial\mathbb{D}$,

$\frac{2^{5\left(1-K^{2}\right)/2}}{\left(K^{2}+K-1\right)^{K}}\leq\left|F^{\prime}(t)\right|\leq
K^{3K}2^{5(K-1/K)/2}.$

In 2020, Chen and Ponnusamy [4] proved the following lemma which characterizes
elliptic mappings.

Lemma G. [4] Let $f\in\mathrm{C}^{1}$ be a sense-preserving mapping. Then $f$
is a elliptic mapping if, and only if, there exist constants $k_{1}\in[0,1)$
and $k_{2}\in[0,\infty)$ such that

$\left|f_{\bar{z}}(z)\right|\leq k_{1}\left|f_{z}(z)\right|+k_{2}$

for $z\in\mathbb{D}.$
In particular, if $f$ is $(K,K^{\prime})$-elliptic mapping then

$k_{1}=\frac{K-1}{K+1},\ \ k_{2}=\frac{\sqrt{K^{\prime}}}{1+K}.$

The main aim of this paper is to establish new versions of Landau-type
theorems for elliptic and $K$-quasiregular harmonic mappings and some
coefficient estimates for $K$-quasiconformal harmonic mappings.

##  2\. Main results

We first prove coefficient estimates for elliptic harmonic mappings and
$K$-quasiregular harmonic mappings.

######  Theorem 2.1.

Let $f=h+\bar{g}$ be $(K,K^{\prime})$-elliptic harmonic mapping defined on
unit disc $\mathbb{D}$, where

(2.2) $h(z)=\sum_{n=1}^{\infty}a_{n}z^{n}\quad\text{and}\quad
g(z)=\sum_{n=1}^{\infty}b_{n}z^{n}.$

If $\lambda_{f}(z)\leq\lambda$ for $z\in\mathbb{D}$, then

(2.3)
$\left|a_{n}\right|+\left|b_{n}\right|\leq\frac{K\lambda+\sqrt{K^{\prime}}}{n}.$

We see that for the case $K=1$ and $K^{\prime}=0$, we get

$\left|a_{n}\right|+\left|b_{n}\right|\leq\frac{\lambda}{n}.$

Using Theorem 2.1, we obtain following corollary.

######  Corollary 2.4.

Let $f=h+\bar{g}$ be a harmonic $K$-quasiregular mapping defined on unit disk
$\mathbb{D}$, where

$h(z)=\sum_{n=1}^{\infty}a_{n}z^{n}\quad\text{and}\quad
g(z)=\sum_{n=1}^{\infty}b_{n}z^{n}.$

If $\lambda_{f}(z)\leq\lambda$ for $z\in\mathbb{D}$, then

$\left|a_{n}\right|+\left|b_{n}\right|\leq\frac{K\lambda}{n}.$

Next, we state Landau-type theorems for elliptic harmonic mappings and
$K$-quasiregular mappings.

######  Theorem 2.5.

Let $f=h+\bar{g}$ be $(K,K^{\prime})$-elliptic harmonic mapping defined on
unit disc $\mathbb{D}$ such that $f(0)=0$, $\lambda_{f}(0)=1$ and
$\lambda_{f}(z)\leq\lambda$ for all $z\in\mathbb{D}$. Then $f$ is univalent on
a disk $\mathbb{D}_{\rho_{1}}$ with

$\rho_{1}=\rho_{1}(\lambda)=\frac{1}{1+K\lambda+\sqrt{K^{\prime}}}$

and $f\left(\mathbb{D}_{\rho_{1}}\right)$ contains a schlicht disc with radius

$R_{1}=\rho_{1}(\lambda)+(K\lambda+\sqrt{K^{\prime}})(\rho_{1}(\lambda)+\ln\left(\left(K\lambda+\sqrt{K^{\prime}}\right)\rho_{1}(\lambda)\right).$

######  Corollary 2.6.

Let $f=h+\bar{g}$ be $K$-quasiregular harmonic mapping defined on unit disk
$\mathbb{D}$ such that $f(0)=0$, $\lambda_{f}(0)=1$ and
$\lambda_{f}(z)\leq\lambda$ for all $z\in\mathbb{D}$. Then $f$ is univalent on
a disc $\mathbb{D}_{\rho_{1}}$ with

$\rho_{1}=\rho_{1}(\lambda)=\frac{1}{1+K\lambda}$

and $f\left(\mathbb{D}_{\rho_{1}}\right)$ contains a schlicht disc with radius

$R_{1}=\rho_{1}(\lambda)+K\lambda\left(\rho_{1}(\lambda)+\ln\left(K\lambda\
\rho_{1}(\lambda)\right)\right).$

Using Corollary 2.4, we proved the following theorem.

######  Theorem 2.7.

Let $f=h+\bar{g}$ be a harmonic K-quasiregular mapping defined on unit disc
$\mathbb{D}$ such that $f(0)=0,J_{f}(0)=1$ and $\lambda_{f}{(z)}\leq\lambda$
for all $z\in\mathbb{D}$. Then $f$ is univalent on a disc $\mathbb{D}_{\rho}$
with

$\rho(\lambda)=\frac{1}{1+\lambda K^{3/2}}$

and $f\left(\mathbb{D}_{\rho}\right)$ contains a schlicht disc with radius

$R(\lambda)=\frac{\rho(\lambda)}{\sqrt{K}}+K\lambda\left(\rho(\lambda)+\ln(\lambda
K^{3/2}\rho(\lambda)\right).$

Finally, using Theorem G, we prove the coefficient bounds for
$K$-quasiconformal self mappings defined on the unit disc $\mathbb{D}$.

######  Theorem 2.8.

Let $f=h+\bar{g}$ be a $K$-quasiconformal harmonic self mapping on the unit
disk $\mathbb{D}$, where $h$ and $g$ are analytic in $\mathbb{D}$ and given by

$h(z)=\sum_{n=1}^{\infty}a_{n}z^{n}\quad\text{and}\quad
g(z)=\sum_{n=1}^{\infty}b_{n}z^{n}$

for $z\in\mathbb{D}$. Then,

$\left|a_{n}\right|,\left|b_{n}\right|\leq\frac{K^{3K}2^{5\left(K-\frac{1}{K}\right)/2}}{n}.$

##  3\. Proof of the main results

Proof of Theorem 2.1. Let $f=h+\bar{g}$ be $(K,K^{\prime})$\- elliptic
harmonic mapping defined on the unit disc $\mathbb{D}$, where $h$ and $g$ are
given by (2.2). For fixed $r\in(0,1)$ and $z=re^{i\theta}$ for some
$\theta\in[0,2\pi]$, we have

(3.1) $\frac{\partial h}{\partial
r}(re^{i\theta})=\sum_{n=1}^{\infty}na_{n}r^{n-1}e^{in\theta}\ \text{ and }\
\frac{\partial g}{\partial
r}(re^{i\theta})=\sum_{n=1}^{\infty}nb_{n}r^{n-1}e^{in\theta}.$

Integrating both sides of (3.1) with respect to $r$ from $0$ to $2\pi$ after
multiplying $e^{-in\theta}$ to the left equation and multiplying
$e^{in\theta}$ to the right equation, we obtain

$na_{n}r^{n-1}=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\partial
h(re^{i\theta})}{\partial r}e^{-in\theta}d\theta$

and

$nb_{n}r^{n-1}=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\partial
g(re^{i\theta})}{\partial r}e^{in\theta}d\theta.$

for $n=1,2...$. Therefore, we have

$n\left|a_{n}\right|r^{n-1}\leq\frac{1}{2\pi}\int_{0}^{2\pi}\left|{\frac{\partial}{\partial
r}}h(re^{i\theta})\right|d\theta$

and

$n\left|b_{n}\right|r^{n-1}\leq\frac{1}{2\pi}\int_{0}^{2\pi}\left|{\frac{\partial}{\partial
r}}g(re^{i\theta})\right|d\theta.$

Hence, we get

$n\left(\left|a_{n}\right|+\left|b_{n}\right|\right)r^{n-1}\leq\frac{1}{2\pi}\int_{0}^{2\pi}\left(\left|\frac{\partial}{\partial
r}h(re^{i\theta})\right|+\left|\frac{\partial}{\partial
r}g(re^{i\theta})\right|\right)d\theta$

It is easy to see that

$\displaystyle\frac{\partial}{\partial r}h(re^{i\theta})$
$\displaystyle=\frac{\partial h}{\partial z}\frac{\partial z}{\partial
r}+\frac{\partial h}{\partial\bar{z}}\frac{\partial\bar{z}}{\partial r}$
$\displaystyle=h_{z}e^{i\theta}+0\quad\quad\ \ \left(\because
h_{\bar{z}}=0\right)$ $\displaystyle=h_{z}e^{i\theta}$

and similarly,

$\frac{\partial}{\partial r}g(re^{i\theta})=h_{z}e^{i\theta}.$

We also see that

(3.2) $\left|f_{z}\right|=\left|h_{z}\right|\text{ and
}\left|f_{\bar{z}}\right|=\left|g_{z}\right|.$

In view of (3.2), we have

(3.3)
$n\left(\left|a_{n}\right|+\left|b_{n}\right|\right)r^{n-1}\leq\frac{1}{2\pi}\int_{0}^{2\pi}\left(\left|h_{z}e^{i\theta}\right|+\left|g_{z}e^{i\theta}\right|\right)d\theta=\frac{1}{2\pi}\int_{0}^{2\pi}\left(\left|f_{z}\right|+\left|f_{\bar{z}}\right|\right)d\theta.$

Since $\lambda_{f}(z)\leq\lambda$ and $f$ is sense preserving, we have

$\left|f_{z}\right|\leq\lambda+\left|f_{\bar{z}}\right|.$

Using Lemma H, we obtain

$\left|f_{z}\right|\leq\lambda+k_{1}\left|f_{z}\right|+k_{2},$

which implies that

$\left|f_{z}\right|\leq\frac{\lambda+k_{2}}{1-k_{1}}.$

Therefore, we have

$\displaystyle|f_{z}|+|f_{\bar{z}}|$
$\displaystyle\leq|f_{z}|+k_{1}|f_{z}|+k_{2}$
$\displaystyle=\left(1+k_{1}\right)\left|f_{z}\right|+k_{2}$
$\displaystyle\leq\left(1+k_{1}\right)\frac{\lambda+k_{2}}{1-k_{1}}+k_{2}.$

In view of Lemma H, we have

(3.4) $|f_{z}|+|f_{\bar{z}}|\leq K\lambda+\sqrt{K^{\prime}}.$

By equation (3.3) and equation (3.4), we obtain

$\displaystyle n\left(\left|a_{n}\right|+\left|b_{n}\right|\right)r^{n-1}$
$\displaystyle\leq\frac{1}{2\pi}\int_{0}^{2\pi}\left(K\lambda+\sqrt{K^{\prime}}\right)d\theta$
$\displaystyle=K\lambda+\sqrt{K^{\prime}}.$

Now by taking $r\rightarrow 1^{-}$, we obtain

$\left|a_{n}\right|+\left|b_{n}\right|\leq\frac{K\lambda+\sqrt{K^{\prime}}}{n}.$

This completes the theorem. ∎

###### Proof of Corollary 2.4.

Let $f$ is $K$-quasiregular harmonic mappiing, then $f$ is a $(K,0)$-elliptic
harmonic mapping. By substituting $K^{\prime}=0$ in (2.3) we obtain

$\left|a_{n}\right|+\left|b_{n}\right|\leq\frac{K\lambda}{n}.$

∎

###### Proof of Theorem 2.5.

Let $f=h+\bar{g}$ be $(K,K^{\prime})$\- elliptic harmonic mapping defined on
the unit disc $\mathbb{D}$, where $h$ and $g$ are given by (2.2). Then

$f(z)=\sum_{n=1}^{\infty}a_{n}z^{n}+\sum_{n=1}^{\infty}\bar{b}_{n}\bar{z}^{n}\quad\text{for}\
z\in\mathbb{D}.$

Fix $r\in(0,1)$ and $z_{1},z_{2}\in\overline{\mathbb{D}}_{r}$ with $z_{1}\neq
z_{2}$. It is easy to see that

$f\left(z_{2}\right)-f\left(z_{1}\right)=h\left(z_{2}\right)-h\left(z_{1}\right)+\overline{g\left(z_{2}\right)-g\left(z_{1}\right)}.$

Therefore, we have

$|f(z_{2})-f(z_{1})|=\left|\int_{[z_{1},z_{2}]}h^{\prime}(z)dz+\overline{\left(\int_{[z_{1},z_{2}]}g^{\prime}(z)dz\right)}\right|,$

where $[z_{1},z_{2}]$ is the line segment joining $z_{1}$ and $z_{2}$. It is
easy to see that

$\displaystyle\begin{split}\left|f\left(z_{2}\right)-f\left(z_{1}\right)\right|&\geq\left|\int_{\left[z_{1},z_{2}\right]}\left(h^{\prime}(z)-h^{\prime}(0)+h^{\prime}(0)\right)dz+\overline{\int_{\left[z_{1},z_{2}\right]}\left(g^{\prime}(z)-g^{\prime}(0)+g^{\prime}(0)\right)dz}\right|\\\
&\geq\left|\int_{[z_{1},z_{2}]}\left(h^{\prime}(z)-h^{\prime}(0)\right)dz+\int_{\left[z_{1},z_{2}\right]}h^{\prime}(0)dz+\overline{\int_{[z_{1},z_{2}]}g^{\prime}(0)dz}\right.\\\
&\ \
\quad+\left.\overline{\int_{\left[z_{1},z_{2}\right]}\left(g^{\prime}(z)-g^{\prime}(0)\right)dz}\right|\\\
&\geq\left|\int_{\left[z_{1},z_{2}\right]}h^{\prime}(0)dz+\overline{\int_{z_{1},z_{2}]}g^{\prime}(0)dz}\right|\\\
&\ \
\quad-\left|\int_{\left[z_{1},z_{2}\right]}\left(h^{\prime}(z)-h^{\prime}(0)\right)dz+\overline{\int_{[z_{1},z_{2}]}\left(g^{\prime}(z)-g^{\prime}(0)\right)}dz\right|\end{split}$

Further,we have

$\displaystyle\left|f\left(z_{2}\right)-f\left(z_{1}\right)\right|$
$\displaystyle\geq\lambda_{f}(0)\left|z_{2}-z_{1}\right|-\left|\int_{\left[z_{1},z_{2}\right]}h^{\prime}\left(z\right)dz-\int_{\left[z_{1},z_{1}\right]}h^{\prime}(0)dz\right|$
$\displaystyle\ \
-\left|\int_{\left[z_{1},z_{2}\right]}g^{\prime}(z)dz-\int_{\left[z_{1},z_{2}\right]}g^{\prime}(0)dz\right|$
$\displaystyle=\lambda_{f}(0)\left|z_{2}-z_{1}\right|-\left|h\left(z_{2}\right)-h\left(z_{1}\right)-h^{\prime}(0)\left(z_{2}-z_{1}\right)\right|$
$\displaystyle\ \ \
-\left|g\left(z_{2}\right)-g\left(z_{1}\right)-g^{\prime}(0)\left(z_{2}-z_{1}\right)\right|$
$\displaystyle=|z_{2}-z_{1}|\left(1-\left|\frac{h\left(z_{2}\right)-h\left(z_{1}\right)}{z_{2}-z_{1}}-h^{\prime}(0)\right|-\left|\frac{g\left(z_{2}\right)-g\left(z_{1}\right)}{z_{2}-z_{1}}-g^{\prime}(0)\right|\right).$

Since $z_{1},z_{2}\in\mathbb{\overline{D}}_{r}$, we obtain

$\left|\frac{h\left(z_{2}\right)-h\left(z_{1}\right)}{z_{2}-z_{1}}-h^{\prime}(0)\right|\leq\sum_{n=2}^{\infty}\left|a_{n}\right|nr^{n-1}$

and

$\left|\frac{g\left(z_{2}\right)-g\left(z_{1}\right)}{z_{2}-z_{1}}-g^{\prime}(0)\right|\leq\sum_{n=2}^{\infty}\left|b_{n}\right|nr^{n-1}.$

Hence

(3.5)
$|f(z_{2})-f(z_{1})|\geq\left|z_{2}-z_{1}\right|\left(1-\sum_{n=2}^{\infty}\left(\left|a_{n}\right|+\left|b_{n}\right|\right)nr^{n-1}\right).$

Using (3.5) and Theorem 2.1, we obtain

$\left|\frac{f(z_{2})-f(z_{1})}{z_{2}-z_{1}}\right|\geq
1-\left(K\lambda+\sqrt{K^{\prime}}\right)\sum_{n=2}^{\infty}r^{n-1}=1-\left(K\lambda+\sqrt{K^{\prime}}\right)\frac{r}{1-r}.$

Since $z_{2}\neq z_{1}$, so $f\left(z_{2}\right)\neq f\left(z_{1}\right)$ only
if

$1-\left(K\lambda+\sqrt{K^{\prime}}\right)\frac{r}{1-r}>0$

which implies that

$r<\frac{1}{1+K\lambda+\sqrt{K^{\prime}}}=\rho_{1}.$

Hence, $f$ is univalent in the disc $\mathbb{D}_{\rho_{1}}.$ This completes
the proof of the first part of the theorem.

To prove the second part of the theorem, let $|z|=\rho_{1}$, by Theorem 2.1,
we obtain

$\displaystyle|f(z)|$
$\displaystyle=\left|\sum_{n=1}^{\infty}a_{n}z^{n}+\sum_{n=1}^{\infty}\bar{b}_{n}\bar{z}^{n}\right|$
$\displaystyle=\left|a_{1}z+\bar{b}_{1}\bar{z}^{\cdot}+\sum_{n=2}^{\infty}a_{n}z^{n}+\sum_{n=2}^{\infty}\bar{b}_{n}\bar{z}^{n}\right|$
$\displaystyle\geq\left|a_{1}z+\bar{b}_{1}\bar{z}\right|-\left|\sum_{n=2}^{\infty}a_{n}z^{n}+\sum_{n=2}^{\infty}\bar{b}_{n}\bar{z}^{n}\right|$
$\displaystyle\geq\left|\left|a_{1}\right|-\left|b_{1}\right|\right||z|-\sum_{n=2}^{\infty}\left(\left|a_{n}\right|+\left|b_{n}\right|\right)\left|z^{n}\right|$
$\displaystyle=\lambda_{f}(0)|z|-\sum_{n=2}^{\infty}\left(\left|a_{n}\right|+\left|b_{n}\right|\right)\left|z^{n}\right|$
$\displaystyle\geq\rho_{1}-(K\lambda+\sqrt{K^{\prime}})\sum_{n=2}^{\infty}\frac{\rho_{1}^{n}}{n}$
$\displaystyle=\rho_{1}+\left(K\lambda+\sqrt{K^{\prime}}\right)\left(\rho_{1}+\ln\left(1-\rho_{1}\right)\right):=R_{1}.$

Hence $f\left(\mathbb{D}_{\rho_{1}}\right)$ contains a schlicht disc of radius
$R_{1}$. This completes the proof.

∎

By taking $K^{\prime}=0$ in Theorem 2.5, we obtain the proof of Corollary 2.6.

###### Proof of Theorem 2.7.

Since $J_{f}(0)=1$, we have

$\displaystyle 1=|f_{z}(0)|^{2}-|f_{\bar{z}}(0)|^{2}$
$\displaystyle=(|f_{z}(0)|+|f_{\bar{z}}(0)|)(|f_{z}(0)|-|f_{\bar{z}}(0)|)$
$\displaystyle\leq K\left(\lambda_{f}(0)\right)^{2}$

which implies that

$\lambda_{f}(0)\geq\frac{1}{\sqrt{K}}.$

Fix $r\in(0,1)$ and $z_{1},z_{2}\in\mathbb{\overline{D}}_{r}$ with $z_{1}\neq
z_{2}$. Then by proof of Theorem 2.5 and Corollary 2.4 we obtain

$\displaystyle\left|\frac{f(z_{2})-f(z_{1})}{z_{2}-z_{1}}\right|$
$\displaystyle\geq\lambda_{f}(0)-\sum_{n=2}^{\infty}\left(|a_{n}|+|b_{n}|\right)nr^{n-1}$
$\displaystyle\geq\frac{1}{\sqrt{K}}-K\lambda\sum_{n=2}^{\infty}r^{n-1}$
$\displaystyle=\frac{1}{\sqrt{K}}-\frac{K\lambda r}{1-r}.$

Since $z_{1}\neq z_{2}$, so $f\left(z_{2}\right)\neq f\left(z_{2}\right)$ only
if

$\frac{1}{\sqrt{k}}-\frac{k\lambda r}{1-r}>0,$

which implies that

$r<\frac{1}{1+\lambda k^{3/2}}=\rho(\lambda).$

Hence, $f\left(z_{2}\right)\neq f\left(z_{1}\right)$ for $r\in(0,\rho)$. This
shows that $f$ is univalent in $\mathbb{D}\rho$.
For $|z|=\rho$, by using the proof of the Theorem 2.5, we have

$\displaystyle|f(z)|$
$\displaystyle\geq\lambda_{f}(0)|z|-\sum_{n=2}^{\infty}\left(|a_{n}|+|b_{n}|\right)|z^{n}|$
$\displaystyle\geq\frac{\rho}{\sqrt{K}}-\sum_{n=2}^{\infty}\frac{K\lambda}{n}\rho^{n}$
$\displaystyle=\frac{\rho}{\sqrt{K}}+K\lambda(\rho+\ln(1-\rho))$
$\displaystyle=\frac{\rho}{\sqrt{K}}+K\lambda\left(\rho+\ln(\lambda
K^{3/2}\rho)\right)=R(\lambda).$

Hence, $f(\mathbb{D}_{\rho})$ contains a schlicht disc of radius $R(\lambda).$
This completes the proof. ∎

###### Proof of Theorem 2.8.

Let $f=h+\bar{g}$ be a $K$-quasiconformal harmonic self mapping on the unit
disk $\mathbb{D}$, where $h$ and $g$ are analytic in $\mathbb{D}$ and given by

$h(z)=\sum_{n=1}^{\infty}a_{n}z^{n}\quad\text{and}\quad
g(z)=\sum_{n=1}^{\infty}b_{n}z^{n}.$

Fix $r\in(0,1)$, then for $z=re^{i\theta}$, we have

$f(re^{i\theta})=\sum_{n=1}^{\infty}a_{n}r^{n}e^{in\theta}+\sum_{n=1}^{\infty}\bar{b}_{n}r^{n}e^{-in\theta}.$

Multiplying both the sides by $e^{-in\theta}$ and integrating with respect to
$\theta$ from 0 to $2\pi$, we get

(3.6)
$a_{n}r^{n}=\frac{1}{2\pi}\int_{0}^{2\pi}f\left(re^{i\theta}\right)e^{-in\theta}d\theta.$

Similarly,

(3.7)
$\bar{b}_{n}r^{n}=\frac{1}{2\pi}\int_{0}^{2\pi}f\left(re^{i\theta}\right)e^{in\theta}d\theta.$

It is easy to see that an integration by part of (3.6) gives

$a_{n}r^{n}=\frac{1}{2\pi
in}\int_{0}^{2\pi}f_{\theta}\left(re^{i\theta}\right)e^{-in\theta}d\theta.$

which implies that

(3.8) $|a_{n}|r^{n}\leq\frac{1}{2\pi
n}\int_{0}^{2\pi}\left|f_{\theta}(re^{i\theta})\right|d\theta.$

Similarly, we see that an integration by part of (3.7) gives

(3.9) $|b_{n}|r^{n}\leq\frac{1}{2\pi
n}\int_{0}^{2\pi}\left|f_{\theta}(re^{i\theta})\right|d\theta.$

In 2008, Kalaj [10] proved that the radial limits of $f_{\theta}$ and $f_{r}$
exists almost everywhere and

$\lim_{r\to 1^{-}}f_{\theta}\left(re^{i\theta}\right)=f^{\prime}(\theta).$

Therefore taking $r\to 1^{-}$ and by Theorem F in (3.8), we get

$\displaystyle|a_{n}|$
$\displaystyle\leq\frac{1}{2\pi}\int_{0}^{2\pi}\left|f^{\prime}(\theta)\right|d\theta$
$\displaystyle\leq\frac{K^{3K}2^{5(K-1/K)/2}}{n}.$

Similarly, by letting $r\to 1^{-}$ and by Theorem F in (3.9), we obtain

$|b_{n}|\leq\frac{K^{3K}2^{5(K-1/K)/2}}{n}.$

This completes the proof. ∎

Acknowledgement: The first author thank SERB-MATRICS and the second author
thank CSIR for their support.

## References

  * [1] D. Bshouty, S. B. Joshi and S. S. Joshi, On close-to-convex harmonic mappings, Complex Var. Elliptic Equ. 58 (2013), 1195–1199.
  * [2] D. Bshouty and A. Lyzzaik, Close-to-convexity criteria for planar harmonic mappings, Complex Anal. Oper. Theory 5 (2011), 767–774.
  * [3] S. Chen, P.M. Gauthier and W. Hengartner, Bloch constants for planer harmonic mappings, Proc. Amer. Math. Soc. 128 (2000), 3231–3240.
  * [4] S. Chen  and S. Ponnusamy,On certain quasiconformal and elliptic mappings, J. Math. Anal. Appl. 486 (2020), 1–16.
  * [5] S. Chen, S. Ponnusamy and X. Wang, Remarks on Norm Estimates of the Partial Derivatives for Harmonic Mappings and Harmonic Quasiregular Mappings, J. Goem. Anal. 31 (2021), 11051–11060.
  * [6] J.G. Clunie and T. Sheil-Small, Harmonic univalent functions, Ann. Acad. Sci. Fenn. Ser. A. I. 9 (1984), 3–25.
  * [7] P.L. Duren, Harmonic mapping in the plan, Cambridge University Press, (2004).
  * [8] R. Finn and J. Serrin, On the Hölder continuity of quasiconformal and elliptic mappings, Trans. Amer. Math. Soc. 89 (1958), 1–15.
  * [9] A. Grigoyan, Landau and Bloch theorems for harmonic mappings, Complex Variable Theory Appl. 51(1) (2006), 81–87.
  * [10] D. Kalaj, Quasiconformal and harmonic mappings between Jordan domains, Math. Z. 260 (2008), 237–252.
  * [11] D. Kalaj, S. Ponnusamy and M. Vuorinen, Radius of close-to-convexity and full starlikeness of harmonic mappings, Complex Var. Elliptic Equ. 59 (2014), 539–552.
  * [12] E. Landau, Der Picard-Schottkysche Satz und die Blochsche Konstanten, Sitzungsber. Preuss. Akad. Wiss. Berlin Phys.-Math. Kl. (1926), 467–474.
  * [13] H. Lewy, On the non-vanishing of the Jacobian in certain one-to-one mappings, Bull. Amer. Math. Soc. 42 (1936), 689–692.
  * [14] L. Nirenberg, On nonlinear elliptic partial differential equations and Hölder continuity, Commun. Pure. Appl. Math. 6 (1953), 103–156.
  * [15] D. Partyka and K. Sakan, On bi-Lipschitz type inequalities for quasiconformal harmonic mappings, Ann. Acad. Sci. Fenn. Math. 32 (2007), 579–594.
  * [16] X-T. Wang and X-Q. Liang, Precise coefficient estimates for close-to-convex harmonic univalent mappings, J. Math. Anal. Appl 263 (2001), 501–509.


#  Multipolar nematic state of nonmagnetic FeSe based on the DFT+$U$ method

Takemi Yamada<EMAIL_ADDRESS>Takami Tohyama  Department of Applied
Physics, Facility of Science, Tokyo University of Science, 6-3-1 Niijuku,
Katsushika, Tokyo 125-8585, Japan

###### Abstract

Clarifying the origin of nematic state in FeSe is one of urgent problems in
the field of iron-based superconductivity. Motivated by the discovery of a
nematic solution in the density-functional theory implemented by on-site
Coulomb interaction (DFT+$U$) [X. Long et al., npj Quantum Mater. 5, 50
(2020)], we reexamine the $U$ dependence of electronic states in the
nonmagnetic normal state of FeSe and perform full multipolar analyses for the
nematic state. We find that with increasing $U$ the normal state experiences a
topological change of the Fermi surfaces before the emergence of a nematic
ground state. The resulting nematic ground state is a multipolar state having
both antiferro-hexadecapoles in the $E$-representation and ferro-multipoles in
the $B_{2}$-representation on each Fe site. Cooperative coupling between the
$E$ and $B_{2}$ multipoles in local coordinate with the $D_{2d}$ point group
will play an important role in the formation of the $d_{xz},d_{yz}$ orbital-
splitting nematic state not only in FeSe but also in other iron pnictides.

FeSe [1, 2] is one of the most intensively studied iron-based superconductors
[3, 4, 5, 6, 7] because of its complex and versatile ground state under
pressure $P$ [8, 9, 10] and substitution of Se [11, 12]. Below the tetragonal-
orthorhombic structural phase transition at $T_{\mathrm{S}}$=90 K, the
electronic state of FeSe shows a behavior of the nematic state breaking the
$C_{4}$ rotational symmetry while keeping translational symmetry without any
magnetic ordering unlike other iron-based superconductors. In spite of a tiny
orthorhombicity [13], the system exhibits a large band splitting associated
with the orbital differentiation of $d_{xz},d_{yz}$, $|E_{yz}-E_{xz}|$=50 meV
[14], which is too large to attribute to the lattice deformation. This
strongly supports the electronic origin of the nematic state, being consistent
with the enhancement of the nematic fluctuation when $T$ approaches
$T_{\mathrm{S}}$ as observed in the nematic susceptibility [15, 11] and the
electronic Raman response [16, 17]. Therefore, the clarification of the
nematic state is significant for the microscopic understanding of recent
interesting phenomena such as the orbital-selective correlation effect [18,
19, 20, 21] and the BCS-BEC crossover [22, 23, 24, 25].

The Fermi surface (FS) and low-energy bandstructures of FeSe have extensively
been investigated [26, 27, 28, 29, 14, 30, 31, 32, 33, 34, 35, 36, 37, 38],
where the multiorbital compensated metal with Fe-$d$ orbitals is confirmed [6,
7]. In normal (nonnematic) state, two hole FSs (2h-FSs) around $\Gamma$ and
two electron FSs (2e-FSs) around M have commonly been reported but their size
is extremely small only occupied 2-3 % in the Brillouin Zone (BZ). These small
FSs and the low-energy bandstructure have not still been reproduced by the
density-functional theory (DFT) [37], DFT+$U$ [38], the dynamical mean field
theory (DMFT) [18, 19, 20, 39] and quasiparticle self-consistent $GW$ (QS$GW$)
[40]. Several studies based on adjusted models to reproduce the low-energy
bands of angle resolved photoemission spectroscopy (ARPES) [31] can explain
the enhancement of orbital and magnetic fluctuations in $T$-$P$ phase [41, 42,
43].

As for the nematic state, several FSs have been reported by the Shubnikov-de
Haas [26, 27] and ARPES experiments [28, 29, 14, 30, 31, 32, 33, 34, 35, 36],
where a single hole FS (1h-FS) near $\Gamma$ is common while it is still
unsettled whether the electron FS near M is a single (1e-FS) or two. The sign
change of the orbital-splitting at $\Gamma$ and M points in the BZ has been
observed [30, 31], whose origin and mechanism have been discussed [44, 45, 46,
47]. The recent DFT study [48] has provided a new nematic ground state with
the $E_{u}$ irreducible representation of the $D_{4h}$ symmetry, which
contains 1e-FS and additional hybridization between $d_{xy}$ and
$d_{xz},d_{yz}$ orbitals [49, 50]. Although this nematic state seems to
explain the recent experiment [36], it is unclear how the nematic state is
reached from the well-known three hole FSs (3h-FSs) of the DFT normal state
[38, 48]. Therefore, a systematic investigation of the normal state on the
verge of nematic ordering and a detailed multipolar analysis in the nematic
state are highly desirable.

In this Letter, we examine the $U$ dependence of the electronic states of FeSe
by the DFT+$U$ method and find a topological change of FSs before a nematic
order occurs. The resulting nematic ground state is found to be a multipolar
state having both antiferro-hexadecapoles in the $E$-representation and ferro-
multipoles in the $B_{2}$-representation on each Fe site with locally $D_{2d}$
point group. This coexistence indicates that cooperative coupling between the
$E$ and $B_{2}$ multipoles can be a source of the formation of the
$d_{xz},d_{yz}$ orbital-splitting nematic state in FeSe and related materials.

We employ the all-electron first-principles code WIEN2k [51] with the
generalized gradient approximation for the exchange-correlation potential [52]
and the DFT+$U$ method within the self-interaction correction [53, 54]. The
Coulomb interaction $U$ for $d$ electrons in the muffin-tin radius $R_{\rm
MT}^{\nu}$ with atomic sites $\nu$=Fe1,Fe2 in the unit cell is introduced,
where the DFT+$U$ correction energy is given by
$E_{U}=\frac{U}{2}\sum_{\nu}\left(n^{d,\nu}-\sum_{mm^{\prime}}\sum_{\sigma}n_{mm^{\prime}}^{\sigma,\nu}n_{m^{\prime}m}^{\sigma,\nu}\right),$
where $n^{d,\nu}$ is total occupation number of $d$ electrons within $R_{\rm
MT}^{\nu}$ as $n^{d,\nu}=\sum_{m\sigma}n_{mm}^{\sigma,\nu}$ and the density
matrix $n_{mm^{\prime}}^{\sigma,\nu}$. Hereafter, we drop the spin index
$\sigma$ and use
$n_{mm^{\prime}}^{\uparrow,\nu}$=$n_{mm^{\prime}}^{\downarrow,\nu}$=$n_{mm^{\prime}}^{\nu}$
and $n_{mm}^{\nu}$=$n_{m}^{\nu}/2$, since only the nonmagnetic state is
studied.

The Kohn-Sham (KS) equation for valence electrons is written as
$\left(-\frac{1}{2}\bm{\nabla}^{2}+v^{\rm KS}_{{\rm
eff}}(\bm{r})+v_{U}(\bm{r})\right)\psi_{\bm{k}n}(\bm{r})=\varepsilon_{\bm{k}n}\psi_{\bm{k}n}(\bm{r})$,
where $\psi_{\bm{k}n}(\bm{r})$ is the KS eigenfunction with wavevector
$\bm{k}$ and band-index $n$ with the band energy $\varepsilon_{\bm{k}n}$ and
$v^{\rm KS}_{\rm eff}(\bm{r})$ is the effective KS potential. The +$U$
potential $v_{U}(\bm{r})$ is given by [54, 55]

$\displaystyle
v_{U}(\bm{r})=\sum_{\nu}\sum_{mm^{\prime}}v_{mm^{\prime}}^{\nu}\frac{\delta
n_{mm^{\prime}}^{\nu}}{\delta\psi_{\bm{k}n}^{*}(\bm{r})},$ (1) $\displaystyle
v_{mm^{\prime}}^{\nu}=\delta_{mm^{\prime}}\frac{U}{2}(1-n_{m}^{\nu})+(1-\delta_{mm^{\prime}})(-Un_{mm^{\prime}}^{\nu}),$
(2)

where $\delta n_{mm^{\prime}}^{\nu}/\delta\psi_{\bm{k}n}^{*}$ in Eq. (1) is
the projector for $\psi_{\bm{k}n}$ and 1st (2nd) term in Eq. (2) is
proportional to the diagonal (off-diagonal) density matrix $n_{mm}^{\nu}$
($n_{mm^{\prime}}^{\nu}$). In calculation, we use $\bm{k}$-mesh of
30${}^{2}\times$21 and the tetragonal crystallographic parameters are taken
from the experiment [2], and $R_{\rm MT}^{\rm Fe(Se)}$=2.25 (2.14) bohr with
$R_{\rm MT}K_{\rm max}$=8.

Figure 1: (Color online) $U$ dependence of pDOS $\rho_{m}^{d}(\varepsilon_{\rm
F})$. The inset shows the 3h-FSs for $U$=0 eV (left center) and the 2h-FSs for
$U$=3 eV (right top).  Figure 2: (Color online) (a) Energy differences of the
total energy $\Delta E_{\rm tot}$ (left axis), and kinetic and potential
energies $\Delta T_{s}$, $\Delta U_{\rm pot}$ (right axis) and (b)
$n_{m}^{d}$. (c)-(f): The obtained band energy
$\varepsilon_{\bm{k}n}-\varepsilon_{\rm F}$ around $\varepsilon_{\rm F}$ for
$U=3\sim 4$ eV.

First, we investigate the $U$ dependence of the normal electronic state in
FeSe. Figure 1 shows the partial density-of-states (pDOS)
$\rho_{m}^{d}(\varepsilon_{\rm F})$ at the Fermi energy $\varepsilon_{\rm F}$.
When $U=0$ eV, 3h-FSs consisting of the 28th, 29th, and 30th bands are
obtained as shown in the inset of Fig. 1, being similar to the previous DFT
results [37, 38]. With increasing $U$, $\rho_{m}^{d}(\varepsilon_{\rm F})$ of
$d_{xy}$ orbital drops at $U\simeq 2$ eV (=$U_{\rm LT}$), while in turn that
of $d_{xz},d_{yz}$ increases gradually. Since the 28th band constructing the
most inner hole FS originates from the $d_{xy}$ orbital, it falls below
$\varepsilon_{\rm F}$ and the FS vanishes at a Lifshitz transition point
$U_{\rm LT}$, where the 3h-FSs change to the 2h-FSs from the 29th and 30th
bands as shown in the inset of Fig. 1 [56]. This change is caused by
$v_{U}(\bm{r})$ that induces orbital-dependent energy shifts, being
proportional to $1-n_{m}^{\nu}$ as shown in the first term of Eq. (2). We can
find that the occupied number $n_{m}^{d}$ in $d_{xy}$ and $d_{z^{2}}$ at $U$=0
eV, which is larger than in $d_{xz},d_{yz}$ and $d_{x^{2}-y^{2}}$, increases
further with increasing $U$, resulting in negatively large value of
$1-n_{m}^{\nu}$ in $d_{xy}$ [57].

In contrast to the h-FSs, the e-FSs from the 31st and 32nd bands with
$d_{xz},d_{yz}$ orbitals at M and A points do not undergo any topological
change but instead $\varepsilon_{\bm{k}n}$ of these bands approaches
$\varepsilon_{\rm F}$ monotonically with increasing $U$. This manifestation of
the $d_{xz},d_{yz}$ orbitals near $\varepsilon_{\rm F}$ as evidenced by the
2h-FSs and e-FSs will triggers the formation of a nematic state mentioned
below.

Next we calculate a nematic solution by preconditioning the initial charge
density [58] as was done in the previous pseudopotential calculation [48].
Figure 2 (a) shows the total energy difference between the normal and nematic
states $\Delta E_{\rm tot}=E_{\rm tot}^{\rm nem}-E_{\rm tot}^{\rm normal}$ as
a function of $U$, together with the kinetic energy and potential energy
differences $\Delta T_{s}$ and $\Delta U_{\rm pot}$, where $E_{\rm
tot}=T_{s}+U_{\rm pot}$ is calculated by the total energy formula of the all-
electron method [59] under the Virial theorem $2T_{s}+U_{\rm pot}=0$. With
increasing $U$ more than $U=3.4$ eV, $\Delta E_{\rm tot}<0$ is realized with
$\Delta T_{s}>0$ and $\Delta U_{\rm pot}<0$ together with the occupied number
splitting between $d_{xz}$ and $d_{yz}$ orbitals as shown in Fig.2(b). The
energy gain is of order of $O(10^{1\sim 2}{\rm meV})$, which basically agrees
with the previous pseudopotential DFT result [48].

The $U$ dependence of $\varepsilon_{\bm{k}n}-\varepsilon_{\rm F}$ near
$\varepsilon_{\rm F}$ at $\Gamma$ and Z [M and A] points is shown in Figs.
2(c) and 2(d) [Figs. 2(e) and 2(f)], respectively. The $d_{xz}$\- and
$d_{yz}$-orbital bands, whose number is denoted in the figures, split due to
the nematic state transition. Since the 32nd band in M and A points rises
above $\varepsilon_{\rm F}$, 2e-FSs in the nonnematic state change to 1e-FS in
the nematic state. Combining with the change from 2h-FSs to 1h-FS at $\Gamma$
point above $U$=3.7 eV, we find the number of FSs consistent with the
experiment [36] and the previous DFT result [48].

Table 1:  The definition and notations of multipole $X_{\alpha}$ denoted by
quadrupole $O_{\Gamma_{\gamma}}$ and hexadecapole $Q_{4_{\Gamma_{\gamma}}}$
together with the irreducible representations (IRRs) in $D_{2d}$ symmetry [60,
61] at each Fe site. $T_{kq}^{(c,s)}$ ($T_{kq}$) is the tesseral (spherical)
tensor operator with a relation
$T_{kq}^{(c,s)}=\frac{(-1)^{q}}{\sqrt{2}}(T_{kq}^{\dagger}\pm T_{kq})$.

IRR | $X_{\alpha}$ | tesseral tensor-rep. |  $(x,y,z)$-rep.
---|---|---|---
$A_{1}$ | $O_{u}$ | $T_{20}$ | $\frac{1}{2}(3z^{2}-r^{2})$
| $Q_{4}$ | $\sqrt{\frac{5}{12}}T_{44}^{(c)}+\sqrt{\frac{7}{12}}T_{40}$ | $\frac{5\sqrt{21}(x^{4}+y^{4}+z^{4}-\frac{3}{5}r^{4})}{12}$
| $Q_{4u}$ | $-\sqrt{\frac{7}{12}}T_{44}^{(c)}+\sqrt{\frac{5}{12}}T_{40}$ | $\frac{7\sqrt{15}[2z^{4}-x^{4}-y^{4}-\frac{6}{7}r^{2}(3z^{2}-r^{2})]}{12}$
$A_{2}$ | $Q_{4\alpha,z}$ | $-T_{44}^{(s)}$ | $\frac{\sqrt{35}}{2}xy(x^{2}-y^{2})$
$B_{1}$ | $O_{v}$ | $T_{22}^{(c)}$ | $\frac{\sqrt{3}}{2}(x^{2}-y^{2})$
| $Q_{4v}$ | $-T_{42}^{(c)}$ | $\frac{7\sqrt{5}[x^{4}-y^{4}-\frac{6}{7}r^{2}(x^{2}-y^{2})]}{4}$
$B_{2}$ | $O_{xy}$ | $T_{22}^{(s)}$ | $\sqrt{3}xy$
| $Q_{4\beta,z}$ | $T_{42}^{(s)}$ | $\frac{\sqrt{5}}{2}xy(7z^{2}-r^{2})$
$E$ | $O_{zx},O_{yz}$ | $T_{21}^{(c)},T_{21}^{(s)}$ | $\sqrt{3}zx,\sqrt{3}yz$
| $Q_{4\alpha,x}$ | $-\sqrt{\frac{1}{8}}T_{43}^{(s)}-\sqrt{\frac{7}{8}}T_{41}^{(s)}$ | $\frac{\sqrt{35}}{2}yz(y^{2}-z^{2})$
| $Q_{4\alpha,y}$ | $-\sqrt{\frac{1}{8}}T_{43}^{(c)}+\sqrt{\frac{7}{8}}T_{41}^{(c)}$ | $\frac{\sqrt{35}}{2}zx(z^{2}-x^{2})$
| $Q_{4\beta,x}$ | $\sqrt{\frac{7}{8}}T_{43}^{(s)}-\sqrt{\frac{1}{8}}T_{41}^{(s)}$ | $\frac{\sqrt{5}}{2}yz(7x^{2}-r^{2})$
| $Q_{4\beta,y}$ | $-\sqrt{\frac{7}{8}}T_{43}^{(c)}-\sqrt{\frac{1}{8}}T_{41}^{(c)}$ | $\frac{\sqrt{5}}{2}zx(7y^{2}-r^{2})$
|  |  |

Figure 3: (Color online) The quadrupole and hexadecapole moments (a), (b)
$O_{\Gamma_{\gamma}}$ and (c), (d) $Q_{4\Gamma_{\gamma}}$ on Fe1 [(a), (c)]
and Fe2 [(b), (d)] sites as a function of $U$. In inset of (a) and (b), atom
sites in the unit-cell are depicted with $x,y$-axis.

The order parameters of the nematic state obtained by DFT+$U$ have been
discussed in previous study [48], where only finite off-diagonal density
matrix elements are taken into account. Here, we derive all the active
multipole moments in the present system more generally [60, 61]. Without the
spin-orbit interaction, the multipole operator can be regarded as a power
series expansion of the rank $k$ of orbital angular momentum operator
$\bm{\ell}=(\ell_{x},\ell_{y},\ell_{z})$. Only the even rank multipoles, i.e.,
quadrupoles ($k=2$) and hexadecapoles ($k=4$), become finite for $d$ electron
basis in each Fe site due to the time reversal symmetry. These multipoles are
classified by the irreducible representations {$A_{1}$, $A_{2}$, $B_{1}$,
$B_{2}$, $E$} at Fe site with $D_{2d}$ symmetry, which are summarized in Table
1.

We calculate all the multipole moments listed in Table 1 from the density
matrix $n_{mm^{\prime}}^{\nu}$, where all the multipole operators are
normalized as ${\rm Tr}[X_{\alpha}X_{\beta}]=\delta_{\alpha\beta}$. We note
that a similar approach has been performed on the magnetic multipole order in
the actinide dioxides [62, 63]. We find that the multipoles directly related
to the nematic order are those in the $B_{2}$ and $E$ representations. Figures
3(a) and 3(c) [3(b) and 3(d)] show the $U$-dependence of the quadrupoles
$O_{xz}$, $O_{yz}$, $O_{xy}$ and the hexadecapoles $Q_{4\beta,z}$,
$Q_{4\alpha,x}$, $Q_{4\alpha,y}$, $Q_{4\beta,x}$, $Q_{4\beta,y}$ at the Fe1
[Fe2] site, respectively. Above $U$=3.4 eV, the $B_{2}$ quadrupole moment for
$O_{xy}$ at the two Fe sites becomes negative as seen in Figs. 3(a) and 3(b)
and the similar behavior is obtained for the $B_{2}$ hexadecapole
$Q_{4\beta,z}$ as seen in Figs. 3(c) and 3(d). This behavior corresponds to
the emergence of a ferro nematic order associated with the orbital
differentiation between $d_{xz}$ and $d_{yz}$. On the other hand, the $E$
hexadecapole moments for $Q_{4\alpha,x}$, $Q_{4\alpha,y}$, $Q_{4\beta,x}$, and
$Q_{4\beta,y}$ are more than ten times larger than the $E$ quadrupole moments
for $O_{xz},O_{yz}$. It is also interesting to notice that $Q_{4\beta,x}>0$,
$Q_{4\beta,y}<0$, and $Q_{4\alpha,x(y)}>0$ at the Fe1 but opposite signs for
the Fe2 site as shown in Figs .3(c) and 3(d). Namely, the $E$-type order
parameter at Fe1 (Fe2) is written as
$\pm(aQ_{4\alpha}^{E}+bQ_{4\beta}^{E}+cO^{E})$, where
$Q_{4\alpha(\beta)}^{E}=\frac{1}{\sqrt{2}}(Q_{4\alpha(\beta),x}\pm
Q_{4\alpha(\beta),y})$ and $O^{E}$=$\frac{1}{\sqrt{2}}(O_{zx}-O_{yz})$ with
$a^{2}$+$b^{2}$+$c^{2}$=1. This result indicates that the antiferro ordering
of the $E$ multipoles with opposite values at two Fe sites coexists with the
ferro ordering of the $B_{2}$ multipoles with the same values at two sites.

Figure 4: (Color online) Bandstructures for (a)-(c) normal state ($U$=3.3 eV)
and (d)-(f) nematic state ($U$=3.4 eV) with the orbital-weight of (a),(d)
$d_{xy}$, (b),(e) $d_{xz},d_{yz}$ and (c),(f) $d_{x^{2}-y^{2}}$ along the high
symmetry line in BZ where ${\rm M}_{x}$[${\rm
M}_{y}$]=$(\frac{\pi}{a},\frac{\pi}{a},0)[(\frac{\pi}{a},-\frac{\pi}{a},0)]$
and ${\rm A}_{x}$[${\rm
A}_{y}$]=$(\frac{\pi}{a},\frac{\pi}{a},\frac{\pi}{c})[(\frac{\pi}{a},-\frac{\pi}{a},\frac{\pi}{c})]$
and $a,c$ are the lattice constants.

Such the coexistence of $E$\- and $B_{2}$-type multipoles can be understood
from the phenomenological intermultipole coupling theory [60, 64] at finite
temperature $T$, where the Ginzburg-Landau free energy with the mean-field
(MF) approximation can be expanded by the multipole moment $X_{\alpha}$ around
the nematic transition as given by $F_{\rm MF}=F_{\rm MF}^{\rm(2nd)}+F_{\rm
MF}^{\rm(3rd)}+\cdots$. Here, $F_{\rm
MF}^{\rm(3rd)}=-\frac{T}{3!}\sum_{\alpha\beta\gamma}g_{\alpha\beta\gamma}X_{\alpha}X_{\beta}X_{\gamma}$,
where $X_{\alpha}$=$O_{\Gamma_{\gamma}}$ or $Q_{4\Gamma_{\gamma}}$ and
$g_{\alpha\beta\gamma}$ is the symmetric constant defined as
$g_{\alpha\beta\gamma}=\frac{1}{2d}{\rm
Tr}\left[(X_{\alpha}X_{\beta}+X_{\beta}X_{\alpha})X_{\gamma}\right]$ [65] with
matrix dimension $d$=5. The coupling terms among $O_{xy},Q_{4\beta,z}$ and $E$
multipoles are explicitly given by $F_{\rm MF}^{\rm(3rd)}$=
$-Tc_{1}O_{xy}\left(-\frac{4}{11}O_{xz}O_{yz}-\frac{7}{11}Q_{4\alpha,x}Q_{4\alpha,y}+Q_{4\beta,x}Q_{4\beta,y}\right)$
$-Tc_{2}Q_{4\beta,z}\left(O_{xz}O_{yz}-\frac{7}{8}Q_{4\alpha,x}Q_{4\alpha,y}+Q_{4\beta,x}Q_{4\beta,y}\right)$,
where
$c_{1}(c_{2})$=$\tfrac{11}{84}\sqrt{\tfrac{15}{14}}(\frac{1}{21}\sqrt{\frac{10}{7}})$
and 2nd and 3rd terms in both parentheses correspond to the coupling between
the $E$ multipoles and are negative while the coefficients of both parentheses
including the $B_{2}$ multipoles are positive as shown in Fig. 3. Therefore,
$F_{\rm MF}^{\rm(3rd)}$ becomes negative as a whole stabilizing the
coexistence state.

Finally, we discuss the band structure and orbital components of the normal
($U$=3.3eV) and nematic ($U$=3.4eV) states as shown in Figs. 4(a)-4(c) and
4(d)-4(f), respectively. It is clearly observed that the degenerated $d_{xz}$
and $d_{yz}$ bands near 0.25 eV (0.5 eV) at $\Gamma$ (Z) point split in the
nematic state, corresponding to the ordering of $O_{xy}$ and $Q_{4\beta,z}$ as
seen in Fig. 3. Band splitting due to the same mechanism is also realized in
the 31st and 32nd electron bands near 0 eV along the M-A direction, which
consist of linear combinations of $d_{xz}$ and $d_{yz}$ orbitals without
$k_{z}$ dependence. This leads to a splitting of the peak structure in pDOS
$\rho_{m}^{d}(\varepsilon_{\rm F})$ of the $d_{xz}$ and $d_{yz}$ orbitals. In
the nematic state as shown in Fig. 3, the $E$ hexadecapoles become finite,
which gives rise to the significant mixings among the $|\ell_{z}|$=2 orbitals
($d_{xy}$ and $d_{x^{2}-y^{2}}$) and $|\ell_{z}|$=1 orbitals ($d_{xz}$ and
$d_{yz}$)[66]. The mixing with $d_{xy}(d_{x^{2}-y^{2}})$ having a large weight
below (above) $\varepsilon_{\rm F}$ pushes down (up) the energy level of the
mixed partner $d_{yz}(d_{xz})$. Such a mixing effect is remarkable in the
My-$\Gamma$ and Ay-Z directions as seen in Figs. 4(a) and 4(d): a mixing gap
is formed near $-0.1$ eV in their directions. On the other hand, the band for
the $d_{xy}$ orbital near $-0.2$ eV in the Mx-$\Gamma$ and Ax-Z directions
remains almost unchanged. This one-side gap opening along the Mx-$\Gamma$-My
direction, which is consistent with the experiment [36], is critically
important for the origin of the nematic state, because it inherently requires
the presence of the $E$ multipoles as pointed out in this paper.

In summary, we have studied the nonmagnetic normal and nematic states of the
iron-based superconductor FeSe by using the DFT+$U$ method with the multipole
analyses. The effect of $U$ on the normal state generates a topological change
of FSs from 3h-FSs to 2h-FSs, leading to a change of dominant orbital near
$\varepsilon_{\rm F}$ from $d_{xy}$ to $d_{xz},d_{yz}$. As a result, the
multipolar nematic state with the $E$ antiferro-hexadecapoles accompanying the
$B_{2}$ ferro-multipoles has been obtained without any assumption of the order
parameters, giving rise to both of the $d_{xz}$-$d_{yz}$ orbital splitting at
$\Gamma$ and the $d_{xy}$-$(d_{xz},d_{yz})$ orbital mixing around M and A
points. From phenomenological analysis, we have found that the intermultipole
coupling of $B_{2}$ and $E$ multipoles on each Fe site can explain the energy
gain larger for the coexisting order than for the quadupole $O_{xy}$ order
alone. This multipolar mechanism for the formation of nematic state will be
applicable not only to FeSe but also to other iron pnictides where the
degenerated $d_{xz}$,$d_{yz}$ orbitals play a crucial role.

###### Acknowledgements.

This work was supported by the “Quantum Liquid Crystals” No. JP19H05825
KAKENHI on Innovative Areas from JSPS of Japan.

## References

  * Hsu  _et al._ [2008] F.-C. Hsu, J.-Y. Luo, K.-W. Yeh, T.-K. Chen, T.-W. Huang, P. M. Wu, Y.-C. Lee, Y.-L. Huang, Y.-Y. Chu, D.-C. Yan, and M.-K. Wu, Proc. Natl. Acad. Sci. U.S.A. 105, 14262 (2008).
  * McQueen  _et al._ [2009] T. M. McQueen, A. J. Williams, P. W. Stephens, J. Tao, Y. Zhu, V. Ksenofontov, F. Casper, C. Felser, and R. J. Cava, Phys. Rev. Lett. 103, 057002 (2009).
  * Mizuguchi and Takano [2010] Y. Mizuguchi and Y. Takano, J. Phys. Soc. Jpn. 79, 102001 (2010).
  * Liu  _et al._ [2015] X. Liu, L. Zhao, S. He, J. He, D. Liu, D. Mou, B. Shen, Y. Hu, J. Huang, and X. J. Zhou, J. Phys.: Condens. Mater. 27, 183201 (2015).
  * B$\ddot{\rm o}$hmer and Kreisel [2018] A. E. B$\ddot{\rm o}$hmer and A. Kreisel, J. Phys.: Condens. Mater. 30, 023001 (2018).
  * Shibauchi  _et al._ [2020] T. Shibauchi, T. Hanaguri, and Y. Matsuda, J. Phys. Soc. Jpn. 89, 102002 (2020).
  * Kreisel  _et al._ [2020] A. Kreisel, P. J. Hirschfeld, and B. M. Andersen, Symmetry 12, 1402 (2020).
  * Sun  _et al._ [2016] J. P. Sun, K. Matsuura, G. Z. Ye, Y. Mizukami, M. Shimozawa, K. Matsubayashi, M. Yamashita, T. Watashige, S. Kasahara, Y. Matsuda, J.-Q. Yan, B. C. Sales, Y. Uwatoko, J.-G. Cheng, and T. Shibauchi, Nat. Commun. 7, 12146 (2016).
  * Terashima  _et al._ [2015] T. Terashima, N. Kikugawa, S. Kasahara, T. Watashige, T. Shibauchi, Y. Matsuda, T. Wolf, A. E. B$\ddot{\rm o}$hmer, F. Hardy, C. Meingast, H. v. L$\ddot{\rm o}$hneysen, and S. Uji, J. Phys. Soc. Jpn. 84, 063701 (2015).
  * Terashima  _et al._ [2016a] T. Terashima, N. Kikugawa, S. Kasahara, T. Watashige, Y. Matsuda, T. Shibauchi, and S. Uji, Phys. Rev. B 93, 180503(R) (2016a).
  * Hosoi  _et al._ [2016] S. Hosoi, K. Matsuura, K. Ishida, H. Wang, Y. Mizukami, T. Watashige, S. Kasahara, Y. Matsuda, and T. Shibauchi, Proc. Natl. Acad. Sci. U.S.A. 113, 8139 (2016).
  * Coldea  _et al._ [2019] A. I. Coldea, S. F. Blake, S. Kasahara, A. A. Haghighirad, M. D. Watson, W. Knafo, E. S. Choi, A. McCollam, P. Reiss, T. Yamashita, M. Bruma, S. C. Speller, Y. Matsuda, T. Wolf, T. Shibauchi, and A. J. Schofield, npj Quantm Mater. 4, 2 (2019).
  * B$\ddot{\rm o}$hmer  _et al._ [2013] A. E. B$\ddot{\rm o}$hmer, F. Hardy, F. Eilers, D. Ernst, P. Adelmann, P. Schweiss, T. Wolf, and C. Meingast, Phys. Rev. B 87, 180505(R) (2013).
  * Shimojima  _et al._ [2014] T. Shimojima, Y. Suzuki, T. Sonobe, A. Nakamura, J. O. M. Sakano, K. Yoshioka, M. Kuwata-Gonokami, K. Ono, H. Kumigashira, A. E. B$\ddot{\rm o}$hmer, F. Hardy, T. Wolf, C. Meingast, H. v. L$\ddot{\rm o}$hneysen, H. Ikeda, and K. Ishizaka, Phys. Rev. B 90, 121111(R) (2014).
  * B$\ddot{\rm o}$hmer  _et al._ [2015] A. E. B$\ddot{\rm o}$hmer, T. Arai, F. Hardy, T. Hattori, T. Iye, T. Wolf, H. v. L$\ddot{\rm o}$hneysen, K. Ishida, and C. Meingast, Phys. Rev. Lett. 114, 027001 (2015).
  * Massat  _et al._ [2016] P. Massat, D. Farina, I. Paul, S. Karlsson, P. Strobel, P. Toulemonde, M.-A. M$\check{\rm e}$asson, M. Cazayous, A. Sacuto, S. Kasahara, T. Shibauchi, Y. Matsuda, and Y. Gallais, Proc. Natl. Acad. Sci. U.S.A. 113, 9177 (2016).
  * Zhang  _et al._ [2020] W. Zhang, S. Wu, S. Kasahara, T. Shibauchi, Y. Matsuda, and G. Blumberg, Proc. Natl. Acad. Sci. U.S.A. 118, e2020585118 (2020).
  * Aichhorn  _et al._ [2010] M. Aichhorn, S. Biermann, T. Miyake, A. Georges, and M. Imada, Phys. Rev. B 82, 064504 (2010).
  * Liebsch and Ishida [2010] A. Liebsch and H. Ishida, Phys. Rev. B 82, 155106 (2010).
  * Yin  _et al._ [2011] Z. P. Yin, K. Haule, and G. Kotliar, Nat. Mater. 10, 932 (2011).
  * Lanat$\check{\rm a}$  _et al._ [2013] N. Lanat$\check{\rm a}$, H. U. R. Strand, G. Giovannetti, B. Hellsing, L. de’ Medici, , and M. Capone, Phys. Rev. B 87, 045122 (2013).
  * Kasahara  _et al._ [2014] S. Kasahara, T. Watashige, T. Hanaguri, Y. Kohsaka, T. Yamashita, Y. Shimoyama, Y. Mizukami, R. Endo, H. Ikeda, K. Aoyama, T. Terashima, S. Uji, T. Wolf, H. v. L$\ddot{\rm o}$hneysen, T. Shibauchi, and Y. Matsuda, Proc. Natl. Acad. Sci. U.S.A. 111, 16309 (2014).
  * Kasahara  _et al._ [2016] S. Kasahara, T. Yamashita, A. Shi, R. Kobayashi, Y. Shimoyama, T. Watashige, K. Ishida, T. Terashima, T. Wolf, F. Hardy, C. Meingast, H. v. L$\ddot{\rm o}$hneysen, A. Levchenko, T. Shibauchi, and Y. Matsuda, Nat. Commun. 7, 12843 (2016).
  * Hanaguri  _et al._ [2019] T. Hanaguri, S. Kasahara, I. E. J. B$\ddot{\rm o}$ker, T. Shibauchi, and Y. Matsuda, Phys. Rev. Lett. 122, 077001 (2019).
  * Kasahara  _et al._ [2020] S. Kasahara, Y. Sato, S. Licciardello, M. $\acute{\rm C}$ulo, S. Arsenijevi$\check{\rm c}$, T. Ottenbros, T. Tominaga, J. B$\ddot{\rm o}$ker, I. Eremin, T. Shibauchi, J. Wosnitza, N. E. Hussey, and Y. Matsuda, Phys. Rev. Lett. 124, 107001 (2020).
  * Terashima  _et al._ [2014] T. Terashima, N. Kikugawa, A. Kiswandhi, E.-S. Choi, J. S. Brooks, S. Kasahara, T. Watashige, H. Ikeda, T. Shibauchi, Y. Matsuda, T. Wolf, A. E. B$\ddot{\rm o}$hmer, F. Hardy, C. Meingast, H. v. L$\ddot{\rm o}$hneysen, and S. Uji, Phys. Rev. B 90, 144517 (2014).
  * Terashima  _et al._ [2016b] T. Terashima, N. Kikugawa, A. Kiswandhi, D. Graf, E.-S. Choi, J. S. Brooks, S. Kasahara, T. Watashige, Y. Matsuda, T. Shibauchi, T. Wolf, A. E. B$\ddot{\rm o}$hmer, F. Hardy, C. Meingast, H. v. L$\ddot{\rm o}$hneysen, and S. Uji, Phys. Rev. B 93, 094505 (2016b).
  * Nakayama  _et al._ [2014] K. Nakayama, Y. Miyata, G. N. Phan, T. Sato, Y. Tanabe, T. Urata, K. Tanigaki, and T. Takahashi, Phys. Rev. Lett. 113, 237001 (2014).
  * Maletz  _et al._ [2014] J. Maletz, V. B. Zabolotnyy, D. V. Evtushinsky, S. Thirupathaiah, A. U. B. Wolter, L. Harnagea, A. N. Yaresko, A. N. Vasiliev, D. A. Chareev, A. E. B$\ddot{\rm o}$hmer, F. Hardy, T. Wolf, C. Meingast, E. D. L. Rienks, B. B$\ddot{\rm u}$chner, and S. V. Borisenko, Phys. Rev. B 89, 220506(R) (2014).
  * Watson  _et al._ [2015] M. D. Watson, T. K. Kim, A. A. Haghighirad, N. R. Davies, A. McCollam, A. Narayanan, S. F. Blake, Y. L. Chen, S. Ghannadzadeh, A. J. Schofield, M. Hoesch, C. Meingast, T. Wolf, and A. I. Coldea, Phys. Rev. B 91, 155106 (2015).
  * Suzuki  _et al._ [2015] Y. Suzuki, T. Shimojima, T. Sonobe, A. Nakamura, M. Sakano, H. Tsuji, J. Omachi, K. Yoshioka, M. Kuwata-Gonokami, T. Watashige, R. Kobayashi, S. Kasahara, T. Shibauchi, Y. Matsuda, Y. Yamakawa, H. Kontani, and K. Ishizaka, Phys. Rev. B 92, 205117 (2015).
  * Zhang  _et al._ [2015] P. Zhang, T. Qian, P. Richard, X. P. Wang, H. Miao, B. Q. Lv, B. B. Fu, T. Wolf, C. Meingast, X. X. Wu, Z. Q. Wang, J. P. Hu, and H. Ding, Phys. Rev. B 91, 214503 (2015).
  * Watson  _et al._ [2016] M. D. Watson, T. K. Kim, L. C. Rhodes, M. Eschrig, M. Hoesch, A. A. Haghighirad, and A. I. Coldea, Phys. Rev. B 94, 201107(R) (2016).
  * Watson  _et al._ [2017] M. D. Watson, A. A. Haghighirad, L. C. Rhodes, M. Hoesch, and T. K. Kim, New J. Phys. 19, 103021 (2017).
  * Pfau  _et al._ [2019] H. Pfau, S. D. Chen, M. Yi, M. Hashimoto, C. R. Rotundu, J. C. Palmstrom, T. Chen, P.-C. Dai, J. Straquadine, A. Hristov, R. J. Birgeneau, I. R. Fisher, D. Lu, and Z.-X. Shen, Phys. Rev. Lett. 123, 066402 (2019).
  * Yi  _et al._ [2019] M. Yi, H. Pfau, Y. Zhang, Y. He, H. Wu, T. Chen, Z. R. Ye, M. Hashimoto, R. Yu, Q. Si, D.-H. Lee, P.-C. Dai, Z.-X. Shen, D. H. Lu, and R. J. Birgeneau, Phys. Rev. X 9, 041049 (2019).
  * Subedi  _et al._ [2008] A. Subedi, L. Zhang, D. J. Singh, and M. H. Du, Phys. Rev. B 78, 134514 (2008).
  * Lohani  _et al._ [2015] H. Lohani, P. Mishra, and B. Sekhar, Physica C 512, 54 (2015).
  * Mandal  _et al._ [2014] S. Mandal, R. E. Cohen, and K. Haule, Phys. Rev. B 89, 220502(R) (2014).
  * Tomczak  _et al._ [2012] J. M. Tomczak, M. van Schilfgaarde, and G. Kotliar, Phys. Rev. Lett. 109, 237010 (2012).
  * Yamakawa  _et al._ [2016] Y. Yamakawa, S. Onari, and H. Kontani, Phys. Rev. X 6, 021032 (2016).
  * Yamakawa and Kontani [2017] Y. Yamakawa and H. Kontani, Phys. Rev. B 96, 144509 (2017).
  * Ishizuka  _et al._ [2018] J. Ishizuka, T. Yamada, Y. Yanagi, and Y. $\bar{\rm O}$no, J. Phys. Soc. Jpn. 87, 014705 (2018).
  * Su  _et al._ [2015] Y. Su, H. Lian, and T. Li, J. Phys.: Condens. Mater. 27, 105702 (2015).
  * Onari  _et al._ [2016] S. Onari, Y. Yamakawa, and H. Kontani, Phys. Rev. Lett. 116, 227001 (2016).
  * Xing  _et al._ [2018] R.-Q. Xing, L. Classen, and A. V. Chubukov, Phys. Rev. B 98, 041108(R) (2018).
  * Kang  _et al._ [2018] J. Kang, R. M. Fernandes, and A. Chubukov, Phys. Rev. Lett. 120, 267001 (2018).
  * Long  _et al._ [2020] X. Long, S. Zhang, F. Wang, and Z. Liu, npj Quantm Mater. 5, 50 (2020).
  * Steffensen  _et al._ [2021] D. Steffensen, A. Kreisel, P. J. Hirschfeld, and B. M. Andersen, Phys. Rev. B 103, 054505 (2021).
  * Rhodes  _et al._ [2021] L. C. Rhodes, J. B$\ddot{\rm o}$ker, M. A. M$\ddot{\rm u}$ller, M. Eschrig, and I. M. Eremin, npj Quantm Mater. 6, 45 (2021).
  * Blaha  _et al._ [2020] P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen, and L. D. Marks, J. Chem. Phys. 152, 074101 (2020).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Anisimov  _et al._ [1993] V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, Phys. Rev. B 48, 16929 (1993).
  * Shick  _et al._ [1999] A. B. Shick, A. I. Liechtenstein, and W. E. Pickett, Phys. Rev. B 60, 10763 (1999).
  * Tran  _et al._ [2008] F. Tran, J. Kune$\check{\rm s}$, P. Nov$\acute{\rm a}$k, P. Blaha, L. D. Marks, and K. Schwarz, Comput. Phys. Commun. 179, 784 (2008).
  * [56] Here we also check $\varepsilon_{\bm{k}n}-\varepsilon_{\rm F}$ around $\varepsilon_{\rm F}$ for the hole (electron) bands at $\Gamma$ and Z (M and A) points as a function $U$, and we obtain a continuous tracking of the connection from the 3h-FSs and 2h-FSs together with a critical interaction of the Lifshitz transition $U_{\rm LT}$, which is a one of the new results for the present study in contrast to the previous DFT+$U$ study [38] where the 2h-FSs were obtained only for several large values of $U$.
  * [57] Throughout the calculation, total $d$ occupation number $n^{d,\nu}$ is almost unchanged from $n^{d,\nu}$=6.17$\sim$6.19 and the total valence electron number $n_{\rm v}$ is strictly fixed $n_{\rm v}$=60 per 2(FeSe) in the unit-cell, where $28$ electrons are attributed to Fe-$3p,3d,4s$ and $32$ electrons to Se-$3d,4s,4p$.
  * [58] In ordering a symmetry broken nonmagnetic solution with WIEN2k, it is needed to prepare an appropriate initial charge density obtained from the SCF converged result in the distorted lattice as $(x_{\rm Fe},y_{\rm Fe})\rightarrow(x_{\rm Fe}+\delta,y_{\rm Fe}+\delta)$ with $\delta=0.03$. In addition to this, the following points should be noted: select the space group $P1$ and directly input the coordinates of each atom in the unit cell in the structure file (case.struct), leaving the lattice constants and $R_{\rm MT}$ the same as in the original structure file. Set the local rotation matrix in the same file to the unit matrix, and set NUMBER OF SYMMETRY OPERATIONS to 1 to give the unit matrix. Skip the x symmetry and x group processes which are the default initial charge calculation.
  * Weinert  _et al._ [1982] M. Weinert, E. Wimmer, and A. J. Freeman, Phys. Rev. B 26, 4571 (1982).
  * Kusunose [2008] H. Kusunose, J. Phys. Soc. Jpn. 77, 064701 (2008).
  * Hayami  _et al._ [2018] S. Hayami, M. Yatsushiro, Y. Yanagi, and H. Kusunose, J. Phys. Soc. Jpn. 98, 165110 (2018).
  * Suzuki  _et al._ [2013] M.-T. Suzuki, N. Magnani, and P. M. Oppeneer, Phys. Rev. B 88, 195146 (2013).
  * Suzuki  _et al._ [2018] M.-T. Suzuki, H. Ikeda, and P. M. Oppeneer, 87, 041008 (2018).
  * Kusunose and Kuramoto [2001] H. Kusunose and Y. Kuramoto, J. Phys. Soc. Jpn. 70, 1751 (2001).
  * [65] Here we normalize the multipoles as $\frac{1}{d}{\rm Tr}[X_{\alpha}X_{\beta}]=\delta_{\alpha\beta}$ with $d=5$ unlike the calculation of the density matrix $n_{mm^{\prime}}^{\nu}$.
  * [66] Here we note that the $E$ quadrupoles $O_{zx},O_{yz}$ also induce the similar mixings but their hybridizations are partial and incomplete, since in $|m\rangle$ orbital basis with $m$=-2,$\cdots$,+2, $O_{\Gamma_{\gamma}}$ with $\Gamma_{\gamma}$=$zx,yz$ has only finite elements as $\langle m|O_{\Gamma_{\gamma}}|m\pm 1\rangle$, while the hexadecapoles $Q_{4\Gamma_{\gamma}}$ with $\Gamma_{\gamma}$=$\alpha,x(y)$ and $\beta,x(y)$ have finite elements as $\langle m|Q_{4\Gamma_{\gamma}}|m\pm 3\rangle$ in addition to $\langle m|Q_{4\Gamma_{\gamma}}|m\pm 1\rangle$.