# Data-Driven Flow and Injection Estimation in PMU-Unobservable Transmission
Systems

Satyaprajna Sahoo, Student Member, IEEE, Anwarul Islam Sifat, Member, IEEE,
and Anamitra Pal, Senior Member, IEEE This work was supported in part by the
National Science Foundation (NSF) grant under Award ECCS-2132904. The authors
are with the School of Electrical, Computer, and Energy Engineering (ECEE) at
Arizona State University (ASU).

###### Abstract

Fast and accurate knowledge of power flows and power injections is needed for
a variety of applications in the electric grid. Phasor measurement units
(PMUs) can be used to directly compute them at high speeds; however, a large
number of PMUs will be needed for computing all the flows and injections.
Similarly, if they are calculated from the outputs of a linear state
estimator, then their accuracy will deteriorate due to the quadratic
relationship between voltage and power. This paper employs machine learning to
perform fast and accurate flow and injection estimation in power systems that
are sparsely observed by PMUs. We train a deep neural network (DNN) to learn
the mapping function between PMU measurements and power flows/injections. The
relation between power flows and injections is incorporated into the DNN by
adding a linear constraint to its loss function. The results obtained using
the IEEE 118-bus system indicate that the proposed approach performs more
accurate flow/injection estimation in severely unobservable power systems
compared to other data-driven methods.

###### Index Terms:

Flow and Injection estimation, Machine learning (ML), Phasor measurement unit
(PMU), Unobservability.

##  I Introduction

Knowledge of active power flows and power injections is fundamental for the
reliable, resilient, and economic operation of the electric grid.
Traditionally, their knowledge had been used to determine the cost for
buying/selling energy as well as for performing vulnerability/security
assessment [1, 2]. More recently, with the growing frequency and intensity of
extreme weather events as well as increasing penetration of renewable energy
resources, a need has been felt for high-speed tracking of the power flowing
through critical equipment, such as transformers, to determine their health
[3].

Power flows and injections can be directly computed at high speeds from the
outputs of a phasor measurement unit (PMU). However, a large number of PMUs
will be needed to estimate all of them. They can be calculated indirectly from
the outputs of a linear state estimator. However, the quadratic relationship
between voltage and power deteriorates the quality of flow/injection estimates
considerably (see Table I in Section II for a sample illustration). This paper
presents a physics-inspired machine learning (ML) formulation to perform data-
driven power flow and injection estimation in PMU-unobservable power systems.
We focus on the transmission system since most of them already have some PMUs
installed. However, the proposed methodology is generic enough that it can be
applied to distribution systems as well.

Literature Review: In [4], a minimum variance unbiased estimator was developed
to calculate the power flows using the state estimation approach. However, the
analysis was restricted to DC power flows. A linearly constrained least-
squares optimization problem was formulated in [5] to estimate nodal power
injections that did not violate power flow constraints. However, no
information was provided regarding how the formulation could perform high
speed estimation in unobservable power systems. A methodology to learn the
topology and estimate the injection statistics in distribution systems with
unobservable nodes at PMU timescales was developed in [6]. However, it
required the unobservable nodes to be non-adjacent, which may not always be
the case.

Main Contributions: This paper proposes the use of ML to simultaneously
estimate all power flows and injections in a system that is sparsely observed
by PMUs. The ML model is built using deep neural networks (DNNs). The flows
and injections in a power system are related by the law of conservation of
energy. This knowledge is embedded into the DNNs as a linear constraint; the
resulting ML model is referred to as a physics-inspired constrained-DNN (PIC-
DNN). The offline training is performed using slow timescale historical data
obtained from the supervisory control and data acquisition (SCADA) system. The
online implementation only uses PMU data, ensuring high speed of estimation.
The superior performance of the proposed model over classical as well as other
ML models is demonstrated for different numbers of PMUs placed in the IEEE
118-bus system.

##  II Sources of Error in Power calculated from State estimates

The most traditional way of determining electrical quantities in a power
system is through state estimation. The widely used static state estimator
provides an estimate of the voltage phasor (magnitude and angle) of every bus
of the system. When a system is completely observed by PMUs, then by only
employing PMU data one can perform linear state estimation (LSE); LSE is
faster and more accurate than SCADA-based state estimation [7]. However, there
are two issues with computing power flows/injections from the outputs of LSE:
(a) quadratic relationship between voltage and power, and (b) non-Gaussian
noise in PMU measurements. These two issues are elaborated below.

Quadratic Relationship: Power flowing through a branch located between bus $i$
and bus $j$, denoted by $p_{ij}$, is related to the voltages of the two buses
by the equation shown below:

$p_{ij}=g_{ij}({v_{i}}^{2}-v_{i}v_{j}\cos{\delta_{ij}})-b_{ij}(v_{i}v_{j}\sin{\delta_{ij}})$
(1)

In (1), $v_{i}$ and $v_{j}$ denote the voltage magnitudes of buses $i$ and
$j$, $\delta_{ij}$ denotes the voltage angle difference between the two buses,
and $g_{ij}$ and $b_{ij}$ are the real and imaginary components of the
admittance matrix, respectively. Similarly, the power injection at bus $i$,
denoted by $p_{i}$, can be expressed as,

$p_{i}=\sum_{j=1}^{br_{i}}p_{ij}=\sum_{k=1}^{bu}v_{i}v_{k}(g_{ik}\cos{\delta_{ik}}+b_{ik}\sin{\delta_{ik}})$
(2)

where $br_{i}$ denotes the total number of branches incident on bus $i$, and
$bu$ denotes the total number of buses in the system. It is clear from (1) and
(2) that power is proportional to the square of the voltage. Therefore, even a
small error in the voltage estimates will get amplified when used for
calculating the power flows and power injections.

Non-Gaussian Noise: LSE is most commonly performed using the least squares
method. However, least squares is the solution to the maximum likelihood
estimation problem under Gaussian noise environments [8]. Recently, it has
been demonstrated that PMU data has non-Gaussian noise [9]. This, in
combination with the previously-mentioned source of error, implies that an
indirect estimation of power flows/injections (by finding the voltages first
using state estimation) can result in lower accuracy. This implication is
further reinforced through a numerical analysis performed on the IEEE 118-bus
system. The results are provided in Table I below.

Table I shows the estimation performance (in terms of root mean square error)
when the output of a linear state estimator is used to calculate the power
flows and injections in the presence of no noise, 1% total vector error (TVE)
Gaussian noise, and 1% TVE non-Gaussian noise. The non-Gaussian noise is
described by a two-component Gaussian mixture model (GMM); its parameters are
given in the last paragraph of Section IV-A. It is observed from the table
that the estimation error increased by at least $60\%$ when the distribution
of the noise changed from Gaussian to non-Gaussian. Furthermore, a purely PMU-
based LSE requires complete observability of the system by PMUs (for obtaining
the results shown in Table I, PMUs were placed at 32 optimal locations that
completely observed the system [10]). This requirement of complete
observability is another concern of the LSE-based approach.

TABLE I: LSE results for voltage and power estimates in IEEE 118-bus system | No noise |  | Gaussian
---
noise (1% TVE)
| Non-Gaussian
---
noise (1% TVE)
| Voltage
---
magnitude (pu)
4.36 e-13 | 0.0036 | 0.0060
| Voltage angle
---
(degrees)
1.88 e-11 | 0.1764 | 0.2890
| Power
---
flow (MW)
1.46 e-10 | 1.3382 | 2.2273
| Power
---
injection (MW)
2.69 e-10 | 1.3613 | 2.2912

Note that a PMU measures voltage and current phasors at the location where it
is placed (subject to its measurement channel limitations [11]). This means
that the outputs of a PMU can be used to directly estimate all the power flows
(and subsequently the power injections using (2)). However, this translates to
the well-known minimum vertex cover problem [12], which requires placing even
more PMUs than those required for complete observability for LSE (which is the
solution to the minimum dominating set problem). For example, PMUs must be
placed at 61 locations in the IEEE 118-bus system to estimate all the power
flows using only PMU data. In this paper, we employ ML to directly estimate
all the power flows and injections at PMU timescales and with reasonable
accuracy while placing significantly fewer PMUs.

##  III Physics-Inspired ML for Power estimation in PMU-unobservable Systems

Data-driven approaches have been shown to attain considerable success for a
variety of power system estimation problems [13]. Some popular data-driven
approaches include linear regression (LR), support vector regression (SVR),
and DNN-based regression. As the linear relation between the input (PMU
measurements) and output (power flows and power injections) for transmission
systems that are sparsely observed by PMUs is not guaranteed, LR may not be a
good fit for the problem considered in this paper. Similarly, the time
complexity of SVR is quadratic w.r.t. the number of training samples [14],
which makes it challenging to implement it in a large network with wide
variations in features. Therefore, we employ DNNs to perform direct estimation
of power flows and injections from PMU measurements.

In Section III-A, we show how DNN-based regression can be used to overcome the
unobservability problem associated with purely PMU-based estimation. In
Section III-B, we embed the physical law that relates flows and injections
into the DNN framework to create the proposed PIC-DNN.

###  III-A DNN Regression to Overcome Unobservability

Recently, we have demonstrated the ability of DNNs to perform high-speed time-
synchronized static state estimation in distribution systems that are
sparingly observed by micro-PMUs [15, 16]. We did this by using the slow
timescale historical smart meter data to create a mapping function between the
fast timescale micro-PMU measurements and the states. The mapping function was
learned using a DNN because it has excellent approximation capabilities. Here,
we leverage this concept to perform power flow and injection estimation from
sparsely-placed PMUs in transmission systems.

We start by making two assumptions: (1) we have access to historical SCADA
data and system information (e.g., topology) for a sufficient time-period
(say, a few months), and (2) PMUs are already placed on select buses of the
system (say, highest voltage buses). Using historical SCADA data and system
information, we solve the power flow problem to produce voltage and current
phasor measurements corresponding to the locations where PMUs are placed. The
power flow problem also generates the flow and injection information that
matches the PMU measurements. Then, we train a DNN whose inputs are the PMU
measurements and outputs are the flows and injections. Finally, during online
implementation, streaming data from a select few PMUs is fed into the trained
DNN to estimate all the flows and injections at PMU timescales. In this way, a
DNN can perform purely PMU-based estimation without needing PMUs to completely
observe the system. This DNN model is henceforth referred to as a Direct DNN.

