# Safe Evaluation For Offline Learning:
Are We Ready To Deploy?

Hager Radi
University of Alberta
Alberta, Canada
<EMAIL_ADDRESS>
&Josiah P. Hanna
University of Wisconsin – Madison
Wisconsin, USA
<EMAIL_ADDRESS>
Peter Stone
The University of Texas at Austin &
Sony AI
Texas, USA
<EMAIL_ADDRESS>
&Matthew E. Taylor
University of Alberta &
Alberta Machine Intelligence Institute
Alberta, Canada
<EMAIL_ADDRESS>

###### Abstract

The world currently offers an abundance of data in multiple domains, from
which we can learn reinforcement learning (RL) policies without further
interaction with the environment. RL agents learning offline from such data is
possible but deploying them while learning might be dangerous in domains where
safety is critical. Therefore, it is essential to find a way to estimate how a
newly-learned agent will perform if deployed in the target environment before
actually deploying it and without the risk of overestimating its true
performance. To achieve this, we introduce a framework for safe evaluation of
offline learning using approximate high-confidence off-policy evaluation
(HCOPE) to estimate the performance of offline policies during learning. In
our setting, we assume a source of data, which we split into a train-set, to
learn an offline policy, and a test-set, to estimate a lower-bound on the
offline policy using off-policy evaluation with bootstrapping. A lower-bound
estimate tells us how good a newly-learned target policy would perform before
it is deployed in the real environment, and therefore allows us to decide when
to deploy our learned policy.

##  1 Introduction

Figure 1: A framework for continual safety-evaluation of offline learning

Suppose someone else is controlling a sequential decision making task for you.
This could be a person trading stocks for you, a hand-coded controller for a
chemical plant, or even a PID for temperature regulation. Off-policy
reinforcement learning allows us to learn policies from fixed data. But when
would you want to switch from the existing controller to your learned policy?
This decision may depend on the cost for continued data collection from the
existing controller (e.g., you pay someone to manage your stock portfolio),
your risk appetite, and your confidence in the performance of the policy you
have learned. This paper takes a critical first step towards the last
question: how to learn off-policy and simultaneously evaluate that policy with
confidence when we have no access to the environment nor the policy generating
the data?

Offline reinforcement learning is a way of training off-policy algorithms
offline using already existing data, generated by humans or other controllers.
Learning with offline data is challenging because of the distribution mismatch
between data collected by the behavior agent and the offline agent (Levine et
al. 2020). What is even more challenging is evaluating offline agents in the
offline RL setting if we assume no access to the environment; in some domains,
we cannot execute our learned policy until it is good enough. This limitation
raises the question about the possibility of using off-policy policy
evaluation (OPE) methods, where we can estimate the value of a policy using
trajectories from another policy, to predict what the performance of the
offline agent is, at any point of time during learning. Further, we
investigate if we can use high-confidence OPE (HCOPE) to control the risk of
overestimating the policy’s performance.

We present a framework for safe evaluation of RL agents learning offline. Safe
refers to defining a lower bound to our policy estimates, a point where we
think the new policy is worse than it is in practice, which is important for
safety-critical applications. The paper tackles the setup where we combine
offline policy improvement to learn a policy from existing data, and off-
policy policy evaluation with bootstrapping confidence intervals to provide a
lower confidence bound estimate of such a policy, simultaneously. We believe
the setup we are tackling is under-studied in the literature. We are aiming to
learn a target policy purely out of a data buffer, and evaluate its
performance while learning so that we can tell when it is ready for
deployment. Our framework can be summarized in Figure 1: we have a source of
data, from which we utilize samples. At each step, we split the data into
training and testing. After each training step, we test the policy using
approximate 111We refer to HCOPE as approximate because bootstrapping lower
bounds may have error rates larger than $\delta$ but they provide a practical
alternative to guaranteed bounds that are too loose to use. high-confidence
off-policy evaluation (HCOPE). We continue the process of training/testing for
a few iterations until the testing shows the policy can outperform the data
with a confidence level $\delta$. We dynamically receive samples, continuously
performs RL updates, and continuously monitors a confidence interval on the
changing policy until it reaches a sufficient level. Our contributions are:

  * •

A framework combining offline RL and approximate high-confidence OPE.

  * •

Investigating the feasibility of HCOPE methods given constantly-improving
target policy (as opposed to a fixed policy previously studied in the
literature)

##  2 Motivation

Offline RL offers an opportunity for learning data-driven policies without
environment interaction. In safety-critical applications found in healthcare
or autonomous driving, there are plenty of data that we can use to learn RL
policies and hence use for decision making (Gottesman et al. 2019). Also,
there are domains where data efficiency is essential as the data collection
process is either expensive or dangerous. If we want to learn policies in such
domains, we need to find a way to tell how good the performance of a policy is
before actually deploying it to the real-world environment. However, the
execution of a new policy can be costly or dangerous if it performs worse than
the policy that is currently being used. Hence, we focus on safety when
evaluating offline agents such that the probability that the performance of
our agent below a baseline is at most $\delta$, where $\delta$ specifies how
much risk is reasonable for a certain domain. This will allow us to know when
we can trust a policy to take control in the real world. If we are given
access to trajectories generated by an unknown behavior policy, we assume an
iterative setting where at each iteration we can either request another batch
of trajectories from the source or deploy our own policy. Our objective is to
only deploy our own policy if a $1-\delta$ lower bound on its expected return
is greater than the expected return of the unknown behavior policy.

##  3 Related work

The current work intersects with the literature of two camps, which are
offline RL and batch RL. Offline RL is about improving a policy from
historical data for control, not just evaluation. A survey paper (Levine et
al. 2020) discusses how to categorize model-free offline RL methods into
policy constraints that constrain the learned policy to be close to the
behavior policy such as the work by Fujimoto, Meger, and Precup (2018), and
uncertainty-based methods that attempt to estimate the epistemic uncertainty
of Q-values to reduce distributional shift. Non-constrained methods include
the traditional Q-learning (Watkins and Dayan 1992), Double DQN Hasselt, Guez,
and Silver (2016), or Soft-Actor Critic (Haarnoja et al. 2018) but they are
not always reliable in the fully offline setup. Imitation learning, more
specifically behavioral cloning (Bain and Sammut 1999), is another way for
learning offline policies from historical data.

Batch RL refers to learning off-policy estimates from historical data, without
environment interaction. It has an important property: given data generated by
a behavior policy, it will estimate a new evaluation/target policy, guaranteed
with high confidence that its performance is not worse than the behavior
policy $\pi_{b}$. This camp includes three sub-directions. Off-policy policy
evaluation (OPE) (Voloshin et al. 2019), high confidence off-policy policy
evaluation (HCOPE) (Thomas, Theocharous, and Ghavamzadeh 2015a), and safe
policy improvement (SPI) (Thomas, Theocharous, and Ghavamzadeh 2015b). An
empirical study of OPE methods (Voloshin et al. 2019) discussed the
applicability of each method and presented method selection guidelines
depending on the environment parameters and the mismatch between
$\pi_{\theta}$ and $\pi_{b}$. This study categorized OPE methods into
importance sampling methods, directed methods, and hybrid methods that combine
aspects from the two worlds. Importance Sampling (IS) (Precup, Sutton, and
Singh 2000) is one of the widely used methods for off-policy evaluation where
rewards are re-weighted by the ratio between $\pi_{\theta}$ and $\pi_{b}$.
There are later versions such as Weighted Importance Sampling (Mahmood, van
Hasselt, and Sutton 2014), Per-Decision Importance Sampling, and Per-Decision
Weighted Importance Sampling (Precup, Sutton, and Singh 2000) that are biased,
but offer lower-variance estimates. Then, we have the direct methods which
rely on regression techniques to directly estimate the value function of the
policy. This category includes model-based methods where the transition
dynamics and reward are estimated from historical data via a model. Then, the
off-policy estimate is computed with Monte-Carlo policy evaluation (Hanna,
Stone, and Niekum 2017). Another direct method is Fitted Q-Evaluation (Le,
Voloshin, and Yue 2019), which is a model-free approach and acts as the policy
evaluation counter-part to batch Q-learning or FQI (Riedmiller 2005). The
third category is the hybrid methods that combine different features from IS
and direct methods. This category mainly involves doubly-robust methods (Jiang
and Li 2015) that uses a direct model to lower the variance of IS.

Given the previous work, none of them discussed the feasibility of evaluating
offline RL agents during learning, and how much we can trust high-confidence
OPE as an approach for testing. The setup we are studying is quite different
from the current literature because previous work assumed a behavior policy
$\pi_{b}$ and a target policy $\pi_{\theta}$, where both policies are fixed
and may be related. As an example, in a study for HCOPE methods by Thomas,
Theocharous, and Ghavamzadeh (2015a), the target policy is initialised as a
subset of the behavior policy such that they are close to each other. In
another study for safe improvement (Thomas, Theocharous, and Ghavamzadeh
2015b), the Daedulus algorithm learns a safe target policy as a continuous
improvement over the behavior policy. This is quite relevant to what we have
here but Daedulus uses data from an older version of the policy it currently
improves to perform the improvement step and safety tests. With the Doubly-
robust estimator, authors present the results with different versions of
$\pi_{\theta}$ such that $\pi_{\theta}$ is always a mixture of $\pi_{b}$ and
$\pi_{\theta}$ with different degrees (Jiang and Li 2015). To the best of our
knowledge, none of the previous work showed how OPE or HCOPE methods would
perform if the target policy is improved independently from the behavior
policy, which is the case for offline learning.

##  4 Methodology

