# Virtual Reality as a Tool for Studying Diversity and Inclusion in Human-
Robot Interaction: Advantages and Challenges

André Helgert1, Sabrina C. Eimler1, Carolin Straßmann1 1Computer Science
Institute, University of Applied Sciences Ruhr West, Bottrop, Germany;
{andre.helgert, sabrina.eimler<EMAIL_ADDRESS>
###### Abstract

This paper investigates the potential of Virtual Reality (VR) as a research
tool for studying diversity and inclusion characteristics in the context of
human-robot interactions (HRI). Some exclusive advantages of using VR in HRI
are discussed, such as a controllable environment, the possibility to
manipulate the variables related to the robot and the human-robot interaction,
flexibility in the design of the robot and the environment, and advanced
measurement methods related e.g. to eye tracking and physiological data. At
the same time, the challenges of researching diversity and inclusion in HRI
are described, especially in accessibility, cyber sickness and bias when
developing VR-environments. Furthermore, solutions to these challenges are
being discussed to fully harness the benefits of VR for the studying of
diversity and inclusion.

##  I Introduction

Virtual Reality (VR) has established itself as a powerful research tool in
human-computer interaction through the use of unique capabilities and
immersive experiences [1]. This is not only due to the fact that VR can be
used to create dynamic and controllable environments, but also because
experimental studies have shown that VR-based assessments are just as
effective compared to traditional surveys and represent human behavior and
physiological patterns just as they do in real-world scenarios [2]. For this
reason, VR is also increasingly used as a research tool for human-robot
interactions (HRI) [3, 4]. Since communication between humans and robots can
be very different and individual due to various diversity characteristics, it
is important that robots can respond to them appropriately. As society becomes
more diverse, diversity-related challenges arise that increasingly affect us
in the private and public places. This development is not only challenging in
human-human interaction, but also affects all computer systems with which
humans interact directly, including robots. Therefore, robotic systems must be
able to respond to and interact with each individual, taking into account
individual diversity factors. In this paper, we will discuss the capabilities
of VR in the study of HRI, especially in terms of diversity and inclusion
research, with advantages in the field of controllability, manipulability,
flexibility and extended measuring methods as well as its challenges and
possible solutions.

##  II VR as a Research Tool

VR is an important subset of immersive technologies that is increasingly being
used in HCI and psychological studies. One reason is the possibility to
conduct studies which are difficult to perform or control in the real world.
Thus, in today’s world, VR appears as one of the most beneficial tools to
achieve effective results in the field of therapy and rehabilitation in
patients [5], such as in anxiety disorders [6] or in developing skills to deal
with pain [7]. Furthermore, the use of VR in the fields of educational [8],
social [9] and experimental psychology [10], among others, has produced
advantages in data collection methods and measurement data that are difficult
to implement in real-world scenarios.

VR is also already being used extensively in HRI, providing researchers with
new ways to explore and understand the complexity of human-robot interaction.
For example, Shariati et al. designed a VR-robot based on a real social robot,
which was specifically designed to improve learning, clinical therapy for
children with chronic diseases, and education [11]. The virtual robot was
compared with a real robot in an experimental study and the results indicate
that the acceptance of a VR-robot is the same as that of a real robot, and
that the virtual robot did not present a significant difference in performance
from the real robot. Shariati et al. conclude that the VR-platform has the
potential to have an important auxiliary solution for social robot research.
Another controlled study was conducted to investigate whether VR can be a
suitable platform for exploring social interactions between humans and robots
[12]. The quantitative results suggest that the core aspects of human-robot
social interactions are preserved in a VR-simulation. These examples show that
the medium VR is already used in HRI to produce valid research results. It is
therefore natural to explore VR in relation to diversity and inclusion in HRI.
Common advantages of VR applications can be applied to the diversity topic,
which will be explained in detail in the next chapter.

##  III Advantages in Exploring Diversity and Inclusion in VR

In addition to the countless applications opportunities of VR, it can also be
used to explore diversity and inclusion by providing exclusive characteristics
of VR, but also by extending common measurement methods. Since VR is used for
a wide variety of research purposes, the following characteristics are
presented in the context of HRI-studies.

###  III-A Controllable Environment

The use of VR as a research tool in HRI provides a significant advantage by
allowing studies to be conducted that are difficult to control in the real
world [13, 14]. Conducting studies in a controlled virtual environment helps
minimize or even eliminate confounding factors that can occur in field
studies. This has the advantage that research results can be generated without
the influence of external factors. In addition, every variable can be
controlled, from visual to auditory stimuli and even haptic impressions. This
level of controllability allows researchers to simulate and manipulate a
variety of scenarios to investigate specific questions. Specifically in
relation to diversity and inclusion, for example, social interactions can be
controlled in their nature of communications and adaptivity to diversity
factors. Investigations can be limited to specific situations, which may not
be as easy to conduct due to a non-controllable environment. This circumstance
ensures a high degree of standardization and reproducibility, which are
crucial aspects of controlled studies.

###  III-B Manipulability

A key advantage of using VR in HRI research is the high degree of
manipulability it offers. Unlike physical robots, VR-simulations can encompass
a wide range of robot variations, allowing researchers to explore various
factors such as robot appearance, behavior, and capabilities [15]. Researchers
can manipulate these variables to study how different robot attributes affect
human perception, attitude, and behavior. For example, they can study how
variations in the robot’s gender, ethnicity, or voice affect users’ trust,
engagement, and willingness to interact. Of course, this manipulability can
also be used to study diversity and inclusion characteristics, as different
people have different expectations and perceptions based on their background,
prior knowledge and so on. However, in addition to manipulating the robot, the
environment itself can also be dynamically adapted. Thus, the robot and the
interaction with humans can be tested and evaluated in different scenarios and
use cases. This allows possible conclusions to be drawn about a wide variety
of variables in the environment and an interaction with respect to the
perception of the HRI.

###  III-C Flexiblity

VR-environments provide researchers with the flexibility to tailor
interactions e.g. based on specific diversity factors [16] which are difficult
or impossible to implement in the real world. By integrating features such as
language preferences, cultural backgrounds, or even physical impairments,
researchers can create personalized experiences for participants. This
customization promotes inclusivity and enables the examination of human-robot
interactions across various demographic groups. This may help to counteract
discriminatory bias and minimize the reproduction of stereotypes by the robot.
In the area of accessibility, the use of VR offers a wide range of
possibilities, which can be explored and tested. Features related to
accessibility could be, for example, speech or gesture recognition or adaptive
interfaces for people with physical disabilities [17]. Understanding how
different individuals or communities interact with robots can inform the
design of more inclusive and culturally sensitive robotic systems. For
example, needs of people with physical limitations or other physical
characteristics such as different body sizes can be captured by simulating
these situations on their own game character in VR. Through these simple
technical adjustments, different diversity characteristics and their needs can
be identified, which could provide for a more inclusive experience with the
robot in future iterations.

###  III-D Extended Measuring Methods

VR is commonly used as a research method to study human behavior and cognition
because, when properly utilized, it provides an expanded data foundation
through specialized measurement methods [1]. There are several measurement
methods for researching human-robot interaction in VR. Common methods include
eye tracking [18] and physiological measurements [19]. For example, eye-
tracking technology can be used to collect data on a person’s eye movements
and what they are looking at and for how long. From this, behavior patterns
and focal points of attention can be identified. Physiological surveys could
unobtrusively collect data [20] on the person’s emotional state or measure
stress in interaction with the robot and cognitive load. With regard to
diversity factors, conclusions could be drawn from people with different
characteristics on the perception of the robot and the environment, which can
subsequently help to make the human-robot interaction more diversity-friendly.
Both methods can be an extension of conventional measurement methods and
reveal implications for HRI.

##  IV Challenges in Exploring Diversity and Inclusion in VR

In addition to the advantages of VR, which can represent a more diversified
exploration and design of different scenarios, the technology can also have an
exclusive or even discriminatory effect. In the following sub-chapters, these
challenges are discussed, in particular the accessibility problems, the cyber
sickness phenomenon and the bias that developers bring to VR-systems.

###  IV-A Accessibility

The VR-technology may not be accessible for all individuals, especially for
people with disabilities. Since conventional VR-devices consist of a headset
and two controllers, the use may work better or worse for different
limitations. For example, visual impairments affect the perception of the
virtual world and, depending on their severity, may be a criterion for
exclusion from VR. Hearing impairments can quickly make the virtual world more
difficult to use if the design is not inclusive [21]. Furthermore, limitations
in the mobility of the person can be a problem, since common VR-applications
provide for two controllers to move around in the virtual world [22].

###  IV-B Cyber Sickness

Cyber Sickness (CS) is a syndrome that can occur when using VR, when there is
a discrepancy between the visual information projected through VR-goggles and
the sensory information perceived by the body. This can lead to nausea and
dizziness and affects the VR-experience to a great extent. CS is considered a
major barrier to the acceptance and adoption of VR, as there is not yet a
fully comprehensive solution for the non-occurrence of the syndrome [23].
There are also gender differences in the occurrence of CS. For example, women
are more likely to experience cyber sickness when using VR-applications
because they have a different average IPD (the distance between the pupils of
both eyes) than men [24] and they cannot adjust their IPD on VR-goggles
because they only support the mean IPD of men [25].

###  IV-C Bias in VR-environments

Researching diversity and inclusion issues can be challenging because, for
example, prejudiced perceptions of developers and researchers and a non-
sensitive perception on other social groups can influence the VR-environment
itself. These prejudices could possibly be (also less) obviously reflected in
the development and design process [26]. Not only can this be potentially
discriminatory, but it can also stigmatize and exclude different social
groups. Algorithms based on artificial intelligence (e.g. robotic systems),
which are needed for a huge amount of projects in VR can be discriminatory
towards certain social groups. For example, a hiring algorithm from Amazon
systematically downgraded applications and resumes of female applicants based
on biased training data [27]. The training set was primarily male centered.
Similar discriminatory occurrences exist, for example, in facial recognition
[28] and in social communication with people [29].

##  V Discussion

VR offers various ways for researching diversity and inclusion in different
settings and szenarios, including HRI. Although we are provided with many new
and effective ways to recognize and implement diversity factors, the medium
itself can discriminate and could have an exclusive effect. As already
discussed, accessiblity, cyber sickness and a bias in VR-environments are a
few relevant variables. The question is how we can overcome these challenges
in order to exploit the opportunities of VR in terms of diversity and
inclusion-friendly implementations.

When looking at accessibility challenges, there is a very diverse and broad
issue. Restrictions on the person are as individual as the person himself.
However, special implementations can help to make VR-applications inclusive.
For example, there are already solutions in the field of vision impairments,
such as those of Zhao et al. who have developed a VR-toolkit, with which
various low vision limitations can be improved by the use of e.g. edge
measurement and depth enhancement [30]. For hearing impairments, the design of
the VR-environments could be adapted to have possible implemented speech
outputs put onto the user’s glasses via text. For limited mobility or motor
impairments it would be possible to use an adaptive controller as a
replacement for the two conventional controllers. An adaptive controller such
as the Xbox Adaptive Controller by Microsoft [31] can be adjusted and
individualized to various constraints. This allows people who are unable to
operate a conventional game controller to explore VR-environments and interact
with them. Only the virtual representation of the VR-controllers and thus the
hand tracking are not (yet) possible. However, this limitation can be
minimized by an inclusive design of the VR-environment.

Cyber sickness can be mitigated by ensuring a minimal latency between physical
movement in the real world and the corresponding virtual movement in the VR-
world. Additionally, optimizing movement settings to avoid abrupt movements
and make them feel more natural, like in real life, can help. For example, a
jerky motion when climbing stairs in unoptimized VR-applications. Approaches
to reduce cyber sickness are currently being investigated. Tian et al., for
instance, examined the seating position and whether the virtual coordinate
system contradicts the received real-world coordinates through our vestibular
system [32]. Promising results were achieved in reducing cyber sickness when
the real vertical axis aligned with the virtual one.

As many of the previous challenges can be attributed to a certain bias among
developers, it is also important to generate a certain awareness among
developers. This could not only address existing challenges in the area of
diversity and inclusion but also avoid bias in future systems. For this
reason, it is important to ensure a certain sensitivity to diversity factors
among researchers and developers, but also to constantly question one’s own
knowledge with regard to diversity.

##  VI Conclusion

VR has established itself as an effective research tool in the field of human-
computer interaction. In particular, VR shows its unique capabilities in
human-robot interaction research focusing on diversity and inclusion. Through
controllable environments, flexibility in customization of interactions and
manipulability of robot attributes, VR enables detailed investigation of the
complex interactions between humans and robots. Another advantage of VR lies
in the advanced measurement methods it offers. Eye-tracking and physiological
measurements can be used to capture behavioral patterns and cognitive
processes of humans during HRI studies. This allows for analysis of individual
perceptions and reactions to the robot and the environment. This helps to make
human-robot interactions more diverse and inclusive.

Nevertheless, there are also challenges in exploring diversity and inclusion
in HRI with VR. Biases and prejudices of developers and researchers can
negatively impact the design and development process and lead to
discrimination and exclusion of certain social groups. Therefore, it is
important for researchers and developers to be sensitive to diversity factors
and continuously challenge their steps regarding diversity. Furthermore, the
topics of accessibility and cyber sickness are still challenges today, where
there are already initial implications or solutions, which nevertheless
represent a hurdle. Addressing these challenges and promoting diversity and
inclusion in VR-based HRI research can lead to more equitable and unbiased
human-robot interactions. By harnessing the potential of VR, researchers can
help develop robotic systems that address individual diversity factors,
promote inclusion, and minimize discriminatory biases. This will ultimately
contribute to improved quality and effectiveness of human-robot interactions.

## References

  * [1] A. J. Martingano and S. Persky, “Virtual reality expands the toolkit for conducting health psychology research,” _Social and personality psychology compass_ , vol. 15, no. 7, 2021.
  * [2] A. C. Roberts, Y. W. Yeap, H. S. Seah, E. Chan, C.-K. Soh, and G. I. Christopoulos, “Assessing the suitability of virtual reality for psychological testing,” _Psychological assessment_ , vol. 31, no. 3, pp. 318–328, 2019.
  * [3] A. Arntz, S. C. Eimler, and H. U. Hoppe, “A virtual sandbox approach to studying the effect of augmented communication on human-robot collaboration,” _Frontiers in robotics and AI_ , vol. 8, p. 728961, 2021\.
  * [4] O. Liu, D. Rakita, B. Mutlu, and M. Gleicher, “Understanding human-robot interaction in virtual reality,” in _2017 26th IEEE International Symposium on Robot and Human Interactive Communication (RO-MAN)_. IEEE, 2017, pp. 751–757.
  * [5] A. A. Rizzo, M. Schultheis, K. A. Kerns, and C. Mateer, “Analysis of assets for virtual reality applications in neuropsychology,” _Neuropsychological Rehabilitation_ , vol. 14, no. 1-2, pp. 207–239, 2004\.
  * [6] J. L. Maples-Keller, B. E. Bunnell, S.-J. Kim, and B. O. Rothbaum, “The use of virtual reality technology in the treatment of anxiety and other psychiatric disorders,” _Harvard review of psychiatry_ , vol. 25, no. 3, pp. 103–113, 2017.
  * [7] E. Chan, S. Foster, R. Sambell, and P. Leong, “Clinical efficacy of virtual reality for acute procedural pain management: A systematic review and meta-analysis,” _PloS one_ , vol. 13, no. 7, p. e0200987, 2018.
  * [8] J. Parong and R. E. Mayer, “Learning science in immersive virtual reality,” _Journal of Educational Psychology_ , vol. 110, no. 6, pp. 785–797, 2018\.
  * [9] G. Carlo and B. A. Randall, “The development of a measure of prosocial behaviors for late adolescents,” _Journal of Youth and Adolescence_ , vol. 31, no. 1, pp. 31–44, 2002.
  * [10] M. Vasser, M. Kangsepp, K. Kilvits, T. Kivisik, and J. Aru, “Virtual reality toolbox for experimental psychology — research demo,” in _2015 IEEE Virtual Reality (VR)_. IEEE, 2015, pp. 361–362.
  * [11] A. Shariati, M. Shahab, A. Meghdari, A. Amoozandeh Nobaveh, R. Rafatnejad, and B. Mozafari, “Virtual reality social robot platform: A case study on arash social robot,” in _Social Robotics_ , ser. Lecture Notes in Computer Science, S. S. Ge, J.-J. Cabibihan, M. A. Salichs, E. Broadbent, H. He, A. R. Wagner, and Á. Castro-González, Eds. Cham: Springer International Publishing, 2018, vol. 11357, pp. 551–560.
  * [12] O. Sadka, J. Giron, D. Friedman, O. Zuckerman, and H. Erel, “Virtual-reality as a simulation tool for non-humanoid social robots,” in _Extended Abstracts of the 2020 CHI Conference on Human Factors in Computing Systems_ , R. Bernhaupt, F. F. Mueller, D. Verweij, J. Andres, J. McGrenere, A. Cockburn, I. Avellino, A. Goguey, P. Bjørn, S. Zhao, B. P. Samson, and R. Kocielnik, Eds. New York, NY, USA: ACM, 2020, pp. 1–9.
  * [13] V. Villani, B. Capelli, and L. Sabattini, “Use of virtual reality for the evaluation of human-robot interaction systems in complex scenarios,” in _2018 27th IEEE International Symposium on Robot and Human Interactive Communication (RO-MAN)_. IEEE, 2018, pp. 422–427.
  * [14] A. Todorov, C. P. Said, A. D. Engell, and N. N. Oosterhof, “Understanding evaluation of faces on social dimensions,” _Trends in cognitive sciences_ , vol. 12, no. 12, pp. 455–460, 2008.
  * [15] S. Negi, T. Arai, K. Inoue, Y. Ujiie, and T. Takubo, “Psychological assessment of humanoid robot appearance and performance using virtual reality,” in _RO-MAN 2008 - The 17th IEEE International Symposium on Robot and Human Interactive Communication_. IEEE, 2008, pp. 719–724.
  * [16] M. Bricken, “Virtual reality learning environments: potentials and challenges,” _ACM SIGGRAPH Computer Graphics_ , vol. 25, no. 3, pp. 178–184, 1991.
  * [17] M. Bhuiyan and R. Picking, “A gesture controlled user interface for inclusive design and evaluative study of its usability,” _Journal of Software Engineering and Applications_ , vol. 04, no. 09, pp. 513–521, 2011.
  * [18] U. J. Pfeiffer, K. Vogeley, and L. Schilbach, “From gaze cueing to dual eye-tracking: novel approaches to investigate the neural correlates of gaze in social interaction,” _Neuroscience and biobehavioral reviews_ , vol. 37, no. 10 Pt 2, pp. 2516–2528, 2013.
  * [19] A. Gaggioli, F. Pallavicini, L. Morganti, S. Serino, C. Scaratti, M. Briguglio, G. Crifaci, N. Vetrano, A. Giulintano, G. Bernava, G. Tartarisco, G. Pioggia, S. Raspelli, P. Cipresso, C. Vigna, A. Grassi, M. Baruffi, B. Wiederhold, and G. Riva, “Experiential virtual scenarios with real-time monitoring (interreality) for the management of psychological stress: a block randomized controlled trial,” _Journal of medical Internet research_ , vol. 16, no. 7, p. e167, 2014.
  * [20] A. Betella, D. Pacheco, R. Zucca, X. Arsiwalla, P. Omedas, A. lanatà, D. Mazzei, A. Tognetti, A. Greco, N. Carbonaro, J. Wagner, F. Lingenfelser, E. Andre, D. de rossi, and P. Verschure, “Interpreting psychophysiological states using unobtrusive wearable sensors in virtual reality,” 2014, pp. 331–336.
  * [21] D. Jain, S. Junuzovic, E. Ofek, M. Sinclair, J. R. Porter, C. Yoon, S. Machanavajhala, and M. Ringel Morris, “Towards sound accessibility in virtual reality,” in _Proceedings of the 2021 International Conference on Multimodal Interaction_ , Z. Hammal, C. Busso, C. Pelachaud, S. Oviatt, A. A. Salah, and G. Zhao, Eds. New York, NY, USA: ACM, 2021, pp. 80–91.
  * [22] M. Mott, J. Tang, S. Kane, E. Cutrell, and M. Ringel Morris, ““i just went into it assuming that i wouldn’t be able to have the full experience”,” in _The 22nd International ACM SIGACCESS Conference on Computers and Accessibility_ , T. Guerreiro, H. Nicolau, and K. Moffatt, Eds. New York, NY, USA: ACM, 2020, pp. 1–13.
  * [23] N. Dużmańska, P. Strojny, and A. Strojny, “Can simulator sickness be avoided? a review on temporal aspects of simulator sickness,” _Frontiers in psychology_ , vol. 9, p. 2132, 2018.
  * [24] H. C. Fledelius and M. Stubgaard, “Changes in eye position during growth and adult life as based on exophthalmometry, interpupillary distance, and orbital distance measurements,” _Acta ophthalmologica_ , vol. 64, no. 5, pp. 481–486, 1986.
  * [25] J. M. Fulvio, M. Ji, and B. Rokers, _Variations in visual sensitivity predict motion sickness in virtual reality_ , 2018.
  * [26] A. Nadeem, O. Marjanovic, and B. Abedin, “Gender bias in ai-based decision-making systems: a systematic literature review,” _Australasian Journal of Information Systems_ , vol. 26, 2022.
  * [27] Jeffrey Dastin, “Amazon scraps secret ai recruiting tool that showed bias against women,” in _Ethics of Data and Analytics_ , 2022.
  * [28] J. Buolamwini and T. Gebru, “Gender shades: Intersectional accuracy disparities in commercial gender classification,” in _Proceedings of the 1st Conference on Fairness, Accountability and Transparency_ , ser. Proceedings of Machine Learning Research, S. A. Friedler and C. Wilson, Eds., vol. 81. PMLR, 23–24 Feb 2018, pp. 77–91. [Online]. Available: https://proceedings.mlr.press/v81/buolamwini18a.html
  * [29] E. Hunt, “Tay, microsoft’s ai chatbot, gets a crash course in racism from twitter,” 2016. [Online]. Available: https://www.theguardian.com/technology/2016/mar/24/tay-microsofts-ai-chatbot-gets-a-crash-course-in-racism-from-twitter
  * [30] Y. Zhao, E. Cutrell, C. Holz, M. R. Morris, E. Ofek, and A. D. Wilson, “Seeingvr,” in _Proceedings of the 2019 CHI Conference on Human Factors in Computing Systems_ , S. Brewster, G. Fitzpatrick, A. Cox, and V. Kostakos, Eds. New York, NY, USA: ACM, 2019, pp. 1–14.
  * [31] C. Godineau, “The new xbox adaptive controller, another step towards digital inclusion.” in _Masters of Media_ , 2018. [Online]. Available: http://mastersofmedia.hum.uva.nl/blog/2018/09/23/the-new-xbox-adaptive-controller-another-step-towards-digital-inclusion/
  * [32] N. Tian, R. Clement, P. Lopes, and R. Boulic, “On the effect of the vertical axis alignment on cybersickness and game experience in a supine posture,” in _2020 IEEE Conference on Games (CoG)_. IEEE, 2020, pp. 359–366.


# Uncertainty aware Multitask Pyramid Vision Transformer for UAV-based Object
Re-Identification

###### Abstract

Object Re-IDentification (ReID), one of the most significant problems in
biometrics and surveillance systems, has been extensively studied by image
processing and computer vision communities in the past decades. Learning a
robust and discriminative feature representation is a crucial challenge for
object ReID. The problem is even more challenging in ReID based on Unmanned
Aerial Vehicle (UAV) as the images are characterized by continuously varying
camera parameters (e.g., view angle, altitude, etc.) of a flying drone. To
address this challenge, multiscale feature representation has been considered
to characterize images captured from UAV flying at different altitudes. In
this work, we propose a multitask learning approach, which employs a new
multiscale architecture without convolution, Pyramid Vision Transformer (PVT),
as the backbone for UAV-based object ReID. By uncertainty modeling of
intraclass variations, our proposed model can be jointly optimized using both
uncertainty-aware object ID and camera ID information. Experimental results
are reported on PRAI and VRAI, two ReID data sets from aerial surveillance, to
verify the effectiveness of our proposed approach.

Index Terms—  UAV-based object ReID, Pyramid Vision Transformer, Uncertainty
Modeling, Multitask Learning

##  1 Introduction

Object Re-IDentification (ReID) [1, 2], the task of matching a particular
object across different camera views, has been widely studied due to its
applications in visual surveillance, especially in the field of person and
vehicle tracking. Most of the existing work on object ReID is mainly focused
on tackling this problem in a normal surveillance domain, e.g. security
cameras installed on the top of a building. With the rapid development in the
UAV industry, visual surveillance using UAV devices has received increasing
attention, such as normal surveillance. However, ReIDs based on drones or UAVs
[3, 4] have remained an under-researched topic. Unlike the normal domain, the
ReID of UAV objects is arguably more challenging because drone-based images
often contain more uncertainty (e.g., view angles, camera distance, and
weather conditions) than standard surveillance images.

In the context of the ReID object based on UAVs, there exist several obstacles
arising from the increased uncertainty. As shown in Figure 1, both the scale
and pose variations are important factors to consider because the altitude
changes of the flying drones can be substantial. As the images captured are by
UAV flying at different altitudes, a ReID model operating on a single scale or
pose cannot create a descriptive feature space that fully characterizes the
entire aerial domain. Therefore, it is natural to develop a ReID model capable
of learning multiscale and pose-invariant features. Moreover, to detect subtle
changes of objects of interest (e.g., a person wearing different colored
dresses, different glasses, shoes, etc.), we target a discriminative feature
space that is generalizable to both coarse- and fine-grained features.

Fig. 1: Example images showing the need for multiscale feature fusion. The
images are taken by drones flying at multiple altitudes. A single scale can
only capture specific parts of objects, lacking the capability to learn fine
scale features. Multiscale features can help a model learn highly
discriminative feature space by fusing information across multiple scales.

Today, the convolutional neural network (CNN) is arguably the most popular
backbone in terms of designing ReID classifiers. However, the latest advances
in deep learning have seen the great success of Transformers [5, 6] in natural
language processing and computer vision. This class of convolution-free models
can be a better fit for the problem of object ReID using self-attention. Most
existing research on object (Persons, Vehicles) ReID has been conducted for
normal surveillance. Limited attention has been paid to object ReID in aerial
surveillance. Most approaches (e.g., part-based convolutional baseline [7],
generative adaptive alignment network [8], multiple granularity network [9])
follow CNN-based architectures. Transformer [10] is a new trend in the
solution of object ReID, establishing it as a strong baseline that beats
current state-of-the-art methods.

Motivated by the recent success of the transformer in normal surveillance
object ReID [1], we adopt an all-attention transformer-based approach for UAV-
based object Re-ID. To address the challenges arising from aerial images, we
propose an uncertainty-aware approach that exploits hierarchical feature maps
and global channel attention gate for object Re-identification in the UAV
domain. Our proposed approach uses a modified version of Pyramid Vision
Transformer (PVT) [11] as the backbone, which is a convolution-free
architecture. Then we apply spatial attention [12] to these multi-resolution
feature maps to put more focus on important features, filtering out irrelevant
ones. To incorporate camera information, we add additional head to the PVT
model and optimize the network using both camera ID, object ID, and center
loss in a joint fashion. Moreover, we normalize the style variance present in
different cameras, incorporating Batch Instance Normalization (BIN) [13] in
our model. Finally, with the help of the channel attention aggregation gate
[14], the model selectively learns the feature maps with higher weights.
Inspired by recent work on modeling feature uncertainty for personal Re-ID
[15], we propose to model aerial uncertainty by predicting the variance of
data as the model output.

PVT-based object recognition with uncertainty estimation is particularly
suitable for UAV-based Re-ID of people and vehicles in aerial surveillance and
long-range biometrics [16]. Our technical contributions can be summarized as
follows: (I) We seek to explore a modified Pyramid Vision Transformer (PVT)
tailored for object re-identification in UAV-based scenarios. The proposed
model utilizes multiscale features to re-identify objects for aerial
surveillance. (II) Spatial attention module helps focus on the relevant
information of the feature maps by filtering out noise. The channel attention
gate follows an adaptive fusion scheme, which dynamically selects the
appropriate feature maps to exploit channel-wise dependencies. (III) We train
our model considering uncertainty in terms of object identity and camera
identity for multitask learning. The model is regularized using Batch Instance
Normalization (BIN) to mitigate style variations in multiple cameras. (IV) Our
proposed framework achieves state-of-the-art performance on two aerial
surveillance datasets, PRAI-1581 [3] and VRAI [4], respectively.

##  2 Methodology

In this section, we present our multiscale approach for object Re-
Identification. An overview of our proposed method is outlined in Figure 2. We
propose a multi-task Pyramid Vision Transformer (PVT) [11], a convolution-free
backbone designed to learn multiscale feature maps. We use two heads for
object and camera ID recognition, respectively. Batch Instance Normalization
(BIN) is incorporated in our model to achieve camera style invariance. To make
feature maps spatially aware of the location of important objects, we apply
spatial attention [12] to feature maps of different resolutions. Finally, we
combine multiscale feature maps in an adaptive way using a Channel Attention
Gate [12, 14]. To make the identifier robust to occlusion, we estimate the
aleatory uncertainty[17] present in the data while computing the loss for the
model.

###  2.1 Multi-task Pyramid Vision Transformer

Fig. 2: The overall architecture of our proposed model. The model uses multi-
task PVT as the backbone. We apply Spatial Attention (SA) on feature maps of
individual scales which helps to focus the network on the most informative
features. The channel attention (CA) is a mini-network with shared wights. The
Identification (ID) head is responsible for classifying the object instances
whereas CamID head gives prediction on camera ID labels considering
uncertainty of prediction in both cases. Batch Instance Normalization (BIN)
layer helps the network to achieve style generalization among different
cameras.

Pyramid Vision Transformer (PVT) [11] generates hierarchical feature maps of
multiple resolutions. This architecture has four blocks. Each block is
responsible for generating a feature map of certain resolution. Each stage
consists of a patch embedding layer and a transformer encoder layer. Each
transformer encoder layer is composed of a modified attention layer named
spatial reduction attention (SRA) and a feedforward layer. SRA is designed to
reduce memory cost so that high-resolution feature maps can be processed.
First, the input image with a size of $H\times W\times 3$ is separated into
$\frac{HW}{4^{2}}$ patches and the patch size is $4\times 4$. Then, each patch
is flattened through a linear projection and passed through a number of
transformer encoders. After that, we get the output feature map of size
$\frac{H}{4}\times\frac{W}{4}\times C_{1}$. Similarly, we obtain the output
feature maps of sizes $\frac{H}{8}\times\frac{W}{8}\times C_{2}$,
$\frac{H}{16}\times\frac{W}{16}\times C_{3}$, and
$\frac{H}{32}\times\frac{W}{32}\times C_{4}$, respectively.

The original PVT uses non-overlapping patches. Inspired by the recent work of
TransReID [1], we have generated overlapping patches using a sliding window or
shifting operator. The head in PVT is used to identify objects. Differently
from the original PVT, we use an additional head for camera ID recognition. By
doing this, we try to take advantage of the camera ID labels to fuse camera
information into the original model. As the generated feature maps are
produced by multiple transformer encoders, the generated feature space is
complex, and it can be confusing for the model to extract meaningful features
for the ReID task. To obtain a more refined feature map so that the model can
learn more informative spatial features, we apply spatial attention to the
feature maps of different resolutions. In spatial attention, both MaxPooling
and AveragePooling are performed on the channel dimension, and the pooled
feature space is concatenated to generate the 2D spatial attention map.

Additionally, we have added Batch Instance Normalization (BIN) [13] to reduce
style variations on multiple cameras. The purpose of using BIN is to normalize
the style-preserving discriminative features of the ReID objectReID. To fully
exploit the functionality of multiscale features, we combine feature maps of
multiple resolutions using a global channel attention gate. To tune the
channel weights of multiscale feature maps in a dynamic fashion, we use a
shared channel attention gate to learn the inter-channel relationship within a
feature map regarding the importance of feature maps. We follow the design
procedure mentioned in [14] to implement this gate as a mini-network composed
of several global average pooling layers (GAP) and a multi-layer perceptron
(MLP) with reduced hidden dimension and one ReLU-activated hidden layer
followed by sigmoid activation.

###  2.2 Loss Function Formulation

To train our model, we use a combination of three uncertainty-aware losses:
identity loss,cameraID and center loss. Identity loss consists of softmax
cross-entropy loss and triplet loss [18, 19]. Camera ID information is learned
through centroid triplet loss [20] and center loss considers the distance from
center information.
Uncertainty-aware ID loss: This loss is computed using a hybrid of softmax
cross-entropy loss and triplet loss [21] taking into account the uncertainty
between classes. The uncertainty-aware softmax cross-entropy loss is
calculated by [18]:

$\mathcal{L}_{softmax}=\frac{1}{2N\sigma(x_{i})^{2}}\sum_{i=1}^{N}log(p_{id}(h_{i}^{id},y_{i}))+\frac{1}{2}{log\sigma(x_{i})^{2}}$
(1)

where $N$ is the number of samples, $y_{i}$ is the ground truth label, and
$\sigma(x_{i})^{2}$ is the variance of the data. We model our uncertainty-
aware soft-margin triplet loss using the following formula [19]:

$\mathcal{L}_{triplet}=\frac{1}{\sigma(x_{a})^{2}}log(1+exp(f(x_{a},x_{n})-f(x_{a},x_{p}))+\\\
\frac{1}{2}{log\sigma(x_{a})^{2}}$ (2)

where a triplet consists of $<x_{a},x_{p},x_{n}>$ in which $x_{a}$ is the
anchor image of a person, $x_{p}$ is the positive anchor image belonging to
the identity of the same person, and $x_{n}$ is the negative anchor image
belonging to a different person. Note that triplet loss cannot measure the
overall spatial distribution of features, while cross-entropy loss does not
have enough discriminant power among features. Therefore, it is better to
combine these two as follows:

$\mathcal{L}_{ua\\_id}={L}_{softmax}+{L}_{triplet}$ (3)

Uncertainty aware camera ID Loss: To tackle intraclass variations arising from
view angle, camera style, distance, etc., we apply a soft-margin version of
centroid triplet loss since the class centroid can be considered as the mean
representation for the retrieval task. Inspired by the unreasonable
effectiveness of centroids in image retrieval [20], we propose to calculate
uncertainty-aware camera ID loss based on centroids as follows:

$\mathcal{L}_{ua\\_camid}=\frac{1}{\sigma(x_{a})^{2}}log(1+exp(f(x_{a},c_{n})-f(x_{a},c_{p}))+\\\
\frac{1}{2}{log\sigma(x_{a})^{2}},$ (4)

where $c_{p}$ and $c_{n}$ are the corresponding centroids of the class for the
positive and negative classes.

Uncertainty-aware Center Loss: We also analyze uncertainty-aware center loss
[22] using the following formula:

$\mathcal{L}_{ua\\_center}=\frac{1}{2\sigma}\sum_{i=1}^{B}{\lvert\lvert{f}_{t_{i}}-{c_{y_{i}}}\rvert\rvert}^{2}_{2},$
(5)

where $y_{i}$ is the label of the $i^{t}h$ image in a mini-batch and B is the
batch size. $c_{y_{i}}$ is the center of deep features in the $y_{i}$th class.

The overall loss can be formulated as follows.

$\mathcal{L}_{total}=\alpha_{1}\mathcal{L}_{ua\\_id}+\alpha_{2}\mathcal{L}_{ua\\_camid}+\alpha_{3}\mathcal{L}_{ua\\_center}$
(6)

Here, $\alpha_{1}$, $\alpha_{2}$, and $\alpha_{3}$ are the regularization
parameters for the corresponding losses.

##  3 Experiments

###  3.1 Datasets

We have conducted our experiments on two aerial surveillance datasets named
Person ReID for Aerial Imagery (PRAI-1581) [3] dataset and Vehicle Re-
identification for Aerial Image (VRAI) [4] dataset. PRAI is a newly released
aerial surveillance dataset which contains 39,461 person images of 1581
classes captured by two UAV drones with a flight altitude ranging from 20 to
60 meters above the ground. VRAI dataset consists around 137,613 images of
13,022 vehicles taken by two UAV drones. This is the largest UAV based vehicle
dataset to date.

###  3.2 Implementation Details

For the PRAI dataset, the training set includes 19,523 images from 782
classes. For the test set, the number of query and gallery images are 4680 and
15258, respectively. For the VRAI dataset, the training set contains 66,113
images with 6,302 classes. For the test set, the query set contains 15,747
images and the gallery set contains 55,753 images, respectively. In our
experiment, we investigated PVT, a multiscale transformer, as the backbone
network. The backbone network is pre-trained on ImageNet 2012 dataset. We
train our model using 4 Titan 1080GTX GPUs. Before training, the images are
resized to $224\times 224$. The batch size is set to 128. ADAM optimizer is
used with a momentum of 0.9 and a weight decay of $1e^{-4}$. The learning rate
is initialized as 0.000015 with a cosine rate decay. For performance
evaluation, we use two metrics: Cumulative Matching Characteristic (CMC) and
mean Average Precision (mAP).

(a) Person

(b) PVT-Large

(c) Ours

Fig. 3: Visualization of feature maps using guided back-propagation. Baseline
PVT-Large [11] fails to retrieve fine features, while ours captures more
discriminative features.

####  3.2.1 Comparison with state-of-the-art

We compare our proposed approach with the latest methods, and the results are
reported in Table 1. Our proposed approach outperforms the previous state-of-
the-art in both the PRAI-1581 and VRAI dataset, respectively. For the VRAI
dataset, our approach achieves 84.47% Rank-1 accuracy and an mAP of 82.86%.
For the PRAI dataset, the accuracy of Rank-1 is 59.18% and the mAP is 51.45%.
We report results based on single-query settings for both datasets. It is
worth mentioning that the gain over previous SOTA methods is consistent across
object domains (person vs. vehicle) and performance metrics (Rank-1 vs. mAP).
The performances of our model for mAP and different rank scores are presented
in Figure 4 for the PRAI-1581 and VRAI data set. It can be observed that PRAI
is a lot more challenging than VRAI because of people’s relatively smaller
size, large pose variations, and deformable motion.

(a) CMC curve

(b) mAP

Fig. 4: Performance of our model for PRAI-1581 and VRAI. Table 1: Performance comparison with state-of-the-art methods for PRAI-1581 and VRAI. Method (Person ReID) | Rank-1 | mAP
---|---|---
PCB [7, 3] | 47.47 | 37.15
SVDNet [23, 3] | 46.10 | 36.70
MGN [9, 3] | 49.64 | 40.86
OSNet [14, 3] | 54.40 | 42.10
TransReID [1] | 56.30 | 49.81
Ours | 59.18 | 51.45
Method (Vehicle ReID) | Rank-1 | mAP
MGN [9, 4] | 67.84 | 69.49
RAM (ResNet-50) [24, 4] | 68.58 | 69.37
RAM (VGG-16) [24, 4] | 72.05 | 57.33
Multi-task+DP [4] | 80.30 | 78.63
TransReID [1] | 82.68 | 81.48
Ours | 84.47 | 82.86

##  4 Conclusion

We have presented an uncertainty-aware multiscale transformer-based approach
for the UAV-based object Re-ID. Our approach captures the information of
instances with different levels of detail by multitasking PVT-based backbone
architecture. The proposed model tries to solve object Re-ID as multitask
learning problem using a unified framework trained with object ID, camera ID,
and center loss. We quantitatively and qualitatively evaluated our proposed
method on two UAV-based aerial surveillance datasets. The experimental results
demonstrate the superiority of the proposed model over the previous state-of-
the-art.

## References

  * [1] Shuting He, Hao Luo, Pichao Wang, Fan Wang, Hao Li, and Wei Jiang,  “Transreid: Transformer-based object re-identification,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2021, pp. 15013–15022.
  * [2] A. Bedagkar-Gala and S. K Shah,  “A survey of approaches and trends in person re-identification,”  Image and vision computing, vol. 32, no. 4, pp. 270–286, 2014.
  * [3] Shizhou Zhang, Qi Zhang, Yifei Yang, Xing Wei, Peng Wang, Bingliang Jiao, and Yanning Zhang,  “Person re-identification in aerial imagery,”  IEEE Trans. on Multimedia, vol. 23, pp. 281–291, 2020.
  * [4] Peng Wang, Bingliang Jiao, Lu Yang, Yifei Yang, Shizhou Zhang, Wei Wei, and Yanning Zhang,  “Vehicle re-identification in aerial imagery: Dataset and approach,”  in ICVV, 2019, pp. 460–469.
  * [5] Ze et al. Liu,  “Swin transformer: Hierarchical vision transformer using shifted windows,”  arXiv preprint arXiv:2103.14030, 2021.
  * [6] Daquan et al. Zhou,  “Deepvit: Towards deeper vision transformer,”  arXiv preprint arXiv:2103.11886, 2021.
  * [7] Yifan Sun, Liang Zheng, Yi Yang, Qi Tian, and Shengjin Wang,  “Beyond part models: Person retrieval with refined part pooling (and a strong convolutional baseline),”  in Proceedings of the European conference on computer vision (ECCV), 2018, pp. 480–496.
  * [8] Guan’an Wang, Tianzhu Zhang, Jian Cheng, Si Liu, Yang Yang, and Zengguang Hou,  “Rgb-infrared cross-modality person re-identification via joint pixel and feature alignment,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 3623–3632.
  * [9] Guanshuo Wang, Yufeng Yuan, Xiong Chen, Jiwei Li, and Xi Zhou,  “Learning discriminative features with multiple granularities for person re-identification,”  in Proceedings of the 26th ACM international conference on Multimedia, 2018, pp. 274–282.
  * [10] Alexey et al. Dosovitskiy,  “An image is worth 16x16 words: Transformers for image recognition at scale,”  arXiv preprint arXiv:2010.11929, 2020.
  * [11] Wenhai et al. Wang,  “Pyramid vision transformer: A versatile backbone for dense prediction without convolutions,”  arXiv preprint arXiv:2102.12122, 2021.
  * [12] Sanghyun Woo, Jongchan Park, Joon-Young Lee, and In So Kweon,  “Cbam: Convolutional block attention module,”  in ECCV, 2018, pp. 3–19.
  * [13] Hyeonseob Nam and Hyo-Eun Kim,  “Batch-instance normalization for adaptively style-invariant neural networks,”  arXiv preprint arXiv:1805.07925, 2018.
  * [14] Kaiyang Zhou, Yongxin Yang, Andrea Cavallaro, and Tao Xiang,  “Omni-scale feature learning for person re-identification,”  in ICCV, 2019, pp. 3702–3712.
  * [15] Tianyuan Yu, Da Li, Yongxin Yang, Timothy M Hospedales, and Tao Xiang,  “Robust person re-identification by modelling feature uncertainty,”  in ICCV, 2019, pp. 552–561.
  * [16] Mang Ye, Jianbing Shen, Gaojie Lin, Tao Xiang, Ling Shao, and Steven CH Hoi,  “Deep learning for person re-identification: A survey and outlook,”  IEEE Trans. on PAMI, 2021.
  * [17] Guotai Wang, Wenqi Li, Michael Aertsen, Jan Deprest, Sébastien Ourselin, and Tom Vercauteren,  “Aleatoric uncertainty estimation with test-time augmentation for medical image segmentation with convolutional neural networks,”  Neurocomputing, vol. 338, pp. 34–45, 2019.
  * [18] Zhizheng Zhang, Cuiling Lan, Wenjun Zeng, Zhibo Chen, and Shif-Fu Chang,  “Rethinking classification loss designs for person re-identification with a unified view,”  arXiv preprint arXiv:2006.04991, 2020.
  * [19] Alexander Hermans, Lucas Beyer, and Bastian Leibe,  “In defense of the triplet loss for person re-identification,”  arXiv preprint arXiv:1703.07737, 2017.
  * [20] Mikolaj Wieczorek, Barbara Rychalska, and Jacek Dabrowski,  “On the unreasonable effectiveness of centroids in image retrieval,”  arXiv preprint arXiv:2104.13643, 2021.
  * [21] Zihao Hu, Huiyan Wu, Shengcai Liao, Hai-Miao Hu, Si Liu, and Bo Li,  “Person re-identification with hybrid loss and hard triplets mining,”  in 2018 IEEE Int. Conf. on Multimedia Big Data (BigMM), 2018, pp. 1–5.
  * [22] Hao Luo, Youzhi Gu, Xingyu Liao, Shenqi Lai, and Wei Jiang,  “Bag of tricks and a strong baseline for deep person re-identification,”  in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition workshops, 2019, pp. 0–0.
  * [23] Yifan Sun, Liang Zheng, Weijian Deng, and Shengjin Wang,  “Svdnet for pedestrian retrieval,”  in Proceedings of the IEEE International Conference on Computer Vision, 2017, pp. 3800–3808.
  * [24] Xiaobin Liu, Shiliang Zhang, Qingming Huang, and Wen Gao,  “Ram: a region-aware deep model for vehicle re-identification,”  in ICME, 2018, pp. 1–6.


# Finite-size effects on the ferroelectricity in rhombohedral HfO2

Francesco Delodovici  Université Paris-Saclay, CentraleSupélec, CNRS,
Laboratoire SPMS, 91190, Gif-sur-Yvette, France  Consiglio Nazionale delle
Ricerche CNR-SPIN, c/o Università degli Studi “G. D’Annunzio”, I-66100 Chieti,
Italy  Paolo Barone  Consiglio Nazionale delle Ricerche CNR-SPIN, Area della
Ricerca di Tor Vergata, Via del Fosso del Cavaliere, 100, I-00133 Rome, Italy
Silvia Picozzi  Consiglio Nazionale delle Ricerche CNR-SPIN, c/o Università
degli Studi “G. D’Annunzio”, I-66100 Chieti, Italy

###### Abstract

In this work we analyze the finite-size effects on the structural properties
and on the polarization of the rhombohedral phase of HfO2 subjected to a
biaxial compressive strain. We show how the presence of surface charges
affects the polarization, leading to a strong reduction with respect to its
bulk value. This reduction can be ascribed to two mechanisms: i) the coupling
between compressive strain and the phase-transition order parameter; ii) the
changes in the ferroelectric distortion. We give two alternative explanations
of this phenomenon: from an atomistic point of view, analyzing the evolution
of the bond lengths, and from a symmetry-analysis point of view, considering
the changes in the amplitude of the symmetry-allowed distortions, when a slab
configuration is considered. These results are independent on the slab-
thickness in the considered range, suggesting the absence of a critical
thickness for ferroelectricity in HfO2, in agreement with the proposed
improper nature of hafnia ferroelectricity.

##  I Introduction

Ferroelectricity in hafnium oxide has been at the center of continuous
interests and deep analysis in the last decade. The stabilization of different
crystalline phases represents a particularly delicate issue. In closer detail,
the detected polar phases [1, 4, 3, 2] of the pristine oxide [5] are only
metastable in bulk configurations at standard thermodynamical conditions, its
ground state rather being a non-polar monoclininc phase [6]. However, the
polar phases become thermodynamically competitive with the monoclinic phase
when hafnia is realized in the form of films [1]. In this respect, a variety
of different factors play a role, such as mechanical strain induced by the
substrate [7, 8, 9, 10], doping [7, 11], defects [13, 14, 12], grains size and
films thickness [15]. Despite many open issues, it is anyway beyond question
that its compatibility with Si-based semiconductors [16], along with the
robustness of the polarization hysteresis at the nanoscale level [17, 18] and
its simple chemistry, make HfO2 extremely appealing for ferroelectric-based
technological applications, when compared to ordinary ferroelectric
perovskite–like films [19].

From the theoretical modelling point of view, an extremely interesting topic
is represented by the investigation of how the order parameter describing
ferroelectricity couples to surfaces and interfaces [20, 21]. In closer
detail, it is particularly interesting to determine whether a minimal
thickness for the development of ferroelectric instabilities in the film
exists or not. The presence of such a critical thickness is generally related
to the screening of surfaces charges associated with the film spontaneous
polarization via depolarizing fields [21, 22, 23]. In proper ferroelectrics,
the occurrence of surfaces can have dramatic consequences on the stability of
the polar phase for films thinner than the critical thickness, which in the
case of prototypical BaTiO3 is on the nanometer order under short-circuit
condition [24, 22]. Ferroelectrics where the onset of polarization emerges as
a secondary effect, being mainly governed by the coupling to a primary non-
polar distortion, may be expected to present robustness against depolarizing
effects. This characteristic makes improper ferroelectrics particularly
appealing for their potentiality in technological applications. Nonetheless,
their behaviour in the ultra-thin limit remains unclear, as well as the
existence of a critical thickness [25, 26].

In this paper, we present a theoretical analysis, based on ab-initio
simulations, of the effects of finite-size on the stability of polar R3m phase
of hafnium oxide. We considered different slab configurations under open-
boundary conditions with clean surfaces in the range of experimentally-
reachable nanometric thicknesses. We stress that these electrostatic boundary
conditions are highly unfavourable to proper ferroelectricity. Indeed, when
polarization can be considered as the leading order parameter of the phase
transition, these conditions correspond to vanishing electric displacement in
the vacuum [21]. Under this configuration, below a critical thickness, the
intense depolarizing field experienced by the slab arising from the unscreened
surface charges is strong enough to completely suppress the polar phase.
Further on, we analyze the evolution of polarization as a function of the
depth in the slab (i.e. distance from the surfaces) by means of Born effective
charges of relaxed bulk phases. Finally, we justify the observed effects in
terms of microscopic modelling and in terms of suppression of the dominant
symmetry-allowed distortions connecting the high and the low-symmetry phase.

##  II Slab configurations

In this paragraph we report some structural information about the simulations
and refer to appendices for further computational details. We orient the
relaxed primitive cell according to the experimentally-detected growth
direction [3]. Thus, the [111] direction of the rhombohedral setting is
parallel to the normal of the film surface: we chose this direction to
coincide with the z axis in our simulations. In this configuration, the
polarization is also directed along z. Hereafter, we refer to calculations
performed in the equivalent hexagonal setting. The relaxed primitive vectors
are: $\rm{\bf a}=(a,0,0)$, $\rm{\bf b}=(-a/2,a\sqrt{3}/2,0)$, $\rm{\bf
c}=(0,0,c)$, where a$=7.134$ Å, c=$8.738$ Å. With this choice, the unit cell
consists of 12 formula units or, equivalently, of 36 atoms.

Previous computational modelling works [3, 27] show that the polarization
undergoes a rapid increase, when biaxial compressive strain is applied
normally to the polarization direction. In closer detail, the polarization of
the relaxed bulk configuration at equilibrium is approximately 0.06
$\mu$C/cm2, as computed with the Berry phase method [28], while higher values,
above 25 $\mu$C/cm2, are obtained when compressing a above 5%. Given this
significant difference in the P value, we expect finite size effects to be
“masked” in equilibrium conditions by the small polarization magnitude;
rather, we expect finite size effects to be more noticeable (and therefore
easier to investigate), upon application of compressive strain. Therefore, we
focus hereafter on the hexagonal cell subjected to a bi-axial 5% in-plane
compression. The bulk primitive cell under this strain condition is defined by
a=6.78 Å and copt=9.66 Å. The optimum copt is obtained by minimizing the total
energy of the primitive cell with fixed (compressed) in-plane lattice constant
and relaxed internal degrees of freedom describing atomic positions. The
polarization of this phase, computed with the Berry-phase method, is 22.8
$\mu$C/cm2, in good agreement with values reported in literature [27, 29]

In order to build the slab configurations, the compressed unit cell is
repeated along its [001] hexagonal direction and “cut” at three different
thicknesses. The resulting asymmetric slabs, before relaxation, are 12.88,
19.32, 38.64 Å thick, corresponding to 4 layers, 6 layers and 12 layers,
respectively. Hereafter, we define a “layer” as the stoichiometric combination
of Hf and O occupying one third of the primitive hexagonal cell volume. As
such, each layer consists of 4 Hf atoms and 8 O atoms. With this definition,
the 6 (12) layers slab consists of a primitive cell repeated 2 (4) times along
the c axis.

##  III RESULTS

###  III.1 Atomic displacements

Figure 1:  Top: the primitive cell of the rhombohedral phase defined by the
vectors aR,bR,cR. Hf (O) atoms shown by a yellow (red) sphere. Center: panel a
and b report a schematic representation of the ferroelectric distortion in
bulk and slab (central layer) configurations, respectively. The atomic
positions of the paraelectric phase are represented by the shaded spheres.
Bottom: structural relaxation along the c primitive vector for the 6-layers
slab (top panel) and the 12-layers slab (bottom panel). The black triangles
and the red circles show, respectively, the displacements along the $z$
direction of Hf and O symmetry-inequivalent types for each layer in the
relaxed slab with respect to the relaxed paraelectric R-3m. The dot-dashed
lines represent the ferroelectric displacements in the bulk configurations.
The atoms are labelled with different Wyckoff positions and multiplicity in
each panel to facilitate the comparison between the panels.

We observe the geometrical relaxation to produce two main effects in all the
considered configurations: i) the contraction of the film thickness; ii) the
relaxation of the atomic positions within the film. To give a quantitative
description of the first effect, we consider, as a marker for the contraction
experienced by the slabs, the thickness of the most central unit cell i.e. the
c axis of an embedded hexagonal cell. In fact, the central layers can be
reasonably assumed to be marginally affected by the surfaces. The comparison
between the Hf-projected DOS (not shown) in these central layers with the bulk
Hf-projected DOS indeed confirms such hypothesis. For the 4-layer-thick slab
the marker is ill-defined; thus, in order to have a value comparable to those
obtained for thicker slabs, we multiply the distance between the two central
layers by 3. The 12, 6 and 4-layer slab experience respectively a contraction
of $\approx$1.4%, $\approx$1.6%, $\approx$1.5%. This contraction can be traced
back to the coupling of the strain acting along c with the different symmetry-
allowed distortions (as will be detailed below, see section IV.2).

