#  OGLE-2018-BLG-0567Lb and OGLE-2018-BLG-0962Lb: Two Microlensing Planets
through Planetary-Caustic Channel

Youn Kil Jung1,2,28, Cheongho Han3,28, Andrzej Udalski4,29, Andrew
Gould1,5,6,28, Jennifer C. Yee7,28,
and
Michael D. Albrow8, Sun-Ju Chung1,2, Kyu-Ha Hwang1, Yoon-Hyun Ryu1, In-Gu
Shin1, Yossi Shvartzvald9, Wei Zhu10, Weicheng Zang11, Sang-Mok Cha1,12, Dong-
Jin Kim1, Hyoun-Woo Kim1, Seung-Lee Kim1,2, Chung-Uk Lee1,2, Dong-Joo Lee1,
Yongseok Lee1,12, Byeong-Gon Park1,2, Richard W. Pogge5
(The KMTNet Collaboration)
Przemek Mróz4,13, Michał K. Szymański4, Jan Skowron4, Radek Poleski4,5, Igor
Soszyński4, Paweł Pietrukowicz4, Szymon Kozłowski4, Krzystof Ulaczyk14,
Krzysztof A. Rybicki4, Patryk Iwanek4, Marcin Wrona4
(The OGLE Collaboration)
1Korea Astronomy and Space Science Institute, Daejon 34055, Republic of Korea
2University of Science and Technology, Korea, 217 Gajeong-ro Yuseong-gu,
Daejeon 34113, Korea  3Department of Physics, Chungbuk National University,
Cheongju 28644, Republic of Korea  4Warsaw University Observatory, Al.
Ujazdowskie 4, 00-478 Warszawa, Poland  5Department of Astronomy, Ohio State
University, 140 W. 18th Ave., Columbus, OH 43210, USA  6Max-Planck-Institute
for Astronomy, K$\rm\ddot{o}$nigstuhl 17, 69117 Heidelberg, Germany  7Center
for Astrophysics $|$ Harvard & Smithsonian, 60 Garden St., Cambridge, MA
02138, USA  8University of Canterbury, Department of Physics and Astronomy,
Private Bag 4800, Christchurch 8020, New Zealand  9Department of Particle
Physics and Astrophysics, Weizmann Institute of Science, Rehovot 76100, Israel
10Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St George Street, Toronto, ON M5S 3H8, Canada  11Department of Astronomy,
Tsinghua University, Beijing 100084, China  12School of Space Research, Kyung
Hee University, Yongin 17104, Republic of Korea  13Division of Physics,
Mathematics, and Astronomy, California Institute of Technology, Pasadena, CA
91125, USA  14Department of Physics, University of Warwick, Gibbet Hill Road,
Coventry, CV4 7AL, UK

###### Abstract

We present the analyses of two microlensing events, OGLE-2018-BLG-0567 and
OGLE-2018-BLG-0962. In both events, the short-lasting anomalies were densely
and continuously covered by two high-cadence surveys. The light-curve modeling
indicates that the anomalies are generated by source crossings over the
planetary caustics induced by planetary companions to the hosts. The estimated
planet/host separation (scaled to the angular Einstein radius $\theta_{\rm
E}$) and mass ratio are $(s,q)=(1.81,1.24\times 10^{-3})$ and
$(s,q)=(1.25,2.38\times 10^{-3})$, respectively. From Bayesian analyses, we
estimate the host and planet masses as $(M_{\rm h},M_{\rm
p})=(0.24_{-0.13}^{+0.16}\,M_{\odot},0.32_{-0.16}^{+0.34}\,M_{\rm J})$ and
$(M_{\rm h},M_{\rm
p})=(0.55_{-0.29}^{+0.32}\,M_{\odot},1.37_{-0.72}^{+0.80}\,M_{\rm J})$,
respectively. These planetary systems are located at a distance of
$7.07_{-1.15}^{+0.93}\,{\rm kpc}$ for OGLE-2018-BLG-0567 and
$6.47_{-1.73}^{+1.04}\,{\rm kpc}$ for OGLE-2018-BLG-0962, suggesting that they
are likely to be near the Galactic bulge. The two events prove the capability
of current high-cadence surveys for finding planets through the planetary-
caustic channel. We find that most published planetary-caustic planets are
found in Hollywood events in which the source size strongly contributes to the
anomaly cross section relative to the size of the caustic.

gravitational lensing: micro – planetary systems

2828footnotetext: The KMTNet Collaboration.2929footnotetext: The OGLE
Collaboration.

##  1 Introduction

The signature of a microlening planet is almost always a short-lasting anomaly
in the smooth and symmetric lensing light curve produced by the host of the
planet. In principle, the signature can appear at any position of the lensing
light curve (Gaudi, 2012). In reality, however, the signatures of planets
detected in the earlier phase of lensing experiments appeared mainly near the
peak of lensing light curves.

The bias toward central anomalies is mostly attributed to the limitation of
early lensing surveys. With a roughly 1 day cadence of the first-generation
survey experiments, e.g., the MACHO (Alcock et al., 1995), OGLE-I (Udalski et
al., 1992), MOA-I (Bond et al., 2001) surveys, it was difficult to detect
planetary signals lasting of order $1$ day or less by the survey experiments.
To meet the cadence requirement for planet detections, Gould & Loeb (1992)
proposed an observational mode, in which wide-field surveys with a low cadence
monitor a large area of sky mainly to detect lensing events, and followup
experiments conduct high-cadence observations for a small number of lensing
events detected by the surveys using a network of multiple narrow-field
telescopes. However, this mode of observations had the drawback that only a
handful of lensing events could be monitored by followup observations.
Combined with the low probability of planetary perturbations, this implied a
low planet detection rate for this phase of the experiments. In fact, for the
first several years, there were no securely detected planets using this mode,
although there was one tentative detection in the event MACHO 98-BLG-35 (Rhie
et al., 2000).

The first three microlensing detections were found using the survey+followup
strategy. In the first event, OGLE 2003-BLG-235/MOA 2003-BLG-53 (Bond et al.,
2004), the planet was found by the surveys, but the MOA survey carried out
additional followup observations in response to the planetary anomaly. The
next two planets, OGLE-2005-BLG-071Lb (Udalski et al., 2005) and
OGLE-2005-BLG-390Lb (Beaulieu et al., 2006), were both found through extensive
followup observations of known microlensing events that were initiated before
the planetary anomaly began. The discovery of OGLE-2005-BLG-071Lb provided a
practical lesson in the value of high-magnification events (Griest &
Safizadeh, 1998) for detecting planets through followup observations.

In the following years, the planet detection rate using the survey+followup
mode was substantially increased by focusing on events with very high
magnifications. Several factors contributed to the increase of the detection
rate. First, the planet detection efficiency for high-magnification events is
high. This is because a planet located in the lensing zone of its host always
induces a small central caustic near the position of the host, and the
trajectory of a high-magnification event passes close to the central caustic.
This yields a high probability that a planet will produce a perturbation and
also confine that perturbation to a short duration of time while the event is
highly magnified, not throughout the whole event. As a result, the time of the
planetary signal, i.e., the peak of the light curve, can be predicted in
advance and enable one to efficiently use resources for followup observations.
By contrast, predicting the time of a planetary signal through other channels
is difficult. Finally, highly magnified source stars are bright enough to be
observed with small-aperture telescopes, down to sub-meter amateur-class
telescopes, and this enables one to maximize available telescopes for followup
observations, e.g., OGLE-2005-BLG-071 (Udalski et al., 2005; Dong et al.,
2009). Thus, the planets detected from the survey+followup experiments were
detected mainly through the high-magnification channel, and this led to the
bias toward central caustic perturbations.

The current planetary lensing experiments are in the second phase, in which
lensing events are observed by high-cadence surveys. The observational cadence
of the lensing surveys in this phase has greatly increased with the employment
of large-format cameras yielding very wide fields of view. The MOA experiment
entered a new phase (MOA-II) by upgrading its instrument with a new wide-field
camera composed of ten $2{\rm k}\times 4{\rm k}$ chips yielding a $2.2\,{\rm
deg}^{2}$ field of view (Sumi et al., 2013). The OGLE survey is in its fourth
phase (OGLE-IV) using $1.4\,{\rm deg}^{2}$ camera composed of 32 $2{\rm
k}\times 4{\rm k}$ chips (Udalski et al., 2015). The KMTNet survey, which
commenced its full operation in 2016 (Kim et al., 2016), utilizes three
globally distributed telescopes, each of which has a camera with a $4.0\,{\rm
deg}^{2}$ field of view. Being able to cover a large area of sky from a single
exposure, the observational cadence of the current survey experiments now
reaches $\Gamma\sim 4\,{\rm hr}^{-1}$ toward the dense bulge fields. This
enables planet detections without additional followup observations.

With the operation of the high-cadence surveys, the detection rate of planets
is rapidly increasing. One important reason for the rapid increase of the
detection rate is that planets can be detected not only through the central-
caustic channel but also through the additional planetary-caustic channel.
Planets are detected through the planetary-caustic channel as anomalies
produced by the source’s approach close to the “planetary caustic”, which
denotes one of the two sets of planet-induced caustics lying away from the
host. The planetary caustic lies at a position with a separation from the host
of $s-1/s$, and thus planetary signals produced by this caustic can appear at
any part of the lensing light curve depending on the planet-host separation
$s$ (normalized to the angular Einstein radius $\theta_{\rm E}$). The
planetary caustic is substantially larger than the central caustic, and thus
the probability of a planetary perturbation is higher. Another importance of
detecting planets through the planetary-caustic channel is that interpreting
the planetary signal is usually not subject to the close-wide degeneracy
(Griest & Safizadeh, 1998; Dominik, 1999), which causes ambiguity in
estimating the planet-host separations for most planets detected through the
central-caustic channel.

In this paper, we present the analysis of two planetary microlensing events
OGLE-2018-BLG-0567 and OGLE-2018-BLG-0962, for which planets are both detected
through a planetary-caustic channel. For both events, the signatures of the
planets were densely and continuously covered by two high-cadence lensing
surveys, and this leads us to unambiguously interpret the planetary signals.

##  2 Observation

The two planetary events were observed by the two lensing surveys conducted by
the OGLE and KMTNet groups. The OGLE survey uses the $1.3\,{\rm m}$ telescope
that is located at the Las Campanas Observatory in Chile. The KMTNet survey
utilizes three $1.6\,{\rm m}$ telescopes that are located at the Siding Spring
Observatory in Australia (KMTA), the Cerro Tololo Interamerican Observatory in
Chile (KMTC), and the South African Astronomical Observatory in South Africa
(KMTS). The global distribution of the KMTNet telescopes makes it possible to
continuously monitor the events. In both surveys, observations were mainly
conducted in the $I$-band, and a fraction of the images were taken in the
$V$-band to determine the color of the microlensed source stars.

OGLE-2018-BLG-0567, $({\rm RA},{\rm Dec})_{\rm J2000}=$(17:56:04.42,
$-27$:59:13.6), or
$(l,b)=(1^{\circ}\hskip-2.0pt.99,-1^{\circ}\hskip-2.0pt.49)$ in Galactic
coordinates, was discovered on 2018 April 14 by the OGLE Early Warning System
(EWS: Udalski, 2003). The event was independently found by the KMTNet survey
as KMT-2018-BLG-0890 from its event-finding algorithm (Kim et al., 2018). The
observational cadence for the event is $\Gamma=1\,{\rm hr}^{-1}$ for OGLE and
$\Gamma=2\,{\rm hr}^{-1}$ for KMTNet.

OGLE-2018-BLG-0962 was discovered by the OGLE EWS on 2018 June 2. It is
located at $({\rm RA},{\rm Dec})_{\rm J2000}=$(17:52:41.95, $-32$:18:33.3) or
Galactic coordinates of
$(l,b)=(-2^{\circ}\hskip-2.0pt.11,-3^{\circ}\hskip-2.0pt.04)$. The OGLE
cadence for this direction is 3-10 times per night. The KMTNet collaboration
also observed the event, with designation KMT-2018-BLG-2071, located in their
two overlapping fields (BLG41 and BLG22). In combination, the KMTNet
observations have a frequency of $\Gamma=3\,{\rm hr}^{-1}$ for KMTC and
$\Gamma=2.25\,{\rm hr}^{-1}$ for KMTS and KMTA.

For both events, the data sets were reduced based on the image subtraction
methodology (Tomaney & Crotts, 1996; Alard & Lupton, 1998), specifically
Albrow et al. (2009) for KMTNet and Woźniak (2000) for OGLE. The photometric
error bars were then readjusted following the prescription presented in Yee et
al. (2012). We note that for the source color measurement, we additionally
carried out pyDIA (Albrow, 2017) reductions for a subset of the KMTNet data,
which simultaneously returns the light curve and field-star photometry on the
same system.