In this work, we present a framework to close the gap between offline RL and
approximate high-confidence off-policy evaluation as a feasible evaluation
method for offline agents in safety-critical domains. In our framework, we
sample data dynamically at each iteration, perform policy updates, and
continuously monitor its performance with a confidence interval till it
reaches a good performance and is ready to be deployed. The policy we are
learning offline does not interact with the environment unless it passes the
safety test. Specifically, we have $k$ iterations where in each iteration $i$,
the current data size $n$ is denoted as $D_{i}$. $D_{i}$ is split between
training our target policy $\pi_{\theta_{i}}$ and performing the safety test
using approximate high-confidence off-policy evaluation methods. We return
when the learned policy is ready to be deployed with appropriate confidence.
To summarize an interaction between an agent and a data source: 1) an agent
requests a set of data $D_{i}$ without knowledge of the behavior policy, 2) we
improve a policy $\pi_{\theta_{i}}$ offline using $D_{train}$ with any offline
RL algorithm, 3) we perform a safety check using $D_{test}$ with high-
confidence evaluation methods, and 4) our new policy is either ready to be
deployed or we go back to step 1).
Instantiating our framework requires selecting an offline policy improvement
method and a method for computing off-policy lower-bound estimates. For these
two components, we study a variety of methods. We refer to each offline policy
improvement method with $\Phi$, and high-confidence off-policy evaluation
method with $\Psi$. Our safe evaluation framework is further explained in
Algorithm 1.

Algorithm 1 Offline Safe Evaluation Framework

Input: initial $\pi_{\theta}$, dataset $D$ of $n$ trajectories, confidence
level $\delta\in[0,1]$, number of bootstrap estimates $B$, policy improvement
method $\Phi$, HCOPE method $\Psi$
Output: $\pi_{\theta}$, $\hat{v_{\delta}}(\pi_{\theta})$: $1-\delta$ lower-
bound on $v(\pi_{\theta})$

1: Let $i=0$

2: while $\hat{v_{\delta}}(\pi_{\theta})\leq\hat{v}(\pi_{b})$ do

3: Request data-set of trajectories $D_{i}$ with size $n$

4: Split $D_{i}$ into $D_{train}$ and $D_{test}$

5: Improve policy $\pi_{\theta_{i}}=\Phi(D_{train})$

6: Evaluate policy $\hat{v_{\delta}}(\pi_{\theta})=\Psi(D_{test})$

7: end while

8: return $\pi_{\theta}$, $\hat{v_{\delta}}(\pi_{\theta})$

###  4.1 Policy improvement

We investigate different offline learning techniques such as naive off-policy
methods (Category A), imitation learning (Category B), and policy constraint
approaches specific for offline learning (Category C). We use Double DQN
(Hasselt, Guez, and Silver 2016) in Category A, Behavioral Cloning (BC) in
Category B, and Batched Constrained Q-Learning (BCQ) (Fujimoto et al. 2019) in
Category C.

###  4.2 Approximate high-confidence off-policy evaluation with bootstrapping

We use off-policy estimators to evaluate a policy that is continuously
learning. An off-policy estimator is a method for computing an estimate
${\hat{v}(\pi_{\theta})}$ on the true value of the target policy
${v(\pi_{\theta})}$ using trajectories $D$ from another policy $\pi_{b}$,
which is what Off-policy policy evaluation (OPE) methods do. For high-
confidence OPE, this extends to lower-bounding the performance of the target
policy $\pi_{\theta}$ where we combine OPE with bootstrapping confidence
intervals. Bootstrapping is considered semi-safe since it requires the
assumption that the bootstrap distribution can represent the statistic of
interest which is not true for a finite sample, which in result makes our
HCOPE approximate. However, it is safe enough as well as data efficient for
high-risk domains like medical predictions (Morrison et al. 2007). When
evaluating, we use new samples each time, to avoid the multiple comparisons
problem222The problem occurs when conducting multiple statistical tests
simultaneously with reusing data so the confidence bound does not strictly
hold anymore. This also ensures that we do not overfit our OPE estimates or
tune training parameters to reduce the estimate error. The evaluation
approximates a confidence lower bound of $\hat{v_{\delta}}(\pi_{\theta})$ on
$v(\pi_{\theta})$ such that $\hat{v_{\delta}}(\pi_{\theta})\leq
v(\pi_{\theta})$ with probability at least $1-\delta$. There are 3 main
categories of estimators as mentioned in Section 3; we use weighted importance
sampling with bootstrapping (WIS), direct model-based bootstrapping (MB), and
weighted doubly-robust with bootstrapping (WDR) to represent each of the 3
categories of off-policy estimators. Further details can be found in Appendix
A.2.

##  5 Experiments

In our experiments, we focus on discrete control in the classic MountainCar
(Sutton and Barto 2018) with a shortened horizon333Further details in Appendix
A.3. To mimic the setup of Figure 1, we simulate a data source with an online
partially-trained actor-critic that collects data of a medium quality. To
include exploratory behavior in data, an agent takes a random action $30\%$ of
the time instead of following the online policy when collecting trajectories.
We improve a policy offline using 3 different improvement methods (BC, BCQ,
and Double DQN), and evaluate each policy simultaneously using $3$ different
off-policy evaluation methods (WIS, MB, and WDR). We follow Algorithm 1, but
we limit the loop to maximum iterations of $10$ because not all evaluation
methods will be able to achieve the stopping condition. In each iteration, we
sample $300$ trajectories from the data source and split between training and
testing such that training uses $20$ trajectories only since evaluation
requires more data for tighter bounds. For each improvement method, we pass
training data and do $m$ policy updates where $m$ is tuned per algorithm such
that the offline policy converges within the total number of iterations. For
bootstrapping, we use ($\delta=0.05$) to get a $95\%$ confidence lower-bound
using $B=2000$ bootstrap estimates as recommended by practitioners (Efron
1979). For WDR only, we use a value of $B=224$ to avoid the heavy computation;
this can still get us a good approximation as suggested by MAGIC (Thomas and
Brunskill 2016).

Figure 2: Results of safe evaluation using different offline improvement
algorithms on MountainCar. Double DQN and BCQ outperforms the data
distribution while BC reaches $\pi_{b}$ at most. WIS, MB, and WDR refer to the
evaluation methods: weighted importance sampling, model-based, and weighted
doubly-robust respectively.

Results: Figure 2 shows how different offline policy improvement methods
perform given medium-quality data along with HCOPE estimators, indicating
which method can tell when $v(\pi_{\theta})>v(\pi_{b})$. Behavior cloning as
an offline policy improvement method can only perform as well as the data by
$\pi_{b}$ while Double DQN and BCQ were able to outperform $\pi_{b}$. The true
value of a target policy is calculated as the average return when running the
policy in the actual environment (not possible in practice) for $1000$
episodes. All reported results are average of 40 runs while the shaded area
shows the standard error. The value of the behavioral policy
$\hat{v}(\pi_{b})$ is the sum of undiscounted rewards of the data set. Since
$\pi_{b}$ is not known, we first estimate $\hat{\pi_{b}}$ given the test data
(Hanna, Niekum, and Stone 2021).

For the safe evaluation of each offline agent, weighted importance sampling
(WIS) with bootstrapping failed to detect that the offline policy outperforms
$\pi_{b}$ for all improvement methods. The model-based estimator (MB) and
weighted doubly-robust with bootstrapping (WDR) have much less error with the
true value of a policy and were able to inform when an offline agent
outperforms $\pi_{b}$ and hence ready to be deployed. For instance, with 95%
confidence, in case of offline improvement with Double DQN, we were able to
tell that our new target policy is better than the behavior policy at
iteration 4 using two estimators (MB and WDR). Our ability to detect an
improved policy is a function of how much better that policy is and how good
our OPE methods are. It is better to rely on estimators that do not take
$\pi_{b}$ into account (e.g. direct methods or hybrid methods) so that
estimates are less affected by the offline improvement method and its
divergence from $\pi_{b}$ (as the case for WIS). Further analysis can be found
in Appendix A.4. Accordingly, high-confidence off-policy estimators (e.g.
direct methods or hybrid methods) are a safe evaluation method for offline
learning with enough confidence that controls the risk of overestimating the
true performance of an offline policy.

##  6 Conclusion

In this paper, we propose a framework for safe evaluation of offline RL
methods. While dynamically receiving data, we train offline RL agents and run
safety tests to estimate a lower-bound on the value of the target policy and
control the risk of overestimating its true value. This is essential for
safety-critical applications to be able to tell when it is safe to deploy a
new policy. We believe safe evaluation is an important step for offline RL;
offline RL has great potential for control in the actual environment, if they
are good enough, where our proposed framework is applicable. In future work,
we will test this framework on more complex environments, both discrete and
continuous. We will include learning offline from multiple sources of data
with different qualities and explore how the quality of the data affects safe
evaluation. Moreover, we can use safety testing to develop algorithms which
can estimate in advance how much data is needed to learn a good-enough policy
that can take over the real environment with maximum data efficiency.

## Acknowledgments and Disclosure of Funding

This work has taken place in the Intelligent Robot Learning (IRL) Lab at the
University of Alberta, which is supported in part by research grants from the
Alberta Machine Intelligence Institute (Amii); a Canada CIFAR AI Chair, Amii;
NSERC; and Compute Canada. A portion of this work has taken place in the
Learning Agents Research Group (LARG) at UT Austin. LARG research is supported
in part by NSF (CPS-1739964, IIS-1724157, FAIN-2019844), ONR (N00014-18-2243),
ARO(W911NF-19-2-0333), DARPA, Lockheed Martin, GM, Bosch, and UT Austin’s Good
Systems grand challenge. Peter Stone serves as the Executive Director of Sony
AI America and receives financial compensation for this work. The terms of
this arrangement have been reviewed andapproved by the University of Texas at
Austin in accordance with its policy on objectivity in research.