To analyze the relaxations of atomic positions in the film, we compare the
ferroelectric distortion in the bulk phases and in the relaxed slab. To this
purpose, we consider the displacements of each atom with respect to a
reference paraelectric (PE) phase: the $R\bar{3}m$. In a bulk configuration,
the polar displacements from a PE phase characterize the ferroelectric (FE)
phase. From the same quantities, in a slab configuration, one can deduce how
finite-size effects affect ferroelectricity. To decouple the displacements
from the effects of the cell compression, we adjust the reference PE phase to
have the same height of the central cell in the relaxed slab. Fig. 1
illustrates the vertical displacements $\rm\delta z=z_{slab}-z_{PE}$ of each
atom in the HfO2 slab. The top and bottom panel report respectively the
displacement of the 6-layers and of the 12-layers slabs. The dot-dashed lines
represent the displacements characterizing the bulk phases. The R3m bulk phase
has two different Wyckoff positions for hafnium and four for oxygen, thus 6
lines appear in the panels (two oxygen lines at -0.2, in red, are almost
degenerate, while a third one close to 0 overlaps with a Hf line in black) .
Each of this high-symmetry position transforms in a unique way under geometric
relaxation, therefore each layer presents 4 different oxygen displacements and
2 for hafnium. The surface layers experience the largest changes: the absolute
displacements of hafnium at the bottom (first) layer change on average by 0.06
Å , the oxygens closer to vacuum on average by 0.07 Å whereas those toward the
bulk change by 0.08 Å . At the top layer, the hafnium displacements change by
0.16 Å on average, while those of oxygens remain closer to the bulk value,
changing on average by 0.03 Å for O toward the bulk and of 0.04 Å for O toward
the vacuum. We stress that the four hafnium atoms at the top surface in the
relaxed slab have the same $z$ coordinate. This is analogous to what happens
in the PE bulk phase. The changes in the ferroelectric distortion rapidly
converge to constant values in the central layers, with the second layer
acting as an intermediate buffer. In the bulk region the displacements of both
hafnium and oxygen atoms do not match the displacements characterizing the
bulk phases: the discrepancy for oxygen reaches 0.06 Å and 0.03 Å for hafnium.
Based on these observations, it is reasonable to expect some changes in the
electric polarization throughout the slab. Nonetheless, the complexity of the
distortions makes these effects hard to deduce only from the evolution of the
ferroelectric displacements. Therefore, we move to a direct analysis of the
polarization throughout the slab.