###  III-B Incorporating Physics-based Constraints into DNN

The Direct DNN model developed in Section III-A can estimate the power flows
and power injections independently from the phasor measurements coming from a
select few PMUs. However, as seen in (2), the flows and injections are related
by the law of conservation of energy. To account for this physical law, we
modify the DNN architecture by appending a linear constraint to its loss
function; the resulting model, called PIC-DNN, is described below. In the
following, small-case letters indicate vectors, upper-case letters indicate
matrices, and the symbol $\hat{}$ indicates estimates.

For an input, $z$, and output, $y$, a regular DNN with weights, $W$, and
biases, $b$, tries to find a function, $F$, that minimizes the difference
between $y$ and $\hat{y}$, where $\hat{y}=F(z,W,b)$. The Direct DNN model
described in Section III-A performs this minimization using a mean squared
error (MSE) loss function, where $z=[V\;;\;I]_{n\times 1}$ and
$y=[P_{f}\;;\;P_{in}]_{m\times 1}$. To incorporate the law of conservation of
energy into the Direct DNN model, we start by defining the following loss
function:

$\displaystyle\min_{W,b}$ $\displaystyle(y-\hat{y})^{2}$ (3) s.t.
$\displaystyle\hat{P}_{in}=A\times\hat{P}_{f}$ where $\displaystyle
A=\begin{bmatrix}a_{11}&a_{12}&\ldots&a_{1m_{b}}\\\
a_{21}&a_{22}&\ldots&a_{2m_{b}}\\\ \vdots&&\ddots&\vdots\\\
a_{(m-m_{b})1}&\ldots&&a_{(m-m_{b})m_{b}}\end{bmatrix}$

In (3), $A$ is a sparse matrix ($a_{ij}\in\\{0,1\\}^{(m-m_{b})\times m_{b}}$)
which denotes the power flows that must be added to calculate the power
injection of a bus, and $m_{b}$ and $m-m_{b}$ denote the number of branch
power flows and bus power injections in the system, respectively. Now, power
injections can be removed from the output features of the DNN and calculated
from the power flow estimates in the following way:

$\hat{y}=B\times\hat{P}_{f},\;\;\;B=[\mathbb{I}_{m_{b}}\;;\;A]\in\\{0,1\\}^{m\times
m_{b}}$ (4)

The conversion matrix, $B$, is the vertical concatenation of an appropriately
sized identity matrix with $A$, and enables us to get both flow and injection
estimates. Thus, by modifying the loss function of the conventional DNN in the
manner shown in (3) and (4), we are able to simultaneously minimize the error
in the estimates of power flows and power injections while also accounting for
the law that relates the two.

There are several ways to implement the above-mentioned constrained
formulation inside a DNN model. A simple strategy would be to add a convex
optimization layer to the model output. However, the inclusion of such a layer
will result in solving two different optimization problems simultaneously
during training, which can become computationally burdensome. A more intuitive
strategy is to add a static layer after the DNN output layer that calculates
the combined $\hat{y}$ from the estimated $\hat{P}_{f}$ and backpropagates the
net error. The updated loss function is given by,

$\displaystyle\min_{W,b}$
$\displaystyle(B^{T}B)^{-1}B^{T}({y}-B\times\hat{P}_{f})^{2}$ (5) where
$\displaystyle\hat{P}_{f}=F(z,W,b)$

The static weight, $(B^{T}B)^{-1}B^{T}$, is the Moore–Penrose inverse of the
conversion matrix, $B$. Using (5), we are not only able to reduce the number
of DNN output variables to just the power flows but also implicitly include
the linear constraints into the loss function without needing a separate
convex optimization layer.

Next, the training dataset, $y_{train}$, which is only composed of the branch
power flows, is split into multiple bins in accordance with the variations
that occur in the flows. This improves the accuracy of the DNN because it now
has to estimate features that have similar variations. The selection of the
optimal number of bins, $n_{bin}$, involves a trade-off. Increasing the number
of bins homogenizes the output data, thereby improving the accuracy of the
model. However, unrestricted increment in bin count can increase computational
burden while decreasing performance utility.

In summary, the two ways in which we modify the Direct DNN model to create the
PIC-DNN model are: (1) modifying the loss function to include a linear
constraint, and (2) binning $y_{train}$ in accordance with the observed
variations in the power flows in the training database. The overall
implementation of PIC-DNN is shown in Algorithm 1. The results obtained when
the PIC-DNN as well as other ML models are used to estimate power flows and
injections in the IEEE 118-bus system is provided in the next section.

Algorithm 1 Implementation of PIC-DNN

Input: Noisy measurements, $z$
Output: Power flow and injection estimates, $\hat{y}$

1:procedure

2: Define $A$, $B$, $n_{bin}$, $epoch$

3: Split $y_{train}$ into $n_{bin}$ bins

4: for i=1 to $epoch$ do

5: for j=1 to $n_{bin}$ do

6: $\hat{p}_{j}\leftarrow F_{j}(z,W_{j},b_{j})$

7: end for

8: $\hat{P_{f}}\leftarrow[\hat{p}_{1},\hat{p}_{2},...,\hat{p}_{n_{bin}}]$

9: $\hat{y}\leftarrow B\times\hat{P}_{f}$

10: $\delta\leftarrow(B^{T}B)^{-1}B^{T}(y-\hat{y})^{2}$

11: backpropagate $\delta$

12: end for

13:end procedure

##  IV Results

The proposed PIC-DNN model for power flow and injection estimation was applied
to the IEEE 118-bus system. This system has 99 loads, 54 generators, 11 high
voltage (HV) buses, and 186 branches. PMUs were assumed to be placed by
default on the 11 HV buses such that all the branches coming out of these
buses was directly monitored by them. To obtain branch power flows and bus
power injections for this system, we solved an AC optimal power flow (ACOPF)
using MATPOWER [17]. The process is similar to the approach in [16], where a
distribution kernel was fit over historical slow timescale data, and then
multiple samples were drawn from it to generate realistic load variation data.
For the given application, the source of slow timescale data was the SCADA
system. However, SCADA data is not available for the IEEE 118-bus system.
Therefore, we superimposed the variations of similar loads found in the
publicly available 2000-bus Synthetic Texas system [18] onto loads of the IEEE
118-bus system. Doing so ensured that our load variations were realistic.
Afterward, the outputs of the ACOPF were used to train the ML models. The
training and validation database had a size of $20,000\times n$ and
$2,000\times n$, respectively, while the test database had a size of
$6,000\times n$, where $n$ is the number of phasor measurements.

The performance of five ML models was investigated in this study. The first
two models are LR and SVR, with the latter implemented using a radial basis
function kernel. The third model was an Indirect DNN. This model first
performs state estimation from sparsely-placed PMUs using the methodology
proposed in [16], and then computes the flows and injections from the outputs
of the state estimator. The fourth model was the Direct DNN described in
Section III-A, while the last model was the PIC-DNN developed in Section
III-B. The hyperparameters used for the three DNN models are given in Table
II. The loss function for the Direct and Indirect DNN models is the MSE, while
for the PIC-DNN model we used (5), with $n_{bin}=5$. Note that the number of
neurons in each hidden layer is kept flexible (using the scaling factor,
$\eta$) to account for the increase in number of input features due to
increase in number of PMUs (see Section IV-B).

TABLE II: Hyperparameters of DNNs created in this study Parameter | Value
---|---
| Layers
---
3
| Neurons in each hidden layer
---
Number of input features $\times\eta$
| Learning rate
---
1e-3
| Number of epochs
---
200
| Activation function
---
Rectified Linear Unit
| Optimizer
---
Adam
| Batch size
---
64

The simulations done to compare the performance of the five ML models were
conducted on a computer with an Intel Core (TM) i7-11800H CPU @2.3GHz with
16GB of RAM and an RTX 3070Ti GPU. The following two case-studies were
designed. The first case-study is a comparative analysis between PIC-DNN and
other ML models when PMUs are placed on only the 11 HV buses of the test
system. In the second case-study, the performance of all five ML models is
investigated as more PMUs are added to the system.

###  IV-A Case-Study I: Performance Comparison with PMUs placed only on HV
buses

In this case-study, the inputs to the ML models are the voltage and currents
obtained from PMUs placed on the 11 HV buses of the IEEE 118-bus system.
Therefore, the number of phasor measurements is 41 (= 11 voltage phasors + 30
current phasors), i.e., $n=82$, since one phasor has two components, while the
number of active power flows and injections are 490 (= 372 power flows + 118
power injections), i.e., $m=490$. To ensure the consistency of estimates
obtained from the five ML models, the analysis was repeated 100 times with
different random subsets of the test dataset.

The statistics of the error metric (RMSE) for the five ML models are recorded
in Table III. The results show that PIC-DNN does better than the other ML
models in both the mean value of RMSE across all the flows and injections as
well as the variation of RMSE over the 100 trials. For example, the proposed
PIC-DNN outperforms others in terms of the mean by at least $15\%$ and in
terms of standard deviation by at least $40\%$; this denotes a significant
improvement in both accuracy and consistency of power flow and injection
estimates. The poor performance of the Indirect DNN model can be attributed to
the quadratic relationship between voltage and power, as highlighted in
Section II. The relatively poor performance of SVR is due to the input dataset
being restricted to $5,000$ samples. This was done to compensate for its
(longer) training time, indicating the handicap of using SVR for large
datasets as its scalability is a concern [14].

Lastly, note that the results shown in Table III (as well as the last column
of Table I and Fig. 1) were obtained when a two-component GMM noise was added
to the voltage and current phasors. The mean, standard deviation, and weights
of this noise for magnitudes and angles are $(-0.4\%,0.6\%)$ and
$(-0.2\degree,0.3\degree)$, $(0.25\%,0.25\%)$ and $(0.12\degree,0.12\degree)$,
and $(0.4,0.6)$, respectively. These values corresponded to a $1\%$ TVE error
for voltage and current phasors, as mandated in the IEEE/IEC Standard for PMUs
[19].

TABLE III: Performance comparison of ML models for power estimation in terms of their root mean square error (RMSE) ML Model | Mean of RMSE |  | Standard deviation of RMSE
---
| LR
---
5.80 | 0.020
| SVR
---
7.75 | 0.120
| Indirect DNN
---
8.37 | 0.020
| Direct DNN
---
5.91 | 0.024
| PIC-DNN
---
4.92 | 0.012