## References

  * Bain and Sammut (1999) Bain, M.; and Sammut, C. 1999.  A Framework for Behavioural Cloning.  In _Machine Intelligence 15, Intelligent Agents [St. Catherine’s College, Oxford, July 1995]_ , 103–129. GBR: Oxford University.  ISBN 0198538677.
  * Efron (1979) Efron, B. 1979.  Bootstrap Methods: Another Look at the Jackknife.  _Ann. Statist."_ , 7.
  * Fujimoto et al. (2019) Fujimoto, S.; Conti, E.; Ghavamzadeh, M.; and Pineau, J. 2019.  Benchmarking Batch Deep Reinforcement Learning Algorithms.  _CoRR_ , abs/1910.01708.
  * Fujimoto, Meger, and Precup (2018) Fujimoto, S.; Meger, D.; and Precup, D. 2018.  Off-Policy Deep Reinforcement Learning without Exploration.  _CoRR_ , abs/1812.02900.
  * Gottesman et al. (2019) Gottesman, O.; Johansson, F.; Komorowski, M.; Faisal, A.; Sontag, D.; Doshi-Velez, F.; and Celi, L. 2019.  Guidelines for reinforcement learning in healthcare.  _Nature Medicine_ , 25: 16–18.
  * Haarnoja et al. (2018) Haarnoja, T.; Zhou, A.; Abbeel, P.; and Levine, S. 2018.  Soft Actor-Critic: Off-Policy Maximum Entropy Deep Reinforcement Learning with a Stochastic Actor.  In Dy, J.; and Krause, A., eds., _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , 1861–1870. PMLR.
  * Hanna, Niekum, and Stone (2021) Hanna, J. P.; Niekum, S.; and Stone, P. 2021.  Importance Sampling in Reinforcement Learning with an Estimated Behavior Policy.  _Machine Learning (MLJ)_.
  * Hanna, Stone, and Niekum (2017) Hanna, J. P.; Stone, P.; and Niekum, S. 2017.  Bootstrapping with Models: Confidence Intervals for Off-Policy Evaluation.  In _Proceedings of the 16th Conference on Autonomous Agents and MultiAgent Systems_ , AAMAS ’17, 538–546. Richland, SC: International Foundation for Autonomous Agents and Multiagent Systems.
  * Hasselt, Guez, and Silver (2016) Hasselt, H. v.; Guez, A.; and Silver, D. 2016.  Deep Reinforcement Learning with Double Q-Learning.  In _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence_ , AAAI’16, 2094–2100. AAAI Press.
  * Jiang and Li (2015) Jiang, N.; and Li, L. 2015.  Doubly Robust Off-policy Evaluation for Reinforcement Learning.  _CoRR_ , abs/1511.03722.
  * Le, Voloshin, and Yue (2019) Le, H. M.; Voloshin, C.; and Yue, Y. 2019.  Batch Policy Learning under Constraints.  _CoRR_ , abs/1903.08738.
  * Levine et al. (2020) Levine, S.; Kumar, A.; Tucker, G.; and Fu, J. 2020.  Offline Reinforcement Learning: Tutorial, Review, and Perspectives on Open Problems.  _CoRR_ , abs/2005.01643.
  * Mahmood, van Hasselt, and Sutton (2014) Mahmood, A. R.; van Hasselt, H. P.; and Sutton, R. S. 2014.  Weighted importance sampling for off-policy learning with linear function approximation.  In Ghahramani, Z.; Welling, M.; Cortes, C.; Lawrence, N.; and Weinberger, K. Q., eds., _Advances in Neural Information Processing Systems_ , volume 27. Curran Associates, Inc.
  * Morrison et al. (2007) Morrison, A. C.; Bare, L. A.; Chambless, L. E.; Ellis, S. G.; Malloy, M.; Kane, J. P.; Pankow, J. S.; Devlin, J. J.; Willerson, J. T.; and Boerwinkle, E. 2007\.  Prediction of Coronary Heart Disease Risk using a Genetic Risk Score: The Atherosclerosis Risk in Communities Study.  _American Journal of Epidemiology_ , 166(1): 28–35.
  * Precup, Sutton, and Singh (2000) Precup, D.; Sutton, R. S.; and Singh, S. P. 2000.  Eligibility Traces for Off-Policy Policy Evaluation.  In _Proceedings of the Seventeenth International Conference on Machine Learning_ , ICML ’00, 759–766. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.  ISBN 1558607072.
  * Riedmiller (2005) Riedmiller, M. 2005.  Neural Fitted Q Iteration – First Experiences with a Data Efficient Neural Reinforcement Learning Method.  volume 3720, 317–328.  ISBN 978-3-540-29243-2.
  * Sutton and Barto (2018) Sutton, R. S.; and Barto, A. G. 2018.  _Reinforcement Learning: An Introduction_.  The MIT Press, second edition.
  * Thomas (2015) Thomas, P. S. 2015.  _Safe reinforcement learning_.  Ph.D. thesis, University of Massachusetts Libraries.
  * Thomas and Brunskill (2016) Thomas, P. S.; and Brunskill, E. 2016.  Data-Efficient Off-Policy Policy Evaluation for Reinforcement Learning.  _CoRR_ , abs/1604.00923.
  * Thomas, Theocharous, and Ghavamzadeh (2015a) Thomas, P. S.; Theocharous, G.; and Ghavamzadeh, M. 2015a.  High-Confidence Off-Policy Evaluation.  In _AAAI_.
  * Thomas, Theocharous, and Ghavamzadeh (2015b) Thomas, P. S.; Theocharous, G.; and Ghavamzadeh, M. 2015b.  High Confidence Policy Improvement.  In _Proceedings of the 32nd International Conference on International Conference on Machine Learning - Volume 37_ , ICML’15, 2380–2388. JMLR.org.
  * Voloshin et al. (2019) Voloshin, C.; Le, H. M.; Jiang, N.; and Yue, Y. 2019.  Empirical Study of Off-Policy Policy Evaluation for Reinforcement Learning.  _CoRR_ , abs/1911.06854.
  * Watkins and Dayan (1992) Watkins, C. J. C. H.; and Dayan, P. 1992.  Q-learning.  In _Machine Learning_ , 279–292.

##  Appendix A Appendix

###  A.1 Bootstrapping

Consider a sample $D$ of $n$ random variables $H_{j}$ for $j=1,2,...,n$ where
we can sample $H_{j}$ from some $i.i.d.$ distribution of data. From the sample
of data $D$, we can compute a sample estimate $\hat{\theta}$ of a parameter
$\theta$ such that $\hat{\theta}=f(D)$ where $f$ is the function to compute
$\theta$. Given a dataset $D$, we create $B$ resamples with replacement, where
$B$ is the number of bootstrap resamples, and compute $\theta$,
$\hat{\theta}$, on each of these resamples. Bootstrapping (Efron 1979) allows
us to estimate the distribution of $\hat{\theta}$ with confidence intervals.
The estimates computed with different resamples will be used to determine the
$1-\delta$ confidence interval. In our setup, the parameter of interest
$\theta$ is the expected return of a policy $v(\pi)$.

Hence, with a confidence level $\delta\in[0,1]$ and $B$ resamples of the
dataset of trajectories $D$, we use bootstrapping methods to approximate a
confidence lower bound of $v_{\delta}(\pi)$ on $v(\pi)$ such that
$v_{\delta}(\pi)\leq v(\pi)$ with probability at least $1-\delta$.

As the size of data $n\to\infty$, bootstrapping has strong guarantees but it
lacks guarantees for finite samples as it is not an exact method. To use
bootstrapping, we have to assume that the bootstrap distribution is
representative of the distribution of the statistic of interest, which is not
the case for finite samples (Hanna, Stone, and Niekum 2017). As a result,
bootstrapping is considered semi-safe but it is still safe enough for high
risk medical predictions in practice with a known record of producing accurate
confidence intervals (Morrison et al. 2007).

###  A.2 High-confidence off-policy evaluation techniques

#### Importance sampling with bootstrapping

Importance Sampling (IS) (Precup, Sutton, and Singh 2000) is a way for
handling mismatch between distributions and hence presented as a consistent
and unbiased off-policy estimator. For a trajectory $H\sim\pi_{b}$ of length
$L$, as $H=S_{1},A_{1},..,S_{L},A_{L}$, we can define the importance sampling
up to time $t$ for policy $\pi_{\theta}$ as follows:

$IS(\pi_{\theta},\pi_{b},D)=\sum_{i=1}^{m}\rho_{L}^{H}R^{i},\rho_{t}^{H}:=\prod_{j=0}^{t}\frac{\pi_{\theta}(A_{j}|S_{j})}{\pi_{b}(A_{j}|S_{j})}$

In our setup given no access to $\pi_{b}$, we estimate the behavior policy
$\pi_{b}$ from data but it results in a biased estimator (Hanna, Niekum, and
Stone 2021). To produce a high-confidence IS estimate, we first compute the
importance weighted returns then use bootstrapping to get the lower-bound
estimate. A popular method is Bias Corrected and accelerated (BCa)
bootstrapping (Efron 1979). A pseudo code of BCa is well-described as
Algorithm 3 in (Thomas, Theocharous, and Ghavamzadeh 2015b). We use weighted
importance sampling (Mahmood, van Hasselt, and Sutton 2014) for discrete
control along with bootstrapping.

#### Direct model-based with bootstrapping

Model-based estimation is another off-policy estimator that lies under the
direct methods. The model-based off-policy estimator MB computes
$v(\pi_{\theta})$ by building a model using all the available trajectories $D$
to build a model $\hat{M}=(S,A,\hat{P},r,\gamma,\hat{d_{0}})$ where $\hat{P}$
and $\hat{d_{0}}$ are estimated with trajectories sampled from $\pi_{b}$.
Then, MB will compute $\hat{v}(\pi_{\theta})$ as the average return of
trajectories simulated in the estimated model $\hat{M}$ while following
$\pi_{\theta}$. Despite having lower variance than IS methods, MB is an
inconsistent estimator such that as $n\rightarrow\infty$, the model estimates
may converge to a value different from $V(\pi_{\theta})$. This is because it
is dependant on the modeling assumptions we make, whether we assume a linear
or a non-linear model. To get a lower confidence bound on the MB estimate, we
use bootstrap confidence intervals; the exact algorithm used is detailed in
Hanna, Stone, and Niekum (2017). This bootstrapping method is quite similar to
what we mentioned in the previous paragraph but much simpler in implementation
and adaptable to different OPE methods. We rely on the model-based estimator
with bootstrapping in our study as one of OPE direct-methods.