##  3 Light Curve Analysis

Figures 1 and 2 show the light curves of OGLE-2018-BLG-0567 and
OGLE-2018-BLG-0962, respectively. It is found that the two events share
various characteristics in common. First, the apparent peak magnifications of
the baseline single-lens single-source (1L1S) light curves are not high:
$A_{\rm peak}\sim 1.6$ for OGLE-2018-BLG-0567 and $A_{\rm peak}\sim 4.9$ for
OGLE-2018-BLG-0962. Second, the light curves of both events exhibit strong
short-term positive anomalies from the baseline 1L1S curves. Third, the
anomalies appear when the 1L1S-model lensing magnifications are low. All these
characteristics strongly suggest that the anomalies are produced by source
crossings over the planetary caustics induced by planetary companions to the
lenses. We, therefore, start with a binary-lens single-source (2L1S) modeling
of the events under the interpretation that the events were produced by lenses
composed of two masses, $M_{1}$ and $M_{2}$.

The standard 2L1S modeling requires one to include seven fitting parameters to
describe an observed light curve. The first three are the Paczyński (1986)
parameters $(t_{0},u_{0},t_{\rm E})$, respectively the time of closest source
approach to the lens, the impact parameter (scaled to $\theta_{\rm E}$), and
the event timescale. The next three $(s,q,\alpha)$ describe the binary-lens
geometry: the projected binary separation (scaled to $\theta_{\rm E}$), the
binary mass ratio $(q=M_{2}/M_{1})$, and the orientation angle of the binary
axis (relative to the source trajectory), respectively. The last describes the
source radius $\rho=\theta_{*}/\theta_{\rm E}$, where $\theta_{*}$ is the
angular source radius.

The modeling is conducted following the procedure described in Jung et al.
(2015). In the first step, we carry out grid searches for the binary
parameters $(s,q,\alpha)$. At each grid, we fix $(s,q)$ and find the remaining
parameters using Markov Chain Monte Carlo (MCMC) $\chi^{2}$ minimization. In
this modeling, the initial values of the parameters $(t_{0},u_{0},t_{\rm E})$
are given as the values estimated from a 1L1S fit for the data excluding the
anomaly. The initial value of the normalized source radius is estimated from
the caustic-crossing timescale, which is related to the event timescale by
$t_{*}=\rho t_{\rm E}$. Here we use inverse ray shooting (Kayser et al., 1986;
Schneider & Weiss, 1987) to compute the finite-source lensing magnifications.
The flux values from the source, $f_{{\rm S},i}$, and blend, $f_{{\rm B},i}$,
for the data set obtained from the $i$th observatory is estimated by
$f_{i}(t)=f_{{\rm S},i}A(t)+f_{{\rm B},i}$, where $f_{i}$ is the observed
flux. Once local solutions are found from the first-round modeling, we refine
the individual locals by releasing all fitting parameters and allowing them to
be free parameters in an MCMC.

From the modeling, it is found that the observed lensing light curves of both
events are well described by unique 2L1S models, in which the mass ratios
between $M_{1}$ and $M_{2}$ are in the planetary regime. The estimated binary
parameters are $(s,q)=(1.81,1.24\times 10^{-3})$ for OGLE-2018-BLG-0567, and
$(s,q)=(1.25,2.38\times 10^{-3})$ for OGLE-2018-BLG-0962. The full lensing
parameters and their uncertainties are presented in Table 1. The model curves
of the solutions are drawn over the data points in Figures 1 and 2 for
OGLE-2018-BLG-0567 and OGLE-2018-BLG-BLG-0962, respectively.

In Figures 3 and 4, we present the lens-system configurations of the
individual events, showing the source trajectory with respect to the lens
components and resulting caustics. From the configurations, it is found that
the anomalies of both events are produced by the source crossing over the
planetary caustic of the lens system. For OGLE-2018-BLG-0567, the source size
is comparable to the caustic size, and thus the detailed caustic-crossing
features, two caustic spikes and a U-shape trough region between the spikes,
was smeared out by finite-source effects. For OGLE-2018-BLG-0962, on the other
hand, the caustic is much bigger than the source size, and thus detailed
caustic-crossing feature of the anomaly is well delineated. It is found the
first part of the anomaly, centered at ${\rm HJD}^{\prime}(={\rm
HJD}-2,450,000\,{\rm days})\sim 8271.5$, was produced by the source passing
over the two caustic segments that flank the inner cusp (on the binary axis)
of the planetary caustic, and the second part, centered at ${\rm
HJD}^{\prime}\sim 8273.8$, was generated by the source passage over the
adjacent (off-axis) cusp.

We investigate the possibility of other interpretations of the events.
Especially for OGLE-2018-BLG-0567, the short-term perturbation might, in
principle, be produced by a second source (1L2S) with a large flux ratio
(Gaudi, 1998). Hence, we conduct an additional modeling of the event with the
1L2S interpretation (Jung et al., 2017). We search for the solution with eight
fitting parameters: $2\times(t_{0},u_{0},\rho)$ for the two sources, $I$-band
flux ratio $q_{F,I}$, and a shared timescale $t_{\rm E}$. The results are
listed in Table 2. We find that the 1L2S solution is disfavored by
$\Delta\chi^{2}>500$. In addition, the 1L2S model not only provides a worse
fit to the peak of the perturbation, but also fails to recover the decrease in
magnitude seen before and after the peak. See Figure 5. We also investigate
higher-order effects by conducting additional modeling considering the
microlens-parallax and lens-orbital effects (Gould, 1992; Dominik, 1998), and
find that it is difficult to constrain the extra lensing parameters related to
these higher-order effects from the observed light curves.

##  4 Physical Parameters

The results in Table 1 show that for both events, the normalized source radii
are precisely measured. This enables us to estimate $\theta_{\rm
E}=\theta_{*}/\rho$ provided that $\theta_{*}$ is measured. Then, we can use
the estimated $\theta_{\rm E}$ to constrain the lens total mass $M$ and
distance $D_{\rm L}$ as given by

$\theta_{\rm E}^{2}\equiv\kappa M\pi_{\rm rel};\qquad\pi_{\rm rel}={\rm
au}\left({{1\over D_{\rm L}}-{1\over D_{\rm S}}}\right),$ (1)

where $\kappa\equiv{4G/(c^{2}{\rm au})}\simeq 8.14\,{{\rm mas}/M_{\odot}}$,
$\pi_{\rm rel}$ is the lens-source relative parallax, and $D_{\rm S}$ is the
source distance. Hence, we evaluate $\theta_{*}$ using the method of Yoo et
al. (2004).

We first estimate the intrinsic source color $(V-I)_{\rm 0,S}$ and magnitude
$I_{\rm 0,S}$. We do this in the following ways. First, we calibrate the
KMTNet pyDIA photometry to the standard Johnson-Cousins system using the OGLE-
III catalog (Szymański et al., 2011). We next construct a $(V-I,I)$ color-
magnitude diagram (CMD) with field stars around the source. We next find the
source location $(V-I,I)_{\rm S}$ in the CMD from the best-fit model. We also
find the position of the red clump centroid (RCC), i.e., $(V-I,I)_{\rm RCC}$.
Figure 6 shows the source and RCC positions in the CMDs for the individual
events. We then measure the offset $\Delta(V-I,I)=(V-I,I)_{\rm S}-(V-I,I)_{\rm
RCC}$. Finally, we find the intrinsic source position as

$(V-I,I)_{\rm 0,S}=\Delta(V-I,I)+(V-I,I)_{\rm 0,RCC},$ (2)

where $(V-I,I)_{\rm 0,RCC}$ is the intrinsic RCC position measured from
independent observations (Bensby et al., 2013; Nataf et al., 2013). Here we
assume that the source star experiences the same amount of extinction as the
RCC. In Table 3, we list our estimated values of $(V-I,I)_{\rm S}$,
$(V-I,I)_{\rm RCC}$, $(V-I,I)_{\rm 0,RCC}$, and $(V-I,I)_{\rm 0,S}$.

We now derive $\theta_{*}$ based on the estimated intrinsic source position.
For this, we apply $(V-I)_{\rm 0,S}$ to the $VIK$ relation (Bessell & Brett,
1988) to find $(V-K)_{\rm 0,S}$. We then obtain $\theta_{*}$ from the
$(V-K)_{\rm 0,S}-\theta_{*}$ relations, specifically Kervella et al. (2004a)
for OGLE-2018-BLG-0567 and Kervella et al. (2004b) for OGLE-2018-BLG-0962. We
note that we add $5\%$ error to $\theta_{*}$ to consider the uncertainty of
the intrinsic RCC position and the color/surface brightness conversion. From
the measured $t_{\rm E}$ and $\rho$, we then find $\theta_{\rm E}$ and the
lens-source relative proper motion, $\mu_{\rm rel}=\theta_{\rm E}/t_{\rm E}$.
The estimated values of $\theta_{*}$, $\theta_{\rm E}$, and $\mu_{\rm rel}$
are also listed in Table 3.

For both events, we are unable to constrain the microlens parallax vector
$\mbox{\boldmath$\pi$}_{\rm E}$. This implies that we cannot directly derive
$M$ and $D_{\rm L}$ by

$M={\theta_{\rm E}\over\kappa\pi_{\rm E}};\qquad D_{\rm L}={{\rm
au}\over\pi_{\rm E}\theta_{\rm E}+\pi_{\rm S}}$ (3)

from the microlensing data (Gould, 2000). Here $\pi_{\rm S}={\rm au}/D_{\rm
S}$. Hence, we estimate the lens properties based on a Bayesian analysis with
Galactic model priors.

The Galactic model is constructed based on a mass function (MF), a density
profile (DP), and a velocity distribution (VD). For the MF, we use the models
presented in Jung et al. (2018). For the DP, we use the Han & Gould (2003)
model for the bulge and the Robin et al. (2003) model for the disk,
respectively. We note that the former is normalized based on the star counts
results of Holtzman et al. (1998), while the latter is normalized by the local
column density $\Sigma(D_{\rm L}=0)=\Sigma_{\rm disk,0}=36\,{M_{\odot}}\,{{\rm
pc}^{-2}}$ from Han & Gould (2003).

The bulge VD is modeled based on stars in the Gaia catalog (Gaia Collaboration
et al., 2016, 2018). That is, we first find red giant stars in the catalog
within 2 arcmin centered on the event location. We then derive their mean
velocity and its dispersion in Galactocentric Cartesian coordinates $(x,y,z)$
as defined in Han & Gould (1995).

For the disk VD, we adopt Gaussian forms of $f(v_{y},v_{z})=f(v_{y})f(v_{z})$
from Han & Gould (1995), which we then modify to consider the change in the
matter distribution. We do this in two ways. First, we introduce an asymmetric
drift $v_{\rm ad}$ to the $y$-direction velocity $v_{y}$. That is, $v_{\rm
ad}(D_{\rm L})=(0.5\sigma_{xyz,0}^{2}/\bar{v}_{y})[\Sigma(D_{\rm
L})/\Sigma_{\rm disk,0}]^{1/2}$, where $\bar{v}_{y}=220\,{\rm km}\,{\rm
s}^{-1}$ is the mean $y$-direction velocity,
$\sigma_{xyz,0}^{2}=\sigma_{x,0}^{2}+\sigma_{y,0}^{2}+\sigma_{z,0}^{2}$, and
$(\sigma_{x,0},\sigma_{y,0},\sigma_{z,0})=(34,30,20)\,{\rm km}\,{\rm s}^{-1}$
are the mean velocity dispersions along the $(x,y,z)$ directions (in the solar
neighborhood), respectively. The $y$-direction velocity at a given line-of-
sight distance $D_{\rm L}$ is then calculated by $v_{y}(D_{\rm
L})=\bar{v}_{y}-v_{\rm ad}(D_{\rm L})$. Second, we subsequently modify the
velocity dispersion as $\sigma_{y}(D_{\rm L})=\sigma_{y,0}[\Sigma(D_{\rm
L})/\Sigma_{\rm disk,0}]^{1/2}$ for the $y$-direction and $\sigma_{z}(D_{\rm
L})=\sigma_{z,0}[\Sigma(D_{\rm L})/\Sigma_{\rm disk,0}]^{1/2}$ for the
$z$-direction.