###  IV-B Case-Study II: Impact of Increase in Number of PMUs

In this case-study, the number of locations where PMUs must be placed is
increased one at a time for each of the five ML models to determine how their
performance changes with an increase in the number of PMUs in the IEEE 118-bus
system. This is a very practical scenario since, with additional investment in
grid modernization, power utilities will add more PMUs not only in the lower
voltages of the transmission system but also in the distribution system.
Therefore, this case-study indicates the improvements in ML-based power flow
and injection estimation with an increase in PMU coverage.

The results obtained as the number of locations where PMUs must be placed
increased from 11 to 32 is shown in Fig. 1. To ensure fairness in comparison,
the next optimal bus location for placing PMUs is determined independently for
each ML model. This means that the horizontal axis of Fig. 1 denotes the total
number of buses where PMUs are placed; however, the bus names that go with
that number can be different for the five ML models. We stopped at 32 because
by placing PMUs on that many buses, we could directly compare the performance
of the different ML models with the conventional LSE results (shown in Table
I).

From Fig. 1, it can be realized that the proposed PIC-DNN performs better than
the other four ML models for higher degrees of unobservability (fewer PMUs).
With an increase in the number of PMUs, the performance of LR becomes
comparable with the proposed approach (and even better than PIC-DNN for PMUs
placed at 30 or more bus locations). This is expected because linear models
can successfully express the relations between the inputs and the outputs as
the number of sensors increase. Lastly, the superiority of data-driven
approaches is realized from the fact that when the number of bus locations
where PMUs are placed is 32, both LR and PIC-DNN outperform LSE results with
two-component GMM measurement noise by $35\%$ and $28\%$, respectively.

Figure 1: Performance comparison of ML models with increase in number of PMUs

##  V Conclusion

This paper proposes a physics-inspired ML approach that uses DNNs to quickly
and accurately estimate all the power flows and injections directly from PMUs
placed in the transmission system. The proposed PIC-DNN model not only
performs consistent estimation when there are very few PMUs in the system but
also ensures that the law of conservation of energy is always satisfied. The
former is established by intelligently combining inferences drawn from
historical SCADA data with real-time PMU data, while the latter is secured by
adding a linear constraint to the loss function of the DNN. The results
indicate that by placing PMUs on only 11 buses of the IEEE 118-bus system,
PIC-DNN can estimate the flows and injections to within 5 MW of the actual
value, even in the presence of non-Gaussian noise in the PMU measurements.
When the PMU locations increase to 32, the estimation error of the proposed
approach becomes lower than 2 MW.

The focus of this paper has been on the transmission system. In the future, we
plan to extend our analysis to distribution systems where PMUs are gradually
being added and for which fast detection of reverse power flows is crucial. We
will also combine the outcomes of this research to perform a better security
assessment of the power system (e.g., by extending the research of [20]). The
DNN models analyzed here had fully-connected feed-forward architectures. We
are currently investigating the ability of graph neural networks to better
incorporate the physical properties of the power system into ML training and
execution (e.g., for ensuring robustness during topology changes).

## References

  * [1] Y. C. Chen and S. V. Dhople, “Tracing power with circuit theory,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 1, pp. 138–147, 2020.
  * [2] R. S. Biswas, A. Pal, T. Werho, and V. Vittal, “A graph theoretic approach to power system vulnerability identification,” _IEEE Transactions on Power Systems_ , vol. 36, no. 2, pp. 923–935, 2020.
  * [3] A. Azmi, J. Jasni, N. Azis, and M. A. Kadir, “Evolution of transformer health index in the form of mathematical equation,” _Renewable and Sustainable Energy Reviews_ , vol. 76, pp. 687–700, 2017.
  * [4] M. Amini, A. I. Sarwat, S. Iyengar, and I. Guvenc, “Determination of the minimum-variance unbiased estimator for DC power-flow estimation,” in _IECON 2014-40th Annual Conference of the IEEE Industrial Electronics Society_. IEEE, 2014, pp. 114–118.
  * [5] A. Al-Digs, S. V. Dhople, and Y. C. Chen, “Estimating feasible nodal power injections in distribution networks,” in _2016 IEEE Power & Energy Society Innovative Smart Grid Technologies Conference (ISGT)_. IEEE, 2016, pp. 1–5.
  * [6] D. Deka, M. Chertkov, and S. Backhaus, “Joint estimation of topology and injection statistics in distribution grids with missing nodes,” _IEEE Transactions on Control of Network Systems_ , vol. 7, no. 3, pp. 1391–1403, 2020\.
  * [7] P. Chatterjee, A. Pal, J. S. Thorp, and J. De La Ree, “Partitioned linear state estimation,” in _2015 IEEE Power & Energy Society Innovative Smart Grid Technologies Conference (ISGT)_, 2015, pp. 1–5.
  * [8] A. C. Varghese, A. Pal, and G. Dasarathy, “Transmission line parameter estimation under non-Gaussian measurement noise,” _IEEE Transactions on Power Systems_ , pp. 1–16, 2022.
  * [9] S. Wang, J. Zhao, Z. Huang, and R. Diao, “Assessing Gaussian assumption of PMU measurement error using field data,” _IEEE Transactions on Power Delivery_ , vol. 33, no. 6, pp. 3233–3236, 2018.
  * [10] A. Pal, G. A. Sanchez-Ayala, V. A. Centeno, and J. S. Thorp, “A PMU placement scheme ensuring real-time monitoring of critical buses of the network,” _IEEE Transactions on Power Delivery_ , vol. 29, no. 2, pp. 510–517, 2013.
  * [11] A. Pal, A. K. S. Vullikanti, and S. S. Ravi, “A PMU placement scheme considering realistic costs and modern trends in relaying,” _IEEE Transactions on Power Systems_ , vol. 32, no. 1, pp. 552–561, 2017.
  * [12] J. E. Anderson and A. Chakrabortty, “A minimum cover algorithm for PMU placement in power system networks under line observability constraints,” in _2012 IEEE Power and Energy Society General Meeting_ , 2012, pp. 1–7.
  * [13] M. Khodayar, G. Liu, J. Wang, and M. E. Khodayar, “Deep learning in power systems research: A review,” _CSEE Journal of Power and Energy Systems_ , vol. 7, no. 2, pp. 209–220, 2021.
  * [14] A. Abdiansah and R. Wardoyo, “Time complexity analysis of support vector machines (SVM) in LibSVM,” _International Journal Computer and Application_ , vol. 128, no. 3, pp. 28–34, 2015.
  * [15] B. Azimian, R. S. Biswas, A. Pal, and L. Tong, “Time synchronized state estimation for incompletely observed distribution systems using deep learning considering realistic measurement noise,” in _2021 IEEE Power & Energy Society General Meeting_. IEEE, 2021, pp. 1–5.
  * [16] B. Azimian, R. S. Biswas, S. Moshtagh, A. Pal, L. Tong, and G. Dasarathy, “State and topology estimation for unobservable distribution systems using deep neural networks,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 71, pp. 1–14, 2022.
  * [17] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “MATPOWER: Steady-state operations, planning, and analysis tools for power systems research and education,” _IEEE Transactions on Power Systems_ , vol. 26, no. 1, pp. 12–19, 2010.
  * [18] A. B. Birchfield, T. Xu, K. M. Gegner, K. S. Shetye, and T. J. Overbye, “Grid structural characteristics as validation criteria for synthetic networks,” _IEEE Transactions on Power Systems_ , vol. 32, no. 4, pp. 3258–3265, 2017\.
  * [19] _IEEE/IEC International Standard - Measuring relays and protection equipment - Part 118-1: Synchrophasor for power systems - Measurements_ , 2018.
  * [20] R. S. Biswas, A. Pal, T. Werho, and V. Vittal, “Mitigation of saturated cut-sets during multiple outages to enhance power system security,” _IEEE Transactions on Power Systems_ , vol. 36, no. 6, pp. 5734–5745, 2021\.


††thanks: Corresponding author.

# Observation of local atomic displacements intrinsic to the double zigzag
chain structure of 1$T$-$M$Te2 ($M$ = V, Nb, Ta).

N. Katayama<EMAIL_ADDRESS>Department of Applied
Physics, Nagoya University, Aichi 464-8603, Japan  Y. Matsuda  Department of
Applied Physics, Nagoya University, Aichi 464-8603, Japan  K. Kojima
Department of Applied Physics, Nagoya University, Aichi 464-8603, Japan  T.
Hara  Department of Applied Physics, Nagoya University, Aichi 464-8603, Japan
S. Kitou  RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama,
351-0198, Japan  Department of Advanced Materials Science, The University of
Tokyo, Kashiwa, 277-8561, Chiba, Japan.  N. Mitsuishi  RIKEN Center for
Emergent Matter Science (CEMS), Wako, Saitama, 351-0198, Japan  H. Takahashi
Division of Materials Physics and Center for Spintronics Research Network
(CSRN), Graduate School of Engineering Science, Osaka University, Toyonaka,
Osaka 560-8531, Japan  Spintronics Research Network Division, Institute for
Open and Transdisciplinary Research Initiatives, Osaka University, Suita,
Osaka, 565-0871, Japan  S. Ishiwata  Division of Materials Physics and Center
for Spintronics Research Network (CSRN), Graduate School of Engineering
Science, Osaka University, Toyonaka, Osaka 560-8531, Japan  Spintronics
Research Network Division, Institute for Open and Transdisciplinary Research
Initiatives, Osaka University, Suita, Osaka, 565-0871, Japan  K. Ishizaka
RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama, 351-0198,
Japan  Quantum-Phase Electronics Center and Department of Applied Physics, The
University of Tokyo, Bunkyo-ku, Tokyo, 113–8656, Japan  H. Sawa  Department of
Applied Physics, Nagoya University, Aichi 464-8603, Japan

###### Abstract