#### Weighted doubly-robust with bootstrapping

Weighted Doubly-Robust (WDR) (Thomas and Brunskill 2016) is a hybrid method
for off-policy estimation, presented as an extension to the doubly-robust (DR)
method (Jiang and Li 2015). DR is an unbiased estimator of $v(\pi_{\theta})$
that uses an approximate model of the MDP to reduce the variance of importance
sampling (Jiang and Li 2015). Although biased, WDR is based on per-decision
weighted importance sampling (PDWIS) and serves as an improvement over DR
method as it significantly balances the bias-variance trade-off. The
approximate model value functions act as a control variate for PDWIS (Thomas
and Brunskill 2016).

$PDWIS(\pi_{\theta},D,\pi_{b})=\sum_{i=1}^{m}\sum_{t=0}^{L-1}\frac{\rho_{t}^{i}}{\sum_{j=1}^{m}\rho_{t}^{j}}\gamma^{t}R_{t}^{i}$
(1)

$\begin{split}WDR(\pi_{\theta},D,\pi_{b})&=PDWIS(\pi_{\theta},D,\pi_{b})-\\\
&\sum_{i=1}^{n}\sum_{t=0}^{L-1}\gamma^{t}(w_{t}^{i}\hat{q}_{\pi_{\theta}}(S_{t}^{i},A_{t}^{i})-w_{t-1}^{i}\hat{v}_{pi_{\theta}}(S_{t}^{i}))\end{split}$
(2)

In previous methods, we use bootstrapping with WDR to provide lower confidence
bound estimates over $v(\pi_{\theta})$. WDR with bootstrapping is guaranteed
to converge to the correct estimate as $n$ increases given the statistical
consistency of PDWIS (Hanna, Stone, and Niekum 2017). For the approximate
model, a single model is estimated with the available trajectories $D$, and
used to compute the value functions of WDR for each bootstrap data. We choose
the weighted doubly robust estimator to represent OPE hybrid methods in our
study.

###  A.3 Experimental Details

Environment: We use the classic MountainCar (Sutton and Barto 2018) with a
continuous state (velocity and position), and 3 possible discrete actions. At
each time-step, the reward is $-1$ except for in a terminal state when it is
$0$. However, we used the modified version of MountainCar as described by
Thomas (2015). This means we shorten the horizon of the problem by holding an
action $a_{t}$ constant for $4$ updates of the environment state. We also
change the start state such that an episode starts with a random position and
random velocity (Jiang and Li 2015; Thomas 2015).

###  A.4 Analysis of evaluation methods

To better understand how safety testing is affected by the policy improvement
methods, there are multiple measures to tell the similarity between two
probability distributions (as a policy is a distribution over states and
actions). One of the practical measures is the total variation (TV) distance.
TV distance is a way to measure the difference between action probabilities
taken under two policies given the data set in hand. TV would be the sum of
differences in probabilities between the behavior policy $\pi_{b}$ and the
target policy $\pi_{\theta}$ for each state-action pair in the test set of
data. Since $\pi_{b}$ is not known, we estimate $\hat{\pi_{b}}$ given the test
data (Hanna, Niekum, and Stone 2021).

Figure 3: Total Variation Distance between $\hat{\pi}_{b}$ & $\pi_{\theta}$

We analyze how the improvement method is affecting TV distance between
$\hat{\pi_{b}}$ and $\pi_{\theta}$ and hence affecting the high-confidence
off-policy evaluation estimates. TV distance is correlated to KL-Divergence
and hence shows how two policies are different from each other. Figure 3 shows
the total variation distance between the offline policy $\pi_{\theta}$ and the
estimated behavior policy $\hat{\pi}_{b}$ across the different offline
learning methods. It show that behavioral cloning (imitation learning) can
achieve a much lower distance than other improvement methods whether they are
constrained or non-constrained (Double DQN and batch-constrained Q-learning);
this is because behavioral cloning forces its target policy to be close to the
behavior policy while other methods do not. This result also explains why
weighted importance sampling achieves the lowest error between the estimate
and the true value in the case of behavioral cloning; the error grows for
other methods that do not constrain the policy to be close to the data
distribution.


# Interplay between ferroelectricity and metallicity in hexagonal YMnO3

Tara N. Tošić<EMAIL_ADDRESS>Materials Theory, ETH Zürich, Wolfgang-
Pauli-Strasse 27, 8093 Zurich, Switzerland  Yuting Chen  Materials Theory, ETH
Zürich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland  Nicola A. Spaldin
<EMAIL_ADDRESS>Materials Theory, ETH Zürich, Wolfgang-Pauli-Strasse 27, 8093
Zurich, Switzerland

###### Abstract

We use first-principles density functional theory to investigate how the polar
distortion is affected by doping in multiferroic hexagonal yttrium manganite,
h-YMnO3. While the introduction of charge carriers tends to suppress the polar
distortion in conventional ferroelectrics, the behavior in improper geometric
ferroelectrics, of which h-YMnO3 is the prototype, has not been studied to
date. Using both background charge doping and atomic substitution, we find an
approximately linear dependence of the polar distortion on doping
concentration, with hole doping suppressing and electron doping enhancing it.
We show that this behavior is a direct consequence of the improper geometric
nature of the ferroelectricity. In addition to its doping effect, atomic
substitution can further suppress or enhance the polar distortion through
changes in the local chemistry and geometry.

††preprint: APS/123-QED

##  I Introduction

The combination of ferroelectricity and metallicity is seemingly incompatible,
because metallic charge carriers tend to screen the long-range interactions
responsible for ferroelectricity and inhibit polarization switching by an
electric field. In spite of this apparent contra-indication, there have been a
number of recent studies on introducing charge carriers into ferroelectric
(FE) materials, motivated in part by the possible relevance of polar metals
for phenomena such as superconductivity and topology [1]. In a device context,
materials with polarization robust to charge doping, could be used for example
as electrodes to stabilize thin film polarization [2]. The most widely studied
is the prototypical FE BaTiO3 [3, 4, 5, 6, 7, 8], where electron doping
directly modifies the electronic structure of the formally $d^{0}$ Ti4+ site,
hindering its off-centering through the second-order Jahn-Teller effect and
ultimately suppressing the primary FE order parameter [8].

Geometric FEs are promising candidates for persistence of the polar distortion
on doping [9, 10], since their FE distortion is driven by a structural
instability such as a polyhedral rotation that is not caused by the second-
order Jahn-Teller effect. As a result, their polarization is expected to be
less sensitive to the detailed electronic occupation of the atomic orbitals.
Note that, since ferroelectricity is not well defined for a metal, we will
refer to the FE distortion as a polar distortion in doped samples. There have
been a number of studies on the effect of doping in proper geometric FEs, in
which the polar geometric mode is the primary order parameter connecting the
paraelectric (PE) and FE structures. For example, in the AnBnO3n+2 Carpy-Galy
family, a structural distortion consisting largely of rotations of the oxygen
octahedra gives rise to the polarisation. First-principles computations find
the polarisation in La2Ti2O7 and Sr2Nb2O7 to be robust to doping [11],
consistent with the observed coexistence of non-centrosymmetricity and quasi
1D conductivity [12]. On the other hand, ferroelectricity in improper
geometric FEs emerges as a secondary effect, through coupling to a non-polar
structural order parameter which drives the FE phase transition. To our
knowledge, the only studies on the coexistence of itinerant charges and
ferroelectricity in improper FEs have been done on the Ruddlesden-Popper phase
Ca3Ti2O7. In this material, the polar mode couples to two tilt-modes through a
trilinear coupling term [13]. Here, optical second-harmonic imaging and
atomic-resolution scanning transmission electron microscopy have measured
formation of polar domains coexisting with metallicity [14, 15], which was
further confirmed by density functional theory [16, 17].

Figure 1: a) Structure of h-YMnO3 in the PE P63/mmc structure b) The two modes
involved in the FEphase transition. Left: the polar $\Gamma_{2}^{-}$ mode
which creates a polar distortion through relative displacements (shown with
arrows) of the ions. Middle and right: the primary FE zone-boundary $K_{3}$
mode which involves tilting of the MnO5 bipyramids towards their corner
sharing O sites and a buckling of the A- ion planes. The angle $\phi$
characterizes the phase of the bipyramidal tilts. c) Structure of h-YMnO3 in
the FE P63cm structure.

Motivated by the promising studies on doped hybrid geometric FEs, this work
investigates the effects of doping on the prototypical geometric improper FE,
h-YMnO3. h-YMnO3 belongs to the family of hexagonal manganites h-RMnO3, where
R denotes In, Y, Sc or small radius rare-earth elements (Dy-Lu)) [18, 19, 20,
21, 22]. The structure is composed of planes of corner-shared MnO5 trigonal
bipyramids (Mn3+ and O2-), alternating with triangular layers of Y3+ ions. At
high temperature, the structure is centrosymmetric with the P63/mmc space
group, (see Fig. 1a)) [23, 24]. Below $T_{c}\simeq 1270$ K, the MnO5 trigonal
bipyramids tilt towards trimerization centers, inducing a buckling of each Y
layer [25, 19, 20] as shown in the right panel of Fig. 1b). This structural
distortion is associated with the zone-boundary $K_{3}$ mode of the high-
symmetry structure, which triples the unit cell and lowers the symmetry to the
polar P63cm space group (shown in Fig. 1c)), but with no polarization. We
write its associated order parameter as
$\mathbf{Q_{\text{K${}_{3}$}}}=Q_{\text{K${}_{3}$}}(\cos{(\phi}),\sin{(\phi)})$,
where $Q_{\text{K${}_{3}$}}$ describes the mode’s amplitude and
$\phi=n\frac{\pi}{3}$ its phase ($n$ is an integer), as illustrated in Fig.
1b). A secondary zone center $\Gamma_{2}^{-}$ structural distortion in which
anions and cations shift in opposite directions along $c$ couples to K3,
inducing a net polarization described by the scalar order parameter
$Q_{\Gamma_{2}^{-}}$. Note that this two-mode coupling distinguishes this
material from the previously studied hybrid improper Ca3Ti2O7, where a three-
mode coupling term stabilizes the polar displacements [14]. The two modes
couple non-linearly and the Landau free energy in terms of
$\mathbf{Q_{\text{K${}_{3}$}}}$ and $Q_{\Gamma_{2}^{-}}$ up to fourth order is
[25, 18]:

$\displaystyle E(\mathbf{Q_{\text{K${}_{3}$}}},Q_{\Gamma_{2}^{-}})$
$\displaystyle=\alpha Q_{\text{K${}_{3}$}}^{2}+\beta
Q_{\text{K${}_{3}$}}^{4}+\alpha_{P}Q_{\Gamma_{2}^{-}}^{2}+\beta_{P}Q_{\Gamma_{2}^{-}}^{4}$
$\displaystyle+\gamma
Q_{\text{K${}_{3}$}}^{3}Q_{\Gamma_{2}^{-}}\cos{(3\phi)}+\gamma^{\prime}Q_{\text{K${}_{3}$}}^{2}Q_{\Gamma_{2}^{-}}^{2}\quad,$
(1)

where the direction of the K3 displacement, and hence the sign of the
polarization, are set by the value of $\phi$ (see Fig.1b)). Since the polar
distortion arises as a secondary effect through coupling to the K3 primary
order [18, 25, 19], we anticipate that, compared to conventional FEs, it will
again, as in the hybrid improper compound, be robust to the presence of dopant
charges, leading to a possible large region of coexistence between the polar
distortion and metallicity.

##  II Computational details

###  II.1 DFT details

Our calculations are performed within the density functional theory (DFT)
framework using the planewave pseudopotential VASP code [26, 27, 28, 29]. We
use the projector-augmented wave pseudopotentials [30, 31] Y_sv, Mn_sv and O
from the VASP library with (4s)2(4p)6(5s)2(4d)1, (3s)2(4s)1(3p)6(3d)6 and
(2s)2(2p)4 valence configurations, respectively. For the substitution atoms,
we use the Mg_sv, Ca_sv, Hf, Ti_sv and V_sv projector-augmented wave
pseudopotentials, with the following valence shell electrons: (2s)2(3s)2(2p)6,
(3s)2(4s)2(3p)6, (5d)3 (6s)1, (3s)2(4s)1(3p)6(3d)3 and (3s)2(4s)1(3p)6(3d)4.
We treat the exchange correlation functional using the local density
approximation (LDA) + Hubbard U approach [32], with U=5.2 eV on Mn $d$
orbitals and 6.5 eV on V $d$ orbitals, calculated using the linear response
Ansatz [33] (see Supplementary Information). We use a $\Gamma$-centered
k-point grid of $6\times 6\times 3$ for the full structural relaxation and
further relax the internal coordinates with a $12\times 12\times 6$ k-grid
upon introduction of the charges. Densities of states (DOS) are calculated
using an $18\times 18\times 9$ grid. We perform non-constrained collinear
magnetic calculations, with magnetic moments on the Mn ions (initialized to
4µB) parallel to each other within each layer and anti-ferromagnetically
aligned between the consecutive layers, to preserve the structural symmetry.

###  II.2 Charge doping

In the first part of our work, charge doping is simulated by changing the
total number of valence electrons in the system and adding a compensating
background charge to maintain charge neutrality by modifying the NELECT flag
in VASP. Since the relaxation of the lattice parameters is not well-defined
during this process [34], we fix the lattice constants either to their fully
relaxed undoped low symmetry values (FE lattice parameters) or fully relaxed
high symmetry values (PE lattice parameters).

###  II.3 Substitution doping

In the second part of our work, we dope by impurity atom substitution, and
single out three factors that influence $P$: the size of the impurity atom and
the subsequent structural changes, the role of the charge carrier, and the
variations in lattice constants. We quantify their separate contributions as
follows: starting from a relaxed h-YMnO3 cell, we first introduce the impurity
atom but subtract the added charges. For example, in the case of the
quadrivalent Hf dopant substituted on the A1 or A2 site, we subtract one
electron and add a background charge to maintain charge neutrality. In a
second step, we relax the structure containing the dopant atom and do not
modify the charge. The lattice constants are fixed in the first two steps.
Finally, we remove the constraint on the lattice parameters and let the whole
cell relax. By comparing the polar distortion of the undoped cell with that
calculated after each step, we can quantify the different contributions to $P$
stemming from the size of the impurity atom, the role of the charge carrier
and the change in lattice constants, respectively. We name the aforementioned
three separate contributions the impurity atom effect, the charge carrier
effect and the lattice constant effect.

###  II.4 Calculation of polar distortion and COHP

Since we are interested in doped systems where the Berry phase polarization is
undefined, we calculate the size of the polar distortion, $P$, by summing the
displacements ($d_{i}$) of the ions from their reference positions in the
high-symmetry structure, multiplied by their Born Effective Charges (BECs),
$Z^{*}_{i}$:

$P=\frac{e}{\Omega}\sum_{i}Z_{i}^{*}d_{i}\quad,$ (2)

where $e$ is the elementary electronic charge and $\Omega$ the unit cell
volume.

We use the following BEC values calculated for the FE structure in reference
[35]: $Z_{\text{Y}}^{*}=+3.6$, $Z_{\text{Mn}}^{*}=+3.3$,
$Z_{\text{O\textsubscript{ax}}^{*}}=-2.3$ and
$Z_{\text{O\textsubscript{eq}}^{*}}=-2.2$ (where Oax and Oeq refer to the
axial and equatorial oxygens, respectively, see Fig. 1). Note that due to the
geometric nature of the ferroelectricity, the BEC values are close to the
formal charges [36]. Using these values, we obtain $P$=8.06 µC/cm2 for our
fully relaxed undoped structure.

Bond strength is characterized using crystal orbital Hamiltonian populations
(COHPs) and integrated crystal orbital Hamiltonian populations (ICOHPs)
calculated using the LOBSTER package [37, 38, 39], with an integration energy
range between $\rm E_{f}-10$ eV and $\rm E_{f}$, capturing the occupied Mn
$3d$ and O $2p$ states. Our -ICOHP plots indicate the change in bonding
strength as a function of doping. The input WAVECAR file for the COHP
calculations is obtained in VASP without symmetrization of the charge density
(ISYM flag is set to -1) and a k-point grid of $12\times 12\times 6$.

##  III Results and discussion

We begin by presenting the main results of our calculations on the effects of
background charge doping on the polar distortion.

###  III.1 Charge doping

We begin by calculating the structure as a function of background-charge
electron and hole doping concentration at fixed PE and FE lattice parameters
and extracting the size of the polar distortion following Eq. 2. We use our
calculated undoped lattice constants which are slightly underestimated for
both the FE P63cm ($a=6.05$ Å, $c=11.32$ Å) and PE P63/mmc ($a=3.53$ Å,
$c=11.13$ Å) structures compared to measured values ($a=6.14$ Å, $c=11.41$ Å
and $a=3.61$ Å, $c=11.39$ Å, respectively) [40, 23] as well as to the DFT-
calculated values in Ref. [25], which used an LSDA pseudopotential and U=8 eV.
We show the size of our calculated polar distortion as a function of electron
and hole doping concentration in Fig. 2a). In contrast to the conventional FE
BaTiO3, where both charge carriers suppress $P$, in this case $P$ increases
monotonically with electron doping: the added electrons do not hinder the
ferroelectricity and even promote it. Hole doping, on the other hand, lowers
$P$ but does not fully suppress it. We also note that, while the polar
distortion remains slightly larger for the FE lattice parameters than for the
PE lattice parameters at all doping values, both curves follow the same trend.

Figure 2: a) Polar distortion as a function of hole and electron doping
concentration at the lattice parameters of the fully relaxed FE and PE
structures. b) Amplitude ($Q_{\text{K${}_{3}$}}$) of the primary FE distortion
mode K3 as a function of doping concentration at FE and PE lattice parameters.
c) $Q_{\Gamma_{2}^{-}}$ plotted as a function of $Q_{\text{K${}_{3}$}}$
evaluated for the range of hole and electron doping, for FE and PE lattice
parameters. The gray line shows $P\simeq 0.4Q_{\text{K${}_{3}$}}$, which we
use as the reference trend for the rest of our analysis.

Figure 3: Upper panel: orbital-resolved density of states for (left) 3 dopant
holes per 30 atom unit cell, (middle) no doping and (right) 3 dopant electrons
per 30 atom unit cell (left, middle and right columns, respectively) in
h-YMnO3 for the Mn, Oeq and Oax sites (top, middle and bottom rows,
respectively). Vertical dotted line indicates the top of the valence band in
the undoped compound, with the shaded gray areas indicating the change in
occupied states on doping. Lower panel: change in charge density in the
conduction band when doped with electrons (right) and in the valence band and
when doped with holes (left). The isosurface is plotted at 0.005 e/bohr3 in
both cases.