###  III.2 Polarization

Figure 2:  The bulk polarization computed with different methods compared to
the layer-projected value. Three slab thickness are reported; top: four–layers
slab; center: six–layers slab; bottom: twelve–layers slab. The dotted-dashed
lines represent the P computed with the Berry-phase method. The continuous
lines report the P computed using the Born-effective charges. The blue line
corresponds to the P of the bulk phase with compressed c. The triangles report
the value of P per each layer.

In order to investigate the evolution of the polarization in the slab as a
function of thickness, we employed the Born effective charges, computed for
the relaxed and for the compressed bulk phases. We do not employ the Berry-
phase method, due to the presence of metallic states confined at the top
surface and associated with oxygen $p$-orbitals, as revealed by the computed
projected density of states. The polarization difference with respect to the
centrosymmetric phase is by definition the contraction of the tensor Z∗ with
the vector of the atomic displacements from the $R\bar{3}m$ phase, normalized
to the bulk volume. For the layer polarization, the dipole moment is
normalized to the volume of the layer, hereby defined as one third of the
volume of the adjusted bulk. These simplified definitions allow us to
investigate the dependence of the dipole-moment density on the “depth” (i.e.
distance from the surfaces) within the slab.

We decouple the effects contributing to the change in polarization: the slab
contraction and the local atomic distortions. Fig. 2 reports the evolution of
the polarization as a function of the layer for three different slab
thicknesses. The polarization of different bulk configurations obtained with
the BP method is reported as a comparison: 0.07 and 22.8 $\mu$C/cm2
respectively for the relaxed bulk and the bulk subjected to a biaxial (along
a,b primitive vectors) compressive (-5%) strain. In this regard, it should be
noticed that the use of effective charges to compute the polarization results
in an underestimate with respect to the BP value: P computed with effective
charges amounts to 0.06 and 20.03 $\mu$C/cm2, respectively, for the relaxed
equilibrium bulk and for the bulk subjected to a biaxial compressive (-5%)
strain. This discrepancy in the compressed case (i.e. almost 2.8 $\mu$C/cm2)
corresponds to 9% of the BP value, which is still a reasonable margin of error
for our analysis.