We describe the existence of local distortion discovered in the synchrotron
X-ray single crystal structure analysis of layered ditelluride 1$T$-$M$Te2
($M$ = V, Nb, Ta). In 1$T$-TaTe2, the double zigzag chain structure of Ta is
deformed at about 170 K, and heptamer molecules are formed periodically at low
temperatures. We found that some of the Ta atoms that compose the double
zigzag chain structure appearing at high temperatures are locally displaced,
resulting in local dimerization. This tendency weakens when Ta is replaced by
V or Nb. Our results indicate that the local distortion persistently survives
in these ditellurides, where the electronic degrees of freedom, including
orbitals, are weakened. We further discuss the origin of local distortion in
these ditellurides, which is different from many usual material systems where
molecular formation occurs at low temperatures.

††preprint: APS/123-QED

##  I Introduction

Among transition metal compounds with orbital degrees of freedom, there are
many substances whose atoms assemble to form “molecules” at low temperatures
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Examples include dimers in
LiRh2O4 [1] and CuIr2S4 [2], trimers in LiVO2 [3, 4] and LiVS2 [5], and
heptamer (trimer/tetramer pair) in AlV2O4 [6, 14]. Although these molecules
that form spontaneously in crystals have been thought to disappear and regular
lattices are realized at high temperatures, recent local structure studies
have revealed that local lattice distortions appear in various forms in a
precursory manner [14, 15, 16, 17, 18, 19, 20]. For example, in Li2RuO3 [15]
and LiRh2O4 [16], dimers that appear at low temperatures appear as short-range
orders at high temperatures; in CuIr2S4 [18], tetragonal distortions appear
locally at high temperatures; in LiVS2, short-range orders of zigzag chains
that are unrelated to trimer in the low temperature phase appear and slowly
fluctuate at high temperatures [19]. These can be interpreted as local nematic
states in which the spontaneous symmetry lowering of the electronic system is
strongly coupled to the lattice system [18], and the search for various
distortion patterns and the elucidation of their mechanisms are important
research themes that go beyond the category of molecular formation systems and
have a broad impact on physical properties in general.

Layered transition metal ditellurides $M$Te2 provide a unique playground for
such studies. In these material systems, Te-Te covalent bonds derived from
large tellurium ions occur, which cause the formal valence of Te to shift from
2- and transfer additional electrons to the transition metal element $M$.
Depending on the amount of charge transfer, a variety of molecular formation
patterns coupled with charge degrees of freedom appear at low temperatures.
For example, in IrTe2, a charge-ordered stripe state between Ir3+ and Ir4+
occurs at low temperatures, forming an Ir4+-Ir4+ dimer state [21, 22], and
superconductivity appears when this dimer phase is suppressed by Pt doping
[23]. In 1$T$-$M$Te2 ($M$ = V, Nb, Ta), transition metal elements form quasi-
one-dimensional double zigzag chains, named “ribbon chain”, as shown in Fig.
1(a) [24, 25, 26]. This chain is formed from multiple linear trimers, with
each $M$ element offering 2/3 electrons for the formation of one trimer [26,
27]. An interesting feature of this system is that only at 1$T$-TaTe2 the
modulation of the charge changes at low temperatures below $T_{c}$ $\sim$ 170
K and the double zigzag chain changes to Ta heptamers, as shown in Fig. 1(b)
[28]. This transformation is different from the situation in conventional
matter systems where molecular formation occurs from a regular lattice, where
the electron degrees of freedom are expected to be highly degenerate. Do
unique local distortions appear in such ditellurides at high temperatures? In
addition, clarifying whether the nature of local distortion differs between
1$T$-TaTe2, where heptamerization occurs, and 1$T$-NbTe2 and 1$T$-VTe2, where
the double zigzag chain is maintained at low temperatures, will provide
important insights into the background physics that generates local
distortion.

In this article, we report on the structural analysis of 1$T$-$M$Te2 ($M$ = V,
Nb, Ta) single crystals using synchrotron radiation X-rays. The anisotropy of
the atomic displacement parameters (ADPs) obtained from the structural
analysis suggests that the $M$ atom at the center of the double zigzag chain
of these diterulides is locally distorted, in the zigzag chain direction.
Structural analysis using the split-site model shows that local distortion
toward the heptamerization occurs, and this tendency is strongest at
1$T$-TaTe2. This indicates that local distortions appear universally, as in
many molecule-forming systems, even if the average structure of the high-
temperature phase is in a state where the degeneracy of the orbital degrees of
freedom has been resolved.

##  II Results and discussions

###  II.1 Sample Preparation and experimental details

Single crystal samples of 1$T$-$M$Te2 ($M$ = V, Nb, Ta) were synthesized by a
conventional solid phase reaction method. The mixture of the constituent
elements in their amphoteric ratios was vacuum-sealed and sintered at 1000
$\tccentigrade$ for $M$ = Nb and Ta, and at 850 $\tccentigrade$ for $M$ = V
for 15 hours. The obtained samples are basically powders, but some of them
contain tiny single crystals of several tens of micrometers on a side, which
were used for diffraction experiments using synchrotron radiation x-rays.
X-ray diffraction (XRD) experiments were conducted using the BL02B1 beamline
at SPring-8 at an x-ray energy of 40 keV. The typical dimensions of the
1$T$-$M$Te2 ($M$ = V, Nb, Ta) single crystal used for the XRD experiment were
20 × 20 × 20 $\mu$m3. A He-gas-blowing device was employed to cool the sample
to 100 K. A 2D CdTe PILATUS detector was utilized with the diffractometer. The
CrysAlisPro program was used to integrate the diffraction profiles.
Diffraction intensity averaging and refinement of structural parameters were
performed using Jana2006 program [29]. Crystal structure was visualized by
using VESTA [30]. The obtained powder diffraction data were indexed using
Conograph [31], and the analysis was performed using Rietan-FP [32].

Figure 1:  (a-b) The lattice structure of Ta ions at (a) 300 K and (b) 125 K.
(c-d) Single crystal XRD patterns of 1$T$-TaTe2 at (c) 300 K and (d) 125 K.
Both the high and low temperature phases could be refined using the previously
reported structure. Details are shown in the Appendix of this paper. (e)
Geometrical arrangement of merohedral domains sharing the 20$\bar{1}$ plane in
the crystal structure of the high-temperature phase monoclinic $C$2/$m$.

###  II.2 Structural studies of 1$T$-TaTe2

Fig. 1(c) and (d) show single-crystal XRD patterns of 1$T$-TaTe2 in the high
temperature (300 K) and low temperature (125 K) phases. Although the Bragg
peaks remain sharp, the presence of a merohedral domain sharing a 20$\bar{1}$
plane, as shown in Fig. 1(e), is confirmed. In the following, we show the
results of the analysis in which only one domain component is extracted. The
diffraction pattern changes below the phase transition around 170 K, and the
superlattice peak shown by the circle is observed at low temperatures.

The obtained crystal structures are shown in Figs. 1(a) and (b), both of which
are consistent with the previously reported structures [25, 28, 33]. Details
of the structure analysis results are summarized in the Appendix of this
paper. It is noteworthy that the Ta-Ta distance constituting the linear trimer
changes significantly between the high and low temperature phases. In the
high-temperature phase, the linear trimer is composed of equally spaced
Ta2-Ta1-Ta2 arrays, whereas in the low-temperature phase, the displacement of
Ta1b ions is accompanied by a large difference in the distance between
adjacent Ta-Ta inside the Ta2a-Ta1b-Ta2b array. This indicates that the
Ta2a-Ta1b-Ta2b array does not form a linear trimer in the low-temperature
phase, but rather transforms into a Ta2a-Ta1b dimer and isolated Ta2b ion.
Although the contraction of the Ta-Ta distance upon heptamerization is
greatest between Ta1a and Ta1b, the Ta1a-Ta1b distance of $\sim$3.34 Å in the
low temperature phase is still much longer than the Ta2a-Ta1b distance of
$\sim$3.18 Å, indicating that the Ta2a-Ta1b bond is more essential. The large
change in the Ta1a-Ta1b distance associated with the phase transition should
be a side effect of the $b$ axis displacement of Ta1b to form dimers with two
Ta2a simultaneously.

Structural analysis revealed that the atomic positions of the high-temperature
phase are almost the same as those previously reported [25], and the
anisotropic ADP at the Ta1 site also shows an anomalous elongation in the
double zigzag chain direction as shown in the inset of Fig. 2(a), as indicated
in a previous report [28]. Fig. 2(a) shows that the $U_{22}$ parameter has an
unusually large value compared to $U_{11}$ and $U_{33}$. It is important to
note that the $U_{11}$ and $U_{33}$ parameters are zero when extrapolated
toward 0 K, but the $U_{22}$ parameter clearly reaches a finite value. This
indicates that the anomalous increase in $U_{22}$ is not simply anisotropic
strong thermal oscillation, but rather a local distortion at the Ta1 site.
This is consistent with the possibility of dynamic disorder discussed in the
previous report [28].

Figure 2:  (a) Temperature dependence of the atomic displacement parameters
$U_{11}$, $U_{22}$, and $U_{33}$ of the Ta1 site in the high temperature
phase. The equations relating the temperature factor T to these atomic
displacement parameters are given in the Appendix of this paper. $U_{11}$,
$U_{22}$, $U_{33}$ are the mean square amplitudes $\langle$$u^{2}$$\rangle$ in
the reciprocal lattice vector $a^{*}$, $b^{*}$, $c^{*}$ directions. The inset
shows the thermal oscillating ellipsoid showing 99 % probability of the Ta1
site. (b) Temperature dependence of $R$ values at each $r$ obtained by the
split-site model. The inset shows a schematic picture of the split-site model.
(c) Temperature dependence of Ta-Ta distance for average structure and split-
site model, respectively.

In order to clarify the local displacement of the ions at the Ta1 site, we
performed a structural analysis using a split-site model. This is an
analytical method that examines the change in the confidence factor $R$ with
the change in distance $r$, assuming that the Ta1 site ions are not at the
central position, but two ions with occupancy 0.5 exist at positions +$r$ and
-$r$ along the $b$ axis, as shown in the inset of Fig. 2(b). If the atomic
displacement as represented by this model does not actually occur, the $R$
value will show a minimum at $r$ = 0. If it occurs, the $R$ value will be a
minimum at a finite $r$. As shown in Fig. 2(b), the analysis shows that the
$R$ value varies significantly with $r$, reaching a minimum at about $r$ =
0.13 Å. This indicates that there is an intrinsic local atomic displacement in
the high-temperature phase of 1$T$-TaTe2, where the local atomic displacement
$r$ value that minimizes the $R$ value is almost independent of temperature,
and consequently the Ta1-Ta2 bond that splits into two types, as shown in Fig.
2(c). It seems strange that such local distortion is independent of
temperature, but it should be noted that similar behavior has been observed at
high temperatures in AlV2O4 [14] and Li2RuO3 [15] where molecular formation
occurs at low temperature.