To uncover the cause of this change in $P$, we first analyze the structural
changes which occur upon doping by performing a symmetry mode analysis using
the AMPLIMODES program [41, 42]. Fig. 2b) shows the dependence of the primary
FE phase transition order parameter, $Q_{\text{K${}_{3}$}}$, on charge doping,
for FE and PE lattice parameters. We see that the evolution of
$Q_{\text{K${}_{3}$}}$ as a function of charge is similar to that of
$Q_{\Gamma_{2}^{-}}$ in Fig. 2a) for both sets of lattice parameters. We
confirm this in Fig. 2c) where we plot $Q_{\Gamma_{2}^{-}}$ versus
$Q_{\text{K${}_{3}$}}$ and find a linear dependence, consistent with Landau
theory at large $Q_{\text{K${}_{3}$}}$ [25]. In the following, we investigate
the electronic causes for these changes in $Q_{\text{K${}_{3}$}}$ and
$Q_{\Gamma_{2}^{-}}$ upon introduction of charges.

Figure 4: a) -COHP between Oeq $p$ orbitals and Mn $d_{xy}/d_{x^{2}-y^{2}}$,
Mn $d_{xz}/d_{yz}$ and Mn $d_{z^{2}}$ orbitals, for undoped h-YMnO3. Positive
and negative values correspond to bonding and anti-bonding states,
respectively. b) -ICOHP between pairs of Oeq $p$ and Mn $d$ orbitals as a
function of doping.

To understand the increase and decrease of $Q_{\text{K${}_{3}$}}$,
$Q_{\Gamma_{2}^{-}}$ and hence $P$, with electron and hole doping,
respectively, we next calculate the electronic DOS and the charge density as a
function of doping, and show our results in Fig. 3. In the upper panel, we
plot our calculated atom- and orbital-resolved DOSs for undoped h-YMnO3, as
well as for representative electron and hole concentrations of three charge
carriers per 30 atom unit cell. The black-dotted vertical line indicates the
top of the valence band in the undoped compound and the gray shaded rectangle
indicates the shift in $E_{F}$ upon doping, with the edge without the black-
dotted vertical line indicating the $E_{F}$ value in the doped compounds. When
no charges are introduced, the valence band states closest to the top of the
valence band are the hybridized Mn $d_{xy},d_{x^{2}-y^{2}}$ and Oeq
$p_{x},p_{y}$ orbitals, while the Mn $d_{xz}$, dyz states hybridize with the
Oeq $p_{z}$ and Oax $p_{x},p_{y}$ states at 1 to 1.5 eV below the valence bond
edge. Our calculated band is gap is equal to $\simeq 1.3$ meV and the
conduction band consists of Mn $d_{z^{2}}$, Oeq $p_{x},p_{y}$ and the Oax
$p_{z}$ orbitals. In the case of electron doping (right panel of Fig. 3),
$E_{f}$ moves into the conduction band, with electrons first occupying the Mn
$d_{z^{2}}$ states and then, at 1 electron per unit cell, the Oeq
$p_{x}+p_{y}$ orbital; this is clearly seen in the charge density, shown below
the DOS in Fig. 3, and in the occupied DOS Mn $d_{z^{2}}$ and Oeq
$p_{x}+p_{y}$ tails in Fig. S3. As electron doping increases, the occupation
of these states increases further and the dopant electrons also start
occupying the Oax $p_{z}$ orbitals (see the 3 and 5 electron/unit cell DOS and
charge density plots in Fig. 3 and Fig. S3). In the case of hole doping, $\rm
E_{f}$ moves into the valence band (Fig. 3 left panel) and the depletion of
electrons is initially most pronounced from the Mn $d_{x^{2}-y^{2}}$, $d_{xy}$
and Oeq $p_{x}+p_{y}$ orbitals. As hole doping increases, Mn $d_{z}^{2}$ and
Oax $p$ states also start depleting. We see in the charge density plots below
the DOSs that the charge state of the Y remains effectively unchanged.

Two hybridizations are allowed by symmetry between the planar Oeq $p$ orbitals
and the occupied Mn $d$ orbitals in the ideal MnO5 bipyramidal environment:
$p-d_{xz}/d_{yz}$ and $p-d_{xy}/d_{x^{2}-y^{2}}$ (note that, as a result of
the MnO5 tilt in the P63cm structure, small $p_{z}-d_{xy}/d_{x^{2}-y^{2}}$ and
$p_{x,y}-d_{z^{2}}$ hybridizations are also allowed, consistent with optical
spectroscopy measurements [43]). To determine how the hybridization evolves
with doping, we evaluate the Crystal Orbital Hamilton Population (COHP) for
the undoped compound (shown in Fig. 4a)) as well as the Integrated Crystal
Orbital Hamilton Population (ICOHP) up to EF as a function of doping
concentration (shown in Fig. 4b)).

In the undoped COHP in Fig. 4a), the band edges are formed of anti-bonding O
$p-$Mn $d$ states. This indicates that the electron doping will increase O
$p-$Mn $d$ $ab$ character and weaken Mn-O bonds, whereas hole doping will
decrease O $p-$Mn $d$ $ab$ character and strengthen Mn-O bonds. In detail, we
find the $p-d_{z}^{2}$ interactions at the bottom of the conduction band, and
$p-d_{xy}+p-d_{x^{2}-y^{2}}$ at the top of the valence band. Anti-bonding
$p-d_{xz}+p-d_{yz}$ interactions occur at $-2.5$ eV below the top of the
valence band.

On electron doping, electrons start occupying the $p-d_{z^{2}}$ dominated
yellow anti-bonding states at the bottom of the conduction band. This weakens
the planar Mn-Oeq bond, and is reflected in a reduction in the size of the
$p-d_{z}^{2}$ ICOHP (yellow diamonds in Fig. 4b)). This reduction in in-plane
bonding allows for larger $Q_{\text{K${}_{3}$}}$ polyhedral tilting and an in
turn increase in $P$.

On hole doping, the anti-bonding $p-d_{xy}/d_{x^{2}-y^{2}}$ bonds, at the top
of the valence band, are depleted, increasing the corresponding bond strength,
as seen in the increase in the size of the ICOHP (grey circles). This in turn
causes a stronger in-plane Mn-Oeq bond, consistent with a smaller MnO5
bipyramidal K3 tilt and a reduced $P$.

The $p-d_{xz}/d_{yz}$ interactions contribute almost negligibly to the change
in bond strength on either electron or hole doping, as their energy levels are
not close to the band edges. The ICOHP is close to zero for the range of
background doping we considered, and even at high doping value of 5 holes/unit
cell, the Mn $d_{xz}$ and $d_{yz}$ orbitals are only slightly depleted (Fig.
S3).

In summary we have shown that dopant electrons increase $P$, while dopant
holes suppress it. We additionally calculate that $Q_{\text{K${}_{3}$}}$
increases (decreases) with electron (hole) doping and that this $\Delta
Q_{\text{K${}_{3}$}}$ change is a result of the weakening (strengthening) of
the in-plane Mn-O bonding. This in turn explains the response in $\Delta P$
through the linear $Q_{\text{K${}_{3}$}}-Q_{\Gamma_{2}^{-}}$ coupling.

###  III.2 Impurity atoms

Figure 5: a) Contribution to $P$ from introduction of charge carriers and
changes in lattice constants. Faded turquoise level indicates the change in
$P$ due to the corresponding background charge doping concentration shown
previously in Fig. 2). b) Impurity atom effects in substitution compounds
which preserve P63cm symmetry. $Q_{\Gamma_{2}^{-}}$ versus
$Q_{\text{K${}_{3}$}}$ compared to the linear trend calculated for the
background charge method: $Q_{\Gamma_{2}^{-}}=0.4Q_{\Gamma_{2}^{-}}$.

In the second part of the study, we introduce charge carriers by replacing
either one Y3+ (A1 or A2 site atoms) or one Mn3+ (B site atom), by one
impurity atom per 30 atom unit cell. We choose the following substitution
atoms: Ca2+ (1 hole per unit cell), Mg2+ (1 hole per unit cell), Hf4+ (1
electron per unit cell), Ti4+ (1 electron per unit cell) and V5+ (2 electrons
per unit cell).

Substitution doping creates more complex changes in the crystal’s chemistry
and geometry than simple background charge doping. Therefore, in order to
disentangle the different contributions to the changes in $P$, we adopt the
three step approach introduced in Ref. [8] to separate the contributions from
the charge carrier (charge carrier effect), lattice constants (lattice
constant effect) and other factors resulting from the impurity atom (impurity
atom effect), as described in the Computational Details section. Our results
for the charge carrier and lattice constant contributions are shown in Fig.
5a) and we separately show the observed impurity atom effects in Fig. 5b). We
summarize our results from this section in Table 1, listing the ionic radii
and charge carrier contributions and the calculated corresponding lattice
constants, polar distortion $P$ as well as the change in polar distortion
$\Delta P$ from the unsubstituted case, for all the substitution atoms.

First we note that, across all substituted compounds (excluding Mg A2-site
substitution which is a special case that we will discuss later), the charge
carrier effect on the polar distortion (turquoise bars in Fig. 5a)) is
consistent with our results from the background-charge method (turquoise
shading in Fig. 5a)). We see that electron doping tends to promote the polar
distortion while hole doping suppresses it, and the sizes of the changes in
$P$ from the impurity atom charge carriers closely match those of the
background-charge calculations.

Secondly, concerning the lattice constant effect, we find that changes in
lattice constants (yellow bars in Fig. 5a)) contribute negligibly to the
change in polar distortion: they increase $P$ by a small amount for most
substitution compounds apart from the Hf-substituted compound, where we
calculate an associated small decrease in $P$. This correlates with a decrease
in $c$ in the Hf-substituted compound, compared to an increase in $c$ for all
the other substitution elements, as shown in Table 1. Our calculations are
consistent with previous X-ray diffraction experiments, correlating an
increase (decrease) in $c$ with an increase (decrease) in the $K_{3}$ mode
amplitude which in turn, through its coupling to $\Gamma_{2}^{-}$, promotes
(suppresses) the FE polar distortion [44]. Note the contrast to conventional
FEs such as BaTiO3, where lattice contributions are important, and where
elongations of the lattice parameter along the polar direction increase $P$
[8].