We stress once again that we employed the effective charges of the bulk
subjected to a biaxial -5% strain state to compute the effects of finite size
on the polarization. This corresponds to assuming that the local distortions
following the relaxation do not affect significantly the atomic responses to
an external electric field. In order to analyze the effects of the slab
contraction, we compute the change in P arising from the compression along c
of the relaxed R3m bulk to the value of c characterising the central hexagonal
cell embedded in the 6 layers slab. In this way, we observe that the slab
contraction induces a decrease in P by about 33% (6.7 $\mu$C/cm2) with respect
to the c-relaxed bulk, as can be observed in Fig.2: the polarization of the
compressed bulk, 13.33 $\mu$C/cm2, and the one of the c-relaxed bulk, 20.03
$\mu$C/cm2, are respectively represented by the blue and the continuous red
line. Such a reduction of the polarization induced by the slab contraction can
be phenomenologically traced back to the same strain-polarization coupling
which is responsible for the remarkable increase of P following the 5% bi-
axial in-plane compression of the bulk configuration and consequent elongation
of the cell along the out-of-plane direction. As the slab contraction causes a
reduction of the strain deformation along the $c$ axis with respect to the
bulk phase, the strain-polarization coupling is responsible for a
corresponding reduction of P. The introduction of finite-size causing the
local rearrangement of atoms induces a further significant reduction of P. In
fact, when moving from the surfaces where the polarization assumes opposite
sign and different absolute values, the polarization settles to approximately
4.5 and to 3.3 $\mu$C/cm2, respectively, in the 12(6) and 4 layers slabs. The
reduced value of P in the thinnest slab is connected with the lack of a proper
bulk-like volume in the slab, as in this configuration the central layers are
just below surfaces. In all the slabs, the polarization in the central layers
is therefore reduced by more than 70% with respect to the corresponding value
in the uncompressed bulk. Nonetheless, it remains finite. We note that the
bulk Z∗ represent a poor approximation to the electric response at the
surfaces. Therefore, the value of P at the extreme layers is likely biased and
should not be commented.