We now carry out the Bayesian analysis with the constraints $(t_{\rm
E},\theta_{\rm E})$. For this, we follow the procedure of Jung et al. (2018).
The estimated lens properties for the individual events are listed in Table 4.
The corresponding posterior probabilities for $M_{1}$ and $D_{\rm L}$ are
shown in Figure 7. We find that the host mass is
$M_{1}=0.24_{-0.13}^{+0.26}\,M_{\odot}$ for OGLE-2018-BLG-0567 and
$M_{1}=0.55_{-0.29}^{+0.32}\,M_{\odot}$ for OGLE-2018-BLG-0962. The planet
masses $(M_{2}=qM_{1})$ and the projected planet-host separations
$(a_{\perp}=sD_{\rm L}\theta_{\rm E})$ of the individual events are then
estimated to $(M_{2},a_{\perp})=(0.32_{-0.16}^{+0.34}\,M_{\rm
J},2.72_{-0.55}^{+0.49}\,{\rm au})$ and
$(M_{2},a_{\perp})=(1.37_{-0.72}^{+0.80}\,M_{\rm J},3.57_{-1.02}^{+0.68}\,{\rm
au})$, respectively. These estimates suggest that the two planetary systems
are likely composed of M-type dwarfs and giant planets lying beyond the snow
line (Kennedy & Kenyon, 2008). Here, the snow line is the location in the
proto-planetary disk where icy material can condense and where giant planets
are thought to be formed (Ida & Lin, 2004). The distance to the planet is
$D_{\rm L}=7.07_{-1.15}^{+0.93}\,{\rm kpc}$ for OGLE-2018-BLG-0567 and $D_{\rm
L}=6.47_{-1.73}^{+1.04}\,{\rm kpc}$ for OGLE-2018-BLG-0962, indicating that
they are likely to be in or near the Galactic bulge.

##  5 Microlensing Planets in the $({\rm log}\,s,{\rm log}\,q)$ plane

Our two survey-only microlensing planets are detected from the perturbations
caused by the planetary caustics (see Figures 3 and 4). In particular, the
planetary perturbation of OGLE-2018-BLG-0567 was generated by a “Hollywood”
geometry (Gould, 1997), in which the source size contributes strongly to, or
dominates, the anomaly cross section relative to the size of the caustic.
These detections prove the capability of the high-cadence surveys for
detecting planets through the planetary-caustic channel.

Figure 8 illustrates the positions of OGLE-2018-BLG-0567Lb and
OGLE-2018-BLG-0962Lb in the $({\rm log}\,s,{\rm log}\,q)$ plane along with
other microlensing planets111https://exoplanetarchive.ipac.caltech.edu. This
figure shows that OGLE-2018-BLG-0567Lb is located in a previously
underpopulated region. The two green solid lines indicate the boundary between
resonant and non-resonant caustics (Schneider & Weiss, 1986; Dominik, 1999).
For planets just outside of the resonance region, the excess magnification
pattern (from the underlying 1L1S magnification) extends out all the way
between the central and planetary caustics. This implies that the effective
cross-section of these caustics for a planet-induced perturbation and their
lensing behavior are similar to those of resonant caustics. Recently, Yee et
al. (2021) classified such kinds of caustics as the “near-resonant”,
specifically a set of caustics that has a excess magnification contour that
connects the central and planetary caustics and has at least $10\%$ excess
magnification along the entire caustic ridges. The two green dashed lines are
the boundary of near-resonant caustics. Yee et al. (2021) empirically
estimated that at a fixed $q$, the maximum size of a $10\%$ deviation
corresponds to $\sim 3{\rm log}\,{s_{\rm r,c}}$ for $s<1$ (close) and $\sim
1.8{\rm log}\,{s_{\rm r,w}}$ for $s>1$ (wide), where $s_{\rm r,c}$ and $s_{\rm
r,w}$ are the boundary values of $s$ between resonant and non-resonant
caustics (Dominik, 1999). Planets (except for our two planets) with a single
and multiple solutions are coded by black and red (with connected lines)
colors, respectively. OGLE-2018-BLG-0567Lb and OGLE-2018-BLG-0962Lb are marked
by yellow and blue colors, respectively. The shape of the symbols represent
the type of caustics that yield the planetary perturbation: circles for
resonant/near-resonant, squares for central, and triangles for planetary
caustics. The filled triangles are the planets found in the Hollywood events.

Figure 8 shows that the majority of planets are located inside the near-
resonant boundary rather than being due to planetary caustics. The bias toward
resonant caustics mainly comes from the relatively large size of resonant
caustics (scaled as $q^{1/3}$) compared to that of planetary caustics (scaled
as $q^{1/2}$) (Dominik, 1999; Han, 2006; Yee et al., 2021). There is also a
bias against planetary caustics due to the observational strategy of earlier
microlensing experiments that was focused on high-magnification events (Griest
& Safizadeh, 1998; Gould et al., 2010).

Only $24$ planets are placed outside the near-resonant boundary and $18$
planets among them are detected from the perturbations produced by clearly
isolated planetary caustics222 The corresponding planetary-caustic events are
OGLE-2005-BLG-390 (Beaulieu et al., 2006), MOA-bin-1 (Bennett et al., 2012),
OGLE-2006-BLG-109 (Gaudi et al., 2008; Bennett et al., 2010),
OGLE-2008-BLG-092 (Poleski et al., 2014), MOA-2010-BLG-353 (Rattenbury et al.,
2015), MOA-2011-BLG-028 (Skowron et al., 2016), MOA-2012-BLG-006 (Poleski et
al., 2017), OGLE-2012-BLG-0838 (Poleski et al., 2020), OGLE-2013-BLG-0341
(Gould et al., 2014), MOA-2013-BLG-605 (Sumi et al., 2016), OGLE-2014-BLG-1722
(Suzuki et al., 2018), OGLE-2016-BLG-0263 (Han et al., 2017),
OGLE-2016-BLG-1227 (Han et al., 2020), KMT-2016-BLG-1107 (Hwang et al., 2019),
OGLE-2017-BLG-0173 (Hwang et al., 2018), OGLE-2017-BLG-0373 (Skowron et al.,
2018), OGLE-2018-BLG-0596 (Jung et al., 2019), and OGLE-2018-BLG-0962 (this
work).. We find that most of these planetary-caustic planets ($12$ planets)
are found in the Hollywood events and they are located in high-cadence
observational fields of the lensing surveys. This proves the capability of the
Hollywood strategy of following big stars to find planets (Gould, 1997). The
majority of the Hollywood planets are located in the region $s>1$. This is
mainly due to the difference in the size of planetary caustics. For $s>1$,
there is one four-sided planetary caustic. For $s<1$, on the other hand, there
are two triangular planetary caustics and each of which size is much smaller
than that of $s>1$. In addition, the planetary signals from the these smaller
planetary caustics tends to be more significantly diminished by the finite-
source effects (Gould & Gaucherel, 1997). As a result, the wide-planetary
caustic has a larger effective cross section and therefore higher sensitivity
for finding planets.

##  6 Summary and Conclusions

We present the discovery of two cold giant planets orbiting M-dwarfs in two
events OGLE-2018-BLG-0567 and OGLE-2018-BLG-0962. Both events clearly showed
deviations from the 1L1S model, caused by the presence of a companion to the
lens host with precisely measured planet/host mass ratios of $(1.24\pm
0.07)\times 10^{-3}$ and $(2.38\pm 0.08)\times 10^{-3}$, respectively. In both
events, the finite-source effects are clearly detected, but the microlens
parallax effects are not meaningfully constrained. Hence, we constrain the
lens properties using the Bayesian analysis. From this, we estimate planet
masses of $0.32_{-0.16}^{+0.34}\,M_{\rm J}$ for OGLE-2018-BLG-0567Lb and
$1.37_{-0.72}^{+0.80}\,M_{\rm J}$ for OGLE-2018-BLG-0962Lb, and their physical
projected separations of $2.72_{-0.55}^{+0.49}\,{\rm au}$ and
$3.57_{-1.02}^{+0.68}\,{\rm au}$, respectively. These planets likely belong to
a class of giant planets orbiting M-dwarfs outside the snow line. The
detection rate of microlensing planets has rapidly increased with the advent
of high-cadence lensing surveys, and thus, planets presented here and future
detections will expand our understanding of the planet population around
M-dwarfs.

The planet hosts can be precisely constrained by future high-resolution
imaging with adaptive optics (AO) mounted on “$30\,{\rm m}$” class telescopes.
That is, they have only a small probability of being non-luminous (see Figure
7). The Bayesian estimates suggest that the dereddened $H$-band magnitude of
the host is $H_{0}=21.90_{-1.82}^{+1.86}$ for OGLE-2018-BLG-0567 and
$H_{0}=19.47_{-1.52}^{+2.04}$ for OGLE-2018-BLG-0962 (Pecaut & Mamajek, 2013).
In addition, both events have $\mu_{\rm rel}>3.0\,{\rm mas}\,{\rm yr}^{-1}$.
Therefore, in 2030, the hosts will be separated from the microlensed source by
$\Delta\theta\gtrsim 40\,{\rm mas}$.

This research has made use of the KMTNet system operated by the Korea
Astronomy and Space Science Institute (KASI) and the data were obtained at
three host sites of CTIO in Chile, SAAO in South Africa, and SSO in Australia.
Work by CH was supported by the grants of National Research Foundation of
Korea (2017R1A4A1015178 and 2019R1A2C2085965). The OGLE has received funding
from the National Science Centre, Poland, grant MAESTRO 2014/14/A/ST9/00121 to
A.U.