Finally, we analyze the change in $P$ caused by factors other than the charge
and lattice constant effects, referring to them collectively as impurity atom
effects. We find that the A1\- and B-site substituted compounds preserve the
P63cm symmetry of pure YMnO3: we analyze them in terms of the
$Q_{\text{K${}_{3}$}}$ and $Q_{\Gamma_{2}^{-}}$ ($\propto P$) described above.
We quantify the order parameter changes using AMPLIMODES [41], by replacing
the identity of the substitution element in the structure file with the
original Y or Mn ion of the host after the DFT relaxations, and present the
resulting changes in $Q_{\text{K${}_{3}$}}$ and $Q_{\Gamma_{2}^{-}}$ in Fig.
5b).

Fig. 5b) shows the change in $Q_{\Gamma_{2}^{-}}$ as a function of
$Q_{\text{K${}_{3}$}}$, both owing to the impurity atom effect. We observe
that the $(Q_{\text{K${}_{3}$}},P)$ coupling closely follows the linear
relationship predicted by Landau theory. Interestingly, the slope is the same
as obtained in our earlier background charge calculations
($P=0.4Q_{\text{K${}_{3}$}}$, plotted here again for comparison). This
suggests that the proportionality constant $\gamma$ (see Eq. I) is
characteristic for our system and does not change much when substituting on
the B site. We note also that the sizes of the changes in
$Q_{\text{K${}_{3}$}}$ arising from B-site impurity atom effects are
comparable to those stemming from the charge effect, and are related to the
radius of the substituting atom, with $Q_{\text{K${}_{3}$}}$ decreasing
(increasing) for B-site substitution atoms with smaller (larger) radii. For
A1-site substitution, $Q_{\text{K${}_{3}$}}$ is larger for smaller radius
substitution atoms as there is less hindrance to tilting but
$Q_{\Gamma_{2}^{-}}$ is smaller due to additional shifts of the A-site
cations. Our corresponding increases or decreases in $P$ are in agreement with
Ref. [45]

…  Dopant Site Impurity atom Charge carrier Lattice constants $P$ $\Delta P$
$\Delta P$ Ionic radius [Å] $V$ [Å3] $c$ [Å] [µC/cm2] [%] [µC/cm2] Undoped - -
- 360.4 11.33 8.06 0 0 Background - - 1h 360.4 11.33 6.64 $-0.2$ $-0.02$ Mg2+
A1 0.89 $\downarrow$ 1h 349.9 11.35 8.56 $+6.3$ $+0.51$ A2$\star$ 0.89
$\downarrow$ 1h 349.1 11.38 0.08 $-99.0$ $-7.98$ B 0.89 $\uparrow$ 1h 359.2
11.42 7.55 $-6.3$ $-0.51$ Ca2+ A1 1.12 $\uparrow$ 1h 357.9 11.42 7.20 $-10.7$
$-0.86$ A${}_{2}^{\star}$ 1h 357.9 11.43 5.61 $-30.4$ $-2.45$ Background - -
1e 360.4 11.33 8.64 $+0.1$ $+0.01$ Hf4+ A${}_{1}^{\star}$ 0.83 $\downarrow$ 1e
357.3 11.27 6.76 $-16.1$ $-1.30$ A2 1e 357.3 11.26 9.71 $+20.6$ $+1.66$ Ti4+ B
0.51 $\downarrow$ 1e 363.8 11.44 7.79 $-3.3$ $-0.27$ Background - - 2e 360.4
11.33 9.55 $+0.2$ $+0.02$ V5+ B 0.46 $\downarrow$ 2e 365.9 11.50 8.81 $+9.4$
$+0.76$

Table 1: Summary table of the calculated total changes in the substituted
compounds. We list substitution site, Shannon ionic radius [46], number and
type of introduced charge carriers, cell volume, $c$ lattice parameter and $P$
and change in $P$ for all considered doping elements. The radii of Y3+ and
Mn3+ are 1.019 Å and 0.58 Å, respectively [46]. $\star$ indicates the lowest
A-substitution site energy, and $\uparrow$ and $\downarrow$ indicate that the
substitution atom’s radius is larger or smaller than the substituted atom’s
radius, respectively.

A2-site substitution affects the $O_{a}x-A_{2}$ bond above the trimerization
center, disturbing the buckling of A layer and changing the structural
symmetry of the compound. Therefore, we do not analyze the structures in terms
of $Q_{\Gamma_{2}^{-}}$ and $Q_{\text{K${}_{3}$}}$. Instead, we note that, as
the radius of the substitution ion decreases, the substitution atom
increasingly rattles with respect to the other A2 site. $P$ either increases
or decreases, depending on the direction of the substitution atom’s
displacement along $c$. The Ca positive A2 site ion moves down along $c$,
giving a $30\%$ suppression of $P$, whereas A2-site Hf moves up resulting in a
$20\%$ increase in $P$. Finally, Mg A2 substitution is a special case, with
the radius of the Mg2+ ion accidentally resulting in the down-middle-up
A-layer configuration, characteristic of the PE P$\mathrm{bar}{3}$c symmetry
[47, 48, 49, 50].

As a concluding comment, we point out that for A-site impurity substitutions,
we find that the substitution site leading to lower $P$ is energetically more
favorable in each case than the higher $P$ compound; while Ca substitution
energetically favors the A2 site, Hf prefers to substitute on the A1 site, and
Mg substitution favors the PE A2-substituted compound.

##  IV Conclusions

In summary, we have calculated and rationalized the relationship between
doping and FE polar distortion in the improper FE h-YMnO3.

Our main finding is that the polar distortion persists in doped h-YMnO3, with
electron doping increasing $P$, and hole doping reducing it, both linearly
proportional to the doping concentration over the entire range studied. This
result is in stark contrast to the behavior of conventional FEs, where polar
distortion is suppressed for both hole and electron doping even at small
doping concentrations. Our analysis of the changes in chemical bonding on
doping indicates that the planar Mn-O bonds weaken (strengthen) upon electron
(hole) doping which in turn increases (decreases) the magnitude of the primary
order parameter of the FE structure.

This main result holds for both background-charge doping and doping via
explicit atomic substitution. In the case of atomic substitution, we find
additional contributions to the total $P$ from structural changes caused by
the changes in lattice parameters or the impurity atom effect itself, which,
in the latter case, can dominate over the change in electron or hole count.
Our calculations suggest that one should B-site impurity atom substitute with
electron charge carriers, such as Ti4+ or V5+, in order to enhance $P$.

Our findings indicate that doped improper geometric FEs are a promising
candidate for new polar metals. In the case of the multiferroic hexagonal
manganites, the additional magnetic degree of freedom could lead to an
intriguing interplay between metallicity, ferroelectricity and magnetization.
We hope that the results that we present here motivate future experimental and
theoretical work in this direction.

##  V Acknowledgments

We thank Veronica F. Michel and Tobias Esswein for useful discussions
regarding the computational details of this study. This work was funded by the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program project HERO grant agreement No. 810451.
Computational resources were provided by ETH Zürich and the Swiss National
Supercomputing center, project IDs eth3 and s1128.