It is interesting to notice that, despite the open-boundary conditions
employed to simulate the slab, the depolarization field does not completely
suppress the polarization, as it happens for proper-ferroelectrics slabs below
their critical thickness under open-circuit conditions [21]. Since the
polarization of the central layers remains finite as the thickness is lowered,
we may infer that the critical thickness for the $R3m$ phase of HfO2 films is
either smaller than the thickness of the 4-layer slab, or it completely
vanishes. We note that the extreme case of a three layer configuration still
presents a finite value of P in the central layer ($\approx$ 2 $\mu$C/cm2),
but this can hardly be intended as a “proper bulk–like” layer. Therefore, the
incomplete suppression of P may be an indirect sign of the improper character
of ferroelectricity in this polar phase of HfO2.

##  IV Interpretation of polarization suppression

###  IV.1 Atomistic model

Figure 3:  Section of the $R3m$ hexagonal unit cell. The primary and the
secondary (vertical) Hf-O bonds are highlighted respectively by green and blue
rectangles. In the unit cell, these bonds appear with a ratio of 3 to 1. A
shorter vertical bond corresponds to both: above for the primary and below for
the secondary bonds, respectively.

Figure 4:  The trend of the long-short primary (3x) and the long-short
secondary (1x) bonds in the 12-layer slab. The dotted-dashed lines represents
the corresponding bulk value. Each color has to be compared with the closest
black line that represents the corresponding bulk bond lengths. The red line
corresponds to the value of the bonds in the $R-3m$ paraelectric phase. The
values for the two surfaces are not reported, since they reflect a different
behaviour occurring at the surfaces.