## References

  * Alard & Lupton (1998) Alard, C., & Lupton, Robert H. 1998, ApJ, 503, 325
  * Alcock et al. (1995) Alcock, C., Allsman, R. A., Axelrod, T. S., et al. 1995, ApJ, 445, 133
  * Albrow (2017) Albrow, M. D. 2017, MichaelDAlbrow/pyDIA: Initial release on github, V1.0.0, Zenodo, doi:10.5281/zenodo.268049
  * Albrow et al. (2009) Albrow, M. D., Horne, K., Bramich, D. M., et al. 2009, MNRAS, 397, 2099
  * Beaulieu et al. (2006) Beaulieu, J.-P., Bennett, D. P., Fouqué, P., et al. 2006, Nature, 439, 437
  * Bennett et al. (2010) Bennett, D.P., Rhie, S.H., Nikolaev, S. et al. 2010, ApJ, 713, 837
  * Bennett et al. (2012) Bennett, D. P., Sumi, T., Bond, I. A., et al. 2012, ApJ, 757, 119
  * Bensby et al. (2013) Bensby, T., Yee, J. C., Feltzing, S., et al. 2013, A&A, 549, 147
  * Bessell & Brett (1988) Bessell, M. S., & Brett, J. M. 1988, PASP, 100, 1134
  * Bond et al. (2001) Bond, I. A., Abe, F., Dodd, R. J. 2001, MNRAS, 327, 868
  * Bond et al. (2004) Bond, I.A., Udalski, A., Jaroszyński, M. et al. 2004, ApJ, 606, L155
  * Dominik (1998) Dominik, M. 1998, A&A, 329, 361
  * Dominik (1999) Dominik, M. 1999, A&A, 349, 108
  * Dong et al. (2009) Dong, S., Gould, A., Udalski, A. et al. 2009, ApJ, 695, 970
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, 1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaudi (1998) Gaudi, B.S. 1998, ApJ, 506, 533
  * Gaudi (2012) Gaudi, B. S. 2012, ARA&A, 50, 411
  * Gaudi et al. (2008) Gaudi, B.S., Bennett, D.P., Udalski, A. et al. 2008, Science, 319, 927
  * Gould (1992) Gould, A. 1992, ApJ, 392, 442
  * Gould (1997) Gould, A. 1997, in Variables Stars and the Astrophysical Returns of the Microlensing Surveys, ed. R. Ferlet, J.-P. Maillard, & B. Raban (Gif-sur-Yvette, France: Editions Frontieres), 125
  * Gould (2000) Gould, A. 2000, ApJ, 542, 785
  * Gould & Gaucherel (1997) Gould, A. & Gaucherel, C. 1997, ApJ, 477, 580
  * Gould & Loeb (1992) Gould, A., & Loeb, A. 1992, ApJ, 396, 104
  * Gould et al. (2010) Gould, A., Dong, S., Gaudi, B. S., et al. 2010, ApJ, 720, 1073
  * Gould et al. (2014) Gould, A., Udalski, A., Shin, I.-G. et al. 2014, Science, 345, 46
  * Griest & Safizadeh (1998) Griest, K., & Safizadeh, N. 1998, ApJ, 500, 37
  * Han & Gould (1995) Han, C., & Gould, A. 1995, ApJ, 447, 53
  * Han & Gould (2003) Han, C., & Gould, A. 2003, ApJ, 592, 172
  * Han et al. (2017) Han, C., Udalski, A., Gould, A., et al. 2017, AJ, 154, 133
  * Han et al. (2020) Han, C., Udalski, A., Gould, A., et al. 2020, AJ, 159, 91
  * Holtzman et al. (1998) Holtzman, J. A., Watson, A. M., Baum, W. A., et al. 1998, AJ, 115, 1946
  * Hwang et al. (2019) Hwang, K.-H., Ryu, Y.-H., Kim, H.-W., et al. 2019, AJ, 157, 23
  * Hwang et al. (2018) Hwang, K.-H., Udalski, A., Shvartzvald, Y., et al. 2018, AJ, 155, 20
  * Ida & Lin (2004) Ida, S., & Lin, D. N. C. 2004, ApJ, 616, 567
  * Jung et al. (2019) Jung, Y. K., Gould, A., Udalski, A, et al. 2019, AJ, 158, 28
  * Jung et al. (2018) Jung, Y. K., Udalski, A., Gould, A., et al. 2018, AJ, 155, 219
  * Jung et al. (2015) Jung, Y. K., Udalski, A., Sumi, T., et al. 2015, ApJ, 798, 123
  * Jung et al. (2017) Jung, Y. K., Udalski, A., Yee, J. C., et al. 2017, AJ, 153, 129
  * Han (2006) Han, C. 2006, ApJ, 638, 1080
  * Kayser et al. (1986) Kayser, R., Refsdal S., & Stabell, R. 1986, A&A, 166, 36
  * Kennedy & Kenyon (2008) Kennedy, G. M., & Kenyon, S. J. 2008, ApJ, 673, 502
  * Kervella et al. (2004a) Kervella, P., Bersier, D., Mourard, D., et al. 2004a, A&A, 428, 587
  * Kervella et al. (2004b) Kervella, P., Thévenin, F., Di Folco, E., & Ségransan, D. 2004b, A&A, 426, 297
  * Kim et al. (2018) Kim, D.-J., Kim, H.-W., Hwang, K.-H., et al. 2018, AJ, 155, 76
  * Kim et al. (2016) Kim, S.-L., Lee, C.-U., Park, B.-G., et al. 2016, JKAS, 49, 37
  * Nataf et al. (2013) Nataf, D. M., Gould, A., Fouqué, P., et al. 2013, ApJ, 769, 88
  * Paczyński (1986) Paczyński, B. 1986, ApJ, 304, 1
  * Pecaut & Mamajek (2013) Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9
  * Poleski et al. (2014) Poleski, R., Skowron, J., Udalski, A., et al. 2014, ApJ, 795, 42
  * Poleski et al. (2020) Poleski, R., Suzuki, D., Udalski, A., et al. 2020, AJ, 159, 261
  * Poleski et al. (2017) Poleski, R., Udalski, A., Bond, I. A., et al. 2017, A&A, 604A, 103
  * Rattenbury et al. (2015) Rattenbury, N. J., Bennett, D. P., Sumi, T., et al. 2015, MNRAS, 454, 946
  * Rhie et al. (2000) Rhie, S. H., Bennett, D. P., Becker, A. C., et al. 2000, ApJ, 533, 378
  * Robin et al. (2003) Robin, A. C., Reylé, C., Derriére, S., & Picaud, S. 2003, A&A, 409, 523
  * Schneider & Weiss (1986) Schneider, P., & Weiss, A. 1986, A&A, 164, 237
  * Schneider & Weiss (1987) Schneider, P., & Weiss, A. 1987, A&A, 171, 49
  * Skowron et al. (2016) Skowron, J., Udalski, A., Poleski, R., et al. 2016, ApJ, 820, 4
  * Skowron et al. (2018) Skowron, J., Ryu, Y.-H., Hwang, K.-H., et al. 2018, AcA, 68, 43
  * Sumi et al. (2013) Sumi, T., Bennett, D. P., Bond, I. A., et al. 2013, ApJ, 778, 150
  * Sumi et al. (2016) Sumi, T., Udalski, A., Bennett, D. P., et al. 2016, ApJ, 825, 112
  * Suzuki et al. (2018) Suzuki, D., Bennett, D. P., Udalski, A., et al. 2018, AJ, 155, 263
  * Szymański et al. (2011) Szymański, M. K., Udalski, A., Soszyński, I., et al. 2011, AcA, 61, 83
  * Tomaney & Crotts (1996) Tomaney, A.B. & Crotts, A.P.S. 1996, AJ, 112, 2872
  * Udalski (2003) Udalski, A. 2003, Acta Astron., 53, 291
  * Udalski et al. (1992) Udalski, A., Szymański, M., Kałużny, J., Kubiak, M., & Mateo, M. 1992, AcA, 42, 253
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, AcA, 65, 1
  * Udalski et al. (2005) Udalski, A., Jaroszyński, M., Paczyński, B., et al. 2005, ApJ, 628, L109
  * Woźniak (2000) Woźniak, P. R. 2000, AcA, 50, 42
  * Yee et al. (2012) Yee, J. C., Shvartzvald, Y., Gal-Yam, A., et al. 2012, ApJ, 755, 102
  * Yee et al. (2021) Yee, J. C., Zang, W., Udalski, A., et al. 2021, arXiv:2101.04696
  * Yoo et al. (2004) Yoo, J., DePoy, D. L., Gal-Yam, A., et al. 2004, ApJ, 603, 139

Table 1: Lensing Parameters Parameters | OGLE-2018-BLG-0567 | OGLE-2018-BLG-0962
---|---|---
$\chi^{2}_{\rm tot}$/dof | 8677.3/9252 | 6892.1/6833
$t_{0}$ (${\rm HJD^{\prime}}$) | 8244.845$\pm$ 0.025 | 8262.494$\pm$ 0.053
$u_{0}$ | 0.733 $\pm$ 0.026 | 0.207 $\pm$ 0.032
$t_{\rm E}$ (days) | 24.641 $\pm$ 1.064 | 28.739 $\pm$ 0.298
$s$ | 1.806 $\pm$ 0.019 | 1.246 $\pm$ 0.027
$q$ ($10^{-3}$) | 1.240 $\pm$ 0.068 | 2.375 $\pm$ 0.076
$\alpha$ (rad) | 0.623 $\pm$ 0.045 | 0.590 $\pm$ 0.044
$\rho$ ($10^{-3}$) | 17.675 $\pm$ 0.831 | 1.137 $\pm$ 0.048
$f_{\rm S}$ | 0.842 $\pm$ 0.035 | 0.039 $\pm$ 0.007
$f_{\rm B}$ | 1.026 $\pm$ 0.035 | 0.274 $\pm$ 0.007
Table 2: Binary Source Model for OGLE-2018-BLG-0567 Parameters | 1L2S
---|---
$\chi^{2}$/dof | 9223.0/9251
$t_{0,1}$ (HJD′) | 8244.323$\pm$0.043
$u_{0,1}$ | 0.705$\pm$0.039
$t_{0,2}$ (HJD′) | 8670.134$\pm$0.010
$u_{0,2}$ $(10^{-5})$ | 2.603$\pm$35.093
$t_{\rm E}$ (days) | 24.088$\pm$0.868
$\rho_{1}$ | –
$\rho_{2}$ $(10^{-3})$ | 6.295$\pm$0.344
$q_{F,I}$ $(10^{-3})$ | 7.520$\pm$0.461
$f_{{\rm s}}$ | 0.795$\pm$0.082
$f_{{\rm b}}$ | 1.073$\pm$0.082
Table 3: Source Star and Lens Properties Parameters | OGLE-2018-BLG-0567 | OGLE-2018-BLG-0962
---|---|---
$(V-I,I)_{\rm S}$ | (3.59$\pm$0.10, 18.20$\pm$0.01) | (2.37$\pm$0.08, 21.40$\pm$0.03)
$(V-I,I)_{\rm RCC}$ | (3.67$\pm$0.02, 17.24$\pm$0.02) | (2.56$\pm$0.02, 16.59$\pm$0.03)
$(V-I,I)_{0,\rm RCC}$ | (1.06, 14.37) | (1.06, 14.56)
$(V-I,I)_{0,\rm S}$ | (0.98$\pm$0.10, 15.33$\pm$0.03) | (0.87$\pm$0.08, 19.37$\pm$0.04)
$\theta_{*}$ $({\mu}{\rm as})$ | 3.769$\pm$0.419 | 0.504$\pm$0.048
$\theta_{\rm E}$ (mas) | 0.213$\pm$0.026 | 0.443$\pm$0.045
$\mu_{\rm rel}$ $({\rm mas}~{}{\rm yr}^{-1})$ | 3.161$\pm$0.381 | 5.639$\pm$0.572
Table 4: Physical Parameters Parameters | OGLE-2018-BLG-0567 | OGLE-2018-BLG-0962
---|---|---
$M_{1}$ $(M_{\odot})$ | $0.24_{-0.13}^{+0.26}$ | $0.55_{-0.29}^{+0.32}$
$M_{2}$ $(M_{\rm J})$ | $0.32_{-0.16}^{+0.34}$ | $1.37_{-0.72}^{+0.80}$
$a_{\bot}$ (au) | $2.72_{-0.55}^{+0.49}$ | $3.57_{-1.02}^{+0.68}$
$D_{\rm L}$ (kpc) | $7.07_{-1.15}^{+0.93}$ | $6.47_{-1.73}^{+1.04}$
Figure 1:  Light curve of OGLE-2018-BLG-0567. The black solid curve on the
data is the best-fit 2L1S solution. The upper panel shows the enlarged view of
the planet-induced anomaly centered on ${\rm HJD}^{\prime}\sim 8270$. The
second and fourth panels show the residuals from the solution. The lensing
parameters of the solution are listed in Table 1 and the caustic geometry is
shown in Figure 3.  Figure 2:  Light curve of OGLE-2018-BLG-0962. The upper
panels show the close-up views of the regions around ${\rm HJD}^{\prime}\sim
8271.5$ (left) and ${\rm HJD}^{\prime}\sim 8273.8$ (right) when the planet-
induced perturbations occur. The lensing parameters of the 2L1S solution are
listed in Table 1 and the caustic geometry is shown in Figure 4.  Figure 3:
Caustic geometry of OGLE-2018-BLG-0567. The line with an arrow is the source
trajectory relative to the binary axis. The open circles (scaled by the
normalized source radius $\rho$) on the trajectory are the source positions at
the times of observations. The two orange circles are the positions of binary-
lens masses ($M_{1}$ and $M_{2}$). In each panel, the cuspy closed curve drawn
in black color represents the caustic. The upper panel shows the enlarged view
of the planetary caustic. Lengths are scaled to the angular Einstein radius of
the lens system.  Figure 4:  Caustic geometry of OGLE-2018-BLG-0962. Notations
are identical to those of Figure 3.  Figure 5:  Light curve of the 1L2S model
for OGLE-2018-BLG-0567. The dashed gray and solid black lines are the best-fit
models from the 1L2S and 2L1S interpretations, respectively. The lower two
panels show the residuals from the two models.  Figure 6:  Color-magnitude
diagrams of OGLE-2018-BLG-0567 (upper panel) and OGLE-2018-BLG-0962 (lower
panel). In each panel, the CMD is constructed using stars in the
$2^{\prime}\times 2^{\prime}$ field centered on the event location based on
KMTNet pyDIA photometry calibrated to the OGLE-III catalog (Szymański et al.,
2011). The blue and red circles are the positions of source and red clump
centroid, respectively.  Figure 7:  Posterior distributions of $M_{1}$ (left
panels) and $D_{\rm L}$ (right panels) for the individual events. In each
panel, the red and blue distributions are, respectively, the contributions by
the bulge and disk lens populations. The black distribution is the total
contribution of the two lens populations. The median value and its $68\%$
confidence interval are represented by the vertical solid and two dotted
lines, respectively.  Figure 8:  Microlensing planets in the $({\rm
log}\,s,{\rm log}\,q)$ plane, adapted from Figure 9 of Yee et al. (2021).
Planets (except for our two planets) are colored by the number of solutions:
black for one solution and red (with connected line) for degenerate solutions.
The two planets OGLE-BLG-2018-BLG-0567Lb and OGLE-2018-BLG-0962Lb are coded by
yellow and blue colors, respectively. Their shapes indicate the caustic
structure giving rise to the planetary perturbation: circles for
resonant/near-resonant, squares for central, and triangles for planetary
caustics. The filled triangles are the planets from the Hollywood events. The
two green solid and dashed lines are the boundary of resonant and near-
resonant caustics, respectively. We note that for compactness, we compress the
planet names, e.g., OGLE-2018-BLG-0567Lb to OB180567.