## References

  * Rischau  _et al._ [2017] C. Rischau, X. Lin, C. Grams, D. Finck, S. Harms, T. Engelmayer, J. Lorenz, Y. Gallais, B. Fauqué, J. Hemberger, and K. Behnia, A ferroelectric quantum phase transition inside the superconducting dome of Sr1-xCaxTiO3-δ, Nature Physics 13, 643–648 (2017).
  * Puggioni  _et al._ [2018] D. Puggioni, G. Giovannetti, and J. M. Rondinelli, Polar metals as electrodes to suppress the critical-thickness limit in ferroelectric nanocapacitors, J. Appl. Phys. 124, 174102 (2018), https://doi.org/10.1063/1.5049607 .
  * Zhou and Ariando [2020] W. X. Zhou and A. Ariando, Review on ferroelectric/polar metals, Jpn. J. Appl. Phys. 59, 10.35848/1347-4065/ab8bbf (2020).
  * Takahashi  _et al._ [2017] K. S. Takahashi, Y. Matsubara, M. S. Bahramy, N. Ogawa, D. Hashizume, Y. Tokura, and M. Kawasaki, Polar metal phase stabilized in strained La-doped BaTiO3 films, Sci Rep 7, 4631 (2017).
  * Li  _et al._ [2021] G. Li, K. Hu, S. Peng, Q. Yang, Y. Xiong, W. Zhang, Y. Xiao, S. Yan, M. Tang, and Z. Li, Strain-tunable giant anisotropic conductivity in the polar metal of KNbO3/BaTiO3 ferroelectric superlattice, Comput. Mater. Sci. 188, 110235 (2021).
  * Ma  _et al._ [2018] C. Ma, K.-J. Jin, C. Ge, and G.-Z. Yang, Strain-engineering stabilization of BaTiO3-based polar metals, Phys. Rev. B 97, 115103 (2018).
  * Raghavan  _et al._ [2016] S. Raghavan, J. Y. Zhang, O. F. Shoron, and S. Stemmer, Probing the metal-insulator transition in BaTiO3, Phys. Rev. Lett. 117, 037602 (2016).
  * Michel  _et al._ [2021] V. F. Michel, T. Esswein, and N. A. Spaldin, Interplay between ferroelectricity and metallicity in BaTiO3, J. Mater. Chem. C 9, 8640 (2021).
  * Hickox-Young  _et al._ [2023] D. Hickox-Young, D. Puggioni, and J. M. Rondinelli, Polar metals taxonomy for materials classification and discovery, Phys. Rev. Mater. 7, 010301 (2023).
  * Bhowal and Spaldin [2023] S. Bhowal and N. A. Spaldin, Polar metals: Principles and prospects, Annu. Rev. Mater. Res. 53, 10.1146/annurev-matsci-080921-105501 (2023), https://doi.org/10.1146/annurev-matsci-080921-105501 .
  * Zhao  _et al._ [2018] H. J. Zhao, A. Filippetti, C. Escorihuela-Sayalero, P. Delugas, E. Canadell, L. Bellaiche, V. Fiorentini, and J. Íñiguez, Meta-screening and permanence of polar distortion in metallized ferroelectrics, Phys. Rev. B 97, 054107 (2018).
  * Lichtenberg  _et al._ [2001] F. Lichtenberg, A. Herrnberger, K. Wiedenmann, and J. Mannhart, Synthesis of perovskite-related layered AnBnO3n+2 = ABOX type niobates and titanates and study of their structural, electric and magnetic properties, Prog. Solid. State Chem. 29, 1 (2001).
  * Benedek and Fennie [2011] N. A. Benedek and C. J. Fennie, Hybrid improper ferroelectricity: A mechanism for controllable polarization-magnetization coupling, Phys. Rev. Lett. 106, 107204 (2011).
  * Lei  _et al._ [2018] S. Lei, M. Gu, D. Puggioni, G. Stone, J. Peng, J. Ge, Y. Wang, B. Wang, Y. Yuan, K. Wang, Z. Mao, J. M. Rondinelli, and V. Gopalan, Observation of quasi-two-dimensional polar domains and ferroelastic switching in a metal, Ca3Ru2O7, Nano Lett. 18, 3088 (2018).
  * Levanyuk and Sannikov [1974] A. P. Levanyuk and D. G. Sannikov, Improper ferroelectrics, Sov. phys., Usp. , 199 (1974).
  * Zhang  _et al._ [2023] H. Zhang, K. Yuan, H. Tang, P. Zheng, W. Zhou, C. Wang, and W. Liu, The optical and ferroelectric properties of fe-doped hybrid improper ferroelectricity Ca3Ti2O7 via experimental and density functional theory study, J. Solid State Chem. 325, 124180 (2023).
  * Huang  _et al._ [2019] C. Huang, W. Wong-Ng, W. Liu, X. Zhang, Y. Jiang, P. Wu, B. Tong, H. Zhao, and S. Wang, Major improvement of ferroelectric and optical properties in na-doped ruddlesden-popper layered hybrid improper ferroelectric compound, Ca3Ti2O7, J. Alloys Compd. 770, 582 (2019).
  * Artyukhin  _et al._ [2014] S. Artyukhin, K. T. Delaney, N. A. Spaldin, and M. Mostovoy, Landau theory of topological defects in multiferroic hexagonal manganites, Nat. Mater. 13, 42 (2014).
  * Meier  _et al._ [2017] Q. N. Meier, M. Lilienblum, S. M. Griffin, K. Conder, E. Pomjakushina, Z. Yan, E. Bourret, D. Meier, F. Lichtenberg, E. K. H. Salje, N. A. Spaldin, M. Fiebig, and A. Cano, Global formation of topological defects in the multiferroic hexagonal manganites, Phys. Rev. X 7, 041014 (2017).
  * Skjærvø  _et al._ [2019] S. H. Skjærvø, Q. N. Meier, M. Feygenson, N. A. Spaldin, S. J. L. Billinge, E. S. Bozin, and S. M. Selbach, Unconventional continuous structural disorder at the order-disorder phase transition in the hexagonal manganites, Phys. Rev. X 9, 031001 (2019).
  * Meier  _et al._ [2020] Q. N. Meier, A. Stucky, J. Teyssier, S. M. Griffin, D. van der Marel, and N. A. Spaldin, Manifestation of structural Higgs and Goldstone modes in the hexagonal manganites, Phys. Rev. B 102, 014102 (2020).
  * Nordlander  _et al._ [2019] J. Nordlander, M. Campanini, and R. M. et al, The ultrathin limit of improper ferroelectricity, Nat. Commun. 10, 10.1038/s41467-019-13474-x (2019).
  * Łukaszewicz and Karut-Kalicińska [1974] K. Łukaszewicz and J. Karut-Kalicińska, X-ray investigations of the crystal structure and phase transitions of YMnO3, Ferroelectrics 7, 81 (1974).
  * Gibbs  _et al._ [2011] A. S. Gibbs, K. S. Knight, and P. Lightfoot, High-temperature phase transitions of hexagonal YMnO3, Phys. Rev. B 83, 094111 (2011).
  * Fennie and Rabe [2005] C. J. Fennie and K. M. Rabe, Ferroelectric transition in $\mathrm{Y}\mathrm{Mn}{\mathrm{O}}_{3}$ from first principles, Phys. Rev. B 72, 100103 (2005).
  * Dronskowski and Bloechl [1993] R. Dronskowski and P. E. Bloechl, Crystal orbital hamilton populations (COHP): energy-resolved visualization of chemical bonding in solids based on density-functional calculations, J. Phys. Chem. 97, 8617 (1993).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal-amorphous-semiconductor transition in germanium, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996b).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981).
  * Cococcioni and de Gironcoli [2005] M. Cococcioni and S. de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the $\mathrm{LDA}+\mathrm{U}$ method, Phys. Rev. B 71, 035105 (2005).
  * Bruneval  _et al._ [2015] F. Bruneval, C. Varvenne, J.-P. Crocombette, and E. Clouet, Pressure, relaxation volume, and elastic interactions in charged simulation cells, Phys. Rev. B 91, 024107 (2015).
  * Van Aken  _et al._ [2004] B. B. Van Aken, T. Palstra, A. Filippetti, and N. A. Spaldin, The origin of ferroelectricity in magnetoelectric YMnO3, Nat. Mater. 3, 164 (2004).
  * Filippetti and Spaldin [2003] A. Filippetti and N. A. Spaldin, Strong-correlation effects in born effective charges, Phys. Rev. B 68, 045111 (2003).
  * Maintz  _et al._ [2013] S. Maintz, V. L. Deringer, A. L. Tchougréeff, and R. Dronskowski, Analytic projection from plane-wave and PAW wavefunctions and application to chemical-bonding analysis in solids, J. Comput. Chem. 34, 2557 (2013).
  * Maintz  _et al._ [2016] S. Maintz, V. L. Deringer, A. L. Tchougréeff, and R. Dronskowski, LOBSTER: A tool to extract chemical bonding from plane-wave based DFT (2016).
  * Nelson  _et al._ [2020] R. Nelson, C. Ertural, J. George, V. L. Deringer, G. Hautier, and R. Dronskowski, LOBSTER: Local orbital projections, atomic charges, and chemical-bonding analysis from projector-augmented-wave-based density-functional theory, J. Comput. Chem. 41, 1931 (2020).
  * Van Aken  _et al._ [2001] B. B. Van Aken, A. Meetsma, and T. T. M. Palstra, Hexagonal LuMnO3 revisited, Acta Crystallogr. Section E 57, i101 (2001).
  * Orobengoa  _et al._ [2009] D. Orobengoa, C. Capillas, M. I. Aroyo, and J. M. Perez-Mato, Amplimodes: symmetry-mode analysis on the bilbao crystallographic server, J. Appl. Crystallogr. 42, 820 (2009).
  * Perez-Mato  _et al._ [2010] J. Perez-Mato, D. Orobengoa, and M. Aroyo, Mode crystallography of distorted structures, Acta Crystallogr. Section A 66, 558 (2010).
  * Choi  _et al._ [2008] W. S. Choi, D. G. Kim, S. S. A. Seo, S. J. Moon, D. Lee, J. H. Lee, H. S. Lee, D.-Y. Cho, Y. S. Lee, P. Murugavel, J. Yu, and T. W. Noh, Electronic structures of hexagonal RMnO3 (R=Gd, Tb, Dy, and Ho) thin films: Optical spectroscopy and first-principles calculations, Phys. Rev. B 77, 045137 (2008).
  * Katsufuji  _et al._ [2002] T. Katsufuji, M. Masaki, A. Machida, M. Moritomo, K. Kato, E. Nishibori, M. Takata, M. Sakata, K. Ohoyama, K. Kitazawa, and H. Takagi, Crystal structure and magnetic properties of hexagonal RMnO3 (R=Y, Lu, and Sc) and the effect of doping, Phys. Rev. B 66, 134434 (2002).
  * Småbråten  _et al._ [2020] D. R. Småbråten, T. S. Holstad, D. M. Evans, Z. Yan, E. Bourret, D. Meier, and S. M. Selbach, Domain wall mobility and roughening in doped ferroelectric hexagonal manganites, Phys. Rev. Res. 2, 033159 (2020).
  * Shannon [1976] R. D. Shannon, Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides, Acta Crystallogr. Section A 32, 751 (1976).
  * Huang  _et al._ [2013] F.-T. Huang, X. Wang, Y. S. Oh, K. Kurushima, S. Mori, Y. Horibe, and S.-W. Cheong, Delicate balance between ferroelectricity and antiferroelectricity in hexagonal InMnO3, Phys. Rev. B 87, 184109 (2013).
  * Kumagai  _et al._ [2012] Y. Kumagai, A. A. Belik, M. Lilienblum, N. Leo, M. Fiebig, and N. A. Spaldin, Observation of persistent centrosymmetricity in the hexagonal manganite family, Phys. Rev. B 85, 174422 (2012).
  * Huang  _et al._ [2014] F.-T. Huang, X. Wang, S. M. Griffin, Y. Kumagai, O. Gindele, M.-W. Chu, Y. Horibe, N. A. Spaldin, and S.-W. Cheong, Duality of topological defects in hexagonal manganites, Phys. Rev. Lett. 113, 267602 (2014).
  * Griffin  _et al._ [2017] S. M. Griffin, M. Reidulff, S. Selbach, and N. A. Spaldin, Defect chemistry as a crystal structure design parameter: Intrinsic point defects and Ga substitution in InMnO3, Chem. Mater. 29, 2425–2434 (2017).