To explain the suppression of P in the central layers, we refer to the
atomistic model of Zhang et al. reported in [27], where the effects of
compressive strain on the $R3m$ phase are analyzed. The authors show that the
increase of P under certain strain states corresponds to the change in the
contribution to P of two competing Hf-O bonds, hereafter named “primary” and
“secondary”, appearing in the unit cell with the ratio of 3 to 1. The primary
bonds involve the Hf and O with high multiplicity (9b, corresponding to Hf1,
O1 and O2 as reported in Fig. 1). Instead, the secondary bonds involve the Hf
and O with low multiplicity (3a, corresponding to Hf2, O3 and O4 in Fig. 1).
Fig. 3 reports the hexagonal unit cell with a primary bond highlighted by the
green rectangle and the secondary bond by the blue one. For small strains, the
contribution to the polarization of the two bonds compensate one another.
However, above a certain strain threshold, the contribution from the 3 primary
bonds grows and determines the increase of P. Fig. 4 reports the trend of the
two bonds lengths across the 12-layers slab. The secondary and the primary
bonds respectively increase and decrease with respect to their bulk value.
Specifically, the primary ones approach the value of the bonds in the centro-
symmetric phase. Since the two bonds compete in determining the sign of the
polarization, the resulting effect of such a global distortion is to reduce
the polarization in the central portion of the slab.