Figure 3:  (a) Temperature dependence of the $U_{11}$,$U_{22}$ and $U_{33}$
parameters of 1$T$-NbTe2 and 1$T$-VTe2. Inset shows a schematic picture of the
thermal oscillating ellipsoid of 1$T$-NbTe2 and 1$T$-VTe2 compared to that of
1$T$-TaTe2. (b) Temperature dependence of $R$ values at each $r$ obtained by
the split-site model. The arrows in the figure indicate the minimum $R$ value.
(c) Temperature dependence of local distortion $r$ for 1$T$-NbTe2 and
1$T$-VTe2. The $r$ range from the smallest $R$ value to a 0.5 % increase was
defined as the error bar.

###  II.3 Structural studies of 1$T$-NbTe2 and 1$T$-VTe2

Will similar local atomic displacements appear in 1$T$-NbTe2 and 1$T$-VTe2,
where the double zigzag chain structure is maintained down to low
temperatures? The anisotropic ADPs of the high-temperature phases of these two
materials are shown in the inset of Fig. 3(a). Although a larger $U_{22}$ is
realized compared to $U_{11}$ and $U_{33}$, its value is smaller than that of
1$T$-TaTe2. As summarized in Appendix, the values of lattice constants and
bond lengths are close between 1$T$-TaTe2, 1$T$-NbTe2 and 1$T$-VTe2. Thus, the
difference in magnitude of $U_{22}$ indicates a clear difference in local
structure between these telluriums. Also, as shown in Fig. 3(a), the $U_{22}$
parameter decreases in a temperature-dependent manner to zero at 0 K. Fig.
3(b) shows the results of the split-site model analysis for these two
materials: 1$T$-NbTe2 exhibits a minimum $R$ value at finite $r$, but the
value of $r$ that minimizes $R$ at 300 K is about 0.09 Å, which is smaller
than that of 1$T$-TaTe2. The change of $R$ with $r$ is also very small
compared to 1$T$-TaTe2. This trend is more pronounced for 1$T$-VTe2, where the
$R$ value is almost constant over a wide $r$ range. As shown in Fig. 3(c), the
local atomic displacements, if any, of 1$T$-NbTe2 and 1$T$-VTe2 are much
smaller than those of 1$T$-TaTe2 and decrease upon cooling to zero at 0 K.

###  II.4 Discussion based on the single-crystal X-ray diffraction
experimental results

Figure 4:  (a) Energy schemes in the linear trimer state and in the dimer and
isolated-atom states. (b) Schematic diagram of the elemental $M$ dependence of
the relationship between the bonding orbital energy difference $\Delta$ and
the lattice system energy loss $\varepsilon$ in the linear trimer state and
the dimer with isolated atom states, and the temperature dependence of the
local distortion $r$ caused by the difference in $M$.

These experimental results raise two questions. The first question is what is
the origin of these local atomic displacements. In the high-temperature phase
of 1$T$-TaTe2, we can predict that two dimer patterns with equal lattice
energy can be realized, as shown in the right side of Fig. 4(a), which are
thermally fluctuating, resulting in a linear trimer in the average structure,
as shown in the left side of Fig. 4(a). In the three-center, two-electron
bonding state derived from the linear trimer, three orbitals are formed: a
bonding orbital, a non-bonding orbital, and an anti-bonding orbital, with two
electrons stored in the bonding orbitals, which are energetically stabilized,
as shown in the left energy scheme of Fig. 4(a). When local atomic
displacement occurs, the linear trimer is transformed into a dimer-isolated
atom pair, which should result in the formation of a bonding orbital and an
anti-bonding orbital, derived from the dimer, and an isolated atomic orbital,
as shown in the right energy scheme of Fig. 4(a). The three energetically
equally spaced orbitals formed by the linear trimer and dimer-isolated atom
pairs are similar at first glance. However, if the actual interatomic
distances of the dimer are sufficiently closer than those of the linear
trimer, the bonding orbitals should be $\Delta$ lower in energy for the dimer.
If the energy difference $\Delta$ is larger than the energy loss due to
lattice distortion, $\varepsilon$, the local dimer is expected to stabilize.
Since the energies between the two dimer patterns are equal, we expect the two
dimer patterns to appear thermally fluctuating, as shown in Fig. 4(a). The
transition process between the two dimer patterns is via a trimer state, and
the low energy of the bonding orbital in this trimer state should facilitate
the transition between the two patterns. The long-range ordering of this local
distortion at low temperatures leads to the formation of heptamer clusters, as
shown in Fig. 1(b).

Another question is why different behaviors appear for 1$T$-TaTe2, 1$T$-NbTe2,
and 1$T$-VTe2. As shown in Fig. 3(c), local atomic displacements are largest
for Ta(5$d$) and small for Nb(4$d$) and V(3$d$). In 1$T$-TaTe2 with large 5$d$
orbitals, the electronic energy gain $\Delta$ is expected to be the largest
because the orbitals overlap more during dimerization than in 1$T$-NbTe2 and
1$T$-VTe2. If $\varepsilon$ $<$ $\Delta$ is realized at 1$T$-TaTe2, the
overall energy of the system including lattice energy with respect to atomic
positions will be as shown in Fig. 4(b). Since atomic displacement and the
associated local dimerization stabilize the system, and this trend does not
change even when the temperature is lowered, the atomic displacement $r$ does
not show any temperature dependence, and the local distortion is expected to
survive until just above the phase transition. On the other hand, for 4$d$ Nb
and 3$d$ V, the orbital overlap is small and the energy is not very stable due
to local dimerization. Therefore, $\varepsilon$ $\geq$ $\Delta$ is realized
and the energy of the whole system including lattice energy has a flat shape
with respect to $r$. The energy drops to a low central position when the
temperature is lowered. Because of these differences in energy schemes, it is
thought that among the three tellurides, only 1$T$-TaTe2 produces large local
distortions, which are maintained at low temperatures and develop into
heptamers. From the above discussion, it seems likely that in 1$T$-NbTe2 and
1$T$-VTe2, when the distance between adjacent transition metals is closer, the
orbital overlap becomes larger and dimer (heptamer) states such as 1$T$-TaTe2
are formed. It should be noted that it has been reported that applying
pressure to 1$T$-NbTe2 induces a structural evolution from the trimeric to the
dimeric structure [34], confirming the validity of this argument.

An argument similar to the local distortion of 1$T$-TaTe2 might be applied to
other material systems that exhibit molecular formation at low temperatures,
such as the spinel compound AlV2O4. It has been argued that vanadium in AlV2O4
spontaneously forms a heptamer molecule at about 700 K. Based on the Curie
paramagnetic component in the low-temperature phase by magnetization
measurements, the heptamer has been revealed to be composed of 9 bonds by 18
electrons [6]. However, recent structural analysis of the split-site model on
synchrotron XRD data has revealed that the heptamer is actually composed of
trimer/tetramer pairs [14]. The two structures are compared in Fig. 25 of the
paper by D.I. Khomskii and S.V. Streltsov [35]. The nature of the difference
between the two types of molecular structures proposed for AlV2O4 can be
interpreted as whether the three bonds connecting the upper and lower
triangular trimer are linear trimers consisting of three-center two-electron
bonds or pairs of dimers and isolated atoms. This is very similar to the
present case. An interesting difference is that such local distortions appear
in the low-temperature ordered phase in AlV2O4, whereas they appear only in
the high-temperature phase in 1$T$-TaTe2. This indicates that unconventional
orbital ordering states that appear in the low-temperature phase can
universally appear as local distortions in the high-temperature phase, even in
systems like 1$T$-TaTe2, where there is no short-range order in the low-
temperature ordered phase.

The presence of local distortion in the high-temperature phase not only has
important implications for the mechanism of phase transitions, but is also
important for understanding the thermodynamics of the high-temperature phase.
For example, while high thermal conductivity is generally realized in the
metallic state due to itinerant electrons, it has been observed that in
CuIr2S4, which undergoes insulating transition with molecular formation at low
temperatures, the thermal conductivity of the high temperature metallic phase
is lower than that of the insulating low temperature phase. It is argued that
local distortion suppresses phonon thermal conduction [36, 37, 38]. In
addition, the multiple degrees of freedom of electrons are expected to play an
important role in the entropy changes associated with molecular formation [8,
5]. The presence of local distortion will be a new factor, not previously
considered, to the concept of electron degrees of freedom in the high-
temperature phase. In this respect, the finding of differences in local
structure among the systematics of $M$ = V, Nb, and Ta in this study should
correspond to providing an attractive stage for discussing the role of local
distortion on physical properties through comparison.

###### Acknowledgements.

The work leading to these results has received funding from the Grant in Aid
for Scientific Research (Nos. JP20H02604, JP21K18599, JP21J21236) and Research
Foundation for the Electrotechnology of Chubu. This work was carried out under
the Visiting Researcher’s Program of the Institute for Solid State Physics,
the University of Tokyo, and the Collaborative Research Projects of Laboratory
for Materials and Structures, Institute of Innovative Research, Tokyo
Institute of Technology. Single-crystal and powder XRD experiments were
conducted at the BL02B1 and BL02B2 of SPring-8, Hyogo, Japan (Proposals No.
2019B1085, 2021A0070, 2021B1261, 2021B1136, 2021B1261, 2022A0304, 2022B0607,
2022B1570, 2022B1582, and 2022B1862), and at the BL5S2 of Aichi Synchrotron
Radiation Center, Aichi Science and Technology Foundation, Aichi, Japan
(Proposals No. 202202037, 202201033, 202105170, 202104111, 2021L3002,
2021L2002, 2021L1002 and 2020L4002).

##  Appendix A Single crystal X-ray diffraction analysis results

Here, the temperature factor $T$ is expressed as a function of the atomic
displacement parameters $U_{11}$, $U_{22}$, $U_{33}$, $U_{12}$, $U_{13}$ and
$U_{23}$, using the following equation,