# Ferroelectric phases and phase transitions in CsGeBr3 induced by mechanical
load

Joshua Townsend<EMAIL_ADDRESS>Department of Physics, University of South
Florida, Tampa, Florida 33620, USA  Ravi Kashikar  Department of Physics,
University of South Florida, Tampa, Florida 33620, USA  S. Lisenkov  I.
Ponomareva<EMAIL_ADDRESS>Department of Physics, University of South
Florida, Tampa, Florida 33620, USA

###### Abstract

First-principles-based atomistic simulations are used to reveal ferroelectric
phases and phase transitions induced in a semiconductor ferroelectric,
CsGeBr3, by external loads: hydrostatic pressure, uniaxial and biaxial
stresses, and misfit strain. Hydrostatic pressure was found to suppress the
Curie point at the rate -0.45$T_{C}(0)$ K/GPa, where $T_{C}(0)$ is the zero
pressure Curie temperature. Stresses and misfit strains were found to induce
additional ferroelectric phase transitions and phases not available under
normal conditions. We find that tensile load significantly enhances both the
Curie temperature and spontaneous polarization, while compressive load has the
opposite effect but with the difference that the Curie temperature is only
slightly suppressed. The isothermal dependencies of polarization on pressure
and stresses are highly nonlinear, which could result in large nonlinear
piezoelectric responses. The phase diagrams reveal the diversity of the phases
accessible through mechanical load, which include tetragonal, orthorhombic and
monoclinic symmetries in addition to the rhombohedral and cubic ones
realizable under normal conditions. We believe that this work reveals the
potential of Ge-based halide perovskites for applications in energy converting
devices, which is especially significant in the current pursuit of
environmental friendly lead-free technologies.

##  I Introduction

Ferroelectric materials exhibit polarization in the absence of an electric
field. This spontaneous polarization typically appears below a certain
temperature, known as the Curie point. Ferroelectrics gained a significant
attention in the 1940s with the discovery of a group of exceptionally robust
ferroelectrics among oxide perovskites, compounds with a chemical formula
ABO3. The associated ferroelectric phase transitions in these materials give
rise to many scientifically appealing and technologically attractive
properties, such as extraordinary large dielectric constants reaching 10,000
in NBCTOdielectric1 ; dielectric2 , piezoelectricity ($d_{33}>\ 2500\ pC/N$ in
PZN-8% PThighpiezo ), pyroelectricity (up to 3000 to 5000 $\mu C\ m^{-2}\
K^{-1}$ in PMN-PTLargePyro ) , and nonlinear optical
propertiesReviewOxidePerovskites . The diverse properties exhibited by oxide
ferroelectrics have opened the door to many applications, such as mechanical
sensors, actuators, temperature sensing devices, and nonvolatile computer
memory. It is, however, uncommon to find semiconductor ferroelectrics among
oxide perovskites. For oxide perovskites, this is seen due to a large
difference in electronegativity between the metal cation (B) and the oxygen
(O)bandgap2 ; bandgap . Very recently, inorganic lead-free halide perovskite
family CsGeX3 (X=Cl, Br, I) has demonstrated ferroelectricityFerroHalides , as
well as many desirable optoelectronic properties, such as highly tunable
bandgaps in the visible light regimeTunableBandgap ; TunableOpt2 , nonlinear
optical propertiesNonlinearOpt , and high photocurrenthighphotocurrent ;
OpticalGeneral , to name a few. Unlike oxide ferroelectrics, which possess
large band gaps, the band gaps of the inorganic halide perovskites lie between
1.6 to 3.3 eVFerroPerovskites , making them suitable candidates for solar cell
applications, detectors, LED devices, and field effect transistors.
Furthermore, this particular group of inorganic halide perovskites may be able
to provide a lead-free alternative to many applications featuring lead-based
perovskites, providing a path toward environmentally friendly technologies.

This family of materials has been reported to undergo a single phase
transition from cubic Pm$\bar{3}$m phase to rhombohedral R3m phase at
temperatures of 428 K, 511 K, and 550 K for CsGeCl3, CsGeBr3, and CsGeI3,
respectivelythiele . Recent experimental studies have shown the existence of
ferroelectricity in CsGeX3 (X = Cl, Br, I) with spontaneous polarization of
12-15 $\mu$C/cm2 for X=Br and 20 $\mu$C/cm2 for X=IFerroPerovskites .
Currently, no experimental reports have been produced for polarization for
X=Cl. The computational study Ref.PolarAntipolar,  has proposed that a
dynamical antipolar phase may overlap with the cubic phase above the Curie
temperature. One of the exceptional features of ferroelectrics is that they
are also piezoelectrics, often times excellent ones. This suggests that
application of mechanical load may have a profound effect on the ferroelectric
phases and phase transitions in ferroelectrics. Indeed, we find multiple
examples of phase tunability by mechanical load among ferroelectric
oxidesPhenomenal1 ; ExpPress1 ; CompPress1 ; presstune ; straintune ;
straintune2 . At the same time virtually nothing is known about the
piezoelectricity and associated electro-mechanical couplings in Ge-based
halide ferroelectrics. How does the mechanical load affect ferroelectric phase
transitions in such materials? Application of mechanical load is known to
break the crystal symmetry, and therefore, is capable of inducing novel
phases, which often times may have lower symmetry and very different
properties. Lower symmetry phases may cause polarization rotation which in
turn is capable of enhancing piezoelectric response BFOpolarrotate . Moreover,
the Curie point of stressed materials could differ significantly from the
Curie point of unstressed materials and may open a way to further property
tunability through Curie temperature engineering. It should be noted that very
little is known about the response of Ge-based perovskites to mechanical load
at present. Structural phase transitions in CsGeBr3 were reported in response
to hydrostatic pressure, effectively lowering the Curie temperatureExpPress1 .
The bandgap could also be manipulated by the application of hydrostatic
pressureCompPress1 .

Computationally, addressing the aforementioned questions is possible through
application of finite-temperature methodologies that have been developed for
ferroelectrics. The examples are phase-field methodsPhaseField1 ,
phenomenological thermodynamic theoriesPhenomenal1 , first-principles-based
effective HamiltoniansVanderbilt1 ; Vanderbilt2  and bond-order
methodsBondOrder . It should be mentioned that Density Functional Theory (DFT)
is presently not used for finite-temperature simulations of ferroelctrics due
to its extreme computational cost. On the other hand, the effective
Hamiltonian approach Vanderbilt1 ; Vanderbilt2  has been specifically
developed to extend the reach of DFT simulations of ferroelectrics to finite
temperatures. Recently, such an approach has been developed for
CsGeBr3PolarAntipolar , which allows us to pursue the investigation of
ferroelectricity in this material under mechanical load.

The aim of this work is to use a multiscale approach that combines both
effective Hamiltonian and DFT simulations in order to: (i) predict
stress/strain induced stabilization of novel ferroelectric phases in CsGeBr3,
none of which are accessible under unstressed conditions; (ii) report rich and
unusual set of phase diagrams in this material; (iii) to quantify property
enhancement and tunability in this material by mechanical load, which feature
up to a 35.9% increase in Curie point and up to a 5.5% increase in
polarization under 1.0 GPa uniaxial tensile stress.

##  II Computational Methodology

To achieve our goals we take advantage of the effective Hamiltonian approach,
which was originally developed for oxide perovskites Vanderbilt1 ; Vanderbilt2
; ReviewOxidePerovskites  and extended for halides in Ref.PolarAntipolar, . In
such an approach only the degrees of freedom that are relevant for
ferroelectricity are retained. These are the local modes, $\mathbf{u}_{i}$,
which are proportional to the local dipole moment of the perovskite unit cell,
and both homogeneous and inhomogenoeus strain variables, ${\eta}$. The
collective motion of local modes describe the soft zone-center ferroelectric
mode. Strain variables describe long-wavelength acoustic phonons. Each unit
cell contains one local mode vector and one inhomogeneous strain tensor, which
allows capturing complex atomistic orderings of electric dipoles, such as
vortices, nanodomains, and skyrmionsPZTdomains ; nanodomains ; FEdomains .

The homogeneous strain variables allow for supercell deformation as a whole,
while inhomogeneous strain variables allow for the capture of unit cell
deformations. The energy of the effective Hamiltonian is

$E^{total}=E^{self}(\\{\mathbf{u_{i}}\\})\ +\ E^{dpl}(\\{\mathbf{u_{i}}\\})\
+\ E^{short}(\\{\mathbf{u_{i}}\\})\ +\ E^{elas}(\\{\mathbf{\eta_{i}}\\})\ +\
E^{int}(\\{\mathbf{u_{i}}\\},\\{\mathbf{\eta_{i}}\\})\ -\ \sigma_{j}\eta_{j}\
+\ PV$ (1)

and includes local mode self energy, $E^{self}$ (harmonic and anharmonic
contributions), long-range dipole-dipole interaction, $E^{dpl}$, short-range
interaction between local modes, $E^{short}$, elastic energy, $E^{elas}$, the
interaction between the local modes and strains, $E^{int}$, and the energies
that couple external stress, strain, pressure, and volume. The effective
Hamiltonian has been previously used to model ferroelectric phase transitions
in oxide ferroelectrics ReviewOxidePerovskites ; Vanderbilt1 ; Vanderbilt2 ;
RabePT ; BCOphase , and their tunability by the mechanical boundary conditions
Phenomenal1 ; BTOstrain ; McCash ; nanomech . In Ref. PolarAntipolar,  the
effective Hamiltonian was used to study the ferroelectric phase transition in
CsGeBr3. With r2SCAN based parametrization, it predicted the transition from
cubic Pm$\bar{3}$m phase to ferroelectric rhombohedral R3m phase at 340 K,
which underestimates the experimental transition temperature of 511
KPolarAntipolar . The underestimation may be caused by the overestimation of
the energy of the rhombohedral phase by the chosen exchange-correlation
functional. It has been proposed that the SCAN family functionals do not
adequately correct delocalization errorsSCANdelocal ; HSEcgb , which could be
an issue for structures with expressed lone pairs, as the rhombohedral phase
is in our case. Indeed, HSE increases the energy difference between the ground
state and cubic phase by a factor of 1.4, which could bring transition
temperatures predicted by effective Hamiltonian in agreement with experimental
data. However, parametrization with higher-level functionals, like HSE, is
presently computationally prohibitive. It was shown to reproduce experimental
coercive field and polarization. In this study we use the effective
Hamiltonian as parameterized from r2SCAN in the framework of classical
Molecular Dynamics (MD) simulations to model bulk CsGeBr3 at finite
temperatures and under a wide range of mechanical loads. We model a supercell
of 30$\times$30$\times$30 unit cells of CsGeBr3 repeated periodically along
all three Cartesian directions to simulate bulk crystal. Newton’s equations of
motion are solved numerically using a predictor-corrector
integratorUnderstandingSim ; ArtOfSim  and an integration step of 0.5 fs. The
temperature is controlled through the Evans-Hoover thermostatUnderstandingSim
. Equilibrium thermodynamical phases are harvested from annealing simulations,
which are initiated at 700 K and proceeded in steps of 10 K down to 10 K. In
the case of hydrostatic pressure we encountered TC higher than 700 K and those
simulations were redone starting at 800 K. For each temperature we simulated
300,000 MD steps, and used 100,000 steps for equilibration and 200,000 steps
for collecting thermal averages. Annealing simulations have been carried out
for each mechanical load investigated. We studied hydrostatic pressure,
uniaxial and biaxial stresses, and misfit, or epitaxial strain. These are
motivated by their experimental realizability and scientific appeal. For
example, misfit strain is nearly ubiquitous for thin films and is a powerful
tool to control ferroelectric propertiesPhenomenal1 ; ThinFilmStrain1 ;
ThinFilmStrain2 ; ThinFilmStrain3 ; ThinFilmStrain4 . Hydrostatic pressure is
known to have a profound effect on ferroelectric phase transitions. It
destabilizes the ferroelectric phase through weakening long range interactions
relative to short range interactionsLinesGlass  and therefore can be used to
fine tune the Curie point. Uniaxial stress is realizable during piezoelectric
measurements, while biaxial stress occurs in partially relaxed epitaxial thin
films StressThinFilms . Technically, we simulate misfit strain in the range
-5.0% to 5.0% by constraining components of the strain tensor (in Voigt
notations) as follows $\eta_{1}=\eta_{2}=\frac{a-a_{0}}{a_{0}}$ and
$\eta_{6}=0$, where $a_{0}$ is the cubic lattice constant of CsGeBr3. In the
simulated strain interval we consider the points with 0.5% stepsize. For each
value of strain we carry out annealing simulations and harvest the equilibrium
polarization vector and lattice constants. These data are used to analyze the
space group of the structure using ISOTROPY softwareFINDSYM ; FINDSYM2 . Some
representative structures are then subjected to structural relaxation using
DFT as implemented in VASP VASP1 ; VASP2 . In such calculations we keep in-
plane lattice vectors fixed to the same values as in effective Hamiltonian
calculations to maintain the same mechanical boundary conditions, while
allowing the other lattice vectors to relax. This is combined with full ionic
relaxations. For DFT calculations we use the projector-augmented wave (PAW)
pseudopotentialsPAW  with the $r^{2}$SCAN exchange-correlation
functionalr2scan , a 10 $\times$ 10 $\times$ 10 reciprocal space grid (grid
density $\sim$ $0.02~{}\mathring{A}^{-1}$), and an energy cutoff of 800 eV.
For piezoelectric tensor calculations we used the finite difference method
with self-consistent response to finite electric fields as implemented in
VASP.