###  IV.2 Suppression of the main distortion modes

Figure 5:  Symmetry modes analysis connecting the high temperature cubic phase
to the ferroelectric R3m. The top panel reports the ab-initio total energy for
the dominant single-modes as a function of the generalized displacement
coordinate $\lambda$. The central panel reports the energy trend of the
coupled-dominant modes. The bottom panel reports the dependence of the leading
order parameter X2- and of the coupling X2- U $\Gamma_{4-}$ U X5- on the
generalized coordinate, in the case of 0% and -5% compressive strain. In the
bottom panel the configurations at $\lambda=0$ are artificially set to zero to
simplify the comparison: the -5% undistorted phase is more than 220 meV/f.u.
higher in energy with respect to the non-compressed phase.

Figure 6:  The evolution of the structural fingerprint for three dominant
modes in the 6-layer slab. The first panel reports the fingerprint for the
leading order parameter. The angle $\alpha$ is the one connecting the O2 and
O4 with the vertex at Hf in the structure above the panels at left. The
central panel reports the fingerprint for the polarization. The panel at
bottom reports the fingerprint for X5- mode. It represents the z-component of
the distance between Hf atoms with different Wyckoff positions which are
reported in the picture above the panels at right.

An alternative approach relies on the analysis of the main symmetry-allowed
distortion modes connecting the hafnia high-temperature cubic phase ($Fm-3m$)
to the $R3m$ phase. We employed standard crystallographic tools [30, 31] to
identify the symmetry-adapted distortions and found six of them, two of which
have negligible amplitudes. Five symmetry modes are associated with the wave
vector qX at the Brillouin-zone boundary, and the last one to a polar
distortion at $\Gamma$. Under -5% strain condition, the four dominant modes
(in the hexagonal setting), ordered by decreasing amplitude, are: X2-,
X5-,$\Gamma_{4-}$ and X3-. The top panel in Fig. 5 reports the ab-initio total
energy computed along the distortions. The only soft mode is X2- corresponding
to an instability of 8 meV/f.u. On the other hand, the high-symmetry phase is
stable against all the remaining modes. This is a canonical behaviour for an
improper ferroelectric, where the polarization emerges through the coupling to
the leading order parameter, which in this case is a zone-boundary distortion.
When examining the atomic distortions related to each mode, we note that the
leading order mode distorts the angles formed by a Hf atom with the
surrounding oxygen cage and changes the relative bond lengths. The X5- mode
moves apart, along the c direction, the Hf with different Wyckoff positions.
These correspond to the Hf atoms forming the primary and the secondary bonds
in Fig.3. The polar mode $\Gamma_{4-}$ moves apart the oxygen and the hafnium
atoms along the c direction. The X3- mode displaces the Hf forming the primary
bonds along a non-trivial direction and those forming the secondary bonds
simply along c. Going back to Fig.5, the central panel reports the couplings
between the different modes. The coupling of the polarization with X2- does
not significantly affect the soft character of the leading mode: the energy
gain is 7 meV/f.u. . Instead, the X5- significantly enhances the soft
character of the leading mode, the energy gain when coupled to the order
parameter amounting to 54 meV/f.u. . It is worth noticing that including the
$\Gamma_{4-}$ mode further stabilize this coupling, even though only by 9
meV/f.u. . In this sense, we may see the coupling of the X5- mode with the
order parameter as a condition favouring the excitation of the polar
distortion. When all the four dominant modes are coupled, the largest gain in
energy is obtained: 128 meV/f.u. below the high symmetry phase at $\lambda=0$.
The full relaxation, which includes the neglected structural distortions,
leads to a further energy decrease of 20 meV/f.u., i.e. less than one sixth of
the energy gain corresponding to the simultaneous activation of the dominant
modes. The bottom panel in Fig.5 describes the dependence of the energy
dispersion of X2- and of the coupling X2- U $\Gamma_{4-}$ U X5- on the applied
strain. The compression to -5% weakens the character of the leading mode: its
stability is reduced by almost 30 meV/f.u. . Instead, the three-modes-coupled
distortion is strongly stabilized by the compression: when strain is applied,
the energy gain is about 27 meV/f.u. . We stress that this energetically
favourable pairing with the compressive strain is not shown by the X2- U
$\Gamma_{4-}$ distortion. This is a further sign of the special role assumed
by X5- when compressive strain is applied. Based on this analysis, we can
describe the suppression of the polarization in the central portion of the
slab in terms of suppression of the dominant modes. To this purpose, we define
a structural “fingerprint” for each mode, which consists in geometrical
quantities catching unambiguously the distortions of each single mode. From
their trends throughout the slab, we can indeed deduce the variations of the
modes amplitude. The fingerprint describing the order parameter X2- can be
taken as the variation of the angle $\alpha$ (between the oxygen labelled O2
and O4 in Fig.6 with the vertex located at the Hf atom), with respect to the
reference value of $\pi$ (i.e. the $\alpha$ value in the PE phase at
$\lambda=0$). This angle is also affected by $\Gamma_{4-}$, but the change
deriving from this mode is much smaller (only 26%) than the one from X2-. The
other angles in the cage are significantly affected also by the mode X3-: the
distortion coming from this mode is almost 80% of the one coming from X2-.
Therefore, they cannot be taken as fingerprints for the order parameter. The
polar distortion can be easily described in terms of the distance between the
oxygens and the Hf average z-coordinates. The effects of the other dominant
modes on this quantity is negligible. The X5- mode can be well represented by
the finite distance between the z-coordinates of the Hf atoms having different
Wyckoff positions. Specifically, we averaged the z-coordinate of the Hf with
multiplicity 3.

Fig.6 reports the trend of the three structural fingerprints for the different
layers in the 6-layer slab configuration. The first panel reports the
difference of the angle $\alpha-\alpha_{\rm para}$ representing X2- and its
bulk value ($\approx$ 12°) as a function of the layer. The reduction of about
30% with respect to the bulk value, occurring in the central layers, is a
clear indication of the mode suppression. The central panel reports the trend
for the fingerprint of the polar mode. Also in this case we observe a clear
reduction of the fingerprint in the central layers. Given the detected
improper character of ferroelectricity, we note that the suppression of the
polarization should be interpreted as caused by the suppression of the X2-
mode. This dependence can also be inferred by the strong correlation between
the trends of the two fingerprints: their suppression should not be intended
as two independent phenomena. The third panel illustrates that the X5-
distortion is not suppressed as the depth increases. On the contrary, its
amplitude oscillates around the bulk value. This is interesting to notice,
given the considerable contribution to the stabilization of the polar phase
coming from the coupling of this mode with the order parameter. This suggests
that finite-size effects reduce the amplitude of the order parameter in such a
way to lower the polarization in the central volume, while minimizing the
reduction of energetic stability with respect to the pure bulk phase. Given
the connection with the polar mode reported in Fig.5, this result also
explains the incomplete suppression of the polarization. In fact, the finite
size has no effect at all on the X5- mode, which is responsible for
“encouraging” the excitation of $\Gamma_{4-}$. Therefore, the role of X5- when
coupled with the order parameter is two-fold even in finite configuration: it
contributes to the stabilization of the ferroelectric phase and it prevents
the polar mode from completely vanishing.

##  V Conclusions

In this paper we addressed the effects of surfaces on the $R3m$ polar phase of
hafnia. We have shown how the introduction of the finite size leads to two
types of consequences. The first one consists in the contraction of the slab,
resulting from the coupling between strain and the order parameter. This
contraction leads to a reduction of polarization by about one third compared
to the bulk value. This effect should diminish with increasing thickness, but
the thicknesses range that we analyzed does not allow us to observe this trend
clearly. The second consequence is the relaxation of the atomic positions
towards a configuration that further decreases the electric dipole density in
the central layers. The values achieved are finite, but much smaller than the
corresponding bulk values. We explained these effects by means of two
approaches: through an atomistic perspective, looking at bond lengths to
identify polarization, and through a more fundamental approach based on
symmetry-allowed distortion analysis. The trend in bond lengths throughout the
slab clearly shows that the presence of the surfaces tends to balance the
weight of different competing bonds that contribute with opposite sign to
layer-polarization. This balance does not lead to full cancellation though,
thus polarization is not completely suppressed. The analysis of the
distortions allowed us to identify the improper nature of ferroelectricity in
hafnium oxide rhombohedral phase, the leading-order parameter being the mode
at the boundary-zone X2- instead of P. By assigning a structural fingerprint
to each of the leading modes, we showed how the introduction of the surfaces
leads to a suppression of the leading mode and, consequently, of the
polarization in the central layers of the slab. However, even in this case the
suppression is not complete and the amplitude of the distortions is not fully
cancelled. Finally, we found that these results are independent on the slab-
thickness in the considered range (1-3.5 nm). Therefore, we may infer the
critical-thickness for the R3m phase of HfO2 either to be lower than 1 nm or
to completely vanish.

##  VI Acknowledgments

We acknowledge support from the Italian Ministry for Research and Education
through the PRIN-2017 project “TWEET: Towards Ferroelectricity in Two
dimensions” (IT-MIUR Grant No. 2017YCTB59) and computational support from
CINECA through the ISCRA initiative with Grant N. HP10CEI2UQ (FeCoSMO
project). We also acknowledge support from the computational resources at
“Gabriele d’Annunzio” University of Chieti.