$\begin{split}T&=\exp{}\\{-2\pi^{2}(h^{2}a^{*2}U_{11}+k^{2}b^{*2}U_{22}+l^{2}c^{*2}U_{33}\\\
&\quad+2hka^{*}b^{*}U_{12}+2hla^{*}c^{*}U_{13}+2klb^{*}c^{*}U_{23})\\}.\end{split}$

$U_{11}$, $U_{22}$ and $U_{33}$ are the mean square amplitudes
$\langle$$u^{2}$$\rangle$ in the reciplocal lattice vector $a^{*}$, $b^{*}$
and $c^{*}$ directions.

Table 1: Summary of crystallographic data of 1$T$-TaTe2.

Temperature (K) | 110 (low $T$ phase) | 175 (high $T$ phase)
---|---|---
Wavelength (Å) | 0.31007 | 0.31007
Crystal dimension ($\mu$m3) | 20$\times$20$\times$20 | 20$\times$20$\times$20
space group |  $C$2/$m$ |  $C$2/$m$
$a$ (Å) | 14.7669(3) | 14.7408(4)
$b$ (Å) | 10.8702(2) | 3.62940(10)
$c$ (Å) | 9.2926(2) | 9.3287(2)
$\beta$ (∘) | 110.630(8) | 110.864(8)
$V$ (Å3) | 1395.99(5) | 466.36(3)
$Z$ | 18 | 6
$F$(000) | 3186 | 1062
(sin$\theta$/$\lambda$)Max (Å-1) | 1.3513 | 1.3512
$N_{Total,obs}$ | 92313 | 32489
$N_{Unique,obs}$ | 13372 | 4637
Average redundancy | 6.9 | 7.0
Completeness | 0.904 | 0.894
Structural analysis using anisotropic displacement parameters
$R_{1}$ [# of reflections] | 0.066 [12475] | 0.0467 [4259]
$R_{1}$ ($I$ $>$ 4$\sigma$) [# of reflections] | 0.0445 [9365] | 0.0358 [3586]
GOF [# of reflections] | 1.035 [12475] | 1.057 [4259]
Structural analysis without using anisotropic displacement parameters∗
split-site model | Not used | Use | Not used
$R_{1}$ [# of reflections] | 0.0673 [12475] | 0.0468 [4259] | 0.1318 [4259]
$R_{1}$ ($I$ $>$ 4$\sigma$) [# of reflections] | 0.0456 [9365] | 0.0359 [3586] | 0.1152 [3586]
GOF [# of reflections] | 1.037 [12475] | 1.048 [4259] | 1.093 [4259]

* In this analysis, isotropic ADP is used only for the Ta1 site at high temperatures and its derived sites at low temperatures,

and anisotropic ADP is used for the remaining sites. In this table, the Z
value indicates how many molecules of that

molecular composition unit are contained in the crystal lattice, and $F$(000)
is the total number of electrons in the unit cell.

As noted in the text, there are merohedral twinning in the present telluride
that shares certain planes, and some

peaks are composed of contributions from both domains and cannot be separated.
Such peaks were excluded from the

analysis. The resulting completeness, which is the ratio of the number of
experimentally observed Bragg peaks to

the number of Bragg peaks that can be observed in principle, is about 0.9; the
same occurs for Nb and V derivatives.

Table 2: Structural parameters of 1$T$-TaTe2 at 175 K without the split-site
model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$ (Å2)
Ta1 | 2$d$ | 1 | 1/2 | -1/2 | 0 | 0.01177(4)
Ta2 | 4$i$ | 1 | 0.63944(2) | 0 | 0.29083(2) | 0.00437(2)
Te1 | 4$i$ | 1 | 0.79655(2) | 1/2 | 0.37790(2) | 0.00419(3)
Te2 | 4$i$ | 1 | 0.49515(2) | -1/2 | -0.30921(2) | 0.00431(3)
Te3 | 4$i$ | 1 | 0.64883(2) | 0 | 0.01094(2) | 0.00551(3)

Table 3: Anisotropic atomic displacement parameters of 1$T$-TaTe2 at 175 K
without the split-site model.

site |  $U_{11}$ (Å2) |  $U_{22}$ (Å2) |  $U_{33}$ (Å2) |  $U_{12}$ (Å2) |  $U_{13}$ (Å2) |  $U_{23}$ (Å2)
---|---|---|---|---|---|---
Ta1 | 0.00341(6) | 0.02829(13) | 0.00349(7) | 0 | 0.00109(5) | 0
Ta2 | 0.00357(3) | 0.00539(4) | 0.00348(4) | 0 | 0.00044(3) | 0
Te1 | 0.00394(5) | 0.00370(5) | 0.00420(6) | 0 | 0.00054(4) | 0
Te2 | 0.00440(5) | 0.00427(5) | 0.00416(6) | 0 | 0.00139(5) | 0
Te3 | 0.00711(6) | 0.00444(5) | 0.00429(6) | 0 | 0.00119(5) | 0

Table 4: Structural parameters of 1$T$-TaTe2 at 175 K with the split-site
model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$($U_{iso}$) (Å2)
Ta1 | 4$g$ | 0.5 | 1/2 | -0.53555(7) | 0 | 0.00430(3)
Ta2 | 4$i$ | 1 | 0.63944(2) | 0 | 0.29083(2) | 0.00441(2)
Te1 | 4$i$ | 1 | 0.79655(2) | 1/2 | 0.37788(2) | 0.00425(3)
Te2 | 4$i$ | 1 | 0.49516(2) | -1/2 | -0.30919(2) | 0.00432(2)
Te3 | 4$i$ | 1 | 0.64882(2) | 0 | 0.01091(2) | 0.00551(3)

In the split-site model analysis, isotropic ADP is used only for the Ta1 site,
and anisotropic ADP is used for the other

elements. The fact that $U_{iso}$ at the Ta1 site and $U_{eq}$ at other sites
are almost equal indicates the validity of this analysis.

Table 5: Structural parameters of 1$T$-TaTe2 at 110 K without the split-site
model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$ (Å2)
Ta1a | 2$d$ | 1 | 0 | 1/2 | 1/2 | 0.00234(3)
Ta1b | 4$h$ | 1 | 0 | 0.19316(2) | 1/2 | 0.00205(2)
Ta2a | 8$j$ | 1 | 0.85918(2) | 0.33810(2) | 0.21004(2) | 0.00196(2)
Ta2b | 4$i$ | 1 | 0.86247(2) | 0 | 0.20874(2) | 0.00229(2)
Te1a | 8$j$ | 1 | 0.70510(2) | -0.16609(2) | 0.12266(2) | 0.00228(2)
Te1b | 4$i$ | 1 | 0.79909(2) | 0 | -0.12305(3) | 0.00223(3)
Te2a | 4$i$ | 1 | 0.99543(2) | 1/2 | 0.18923(3) | 0.00216(3)
Te2b | 8$j$ | 1 | 0.99419(2) | 0.17096(2) | 0.19070(2) | 0.00226(2)
Te3a | 4$i$ | 1 | 0.86348(2) | 0 | 0.49175(3) | 0.00207(3)
Te3b | 8$j$ | 1 | 0.84603(2) | 0.33834(2) | 0.49008(2) | 0.00215(2)

Table 6: Anisotropic atomic displacement parameters of 1$T$-TaTe2 at 110 K
without the split-site model.

site |  $U_{11}$ (Å2) |  $U_{22}$ (Å2) |  $U_{33}$ (Å2) |  $U_{12}$ (Å2) |  $U_{13}$ (Å2) |  $U_{23}$ (Å2)
---|---|---|---|---|---|---
Ta1a | 0.00143(5) | 0.00343(6) | 0.00215(7) | 0 | 0.00064(5) | 0
Ta1b | 0.00172(4) | 0.00223(4) | 0.00215(5) | 0 | 0.00063(4) | 0
Ta2a | 0.00158(3) | 0.00201(3) | 0.00199(3) | -0.00004(2) | 0.00024(2) | -0.00005(2)
Ta2b | 0.00187(4) | 0.00266(4) | 0.00210(5) | 0 | 0.00040(4) | 0
Te1a | 0.00187(4) | 0.00207(4) | 0.00247(5) | 0.00004(3) | 0.00023(4) | 0.00006(3)
Te1b | 0.00169(6) | 0.00216(6) | 0.00242(7) | 0 | 0.00021(5) | 0
Te2a | 0.00208(6) | 0.00191(5) | 0.00251(7) | 0 | 0.00083(5) | 0
Te2b | 0.00226(4) | 0.00199(4) | 0.00247(5) | -0.00008(3) | 0.00077(4) | 0.00001(3)
Te3a | 0.00197(6) | 0.00182(5) | 0.00246(7) | 0 | 0.00084(5) | 0
Te3b | 0.00186(4) | 0.00195(4) | 0.00258(5) | 0.00001(3) | 0.00072(4) | 0.00001(3)

Ta1a and Ta1b sites are from the high temperature phase Ta1 site, Ta2a and
Ta2b sites are from the high temperature

phase Ta2 site, Te1a and Te1b sites are from the high temperature phase Te1
site, Te2a and Te2b are from the high

temperature phase Te2 site, Te3a and Te3b are from the high temperature phase
Te3 site of the high-temperature phase.

Table 7: Summary of crystallographic data of 1$T$-NbTe2 at 300 K.

Temperature (K) | 300
---|---
Wavelength (Å) | 0.31011
Crystal dimension ($\mu$m3) | 20$\times$20$\times$10
space group |  $C$2/$m$
$a$ (Å) | 14.6619(5)
$b$ (Å) | 3.63760(10)
$c$ (Å) | 9.3144(3)
$\beta$ (∘) | 110.070(8)
$V$ (Å3) | 466.61(3)
$Z$ | 6
$F$(000) | 870
(sin$\theta$/$\lambda$)Max (Å-1) | 1.2499
$N_{Total,obs}$ | 29479
$N_{Unique,obs}$ | 3618
Average redundancy | 8.1
Completeness | 0.878
Structural analysis using anisotropic displacement parameters
$R_{1}$ [# of reflections] | 0.0434 [3494]
$R_{1}$ ($I$ $>$ 4$\sigma$) [# of reflections] | 0.0290 [2778]
GOF [# of reflections] | 0.981 [3494]
Structural analysis without using anisotropic displacement parameters∗
split-site model | Use | Not used
$R_{1}$ [# of reflections] |  0.0441 [3494] | 0.0564 [3494]
$R_{1}$ ($I$ $>$ 4$\sigma$) [# of reflections] | 0.0294 [2778] | 0.0407 [2778]
GOF [# of reflections] | 0.998 [3494] | 1.061 [3494]

* In this analysis, isotropic ADP is used only for the Nb1 site, and anisotropic ADP is used for the remaining sites.

Table 8: Summary of crystallographic data of 1$T$-NbTe2 at 300 K without the
split-site model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$ (Å2)
Nb1 | 2$d$ | 1 | 1/2 | -1/2 | 0 | 0.01070(4)
Nb2 | 4$i$ | 1 | 0.638910(14) | 0 | 0.29045(2) | 0.00829(3)
Te1 | 4$i$ | 1 | 0.797228(11) | 1/2 | 0.378315(18) | 0.00888(3)
Te2 | 4$i$ | 1 | 0.503383(11) | -1/2 | 0.309539(18) | 0.00867(3)
Te3 | 4$i$ | 1 | 0.649411(11) | 0 | 0.009461(18) | 0.00865(3)

Table 9: Anisotropic atomic displacement parameters of 1$T$-NbTe2 at 300 K
without the split-site model.

site |  $U_{11}$ (Å2) |  $U_{22}$ (Å2) |  $U_{33}$ (Å2) |  $U_{12}$ (Å2) |  $U_{13}$ (Å2) |  $U_{23}$ (Å2)
---|---|---|---|---|---|---
Nb1 | 0.00864(8) | 0.01617(12) | 0.00719(9) | 0 | 0.00259(7) | 0
Nb2 | 0.00853(6) | 0.00850(7) | 0.00692(6) | 0 | 0.00147(5) | 0
Te1 | 0.00891(5) | 0.00769(5) | 0.00863(5) | 0 | 0.00119(4) | 0
Te2 | 0.00998(5) | 0.00792(5) | 0.00824(5) | 0 | 0.00327(4) | 0
Te3 | 0.00971(5) | 0.00735(5) | 0.00864(5) | 0 | 0.00284(4) | 0

Table 10: Structural parameters of 1$T$-NbTe2 at 300 K with the split-site
model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$($U_{iso}$) (Å2)
Nb1 | 4$g$ | 0.5 | 1/2 | -0.52321(19) | 0 | 0.00812(5)
Nb2 | 4$i$ | 1 | 0.63892(2) | 0 | 0.29044(2) | 0.00830(3)
Te1 | 4$i$ | 1 | 0.79723(2) | 1/2 | 0.37831(2) | 0.00889(3)
Te2 | 4$i$ | 1 | 0.50338(2) | -1/2 | 0.30954(2) | 0.00868(3)
Te3 | 4$i$ | 1 | 0.64941(2) | 0 | 0.00948(2) | 0.00865(3)

In the split-site model analysis, isotropic ADP is used only for the Nb1 site,
and anisotropic ADP is used for the other

elements. The fact that $U_{iso}$ at the Nb1 site and $U_{eq}$ at other sites
are almost equal indicates the validity of this analysis.

Table 11: Summary of crystallographic data of 1$T$-NbTe2 at 100 K.

Temperature (K) | 100
---|---
Wavelength (Å) | 0.31011
Crystal dimension ($\mu$m3) | 20$\times$20$\times$10
space group |  $C$2/$m$
$a$ (Å) | 14.5770(3)
$b$ (Å) | 3.63410(10)
$c$ (Å) | 9.2961(2)
$\beta$ (∘) | 109.956(8)
$V$ (Å3) | 462.88(3)
$Z$ | 6
$F$(000) | 870
(sin$\theta$/$\lambda$)Max (Å-1) | 1.2500
$N_{Total,obs}$ | 28124
$N_{Unique,obs}$ | 3608
Average redundancy | 7.8
Completeness | 0.882
Structural analysis using anisotropic displacement parameters
$R_{1}$ [# of reflections] | 0.0453 [3504]
$R_{1}$ ($I$ $>$ 4$\sigma$) [# of reflections] | 0.0345 [2842]
GOF [# of reflections] | 0.986 [3504]
Structural analysis without using anisotropic displacement parameters∗
split-site model | Use | Not used
$R_{1}$ [# of reflections] |  0.0459 [3504] | 0.0469 [3504]
$R_{1}$ ($I$ $>$ 4$\sigma$) [# of reflections] | 0.0349 [2842] | 0.0359 [2842]
GOF [# of reflections] | 0.993 [3504] | 1.017 [3504]

* In this analysis, isotropic ADP is used only for the Nb1 site, and anisotropic ADP is used for the remaining sites.

Table 12: Summary of crystallographic data of 1$T$-NbTe2 at 100 K without the
split-site model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$ (Å2)
Nb1 | 2$d$ | 1 | 1/2 | -1/2 | 0 | 0.00464(4)
Nb2 | 4$i$ | 1 | 0.63853(2) | 0 | 0.20955(2) | 0.00385(3)
Te1 | 4$i$ | 1 | 0.79750(2) | 1/2 | 0.37881(2) | 0.00392(3)
Te2 | 4$i$ | 1 | 0.503383(11) | -1/2 | 0.30954(2) | 0.00391(3)
Te3 | 4$i$ | 1 | 0.64991(2) | 0 | 0.00941(2) | 0.00385(3)

Table 13: Anisotropic atomic displacement parameters of 1$T$-NbTe2 at 100 K
without the split-site model.

site |  $U_{11}$ (Å2) |  $U_{22}$ (Å2) |  $U_{33}$ (Å2) |  $U_{12}$ (Å2) |  $U_{13}$ (Å2) |  $U_{23}$ (Å2)
---|---|---|---|---|---|---
Nb1 | 0.00440(9) | 0.00579(10) | 0.00332(9) | 0 | 0.00079(8) | 0
Nb2 | 0.00426(7) | 0.00336(7) | 0.00319(7) | 0 | 0.00030(5) | 0
Te1 | 0.00431(6) | 0.00301(5) | 0.00359(6) | 0 | 0.00025(4) | 0
Te2 | 0.00465(6) | 0.00315(5) | 0.00350(5) | 0 | 0.00086(4) | 0
Te3 | 0.00457(6) | 0.00290(5) | 0.00353(5) | 0 | 0.00069(4) | 0

Table 14: Structural parameters of 1$T$-NbTe2 at 100 K with the split-site
model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$($U_{iso}$) (Å2)
Nb1 | 4$g$ | 0.5 | 1/2 | 0.5113(4) | 0 | 0.00405(6)
Nb2 | 4$i$ | 1 | 0.63854(2) | 0 | 0.29028(2) | 0.00384(4)
Te1 | 4$i$ | 1 | 0.79751(2) | 1/2 | 0.37881(2) | 0.00391(3)
Te2 | 4$i$ | 1 | 0.50338(2) | -1/2 | 0.30954(2) | 0.00390(3)
Te3 | 4$i$ | 1 | 0.64991(2) | 0 | 0.00943(2) | 0.00383(3)

In the split-site model analysis, isotropic ADP is used only for the Nb1 site,
and anisotropic ADP is used for the other

elements. The fact that $U_{iso}$ at the Nb1 site and $U_{eq}$ at other sites
are almost equal indicates the validity of this analysis.

Table 15: Summary of crystallographic data of 1$T$-VTe2 at 100 K.

Temperature (K) | 100
---|---
Wavelength (Å) | 0.311
Crystal dimension ($\mu$m3) | 30$\times$30$\times$30
space group |  $C$2/$m$
$a$ (Å) | 14.3110(2)
$b$ (Å) | 3.59625(4)
$c$ (Å) | 9.09130(10)
$\beta$ (∘) | 109.602(2)
$V$ (Å3) | 440.770(9)
$Z$ | 6
$F$(000) | 762
(sin$\theta$/$\lambda$)Max (Å-1) | 1.7824
$N_{Total,obs}$ | 41843
$N_{Unique,obs}$ | 9196
Average redundancy | 4.6
Completeness | 0.831
Structural analysis using anisotropic displacement parameters
$R_{1}$ [# of reflections] | 0.0708 [8327]
$R_{1}$ ($I$ $>$ 4$\sigma$) [# of reflections] | 0.0607 [7414]
GOF [# of reflections] | 1.323 [8327]
Structural analysis without using anisotropic displacement parameters∗
split-site model | Use | Not used
$R_{1}$ [# of reflections] |  0.0711 [8327] | 0.0720 [8327]
$R_{1}$ ($I$ $>$ 4$\sigma$) [# of reflections] | 0.0610 [7414] | 0.0620 [7414]
GOF [# of reflections] | 1.319 [8327] | 1.316 [8327]

* In this analysis, isotropic ADP is used only for the V1 site, and anisotropic ADP is used for the remaining sites.

Table 16: Summary of crystallographic data of 1$T$-VTe2 at 100 K without the
split-site model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$ (Å2)
V1 | 2$d$ | 1 | 1/2 | -1/2 | 0 | 0.00503(6)
V2 | 4$i$ | 1 | 0.64217(4) | 0 | 0.29701(6) | 0.00431(4)
Te1 | 4$i$ | 1 | 0.79577(2) | 1/2 | 0.37769(2) | 0.00396(2)
Te2 | 4$i$ | 1 | 0.50775(2) | -1/2 | 0.30650(2) | 0.00377(2)
Te3 | 4$i$ | 1 | 0.64568(2) | 0 | 0.01414(2) | 0.00378(2)

Table 17: Anisotropic atomic displacement parameters of 1$T$-VTe2 at 100 K
without the split-site model.

site |  $U_{11}$ (Å2) |  $U_{22}$ (Å2) |  $U_{33}$ (Å2) |  $U_{12}$ (Å2) |  $U_{13}$ (Å2) |  $U_{23}$ (Å2)
---|---|---|---|---|---|---
V1 | 0.00440(13) | 0.00684(16) | 0.00345(13) | 0 | 0.00079(11) | 0
V2 | 0.00447(8) | 0.00465(9) | 0.00319(9) | 0 | 0.00047(7) | 0
Te1 | 0.00407(3) | 0.00373(3) | 0.00334(3) | 0 | 0.00026(3) | 0
Te2 | 0.00420(3) | 0.00365(3) | 0.00297(3) | 0 | 0.00055(3) | 0
Te3 | 0.00403(3) | 0.00353(3) | 0.00327(3) | 0 | 0.00053(3) | 0

Table 18: Structural parameters of 1$T$-VTe2 at 100 K with the split-site
model.

| atomic coordinates |
---|---|---
site | Wyck. | Occ. | $x/a$ | $y/b$ | $z/c$ |  $U_{eq}$($U_{iso}$) (Å2)
V1 | 4$g$ | 0.5 | 1/2 | 0.5135(5) | 0 | 0.00428(7)
V2 | 4$i$ | 1 | 0.64217(4) | 0 | 0.29696(6) | 0.00431(4)
Te1 | 4$i$ | 1 | 0.79577(2) | 1/2 | 0.37769(2) | 0.00396(2)
Te2 | 4$i$ | 1 | 0.50775(2) | -1/2 | 0.30650(2) | 0.00377(2)
Te3 | 4$i$ | 1 | 0.64568(2) | 0 | 0.01414(2) | 0.00378(2)

In the split-site model analysis, isotropic ADP is used only for the V1 site,
and anisotropic ADP is used for the other

elements. The fact that $U_{iso}$ at the V1 site and $U_{eq}$ at other sites
are almost equal indicates the validity of this analysis.

## References

  * Okamoto  _et al._ [2008] Y. Okamoto, S. Niitaka, M. Uchida, T. Waki, M. Takigawa, Y. Nakatsu, A. Sekiyama, S. Suga, R. Arita, and H. Takagi, Phys. Rev. Lett. 101, 086404 (2008).
  * Radaelli  _et al._ [2002] P. G. Radaelli, Y. Horibe, M. J. Gutmann, H. Ishibashi, C. H. Chen, R. M. Ibberson, Y. Koyama, Y. S. Hor, V. Kiryukhin, and S. W. Cheong, Nature (London) 416, 155 (2002).
  * Kojima  _et al._ [2019] K. Kojima, N. Katayama, S. Tamura, M. Shiomi, and H. Sawa, Phys. Rev. B 100, 235120 (2019).
  * Kojima  _et al._ [2023] K. Kojima, N. Katayama, Y. Matsuda, M. Shiomi, R. Ishii, and H. Sawa, Phys. Rev. B 107, L020101 (2023).
  * Katayama  _et al._ [2009] N. Katayama, M. Uchida, D. Hashizume, S. Niitaka, J. Matsuno, D. Matsumura, Y. Nishihata, J. Mizuki, N. Takeshita, A. Gauzzi, M. Nohara, and H. Takagi, Phys. Rev. Lett. 103, 146405 (2009).
  * Horibe  _et al._ [2006] Y. Horibe, M. Shingu, K. Kurushima, H. Ishibashi, N. Ikeda, K. Kato, Y. Motome, N. Furukawa, S. Mori, and T. Katsufuji, Phys. Rev. Lett. 96, 086406 (2006).
  * Katayama  _et al._ [2020] N. Katayama, H. Takeda, T. Yamaguchi, Y. Yamada, K. Iida, M. Takigawa, Y. Ohta, and H. Sawa, Phys. Rev. B 102, 081106(R) (2020).
  * Katayama  _et al._ [2018] N. Katayama, S. Tamura, T. Yamaguchi, K. Sugimoto, K. Iida, T. Matsukawa, A. Hoshikawa, T. Ishigaki, S. Kobayashi, Y. Ohta, and H. Sawa, Phys. Rev. B 98, 081104(R) (2018).
  * Okamoto  _et al._ [2020] Y. Okamoto, H. Amano, N. Katayama, H. Sawa, K. Niki, R. Mitoka, H. Harima, T. Hasegawa, N. Ogita, Y. Tanaka, M. Takigawa, Y. Yokoyama, K. Takehana, Y. Imanaka, Y. Nakamura, H. Kishida, and K. Takenaka, Nat. Commun. 11, 3144 (2020).
  * Miura  _et al._ [2007] Y. Miura, Y. Yasui, M. Sato, N. Igawa, and K. Kakurai, J. Phys. Soc. Jpn. 76, 033705 (2007).
  * Hirai  _et al._ [2022] D. Hirai, K. Kojima, N. Katayama, M. Kawamura, D. Nishio-Hamane, and Z. Hiroi, J. Am. Chem. Soc. 144, 17857 (2022).
  * Kobayashi  _et al._ [2019] S. Kobayashi, N. Katayama, T. Manjo, H. Ueda, C. Michioka, J. Sugiyama, Y. Sassa, O. K. Forslund, M. Månsson, K. Yoshimura, and H. Sawa, Inorg. Chem. 58, 14304 (2019).
  * Schmidt  _et al._ [2004] M. Schmidt, W. Ratcliff, P. G. Radaelli, K. Refson, N. M. Harrison, and S. W. Cheong, Phys. Rev. Lett. 92, 056402 (2004).
  * Browne  _et al._ [2017] A. J. Browne, S. A. J. Kimber, and J. P. Attfield, Phys. Rev. Mater. 1, 052003(R) (2017).
  * Kimber  _et al._ [2014] S. A. J. Kimber, I. I. Mazin, J. Shen, H. O. Jeschke, S. V. Streltsov, D. N. Argyriou, R. Valent$\acute{i}$, and D. I. Khomskii, Phys. Rev. B 89, 081408(R) (2014).
  * Shiomi  _et al._ [2022] M. Shiomi, K. Kojima, N. Katayama, S. Maeda, J. A. Schneeloch, S. Yamamoto, K. Sugimoto, Y. Ohta, D. Louca, Y. Okamoto, and H. Sawa, Phys. Rev. B 105, L041103 (2022).
  * Koch  _et al._ [2022] R. J. Koch, N. Aryal, O. Ivashko, Y. Liu, M. Abeykoon, E. D. Bauer, M. v. Zimmermann, W. Yin, C. Petrovic, and E. S. Bozin, Phys. Rev. B 106, 214516 (2022).
  * Bozin  _et al._ [2019] E. S. Bozin, W. G. Yin, R. J. Koch, M. Abeykoon, Y. S. Hor, H. Zheng, H. C. Lei, C. Petrovic, J. F. Mitchell, and S. J. L. Billinge, Nat. Commun. 10, 3638 (2019).
  * Katayama  _et al._ [2021] N. Katayama, K. Kojima, T. Yamaguchi, S. Hattori, S. Tamura, K. Ohara, S. Kobayashi, K. Sugimoto, Y. Ohta, K. Saitoh, and H. Sawa, npj Quantum Mater. 6, 16 (2021).
  * Torigoe  _et al._ [2018] S. Torigoe, T. Hattori, K. Kodama, T. Honda, H. Sagayama, K. Ikeda, T. Otomo, H. Nitani, H. Abe, H. Murakawa, H. Sakai, and N. Hanasaki, Phys. Rev. B 98, 134443 (2018).
  * Ko  _et al._ [2015] K.-T. Ko, H.-H. Lee, D.-H. Kim, J.-J. Yang, S.-W. Cheong, M. Eom, J. Kim, R. Gammag, K.-S. Kim, H.-S. Kim, T.-H. Kim, H.-W. Yeom, T.-Y. Koo, K. H.-D., and P. J.-H., Nat. Commun. 6, 7342 (2015).
  * Toriyama  _et al._ [2014] T. Toriyama, M. Kobori, T. Konishi, Y. Ohta, K. Sugimoto, J. Kim, A. Fujiwara, S. Pyon, K. Kudo, and M. Nohara, J. Phys. Soc. Jpn. 83, 033701 (2014).
  * Pyon  _et al._ [2012] S. Pyon, K. Kudo, and M. Nohara, J. Phys. Soc. Jpn. 81, 053701 (2012).
  * Bronsema  _et al._ [1984] K. Bronsema, G. Bus, and G. Wiegers, J. Solid State Chem. 53, 415 (1984).
  * Brown [1966] B. Brown, Acta Crystallogr. 20, 264 (1966).
  * Whangbo and Canadell [1992] M.-H. Whangbo and E. Canadell, J. Am. Chem. Soc. 114, 9587 (1992).
  * Mitsuishi  _et al._ [2020] N. Mitsuishi, Y. Sugita, M. Bahramy, M. Kamitani, T. Sonobe, M. Sakano, T. Shimojima, H. Takahashi, H. Sakai, K. Horiba, H. Kumigashira, K. Taguchi, K. Miyamoto, T. Okuda, S. Ishiwata, Y. Motome, and K. Ishizaka, Nat. Commun. 11, 2466 (2020).
  * Sörgel  _et al._ [2006] T. Sörgel, U. Nuss, J. abd Wedig, R. Kremer, and M. Jansen, Mater. Res. Bull. 41, 987 (2006).
  * Petříček  _et al._ [2014] V. Petříček, M. Dušek, L. Palatinus, and K. logr. Z., Cryst. Mater. 229, 345 (2014).
  * Momma and Izumi [2011] K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).
  * Oishi-Tomiyasu [2014] R. Oishi-Tomiyasu, J. Appl. Crystallogr. 47, 593 (2014).
  * Izumi and Momma [2007] F. Izumi and K. Momma, Solid State Phenom. 130, 15 (2007).
  * Petkov  _et al._ [2020] V. Petkov, K. Chapagain, J. Yang, S. Shastri, and Y. Ren, Phys. Rev. B 102, 024111 (2020).
  * Li  _et al._ [2021] S. Li, Q. Dong, J. Feng, Y. Wang, M. Hou, W. Deng, R. A. Susilo, N. Li, H. Dong, S. Wan, C. Gao, and B. Chen, Inorg. Chem. 60, 7857 (2021).
  * Khomskii and Streltsov [2021] D. Khomskii and S. Streltsov, Chemical Reviews 121, 2992 (2021).
  * Ito  _et al._ [2017] M. Ito, K. Sonoda, and S. Nagata, Solid State Commun. 23, 265 (2017).
  * Balcerek  _et al._ [1999] K. Balcerek, C. Marucha, R. Wawryk, T. Tyc, N. Matsumoto, and S. Nagata, Philos. Mag. B. 79, 1021 (1999).
  * Kang  _et al._ [1998] H. Kang, P. Mandal, I. Medvedeva, J. Liebe, G. Rao, K. Bärner, P. A., and G. E., J. Appl. Phys. 83, 6977 (1998).