Similarly, biaxial stress is simulated by constraining components of the
stress tensor as follows $\sigma_{1}=\sigma_{2}=\sigma$ and $\sigma_{6}=0$,
while all other components are allowed to relax. We simulated the values of
$\sigma$ in the range -2.0 GPa to 2.0 GPa with an interval of 0.5 GPa.
Uniaxial stress was modeled by fixing a single component of the stress tensor
to a value from the range -2.0 GPa to 2.0 GPa with an interval of 0.5 GPa.
Likewise, hydrostatic pressure was simulated in the range -2.0 GPa to 2.0 GPa,
with an interval of 0.5 GPa.

##  III Phase diagrams predicted from finite-temperature simulations

Figure 1 gives the temperature evolution of polarization for a few
representative stresses/strains for each type of mechanical load and compares
them to the stress-free case. It reveals that additional phase transitions can
be induced in this material by external load.

For the case of hydrostatic pressure, we find that the positive (negative)
pressure destabilizes (stabilizes) the ferroelectric phase with the Curie
temperature change rate of $-172.3$ K/GPa, or $-0.45T_{C}(0)$ K/GPa. At room
temperature, this yields a structural phase transition at a pressure of 0.94
GPa, which is within the range experimentally measured in Ref. PressureCGBexp,
of 1.0(2) GPa. The material undergoes a single phase transition, just as in
the case of zero pressure. Figure 2(a) shows the associated phase diagram.

Figure 1:  Polarization as a function of temperature for a few representatives
of each type of the mechanical load (a)-(d); Polarization as a function of
hydrostatic pressure (e); uniaxial stress (f); biaxial stress (g); and misfit
strain (h) for a few representative temperatures.

Computational data for uniaxial stress are presented in Fig. 1(b) and reveal
that the material now undergoes two phase transitions, as evident from the
discontinuity in polarization. We find that tensile load increases Curie
temperature, while compressive load decreases it. To the best of our
knowledge, no other perovskite in bulk has shown this type of relationship
before. For example, for bulk BaTiO3, $T_{C}$ increases with both tensile and
compressive stress/strainBTOPhaseDiagram .

The associated phase diagram for uniaxial stress is given in Fig. 2(b). For
tensile uniaxial stress the material undergoes the following sequence of the
phase transitions: paraelectric tetragonal P4/mmm to ferroelectric tetragonal
P4mm to ferroelectric monoclinic Cm. Interestingly, we find that tensile
stress drastically affects the Curie temperature, but the transition
temperature between the two ferroelectric phases (P4mm and Cm) does not
respond to stress. Under compressive uniaxial stress the material undergoes a
phase transition from paraelectric P4/mmm phase to ferroelectric orthorhombic
Amm2 and ferroelectric monoclinic Cm phases. The Curie temperature barely
responds to the compressive stress, while the transition temperature between
the two ferroelectric phases decreases as compressive stress increases in
magnitude. In all cases the supercell remains in a monodomain phase.

Figure 2:  Phase diagrams for CsGeBr3 predicted in computations: phase
transition temperature as a function of hydrostatic pressure (a); uniaxial
stress (b); biaxial stress (c); and misfit strain (d). * indicates phases in
the absence of mechanical load.

Figure 1(c) displays temperature evolution of polarization in the presence of
biaxial stress. We find that biaxial stress induces additional phase
transitions and modifies the Curie temperature of the material. The associated
phase diagram is given in Fig. 2(c). Similar to the case of uniaxial stress,
we find that tensile load increases the Curie temperature, while compressive
load suppresses it, but only very slightly. The phase diagram for biaxial
stress is a mirror image of the one for uniaxial stress, which is rather
intuitive, as compressive uniaxial stress expands the unit cell along the two
perpendicular directions, somewhat mimicking biaxial tensile load.

The last type of the mechanical boundary condition we discuss is the misfit
strain. The temperature evolution of polarization in the presence of misfit
strain is given in Fig. 1(d) and reveals that tensile strain stabilizes the
ferroelectric phase and enhances the polarization, while compressive strain
has an opposite effect. Once again we find that mechanical clamping induces an
additional ferroelectric phase transition. The associated phase diagram is
presented in Fig. 2(d) and is qualitatively similar to the case of biaxial
stress. The main differences are limited to the larger temperature range
associated with stability of orthorhombic and tetragonal phases.

Thus, our simulations predict that most considered mechanical deformations are
capable of inducing tetragonal, orthorhombic and monoclinic phases, none of
which are available under unstressed conditions. The monoclinic phases are
associated with a significant rotation of polarization from the polar axis of
the rhombohedral phase, which could result in enhancement of piezoelectric
response. In particular, we find that the polarization rotates from the
rhombohedral polar axes by 5.64° under 2.0 GPa tensile uniaxial stress and
10.17° under 2.0 GPa compressive uniaxial stress. Under biaxial stress, the
polarization rotates from the rhombohedral polar axes by 6.18° under 2.0 GPa
tensile biaxial stress and by 13.82° under 2.0 GPa compressive biaxial stress.
For misfit strain, the polarization rotates from the rhombohedral polar axes
by 4.13° under 5.0% tensile misfit strain and by 54.73° under 5.0% compressive
misfit strain. It should be noted that all the phases here are single domain.
We did encounter polydomain phases in some of our annealing simulations.
However, they are likely to be metastable as in all investigated cases a
single domain phase can be stabilized using slower annealing rate.

Let us now investigate the origin of the peculiar features of the phase
diagrams presented in Fig 2, namely the decrease in the Curie temperature
under compressive loads. We note that many phase diagrams for oxide
ferroelectrics feature a V-shapeBTOPhaseDiagram ; McCash , where the Curie
temperature is enhanced by both tensile and compressive loads. The effective
Hamiltonian of Eq. 1 can be minimized with respect to its independent
variables following the approach of Ref. McCash, . We focus on uniaxial stress
applied perpendicular to (100) plane, that is $\sigma_{1}$. Note, that the
axes of our Cartesian coordinate system point along [100], [010] and [001]
pseudocubic directions. The local modes associated with the energy minima up
to linear order in stress are

$\displaystyle u_{x}^{2}=u_{0,x}^{2}+a_{x}\sigma_{1}$ (2) $\displaystyle
u_{y}^{2}=u_{0,y}^{2}+a_{y}\sigma_{1}$

where the slopes are

$\begin{aligned}
a_{x}&=\frac{\frac{2\nu_{t}}{3\mu_{t}}(3\alpha^{\prime}+\gamma^{\prime})-\frac{\gamma^{\prime}}{6}\frac{C}{B}}{2\gamma^{\prime}(3\alpha^{\prime}+\gamma^{\prime})}\\\
\end{aligned}\qquad\begin{aligned}
a_{y}&=\frac{-\frac{\nu_{t}}{3\mu_{t}}(3\alpha^{\prime}+\gamma^{\prime})-\frac{\gamma^{\prime}}{6}\frac{C}{B}}{2\gamma^{\prime}(3\alpha^{\prime}+\gamma^{\prime})}\\\
\end{aligned}$ (3)

where

$\begin{aligned} C&=B_{1xx}+2B_{1yy}\\\ \nu_{t}&=\frac{B_{1xx}-B_{1yy}}{2}\\\
\nu_{r}&=B_{5xz}\\\
\alpha^{\prime}&=\alpha-\frac{1}{24}(\frac{C^{2}}{B}+4\frac{\nu_{t}^{2}}{\mu_{t}})\end{aligned}\qquad\begin{aligned}
B&=B_{11}+2B_{12}\\\ \mu_{t}&=\frac{B_{11}-B_{12}}{2}\\\ \mu_{r}&=B_{44}\\\
\gamma^{\prime}&=\gamma+\frac{1}{2}(\frac{\nu_{t}^{2}}{\mu_{t}}-\frac{\nu_{r}^{2}}{\mu_{r}})\end{aligned}\qquad$
(4)

and

$u_{0,x}^{2}=u_{0,y}^{2}=-\frac{\kappa}{6\alpha^{\prime}+2\gamma^{\prime}}$
(5)

Parameters $\kappa,\alpha,\gamma,B_{11},B_{12},B_{44},B_{1xx},B_{1yy}$, and
$B_{5xz}$ are first-principles-based parameters of the effective Hamiltonian
derived from DFT calculations in Ref. PolarAntipolar, . For CsGeBr3$\\\
a_{x}=$ 1260 Bohr5/Ha and $a_{y}=$ 64 Bohr5/Ha, which results in dependencies
given in Fig. 3. The unique feature of these dependencies is the large
difference between the slopes for the two components of the local mode. We
find that the local mode (and therefore polarization) responds very strongly
to the load along the stressed direction, while it barely responds in the
transverse direction. This originates from the interplay between the
$\frac{\nu_{t}}{3\mu_{t}}(3\alpha^{\prime}+\gamma^{\prime})$, which is a
negative quantity, and $\frac{\gamma^{\prime}}{6}\frac{C}{B}$, which is
positive, terms in Eq. 3. These terms compete in $a_{y}$ and cooperate in
$a_{x}$, due to their opposing signs.

Figure 3:  The dependencies of the components of the local mode squared,
$u_{i}^{2}$, on the uniaxial stress obtained for 0 K from the analytical
solution given in Eq. 2 for different perovskite ferroelectrics, as specified
in the titles. Symbols in (a) give the computational data from 10 K
simulations.

This dependence of $u_{x}^{2}$ on stress results in the strong enhancement of
$P_{x}$ under tensile load which leads to the strong enhancement of $T_{C}$.
On the other hand, under compressive load $P_{x}$ is greatly suppressed, while
the $P_{y/z}$ are suppressed only slightly. So the paraelectric to
ferroelectric phase transition is led now by the condensation of the latter
components, which explains why the Curie temperature is slightly suppressed
under compressive load. We note that the analytical predictions from Eqs. 2
closely match results from the simulations at the lowest temperature of 10 K
(see Fig 3(a)).

Our analytical model allows us to compare CsGeBr3 with oxide ferroelectrics
which also undergo phase transitions into rhombohedral phase. Using the
effective Hamiltonian parameters for BaTiO3BTOKNOparams , KNbO3BTOKNOparams
and BiFeO3BFOparams  we obtained the $u_{i}^{2}(\sigma)$ dependencies given in
Fig. 3(b)-(d). Interestingly, for all three of oxide ferroelectrics, the
slopes for $u_{x}^{2}(\sigma)$ and $u_{y/z}^{2}(\sigma)$ have opposite sign,
which results in the V-shape of the upper part of the phase diagram as
previously discussedBTOPhaseDiagram ; McCash . So, the unique feature of
CsGeBr3 is that its $T_{C}$ can be tuned both up and down by stresses and
strains.

##  IV Electromechanical coupling

We now want to find out the isothermal response of polarization to pressure
and stresses. Figure 1(e)-(g) show the dependence of polarization on
pressure/stress for a few representative temperatures. At the lowest
temperature of 150 K, we find from Fig. 1(a) that a pressure of 1.5 GPa
destabilizes the ferroelectric phase. This critical pressure decreases as
temperature increases.

Figure 1(f) shows the isothermal polarization response to uniaxial stress and
reveals its highly nonlinear nature. In the vicinity of the stress-induced
phase transition, the slope of $P(\sigma)$ can reach very large values,
promising outstanding piezoelectric response. Similar trends are found for
biaxial stress [See Fig. 1(g)]. Figure 1(h) shows isothermal dependencies of
polarization on misfit strain, reported for the sake of completeness, although
these might be difficult to realize experimentally.

Figure 4: Polarization components as functions of uniaxial stress computed at
10 K. Solid lines give the linear fit to the data.