##  Appendix A Computational details

We performed density functional theory simulations, as implemented in the
Vienna Ab-initio Simulation Package (VASP) [32] at the level of the revised
Perdew, Burke, and Ernzerhof functional for solids [33, 34]. The bulk phases
were relaxed until the difference in the total energy between two successive
self-consistent steps is smaller than $10^{-7}$ eV and forces are smaller than
$10^{-3}$ eV$\rm/\AA$. The same convergence criteria apply for the 4-layers
and 6-layers slabs. Instead, the 12-layers slab was relaxed until all the
forces are smaller than than $10^{-2}$ eV$\rm/\AA$. The cutoff for the
expansion on the plane waves basis is set to 600 eV. The k-point mesh is set
to 5x5x3 for bulk calculations and to 4x4x1 for the slabs.

We relax the slabs with at least 30 Å of vacuum including dipole corrections
to prevent periodic-images spurious interaction. No external field is applied.

In the distortion modes analysis, we compute the energy of 10 intermediate
structures along the distortion paths, with the same computational parameters
reported above.

## References

  * [1] Böscke, T. S., Müller, J., Bräuhaus, D., Schröder, U. , Böttger, U. , Ferroelectricity in hafnium oxide thin films. _Appl. Phys. Lett._ 99, 102903 (2011).
  * [2] Barabash, S. Prediction of new metastable HfO2 phases:toward understanding ferro-and antiferroelectric films._J. Comput.Electron._ 16, 1227–1235 (2017).
  * [3] Wei, Y., Nukala, P., Salverda, M., Matzen, S., Zhao, H. J.,Momand, J., Everhardt, A. S., Agnus, G., Blake, G. R.,Lecoeur, P.,_et al._. A rhombohedral ferroelectric phase inepitaxially strained Hf${}_{0}.5$Zr${}_{0}.5$O2 thin films. _Naturemater._ 17, 1095–1100 (2018).
  * [4] Lee, H.J., Lee, M., Lee, K., Jo, J., Yang, H., Kim, Y., Chae, S. C.,Waghmare, U. Lee, J. H. , Scale-free ferroelectricity induced by flat phonon bands in HfO2. _Science_ 369, 1343–1347 (2020).
  * [5] Xu, X., Huang, F.T., Qi, Y., Singh, S., Rabe, K. M.,Obeysekera, D., Yang, J., Chu, M.W. Cheong, S.W. , Kinetically-stabilized ferroelectricity in bulk single-crystalline HfO2:Y. _Nature mater._ , (2021)
  * [6] Huan, T. D., Sharma, V., Rossetti, G. A., Ramprasad, R. , Pathways towards ferroelectricity in hafnia. _Phys. Rev. B_ 90, 064111–064116 (2014).
  * [7] Batra, R., Huan, T. D., Jones, J. L., Rossetti, G. , Ramprasad, R. , Factors Favoring Ferroelectricity in Hafnia: A First-Principles Computational Study. _J. Phys. Chem. C_ 121, 4139–4145 (2017).
  * [8] Estandìa, S., Dix, N., Gazquez, J., Fina, I., Lyu, J.,Chisholm, M. F., Fontcuberta, J. Sánchez, F. , Engineering Ferroelectric Hf${}_{0}.5$Zr${}_{0}.5$O2 Thin Films by Epitaxial Stress. _ACS Appl. Electron. Mater._ 1, 1449–1457 (2019).
  * [9] Estandìa, S., Dix, N., Chisholm, M. F., Fina, I. Sánchez, F. , Domain-Matching Epitaxy of Ferroelectric Hf${}_{0}.5$Zr${}_{0}.5$O2(111) on La2/3Sr1/3MnO3(001). _Cryst. Growth Des._ 20, 3801–3806 (2020).
  * [10] Liu, S., Hanrahan, B. M., Effects of growth orientations and epitaxial strains on phase stability of HfO2 thin films. _Phys. Rev. Materials_ 3, 054404–054413 (2019).
  * [11] Schröder, U., Richter, C., Park, M. H., Schenk, T., Pesić, M., Hoffmann, M., Fengler, F. P. G., Pohl, D., Rellinghaus, B., Zhou, C., _et al._. Lanthanum-Doped Hafnium Oxide: A Robust Ferroelectric Material. _Inorg. Chem._ 57, 2752–2765 (2018).
  * [12] Nukala, P. , Ahmadi, M. , Wei, Y., de Graaf, S. , Stylianidis, E., Chakrabortty, T., Matzen, S. , Zandbergen, H.W. , Björling, A. , Mannix , D. , Carbone, D. , Kooi, B., Noheda, B., Reversible oxygen migration and phase transitions in hafnia-based ferroelectric devices. _Science_ , 630-635 (2021)
  * [13] Chouprik, A., Negrov, D. , Tsymbal, E.Y. and Zenkevich, A., Defects in ferroelectric HfO2. _Nanoscale_ , 21, 11635-11678 (2021).
  * [14] Hoffmann, M., Schröder, U., Schenk, T., Shimizu, T.,Funakubo, H., Sakata, O., Pohl, D., Drescher, M., Adelmann, C.,Materlik, R., _et al._. Stabilizing the ferroelectric phase indoped hafnium oxide. _J. Appl. Phys._ 118, 072006 (2015)
  * [15] Polakowski, P. , Müller, J., Ferroelectricity in undoped hafnium oxide. _Appl. Phys. Lett._ 106, 232905 (2015).
  * [16] Chernikova, A., Kozodaev, M., Markeev, A., Negrov, D., Spiridonov, M., Zarubin, S., Zenkevich, A., Ultrathin Hf0.5Zr0.5O2 ferroelectric films on Si. _ACS Appl. Mater. Interfaces_ , 8, 7232-7237, (2016).
  * [17] Yurchuk, E., Müller, J., Knebel, S., Sundqvist, J., Graham, A. P.,Melde, T., Schröder, U. Mikolajick, T. Impact of layer thickness on the ferroelectric behaviour of silicon doped hafnium oxide thin films. _Thin Solid Films_ 533, 88–92 (2013).
  * [18] Dawber, M. , Rabe, K. M., Scott, J. F. , Physics of thin-film ferroelectric oxides. _Rev. Mod. Phys._ , 77, 1083-1130 (2005).
  * [19] Müller, J. ,Polakowski, P. Mueller, S. and Mikolajick, T. , Ferroelectric Hafnium Oxide Based Materials and Devices: Assessment of Current Status and Future Prospects. _ECS J. Solid State Sci. Technol._ 4, N30-N35 (2015)
  * [20] Eliseev, E.A., Kalinin,S.V., and Morozovska A.N. , Finite size effects in ferroelectric-semiconductor thin films under open-circuit electric boundary conditions, J. Appl. Phys, 117, 034102, (2015).
  * [21] , Meyer B., and Vanderbilt D. Ab-initio study of BaTiO3 and PbTiO3 surfaces in external electric fields. _Phys. Rev. B_ 63, 205426 (2001).
  * [22] Junquera, J. Ghosez, P. Critical thickness for ferroelectricity in perovskite ultrathin films. _Nature_ 422, 506–509 (2003).
  * [23] Batra, I. P., Wurfel, P., Silverman, B. D., Phase Transition,Stability, and Depolarization Field in Ferroelectric Thin Films. _Phys. Rev. B_ 8, 3257-3265 (1973).
  * [24] Kim,Y. S. , Kim,D. H. , Kim,J. D. , Chang,Y. J., Noh,T. W. , Kong,J. H. , Char,K. , Park,Y. D. , Bu,S. D. , Yoon,J.-G. and Chung,J.-S., Critical thickness of ultrathin ferroelectric BaTiO3. films, _Appl. Phys. Lett._ , 86, 102907, (2005).
  * [25] Sai, N. , Fennie, C. J. and Demkov, A. A., Absence of Critical Thickness in an Ultrathin Improper Ferroelectric Film. _Phys. Rev. Lett._ , 102, 107601, (2009)
  * [26] Nordlander, J., Campanini, M., Rossell, M.D., Erni, R., Meier, Q.N., Cano, A. Spaldin, N.A., Fiebig, M. , Trassin, M. , The ultrathin limit of improper ferroelectricity. _Nat. Commun._ 10, 5591 (2019).
  * [27] Zhang, Y., Yang, Q., Tao, L., Tsymbal, E. Y., Alexandrov, V. Effects of Strain and Film Thickness on the Stability of the Rhombohedral Phase of Hf${}_{0}.5$Zr${}_{0}.5$O2 . _Phys. Rev. Applied_ 14, 014068-014076 (2020).
  * [28] King-Smith, R.D. , Vanderbilt, D. , Theory of polarization of crystalline solids. _Phys. Rev. B_ 47, 1651-1654 (1993).
  * [29] Nukala, P., Antoja-Lleonart, J., Wei, Y., Yedra, L., Dkhil, B.Noheda, B. , Direct Epitaxial Growth of Polar $(1-x)$HfO2-$(x)$ZrO2 Ultrathin Films on Silicon. _ACS Appl. Electron. Mater._ 1, 2585–2593 (2019).
  * [30] Stokes, H. T. , Hatch, D. M. , Campbell,B. J. , and Tanner, D. E. , ISODISPLACE: A web-based tool for exploring structural distortions, _J. Appl. Cryst._ 39, 607 (2006).
  * [31] Orobengoa, D. , Capillas, C. ,Aroyo, M. I. , and Perez-Mato, J. M. , AMPLIMODES: Symmetry-mode analysis on the Bilbao Crys- tallographic Server, _J. Appl. Cryst._ 42, 820 (2009).
  * [32] Kresse, G. and Joubert, D., From ultrasoft pseudopotentials to the projector augmented-wave method. _Phys. Rev. B_ 59, 1758-1775 (1999).
  * [33] Perdew, J.P. , Ruzsinszky, A. , Csonka, G.I. , Vydrov, O.A. , Scuseria, G.E. , Constantin, L.A. , Zhou, X., Burke, K. , Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. _Phys. Rev. Lett._ 100, 136406-136410 (2008).
  * [34] Csonka, G.I. , Perdew, J.P. , Ruzsinszky, A. , Philipsen, P.H. T. , Lebègue, S. , Paier, J. , Vydrov, O.A. , Ángyán, J.G. , Assessing the performance of recent density functionals for bulk solids. _Phys. Rev. B_ 79, 155107-155121 (2009).