To quantify electromechanical coupling in CsGeBr3 we compute the zero field
slope of $P_{i}(\sigma)$ for the case of uniaxial stress using the data
presented in Fig. 4. The 10 K values (in pseudo cubic setting) are reported in
Table 1. The $d_{11}$ coefficient is 21.7 pC/N, which is lower than in oxide
ferroelectrics but significantly higher than in hybrid organic inorganic
perovskites. For example, the longitudinal piezoelectric coefficient d33 of
typical oxide perovskites include 191 pC/N for BaTiO3BTOpiezo , 117 pC/N for
PbTiO3PropOfMat , and 593 pC/N for PZT-5HPropOfMat ; while typical hybrid
organic inorganic perovskite d33 values include 17.1 pC/N for
(BA)2(MA)2Sn3Br10OrgInorg  and 14 pC/N for MDABCO–NH4I3OrgInorg  and
1.18–21.12 pC/N in formatesFormates . Similarly, the isothermal dependencies
of polarization on hydrostatic pressure can be used to compute the hydrostatic
piezoelectric coefficient, $d_{h}=42.51$ pC/N, which is comparable to poled
BaTiO3, $d_{h}=34.0$ pC/NBTOdh .

Table 1: Piezoelectric stress and strain coefficients ($e_{ij}$ and $d_{ij}$) in $C/m^{2}$ and pC/N, respectively, and elastic constants ($C_{ij}$) in kbar computed from effective Hamiltonian ($H_{eff}$) and/or DFT. Coefficients $d_{ij}$ are reported in both standard (hexagonal) setting as well as the pseudocubic setting ($d_{ij}^{pc})$ to facilitate comparison between effective Hamiltonian and DFT predictions.  $e_{15}$ |  $e_{22}$ | $e_{31}$ |  $e_{33}$ | $C_{11}$ |  $C_{12}$ | $C_{13}$ |  $C_{14}$ | $C_{33}$ |  $C_{44}$ | $C_{65}$ |  $C_{66}$
---|---|---|---|---|---|---|---|---|---|---|---
-0.23 |  0.10 | -0.27 |  -0.52 | 222.15 |  61.30 | 74.44 |  31.51 | 155.29 |  94.14 | 31.50 |  80.01
$d_{15}$ |  $d_{22}$ | $d_{31}$ |  $d_{33}$ | $d_{11}^{pc}$ |  $d_{12}^{pc}$ | $d_{14}^{pc}$ |  $d_{15}^{pc}$ |  $d_{11}^{pc}$ (Heff) |  $d_{12}^{pc}$ (Heff) |  $d_{14}^{pc}$ (Heff) |  $d_{15}^{pc}$ (Heff)
-23.39 |  1.93 | -0.81 |  -33.05 | 16.61 |  1.55 | 5.45 |  16.72 | 21.70 |  1.23 | - |  -

##  V Insights from DFT simulations

In this section we take advantage of DFT simulations to both validate our
effective Hamiltonian predictions and reach additional insights into the
electromehanical properties of CsGeBr3. We begin by validating the phases
predicted by the effective Hamiltonian under misfit strain. For that we pick
one point from each unique area of the temperature-misfit strain phase
diagram, Fig. 2(d), and use the associated structure to initialize DFT
simulations. Technically, the average values of local mode and lattice
constants are used to reconstruct atomic structure from effective Hamiltonian
degrees of freedom. We keep the in-plane lattice vectors constrained to the
values reconstructed from the effective Hamiltonian simulations but allow for
full ionic relaxation and relaxation of the remaining lattice vector. All the
phases predicted by the effective Hamiltonian calculations retained their
space groups. Table 2 summarizes these phases, conditions at which they were
obtained, their structural parameters as well as their energies and
enthalpies. The structural files are also provided in Ref. GITHUB, . The
energies of most phases are significantly higher than the ground state (R3m),
which is due to the mechanical constraints and finite temperatures. Enthalpies
are reported with respect to ground state (R3m) for the phases harvested below
the Curie point and with respect to cubic phase (Pm$\bar{3}$m) for the phases
obtained above the Curie point. In all cases the relative enthalpies are
negative, validating the stability of the predicted phases under chosen
conditions.

Table 2: Structural parameters (a, b, c, $\alpha$, $\beta$, and $\gamma$), relative energies ($\Delta E$) and enthalpies ($\Delta H$) computed from DFT for different phases obtained from the temperature-misfit strain phase diagram points ($\eta_{mis}$, $T$). Enthalpies are reported relative to cubic phase (Pm$\bar{3}$m) for phases above TC and relative to the ground state (R3m) for phases below TC. Structures are symmetrized using the FINDSYM software packageFINDSYM ; FINDSYM2 . Space Group |  $T$ (K) | $\eta_{mis}$ |  a (Å) | b (Å) |  c (Å) | $\alpha$ |  $\beta$ | $\gamma$ |  $\Delta$E (meV/f.u.) |  $\Delta$Hcubic (meV/f.u.) |  $\Delta$HR3m (meV/f.u.)
---|---|---|---|---|---|---|---|---|---|---|---
P4/mmm |  600 | -4.0% |  5.32 | 5.32 |  5.64 | 90.0 |  90.0 | 90.0 |  235 | -146 |  -
P4/mmm |  600 | +2.5% |  5.68 | 5.68 |  5.49 | 90.0 |  90.0 | 90.0 |  149 | -33 |  -
P4mm |  170 | -4.0% |  5.32 | 5.32 |  5.79 | 90.0 |  90.0 | 90.0 |  200 | - |  -220
Amm2 |  500 | +5.0% |  5.53 | 8.23 |  8.23 | 90.0 |  90.0 | 90.0 |  58 | -119 |  -
Cm |  10 | -4.0% |  7.52 | 7.52 |  5.79 | 90.0 |  90.1 | 90.0 |  199 | - |  -69
Cm |  10 | +5.0% |  8.23 | 8.23 |  5.64 | 90.0 |  90.4 | 90.0 |  32 | - |  -32

We next validate and augment our effective Hamiltonian predictions of the
piezoelectric response using DFT calculations. Table 1 reports the
piezoelectric stress and strain tensor components computed from DFT. For the
strain coefficients, we report the values in the hexagonal setting, as well as
the pseudocubic setting to facilitate comparison with the effective
Hamiltonian predictions. Note that we cannot convert the effective Hamiltonian
piezoelectric tensor to hexagonal setting, as we cannot resolve all of its
components from uniaxial stress application. We find good agreement between
effective Hamiltonian and DFT predictions for the piezoelectric response. DFT
predicts the largest piezoelectric coefficient is $d_{33}=33.05$ pC/N, which
is comparable to rhombohedral BiFeO3: $d_{33}=27-45$ pC/NBFOpiezo ; BFOpolar2
.

##  VI Conclusions

In summary, we used a combination of finite-temperature effective Hamiltonian
simulations and zero Kelvin DFT simulations to study ferroelectric phases and
phase transitions in the Ge-based ferroelectric halide perovskite CsGeBr3. We
found that hydrostatic pressure has a significant effect on the Curie
temperature and decreases it at the rate $-0.45$ $T_{C}$(0) K/GPa, where
$T_{C}$(0) is the zero pressure Curie temperature. However, the phase
transition sequence remains the same as for the unstressed case, namely
paraelectric Pm$\bar{3}$m to ferroelectric R3m. Application of uniaxial or
biaxial stress, as well as misfit strain, induces a completely different
sequence of phase transitions. Under biaxial tensile stress/strain, the
material undergoes a phase transition from paraelectric tetragonal phase to
ferroelectric orthorhombic phase, followed by another ferroelectric phase
transition to monoclinic phase. Such tensile load dramatically enhances the
Curie temperature, up to 2 times under the largest investigated strain of 5%.
Biaxial compressive stress/strain slightly suppresses the Curie point. The
material undergoes a sequence of phase transitions from paraelectric
tetragonal phase to ferroelectric tetragonal phase to ferroelectric monoclinic
phase. The ferroelectric-ferroelectric phase transition temperature is highly
sensitive to the load. Uniaxial stress has an opposite effect on the phase
transitions and the associated transition temperatures in comparison with the
biaxial stress, which could be justified by the deformations in the unit cell
produced by such loads. This behavior originates from the interplay between
the material parameters that describe electrostriction, elastic properties,
and the strength of the ferroelectric instability. The zero field
piezoelectric constants were found to be in the range 0.8-33.1 pC/N, which is
comparable to BiFeO3BFOpiezo . Strong nonlinearity in the response of
polarization to stress has been predicted, which could result in outstanding
energy converting properties and find applications in energy converting
devices. We believe that our work lays a foundation for phase and
ferroelectricity manipulation in lead-free perovskite materials, which could
lead to the experimental discoveries of these phases and potential use of Ge-
based ferroelectrics for environmentally friendly lead-free energy converting
technologies.

##  VII Acknowledgement

This work was supported by the U.S. Department of Energy, Office of Basic
Energy Sciences, Division of Materials Sciences and Engineering under Grant
No. DE-SC0005245. Computational support was provided by the National Energy
Research Scientific Computing Center (NERSC), a U.S. Department of Energy,
Office of Science User Facility located at Lawrence Berkeley National
Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award
BES-ERCAP-0025236.

## References

  * (1) Z. Yang, H. Ren, X. Chao, and P. Liang, “High permittivity and low dielectric loss of $na_{0.5}bi_{0.5-x}la_{x}cu_{3}ti_{4}o_{12}$ ceramics,” Materials Research Bulletin, vol. 47, no. 5, pp. 1273–1277, 2012.
  * (2) Y. Su, C. Ba, and F. Wang, “Improved dielectric permittivity of nbcto/epoxy composite films with low dielectric loss,” SN applied sciences, vol. 1, no. 10, pp. 1147–, 2019.
  * (3) S.-E. Park and T. R. Shrout, “Ultrahigh strain and piezoelectric behavior in relaxor based ferroelectric single crystals,” Journal of Applied Physics, vol. 82, pp. 1804–1811, 08 1997.
  * (4) W. Li, G. Tang, G. Zhang, H. M. Jafri, J. Zhou, D. Liu, Y. Liu, J. Wang, K. Jin, Y. Hu, H. Gu, Z. Wang, J. Hong, H. Huang, L.-Q. Chen, S. Jiang, and Q. Wang, “Improper molecular ferroelectrics with simultaneous ultrahigh pyroelectricity and figures of merit,” Science Advances, vol. 7, no. 5, p. eabe3068, 2021.
  * (5) P. Ghosez and J. Junquera, “Modeling of ferroelectric oxide perovskites: From first to second principles,” Annual Review of Condensed Matter Physics, vol. 13, no. 1, pp. 325–364, 2022.
  * (6) F. Wang, I. Grinberg, L. Jiang, S. Young, P. Davies, and A. Rappe, “Materials design of visible-light ferroelectric photovoltaics from first principles,” Ferroelectrics, vol. 483, pp. 1–12, 07 2015.
  * (7) P. Machado, M. Scigaj, J. Gazquez, E. Rueda, A. Sánchez-Díaz, I. Fina, M. Gibert-Roca, T. Puig, X. Obradors, M. Campoy-Quiles, and M. Coll, “Band gap tuning of solution-processed ferroelectric perovskite bife1–xcoxo3 thin films,” Chemistry of materials, vol. 31, no. 3, pp. 947–954, 2019.
  * (8) A. M. Ganose, D. O. Scanlon, A. Walsh, and R. L. Z. Hoye, “The defect challenge of wide-bandgap semiconductors for photovoltaics and beyond,” Nature communications, vol. 13, no. 1, pp. 4715–4715, 2022.
  * (9) L.-K. Gao, Y.-L. Tang, and X.-F. Diao, “First-principles study on the photoelectric properties of csgei3 under hydrostatic pressure,” Applied sciences, vol. 10, no. 15, pp. 5055–, 2020.
  * (10) U. Schwarz, F. Wagner, K. Syassen, and H. Hillebrecht, “Effect of pressure on the optical-absorption edges of csgebr 3 and csgecl 3,” Physical review. B, Condensed matter, vol. 53, no. 19, pp. 12545–12548, 1996.
  * (11) L. C. Tang, J. Y. Huang, C. S. Chang, M. H. Lee, and L. Q. Liu, “New infrared nonlinear optical crystal csgebr3: synthesis, structure and powder second-harmonic generation properties,” Journal of physics. Condensed matter, vol. 17, no. 46, pp. 7275–7286, 2005.
  * (12) T. Krishnamoorthy, H. Ding, C. Yan, W. L. Leong, T. Baikie, Z. Zhang, M. Sherburne, S. Li, M. Asta, N. Mathews, and S. G. Mhaisalkar, “Lead-free germanium iodide perovskite materials for photovoltaic applications,” Journal of materials chemistry. A, Materials for energy and sustainability, vol. 3, no. 47, pp. 23829–23832, 2015.
  * (13) N. Thi Han, V. Khuong Dien, and M.-F. Lin, “Electronic and optical properties of csgex3 (x= cl, br, and i) compounds,” ACS Omega, vol. 7, no. 29, pp. 25210–25218, 2022.
  * (14) Y. Zhang, E. Parsonnet, A. Fernandez, S. M. Griffin, H. Huyan, C.-K. Lin, T. Lei, J. Jin, E. S. Barnard, A. Raja, P. Behera, X. Pan, R. Ramesh, and P. Yang, “Ferroelectricity in a semiconducting all-inorganic halide perovskite,” Science advances, vol. 8, no. 6, pp. eabj5881–eabj5881, 2022\.
  * (15) G. Thiele, H. W. Rotter, and K. Schmidt, “Kristallstrukturen und phasentransformationen von caesiumtrihalogenogermanaten(ii) csgex3 (x = cl, br, i),” Zeitschrift für anorganische und allgemeine Chemie, vol. 545, pp. 148–156, 1987.
  * (16) R. Kashikar, S. Lisenkov, and I. Ponomareva, “Coexistence of polar and antipolar phases in ferroelectric halide perovskite ${\mathrm{csgebr}}_{3}$,” Phys. Rev. B, vol. 109, p. L020101, Jan 2024\.
  * (17) N. A. Pertsev, A. G. Zembilgotov, and A. K. Tagantsev, “Effect of mechanical boundary conditions on phase diagrams of epitaxial ferroelectric thin films,” Physical review letters, vol. 80, no. 9, pp. 1988–1991, 2006.
  * (18) U. Schwarz, H. Hillebrecht, M. Kaupp, K. Syassen, H.-G. von Schnering, and G. Thiele, “Pressure-induced phase transition in csgebr3 studied by x-ray diffraction and raman spectroscopy,” Journal of Solid State Chemistry, vol. 118, no. 1, pp. 20–27, 1995.
  * (19) D.-K. Seo, N. Gupta, M.-H. Whangbo, H. Hillebrecht, and G. Thiele, “Pressure-induced changes in the structure and band gap of csgex3 (x = cl, br) studied by electronic band structure calculations,” Inorganic chemistry, vol. 37, no. 3, pp. 407–410, 1998.
  * (20) Y. Qin, Z. Lv, S. Chen, W. Li, X. Wu, L. Ye, N. Li, and P. Lu, “Tuning pressure-induced phase transitions, amorphization, and excitonic emissions of 2d hybrid perovskites via varying organic amine cations,” Journal of physical chemistry. C, vol. 123, no. 36, pp. 22491–22498, 2019.
  * (21) S. Dhole, A. Chen, W. Nie, B. Park, and Q. Jia, “Strain engineering: A pathway for tunable functionalities of perovskite metal oxide films,” Nanomaterials (Basel, Switzerland), vol. 12, no. 5, pp. 835–, 2022.
  * (22) C. Grote and R. F. Berger, “Strain tuning of tin–halide and lead–halide perovskites: A first-principles atomic and electronic structure study,” Journal of physical chemistry. C, vol. 119, no. 40, pp. 22832–22837, 2015.
  * (23) Z. Chen, X. Zou, W. Ren, L. You, C. Huang, Y. Yang, P. Yang, J. Wang, T. Sritharan, L. Bellaiche, and L. Chen, “Study of strain effect on in-plane polarization in epitaxial bifeo3 thin films using planar electrodes,” Phys. Rev. B, vol. 86, p. 235125, Dec 2012.
  * (24) Y. Su and C. M. Landis, “Continuum thermodynamics of ferroelectric domain evolution: Theory, finite element implementation, and application to domain wall pinning,” Journal of the Mechanics and Physics of Solids, vol. 55, no. 2, pp. 280–305, 2007.
  * (25) W. Zhong, D. Vanderbilt, and K. Rabe, “First-principles theory of ferroelectric phase transitions for perovskites: The case of batio3,” Physical review. B, Condensed matter, vol. 52, no. 9, pp. 6301–6312, 1995.
  * (26) W. Zhong, D. Vanderbilt, and K. Rabe, “Phase transitions in batio3 from first principles,” Physical review letters, vol. 73, no. 13, pp. 1861–1864, 1994\.
  * (27) S. Liu, I. Grinberg, and A. M. Rappe, “Development of a bond-valence based interatomic potential for bifeo3 for accurate molecular dynamics simulations,” Journal of physics. Condensed matter, vol. 25, no. 10, pp. 102202–102202, 2013.
  * (28) I. Ponomareva and L. Bellaiche, “Influence of the growth direction on properties of ferroelectric ultrathin films,” Phys. Rev. B, vol. 74, p. 064102, Aug 2006.
  * (29) I. I. Naumov, L. Bellaiche, and H. Fu, “Unusual phase transitions in ferroelectric nanodisks and nanorods,” Nature, vol. 432, no. 7018, pp. 737–740, 2004.
  * (30) Y. Nahas, S. Prokhorenko, Q. Zhang, V. Govinden, N. Valanoor, and L. Bellaiche, “Topology and control of self-assembled domain patterns in low-dimensional ferroelectrics,” Nature Communications, vol. 11, no. 1, p. 5779, 2020.
  * (31) K. Rabe and U. Waghmare, “Ferroelectric phase transitions from first principles,” Journal of Physics and Chemistry of Solids, vol. 57, no. 10, pp. 1397–1403, 1996.  Proceeding of the 3rd Williamsburg Workshop on Fundamental Experiments on Ferroelectrics.
  * (32) Y. Yang, H. Xiang, and L. Bellaiche, “Revisiting structural phases in some perovskites: The case of ba ce o 3,” Physical review. B, vol. 104, no. 17, 2021.
  * (33) R. Oja, K. Johnston, J. Frantti, and R. M. Nieminen, “Computational study of (111) epitaxially strained ferroelectric perovskites ${\text{batio}}_{3}$ and ${\text{pbtio}}_{3}$,” Phys. Rev. B, vol. 78, p. 094102, Sep 2008.
  * (34) K. McCash, B. K. Mani, C.-M. Chang, and I. Ponomareva, “The role of mechanical boundary conditions in the soft mode dynamics of pbtio3,” Journal of Physics: Condensed Matter, vol. 26, p. 435901, oct 2014.
  * (35) R. Herchig and I. Ponomareva, “Unusual soft mode dynamics in ferroelectric pbtio3 nanowire under different mechanical boundary conditions,” Journal of applied physics, vol. 122, no. 21, 2017.
  * (36) D. Wickramaratne and J. L. Lyons, “Shortcomings of using the scan functional for point defects and polarons in semiconductors,” arXiv preprint arXiv:2311.03634, 2023.
  * (37) M. W. Swift and J. L. Lyons, “Lone-pair stereochemistry induces ferroelectric distortion and the rashba effect in inorganic halide perovskites,” Chemistry of Materials, vol. 35, no. 21, pp. 9370–9377, 2023.
  * (38) “Appendix e - integration schemes,” in Understanding Molecular Simulation (Second Edition) (D. Frenkel and B. Smit, eds.), pp. 533–544, San Diego: Academic Press, second edition ed., 2002.
  * (39) D. C. Rapaport, The Art of Molecular Dynamics Simulation.  Cambridge University Press, 2 ed., 2004.
  * (40) C. Ederer and N. A. Spaldin, “Effect of epitaxial strain on the spontaneous polarization of thin film ferroelectrics,” Phys. Rev. Lett., vol. 95, p. 257601, Dec 2005.
  * (41) H. Meley, Karandeep, L. Oberson, J. de Bruijckere, D. T. L. Alexander, J.-M. Triscone, P. Ghosez, and S. Gariglio, “Structural analysis of LaVO3 thin films under epitaxial strain,” APL Materials, vol. 6, p. 046102, 04 2018\.
  * (42) H. N. Lee, S. M. Nakhmanson, M. F. Chisholm, H. M. Christen, K. M. Rabe, and D. Vanderbilt, “Suppressed dependence of polarization on epitaxial strain in highly polar ferroelectrics,” Phys. Rev. Lett., vol. 98, p. 217602, May 2007.
  * (43) D. Sando, “Strain and orientation engineering in abo3 perovskite oxide thin films,” Journal of Physics: Condensed Matter, vol. 34, p. 153001, feb 2022\.
  * (44) M. E. Lines, Principles and applications of ferroelectrics and related materials.  Oxford classic texts in the physical sciences Principles and applications of ferroelectrics and related materials, Oxford: Clarendon, 2001\.
  * (45) K. N. K.-N. Tu, Electronic thin film science : for electrical engineers and materials scientists.  New York: Macmillan, 1992.
  * (46) H. T. Stokes and D. M. Hatch, “Findsym: program for identifying the space-group symmetry of a crystal,” Journal of applied crystallography, vol. 38, no. 1, pp. 237–238, 2005.
  * (47) D. M. H. H. T. Stokes and B. J. Campbell, “Findsym, isotropy software suite.” https://iso.byu.edu.
  * (48) G. Kresse and J. Furthmüller, “Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set,” Phys. Rev. B, vol. 54, pp. 11169–11186, Oct 1996.
  * (49) G. Kresse and J. Furthmüller, “Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set,” Computational Materials Science, vol. 6, no. 1, pp. 15–50, 1996.
  * (50) P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B, vol. 50, pp. 17953–17979, Dec 1994.
  * (51) J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, “Accurate and numerically efficient r2scan meta-generalized gradient approximation,” The Journal of Physical Chemistry Letters, vol. 11, no. 19, pp. 8208–8215, 2020\.  PMID: 32876454.
  * (52) U. SCHWARZ, H. HILLEBRECHT, M. KAUPP, K. SYASSEN, H.-G. VON SCHNERING, and G. THIELE, “Pressure-induced phase transition in csgebr3 studied by x-ray diffraction and raman spectroscopy,” Journal of solid state chemistry, vol. 118, no. 1, pp. 20–27, 1995.
  * (53) A. Grunebohm, M. Marathe, and C. Ederer, “Ab initio phase diagram of batio3 under epitaxial strain revisited,” Applied Physics Letters, vol. 107, p. 102901, 2015.
  * (54) R. Herchig, C. M. Chang, B. K. Mani, and I. Ponomareva, “Electrocaloric effect in ferroelectric nanowires from atomistic simulations,” Scientific Reports, vol. 5, no. 1, p. 17294, 2015.
  * (55) S. L. C.-M. Chang, B. K. Mani and I. Ponomareva, “Prediction of electromagnons in antiferromagnetic ferroelectrics from first-principles: The case of bifeo3,” Ferroelectrics, vol. 494, no. 1, pp. 68–75, 2016.
  * (56) J. Gao, D. Xue, W. Liu, C. Zhou, and X. Ren, “Recent progress on batio3-based piezoelectric ceramics for actuator applications,” Actuators, vol. 6, p. 24, 07 2017.
  * (57) R. E. Newnham, “Stress and strain,” in Properties of Materials: Anisotropy, Symmetry, Structure, Oxford University Press, 11 2004.
  * (58) T. Zhang, K. Xu, J. Li, L. He, D.-W. Fu, Q. Ye, and R.-G. Xiong, “Ferroelectric hybrid organic–inorganic perovskites and their structural and functional diversity,” National Science Review, vol. 10, p. nwac240, 11 2022.
  * (59) A. Popoola, P. S. Ghosh, M. Kingsland, R. Kashikar, D. DeTellem, Y. Xu, S. Ma, S. Witanachchi, S. Lisenkov, and I. Ponomareva, “First-principles property assessment of hybrid formate perovskites,” The Journal of chemical physics, vol. 159, no. 7, 2023.
  * (60) W. P. . Mason, Physical acoustics : principles and methods : Vol. 1, (part A).  New York and London: Academic Press New York and London, 1964-.
  * (61) “Provided structure files.” https://github.com/USFmatscilab/CGB-structures.git.
  * (62) G. Yuan, S. Or, Y. Wang, Z. Liu, and J. Liu, “Preparation and multi-properties of insulated single-phase bifeo3 ceramics,” Solid State Communications, vol. 138, no. 2, pp. 76–81, 2006.
  * (63) L. Liu, T. Rojac, J. Kimpton, J. Walker, M. Makarovic, J.-F. Li, and J. Daniels, “Poling-induced inverse time-dependent microstrain mechanisms and post-poling relaxation in bismuth ferrite,” Applied Physics Letters, vol. 116, p. 122901, 03 2020.