# Electromechanical enhancement of live jellyfish for ocean exploration

Simon R. Anuszczyk1 and John O. Dabiri1,2∗

1Graduate Aerospace Laboratories, California Institute of Technology,
Pasadena, CA, USA.
2Mechanical and Civil Engineering, California Institute of Technology,
Pasadena, CA, USA.

∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
> The vast majority of the ocean’s volume remains unexplored, in part because
> of limitations on the vertical range and measurement duration of existing
> robotic platforms. In light of the accelerating rate of climate change
> impacts on the physics and biogeochemistry of the ocean, the need for new
> tools that can measure more of the ocean on faster timescales is becoming
> pressing. Robotic platforms inspired or enabled by aquatic organisms have
> the potential to augment conventional technologies for ocean exploration.
> Recent work has demonstrated the feasibility of externally stimulating live
> jellyfish swimming via implanted microelectronics. We present a robotic
> platform that leverages this external electrical swimming control, while
> also using a passive mechanical attachment to streamline the jellyfish shape
> and significantly enhance the payload capacity of the biohybrid robotic
> jellyfish. A six-meter-tall, 13,600-liter saltwater facility was constructed
> to enable testing of the vertical swimming capabilities of the biohybrid
> robotic jellyfish over distances exceeding 35 body diameters. We found that
> the combination of external swimming control and the addition of the
> mechanical forebody resulted in an increase in swimming speeds to 4.5 times
> natural jellyfish locomotion. Moreover, the biohybrid jellyfish were capable
> of carrying a payload volume up to 105% of the jellyfish body volume. The
> added payload decreased the intracycle acceleration of the biohybrid robots
> relative to natural jellyfish, which could also facilitate more precise
> measurements by onboard sensors that depend on consistent platform motion.
> The demonstrated performance of these biohybrid robots suggests an
> opportunity to expand the set of robotic tools for ocean exploration.

One sentence summary: A biohybrid robotic platform that leverages electrical
and mechanical enhancement of live jellyfish increases swimming speeds to 4.5
times natural locomotion, while enabling a payload larger than the animal
itself.

## INTRODUCTION

The ocean contains the majority of the habitable volume on earth, and many of
its habitats and species remain largely unexplored (?). The acceleration of
climate change is leading to increasing ocean warming, acidification, and
deoxygenation and posing a danger to biodiversity and associated ocean
resources (?). Recognizing these impacts, the United Nations has declared this
the ”United Nations Decade of Ocean Science for Sustainable Development” in an
effort to galvanize ocean exploration (?). Less than 80% of the ocean has been
explored, with basic properties such as depth known for only 18% of the ocean
at a resolution of 1 km (?). These knowledge gaps motivate more comprehensive
ocean exploration at timescales relevant to climate change (?).

A variety of robotic technologies for ocean exploration have resulted in
significant advances in ocean science (?, ?). For example, remote sensing
platforms such as uncrewed aerial vehicles and satellites have been used to
produce high resolution maps of sensitive marine habitats (?) and to study
important physical characteristics such as salinity (?) and temperature (?).
However, these technologies are limited in the achievable depth of
observations by water turbidity and solar illumination (?). Deep sea tools
including gliders (?), seafloor mounted buoys (?), and autonomous underwater
vehicles (?) can explore to significant depths but with mission duration often
limited by battery capacity. The sheer size of the ocean, comprising more than
1.3 billion cubic kilometers (?), challenges the cost-effective scaling of
these technologies for comprehensive ocean coverage (?).

Robotic vehicles wholly comprised of mechanical components have increasingly
been inspired by the demonstrated performance of aquatic animals, with the
goal of improving robotic vehicle propulsive efficiency and range. For
example, an acoustically controlled soft robotic fish was developed to observe
aquatic animal behaviour with minimal disruption (?). However, this robot
required an accompanying human diver, thereby limiting operating depths to
regions that are safe for human diving. A benthic walking robot demonstrated
locomotion along the seafloor and geological sediment-sampling using a shovel
(?). While valuable for benthic exploration, this and similar robots are
limited to the seafloor due to their inability to swim. A soft snailfish-
inspired robot demonstrated swimming locomotion at a depth of 10,900 m (?).
This work presented a novel method for electronics hardening by distributing
the electronics in a polymeric matrix, but power consumption limited swimming
to 45 minutes. Each of the aforementioned examples demonstrates the emerging
capabilities of robotic mimics of aquatic organisms; however, the performance
envelope of live aquatic animals in the ocean has not yet been realized by
these approaches.

A different strategy, biologging, embeds sensor tags on live animals such as
pinnipeds (?), sea turtles (?), and jellyfish (?, ?) and records measurements
of the aquatic environment as the animals locomote naturally through the
ocean. Biologgers can be used to study ocean characteristics such as salinity,
temperature, and oxygen distribution using conductivity-temperature-depth
instruments, oxygen probes, hydrophones, and other sensors (?). However, the
timing and locations of the measurements are uncontrolled and limited to the
path of the tagged animal.

Biohybrid robots have the potential to combine the control of fully mechanical
robotic platforms and the swimming performance envelope of live aquatic
animals to achieve the goal of comprehensive ocean exploration. Proof-of-
concept of biohyrid robotic control has been achieved in a variety of
organisms, including cockroaches (?), hawkmoths (?), and tissue engineered
rays (?) and jellyfish (?). By utilizing biological locomotion, biohybrid
robotics can potentially capitalize on the efficiency and adaptability of
existing biological organisms.

Here, we investigate jellyfish as a platform for biohybrid robotic sensors
that can access the full ocean, including its deepest regions. Jellyfish are
naturally occurring in a diverse range of oceanic environments, including
tropical (?) and polar regions (?), and spanning from the surface to depths
exceeding 10,800 meters (?). They are also tolerant to hypoxic conditions
below 1 mg $O_{2}$ per liter (?), in contrast to fish who avoid levels below 2
to 3 mg $O_{2}$ per liter. In terms of their locomotion, jellyfish are the
most energy efficient of all known metazoans, with a cost of transport
3.5-fold lower than even efficiently swimming fish such as salmon (?).

Jellyfish have previously been used as biologgers to investigate their ecology
and behaviour (?, ?). Previous work has also demonstrated that the swimming of
live jellyfish can be controlled by external electrical stimulation from an
implanted microelectronic device (?). When swimming contractions were
stimulated at approximately 0.6 Hz, the jellyfish were observed to swim up to
2.8 times baseline swimming speeds without the swim controller. These
increased speeds required only a twofold increase in energy expenditure,
implying a concomitant increase in propulsive efficiency during external
swimming control.

In consideration of the ethics of the aforementioned robotic manipulations,
the authors of that work collaborated with bioethicists to establish a
framework for responsible implementation of these techniques (?). It is
important to note that the jellyfish genus that is the focus of this work,
Aurelia, lacks a central nervous system and nociceptors, and therefore lack
the capacity to sense pain (?). While they have a known mucosal stress
response, this was not observed to occur in any prior experiments using
external swimming control.

Although previous work demonstrated the feasibility of external control of
jellyfish swimming frequency and quantified associated performance gains, key
gaps remain in development of the biohybrid jellyfish robot as a research
platform for ocean science. First, the prior work did not explicitly
incorporate a payload capacity, which is essential for the carrying of onboard
measurement sensors. Second, assessment of the swimming performance was
limited to tests either in a 2-meter-tall laboratory water tank (?) or a
coastal marine environment (?). Swimming performance over longer distances
consistent with the envisioned deep ocean measurements has not been
characterized. Finally, while natural jellyfish boast a low cost of transport
(?), this performance is achieved despite a body shape that is not
streamlined.

We investigate for the first time the feasibility of simultaneous electrical
and mechanical modification of live jellyfish as biohybrid robotic platforms.
The mechanical modification comprises the addition of a passive forebody to
streamline the natural jellyfish body shape while also facilitating a
significant payload capacity. The performance of the electromechanically
modified jellyfish is tested in a specially constructed six-meter-tall,
13,600-liter saltwater facility that enables the jellyfish to swim distances
up to 35 body diameters. These measurements, combined with a computational
model of the swimming dynamics, enable optimization of the design of the
biohybrid jellyfish robots.

## RESULTS

We modified the shape of the exumbrellar surface of the natural jellyfish by
attaching additively manufactured hemi-ellipsoids of various aspect ratios to
the apex of the animals (see Methods). As in prior work (?), the jellyfish
swim muscles, which are located on the subumbrellar surface, were stimulated
with two electrodes to produce regular contractions at 0.5 Hz with a 3.7-volt,
10-ms square wave. The biohybrid robot is ballasted to create a slight
positive buoyancy and a horizontal restoring torque constraining the animals
to swimming in the vertical direction. This positive ballasting ensured that
the observed downward swimming could be attributed only to the active animal
swimming and not to the passive dynamics of the attached mechanical
forebodies.

Figure 1. Comparison of flow streamlining with and without passive mechanical
forebody.
Jellyfish biohybrid robot with and without attached hemi-ellipsoid forebody.
Cyan arrows indicate downward swimming direction of the animals. Yellow arrows
highlight areas of flow diversion and demonstrate flow disturbance before (A)
and after (B) addition of streamlined passive mechanical forebody (outlined in
red dashed curve). White scale bar in each panel indicates 5 cm. In panel (B),
$H$ is the forebody height and $D$ is the maximum forebody diameter. In both
panels, the flow visualization was created by superimposing several
consecutive images, each shifted spatially so that the jellyfish appears
stationary while particles trace pathlines around the animal.

### Small-scale flow measurements

We first studied the hydrodynamics of the biohybrid robot in a 1.2-m tall, 300
liter saltwater tank to investigate the hydrodynamic impacts of the mechanical
forebodies on the local flow around the jellyfish. Time-lapse images of 13
micron, neutrally buoyant tracer particles showed that the addition of hemi-
ellipsoids of increasing aspect ratios streamline the flow around the
biohybrid robots (Fig. 1; see also Supplementary Movies 1 and 2).
Specifically, comparison of the tracer particle pathlines illustrated the
reduction in flow re-direction at the apex of the biohybrid robots relative to
the natural jellyfish (yellow arrows). While the total lateral diversion of
flow is the same in the two cases, i.e., directing flow around the maximum
body width at the bell margin, the passive forebody achieves this diversion in
a more gradual fashion. It was anticipated that the larger lateral diversion
of flow at the apex of the natural jellyfish shape would lead to greater drag
than on the streamlined forebodies.

Figure 2. Example swimming speed measurements of biohybrid robot.
(A). 30-second swimming speed trace (i.e., swimming distance of 6.2 body
diameters) shows contraction and relaxation cycles of biohybrid robot. (B).
Frequency spectrum analysis confirms animal response to external stimulation
at 0.5 Hz. Semi-transparent band around data curve in panel (A) represents a
tracking error of $\pm 5$ pixels.

### Long-distance swimming performance

To test the impact of the streamlining on swimming performance, the natural
jellyfish and biohybrid robots were tested in a new 6-m tall, 13,600 liter
saltwater facility built for free swimming experiments (see Methods for
facility description). Figure 2A shows a 30-second swimming speed trace (i.e.,
swimming distance of 6.2 body diameters) of the swimming biohybrid robot
showing speed peaks during the jellyfish bell contraction phase and troughs
during bell relaxation. The swimming speed oscillates with each contraction
and relaxation, and it eventually reaches a steady state average swimming
speed. A spectral analysis of the swimming speed trace shows that the animal
responded to external electrical stimulation by contracting at the 0.5 Hz
stimulation frequency (Fig. 2B).

Figure 3. Measurements of swimming speed in six-meter vertical saltwater
facility.
Sample biohybrid robot trajectories at 6.67-s intervals for the biohybrid
robot without mechanical forebody (A) and with mechanical forebody $AR_{3}$
(B). Larger spacing between robot images in panel (B) qualitatively
illustrates faster swimming speed. Zoomed-in image shows middle one-third of
height of facility. Orange dropline is marked every 1.5 m. White scale bar
indicates 20 cm. (C) Normalized biohybrid robot speed increases with
streamlining coefficient (S) demonstrating impact of streamlined forebodies.
Each color represents one animal with diameters listed in the legend;
transparent band represents $\pm$ 1 standard deviation, calculated based on
the three trials at each point. (D) Enhancement factor compares swimming
performance of biohybrid robot to same stimulated jellyfish without mechanical
forebody. Maximum enhancement factor corresponds to 61% increase in swimming
speeds or a combined increase from electrical and mechanical enhancement to
4.5 times typical unstimulated jellyfish locomotion.

We tested the stimulated jellyfish biohybrid robot without a forebody and with
the addition of three different hemi-ellipsoid forebodies of increasing aspect
ratio, $AR=H/D$, where $H$ is the forebody height and $D$ is the maximum
forebody diameter: $AR_{0}=0$ (i.e., no forebody), $AR_{1}=0.31$,
$AR_{2}=0.66$, and $AR_{3}=1.0$. A total of $N=3$ different animals were
tested with bell diameters $d_{1}=16.0\pm 0.5$ cm, $d_{2}=16.5\pm 0.5$ cm, and
$d_{3}=18.0\pm 0.5$ cm. Each forebody configuration was tested 3 times in a
randomized order to ensure that trends in forebody performance were not
affected by any drift in animal performance over time. Sample biohybrid robot
trajectories with and without the addition of a hemi-ellipsoid forebody
$AR_{3}$ are shown in Figure 3A and 3B, respectively (see also Supplementary
Movies 3 and 4). The images in Figure 3 are equally spaced at intervals of
6.67 seconds and show a zoomed-in view of the middle one-third of the total
height of the facility. An orange scaling dropline is marked every 1.5 meters
to provide a reference length scale. The spatial intervals between images of
the animal in Figure 3B are larger than those in Figure 3A, qualitatively
illustrating the faster swimming speed of the biohybrid robot with the
addition of the streamlined mechanical forebody.

The biohybrid robot swimming speed was found to increase with the aspect ratio
of the forebody, with the largest increase occurring between $AR_{0}$ (i.e.,
no forebody) and $AR_{1}$. Figure 3C plots the normalized biohybrid robot
swimming speed (in body diameters per second) versus streamlining coefficient,
$S=-log(C_{D})$, where $C_{D}$ is the forebody drag coefficient measured via
terminal velocity drop tests (see Methods). Swimming speeds increased in
direct proportion with the streamlining coefficient, reaching a plateau for
the most streamlined forebodies. Diameter-normalized swimming speeds were
observed to be larger for the larger jellyfish, likely reflecting a greater
propulsive capacity to translate the mechanical forebodies downward against
their positive buoyancy.

We also compared the swimming performance of each biohybrid robot to the same
jellyfish without a mechanical forebody (Fig. 3D). The maximum enhancement
factor (i.e., performance ratio relative to swimming without a forebody),
which occurred for the mechanical forebody with the largest streamlining
coefficient, i.e., $AR_{3}$, corresponded to an increase in swimming speeds of
up to 61%. Across all cases, we observed an average increase in swimming speed
of $38\%\pm 22\%$. In light of prior work demonstrating an increase of up to
2.8 times the swimming speed of natural jellyfish, the present results suggest
that the combined electromechanical enhancement could enable an overall
increase in swimming speeds to more than 4.5 times natural, unstimulated
jellyfish locomotion.

Figure 4. Intracycle swimming speed range and forebody inertia versus
streamlining coefficient.
(A) Intracycle swimming speed range decreases with streamlining coefficient
$S$ due to increased inertia of the hemi-ellipsoid forebody. Error bars
indicate $\pm$ 1 standard deviation of intracycle speed throughout each trial.
A small, random variability is added to the x-axis to help distinguish each of
the data points. (B) Forebody mass increases with streamlining coefficient
$S$. The largest tested forebody weighed $m_{4}=271.00\pm.05$ grams in air.

The addition of mechanical forebodies adds inertia to the biohybrid robotic
system proportional to the volume of displaced water. The largest forebody
corresponds to an increase in the mass by $m_{4}=271.04\pm.05$ grams of
payload and a volume of $105\%\pm 1\%$ of the animal body volume. Figure 4A
plots the intracycle swimming speed range in body diameters per second vs $S$
of all 24 trials of the biohybrid robot with the three forebodies. With
increasing streamlining coefficient $S$, the intracycle swimming speed range
of the biohybrid robot decreases. Figure 4B illustrates that the forebody mass
increases in direction proportion to the streamlining coefficient, which
explains the concurrent reduction in intracycle swimming speed range.
Specifically, because the biohybrid robot swims in a highly unsteady motion as
shown in Figure 2A, the additional inertia changes the swimming dynamics.

Figure 5. Validation of theoretical model of biohybrid robot performance.
(A) Model prediction (blue) and empirical measurement (orange) of biohybrid
robot swimmming speed for $AR_{3}$ forebody over 5 body contraction cycles.
(B) Model prediction (blue) and empirical measurements (orange) of decrease in
intracycle swimming speed range with increased forebody mass. Error bars
indicate $\pm$ 1 standard deviation of the 3 trials at each point. (C) Model
prediction (blue) and emprical measurements (orange) of average swimming
speeds over all animals and all trials for each forebody configuration. Error
bars indicate $\pm$ 1 standard deviation of all trials.

We developed a theoretical model of the biohybrid robot swimming dynamics to
further investigate the dependence of the swimming performance on the forebody
geometry and properties as well as the animal swimming behavior (see Methods).
Figure 5A compares an example predicted swimming speed trace (blue) and
measurements from swimming tests using the mechanical forebody $AR_{3}$ over 5
body contraction cycles. The model successfully captures the trends in both
swimming speed and intracycle speed range versus streamlining factor. While
the model predicts a decrease in the intracycle swimming speed range with
increasing forebody mass, it underestimates the quantitative decrease observed
between the $AR_{2}$ and $AR_{3}$ cases. We speculate that this might be due a
change in the effective coefficient of drag of the combined biohybrid robotic
system with the $AR_{3}$ forebody, given its smaller absolute drag
coefficient.

Pooling the measurements from all 3 animals tested in these experiments, we
observed a plateau in the average speed between the $AR_{2}$ and $AR_{3}$
cases. This plateau is qualitatively consistent with a scaled version of the
terminal velocity equation, i.e., it indicates that performance at high
swimming speeds is likely limited by the terminal velocity of the mechanical
forebody. In sum, the model successfully captures the biohybrid robot dynamics
and may there be used to extrapolate to larger jellyfish and different
forebody geometries in the future.

## DISCUSSION

We have presented the electromechanical enhancement of a biohybrid jellyfish
robot for ocean exploration through the addition of streamlined, passive
forebodies. We demonstrated the ability to travel at speeds up to 4.5 times
baseline jellyfish swimming speeds while carrying a payload of more than 270
grams in a volume of 105% of the jellyfish volume. The flow around the
biohybrid robot was examined using time-lapse pathline visualization to
demonstrate local streamlining effects, and biohybrid robot swimming
performance was investigated in a custom 6-m tall saltwater facility that
enabled testing of swimming over distances significantly longer than prior
studies. The addition of a payload was found to decrease intracycle
acceleration, which may facilitate more precise onboard measurements in cases
where unsteady platform motion should be minimized. By utilizing live animals
and their ability to extract chemical energy from the water column via
feeding, the presentation approach has the potential to eliminate most of the
battery power typically required for underwater propulsion, instead leveraging
and controlling the jellyfish locomotion This suggests an opportunity to
expand the set of biohybrid robotic tools for ocean exploration.

While this work has focused on maximizing the swimming speed of the biohybrid
jellyfish robots, it has not considered the associated metabolic costs of this
increased performance. Previous work found that increased swimming performance
by biohybrid jellyfish was not penalized by proportional increases in energy
expenditure; in fact, the externally stimulated biohybrid jellyfish swam more
efficiently than their natural counterparts (?). However, energy expenditure
experiments have not been conducted for freely swimming animals. Future work
will aim to address this knowledge gap using a combination of flow
measurements and respirometry. Longer-term experiments on timescales relevant
to ocean exploration (i.e., weeks to months in duration) will be necessary to
confirm that the performance observed here can be sustained while continuing
to ensure the long-term well-being of the jellyfish. Finally, these
experiments were constrained to one-dimensional swimming to focus on vertical
profiling applications in the ocean. Future work will extend the
maneuverability of the biohybrid robotic jellyfish to enable three-dimensional
exploration.

Figure 6. Utilizing biohybrid robotic jellyfish payload capacity for ocean
science.
Demonstrated payload mass and volume capacity (i.e., blue boundaries) can be
utilized for a variety of ocean science and exploration applications. COTS
wildlife trackers with a variety of sensors are shown in red; power sources
with approximate duration powering current electronics in yellow; CTD sensors
shown in purple; and chemical probes in green. Filling the forebody with a
lithium metal battery could power the biohybrid robot for more than 22 weeks.

The demonstrated ability to carry payloads could be utilized to investigate a
variety of ocean science questions. Figure 6 explores some of the ways in
which the mass and volume of the largest $AR_{3}$ hemi-ellipsoid forebody
could be used to further ocean exploration. Several commercial off-the-shelf
(COTS) wildlife trackers are shown in red, which include different
combinations of sensors for parameters such as temperature, depth, light, and
accelerometers for motion tracking (ORI400, Little Leonardo, Tokyo, Japan;
MINIPAT, Wildlife Computers Inc, WA, USA; Daily Diary, Wildbyte Technologies
Ltd, Swansea, United Kingdon). The forebody could instead be utilized as a
large battery (power sources shown in orange) which could power the biohybrid
robot for up to 22 weeks using a disposable 101.6 watt-hour Lithium-metal
battery (?) or from 11 to 22 weeks using other battery types (?, ?). The
relatively high energy density of lithium batteries could be exploited to
preserve additional volume for less dense scientific payloads. For example, an
onboard conductivity-temperature-depth (CTD) sensor, such as those shown in
purple, could provide salinity and temperature measurements to inform ocean
science models. Several COTS CTDs (LeveLine-MINI, AQUAREAD, Broadstairs,
England; NBOSI CT, Neil Brown Ocean Sensors Inc., East Falmouth, MA, USA)
could be accommodated within the forebody volume, as could a custom CTD with
tracker from previous work (?). The biohybrid robot could also be used for
ocean chemistry profiling (shown in green) by embedding a dissolved oxygen
(DO) probe (#Kit-103DX, Atlas Scientific, Long Island City, NY, USA), or to
study ocean chemical plumes using a pH probe (#ENV-20-pH, Atlas Scientific,
Long Island City, NY, USA). While currently too large to fit into the forebody
volumes tested in this work, a CO2 probe to aid in the understanding of carbon
distribution in the ocean may fit in a future forebody design (CO2 probe,
aquams, Maxeville, France).

This biohybrid robot itself can be fabricated for approximately $20 plus the
cost of the jellyfish, which can be procured from local aquaria or collected
from the ocean. This accessibility of the technology can democratize ocean
exploration and make it feasible to deploy swarms of biohybrid robots for
comprehensive ocean exploration on timescales commensurate with a changing
ocean.

## MATERIALS AND METHODS

The objective of this study was to design and test a biohybrid jellyfish robot
with passive mechanical forebodies for streamlining and study its performance
in free swimming conditions.

### Animal husbandry

We conducted large scale swimming experiments with $N=3$ jellyfish of
diameters $d_{1}=16.0\pm 0.5$ cm, $d_{2}=16.5\pm 0.5$ cm, and $d_{3}=18.0\pm
0.5$ cm and small-scale flow experiments with 1 jellyfish of diameter
$d_{4}=19.0\pm 0.5$ cm. Aurelia aurita were obtained from Cabrillo Marine
Aquarium (San Pedro CA, USA) and cared for in a 453-liter psuedokreisel tank
(Jelliquarium 360, Midwater Systems, Thousand Oaks, CA, USA). The tank was
filled with artificial seawater made from sea salt (Instant Ocean Sea Salt,
Spectrum Brands, Blacksburg, VA, USA) and deionized water balanced at 35 parts
per thousand (PPT) and kept at 21∘C. The jellyfish were hand fed twice daily
with live Artemia franciscana brine shrimp (Hatching Shell-Free Brine Shrimp
Eggs E-Z Egg, Brine Shrimp Direct, Ogden, UT, USA). Brine shrimp were hatched
every other day in aerated hatchery cones of artificial seawater heated to
28∘C for 24 hours, then transferred to beakers at 21∘C. Brine shrimp were fed
PHYTO-Feast (Reef Nutrition, Reed Mariculture Inc., Campbell, CA) and enriched
with SELCO (Self-Emulsifying Lipid Concentrate, Brine Shrimp Direct, Ogden,
UT, USA). The animals’ care was in accordance with institutional guidelines.

### System architecture

The biohybrid robot electronics schematic and all components are shown in
Figure 7A. Shown are a 5 cm scale bar and 1. the jellyfish bell, 2. the 3D
printed $2.50\pm 0.03$ cm diameter electronics housing containing the swim
controller and battery, 3. the two electrodes embedded in the jellyfish muscle
tissue, 4. a wooden pin attaching the electronics to the jellyfish, 5. a 3D
printed hemi-ellipsoid with $AR_{1}=0.31$, 6. a hemi-ellipsoid with
$AR_{2}=0.66$, and 7. a hemi-ellipsoid with $AR_{3}=1.0$. Figure B is a sample
image of a jellyfish biohybrid robot with attached $AR_{3}$ forebody swimming
in a lab saltwater tank. This image has been digitally manipulated to remove
reflections and outline the animal body in black.

Fig 7. Hardware schematic.
(A) Biohybrid robot electronics schematic with 5 cm scale bar. 1. Jellyfish
bell. 2. $2.50\pm 0.03$ cm diameter housing containing swim controller and
battery. 3. Electrode embedded in muscle tissue. 4. Wooden pin attaching
electronics to jellyfish. 5. 3D printed PLA forebody with $AR=0.31$. 6.
Forebody with $AR=0.66$. 7. Forebody with $AR=1.0$. (B) Jellyfish biohybrid
robot with attached hemi-ellipsoid with $AR=1.0$. This image has been
digitally manipulated to remove reflections and outline the animal body in
black.

### Components and fabrication

The electronics consist of a TinyLily processor (TinyCircuits, Akron, OH, USA)
and a 15 mAH lithium-polymer battery cell (PGEB181118, PowerStream Technology
Inc., Orem, UT, USA) enclosed in a custom $2.50\pm 0.03$ cm diameter 3D
printed housing printed on a Form 3 resin printer for waterproofing at high
pressure (Formlabs, Somerville, Massachusetts). The housing utilized a male
threaded base and female threaded cap compressing a face sealed soft Buna-N
O-ring with food safe O-ring grease. Silver wire coated in perfluoroalkoxy was
threaded through small holes in the housing and connected to platinum rods
(A-M Systems, Sequim, WA, USA) embedded in the animal muscle tissue to create
two electrodes. The holes were sealed with adhesive and LED’s were connected
in series to provide visual confirmation of stimulation as shown in Figure 7B.
This system was based on previous work (?).

Hemi-ellipsoid forebodies were designed with a density gradient such that the
extremity, the bottom as shown in Figure 7B, is more dense. When the biohybrid
robot turns, the low center of gravity creates a positive restoring torque and
passive stability constraining the biohybrid robot to vertical swimming.
Concavity at the forebody-animal interface creates a close fit against the
exumbrellar bell surface. The forebodies were 3D printed out of PLA
(Polylactic Acid) on an X1 Carbon printer (Bambulab, Shenzhen, China) and
coated with 2-part epoxy resin for waterproofing. After drying, a hole was
drilled into the forebody at the interface with the animal to create a press
fit for the wooden rod as shown in Figure 7A item 4. The forebodies were
ballasted to maintain $0.01\pm 0.01$ grams positive buoyancy in saltwater at
35 PPT.

### Small-scale flow measurement experiments

We conducted the small scale flow measurements in a 1.2 m by 0.5 m by 0.5 m
tank filled with artificial seawater at 35 PPT and seeded with 13 micron
silver coated hollow glass particles (SH400S20, Potters Industries, Carlstadt,
NJ, USA). The biohybrid robot swam down through a laser sheet generated by a
671 nm continuous wave laser (5-Watt LRS0671 DPSS Laser System, Laserglow
Technologies, North York, Ontario, Canada) shone through a 250 mm diameter
condenser lens and a sheet forming optic. We used a high speed camera normal
to the laser sheet recording at 125 frames per second with a shutter speed of
$1/200$ seconds, a resolution of 1024x1024 pixels (FASTCAM SA-Z, Photron USA
Inc, San Diego, CA, USA), and a fixed focal length macro lens (Micro-NIKKOR
105 mm with a 36 mm extension tube, Nikon, Tokyo, Japan). We processed the
images in MATLAB (Mathworks, Natick, MA, USA) by stacking every 3rd frame for
50 frames where in each frame, the biohybrid robot is translated to the
location in the previous frame such that the biohybrid robot is motionless.
This switches the video into the biohybrid robot frame of reference and the
observed particles trace out pathlines as in a long exposure image. This
process was repeated for the biohybrid robot without the forebody and with the
$AR_{3}$ forebody configuration to compare the streamlining effects.

### Large-scale swimming experiments

A 6 m tall 13,600 liter saltwater facility was constructed for testing
jellyfish biohybrid robot swimming performance. This 6 m by 1.6 m by 1.4 m, or
13.4 m3, vertical tank was filled with artificial seawater balanced at 35 PPT
and held at 21∘C (Fig. 8). The facility has a 40 gallons per minute continuous
filtration pump which turns over the entire system water volume every 90
minutes. This pump was turned off at least 2 hours before beginning
experiments to allow residual currents to dissipate. The facility has 3
acrylic windows on each side spanning the width of the tank. A GoPro (HERO9
Black, GoPro Inc, San Mateo, CA) recording at a resolution of 3840x2160 pixels
was used to record videos of the jellyfish biohybrid robot swimming trials.
The camera was centered horizontally at a height of 1.8 m from the base of the
tank and set back 2.1 m from the front. A 100 watt LED floodlight illuminated
the tank from below through a porthole on the bottom of the tank.

While salinity and temperature were carefully matched between the
pseudokreisel and this facility, the jellyfish were also acclimated by slowly
introducing water from this facility to the jellyfish for 30 minutes in a 5
gallon bucket. The jellyfish was then transferred to a glass bowl and equipped
with the electronics by inserting the wooden pin through the jellyfish stomach
and embedding both electrodes into the muscle tissue on opposite sides of the
bell margin. The jellyfish was transferred into the saltwater facility and the
appropriate forebody was mounted to complete the biohybrid robot. Synthetic
0.3 mm diameter line was tied onto the electronics housing to facilitate
recovery. The biohybrid robot was positioned in the center of the tank to
avoid contacting the walls and released to begin swimming. In order to avoid
using swimming data before the animal had reached a constant average swimming
speed and to avoid the impact of the bottom of the tank, swimming data from
the middle 2 m section of the tank was used. The model’s predictions and post
processing of experimental data confirmed that the biohybrid robot reached a
constant average swimming speed for all forebody configurations before the
region of interest. Forebody trial order was randomized and each forebody
configuration was tested 3 times for a total of 12 trials per biohybrid robot.
Between each trial, the biohybrid robot was retrieved using the synthetic line
and the tank was allowed to rest for 10 minutes before beginning the next
trial. Experiments were conducted with all 3 jellyfish over 2 days.

Fig 8. 6 m tall saltwater facility.
6 m by 1.6 m by 1.4 m, or 13,600 liter, saltwater vertical tank facility with
human for scale. The camera was placed at a height of 1.8 m from the base of
the tank and set back 2.1 m to the right of the tank.

Post-processing used a combination of a custom MATLAB program and Blender’s
(Blender Foundation, Amsterdam, Netherlands) motion tracking feature to track
the jellyfish. The MATLAB program assigned a tracking box around the biohybrid
robot and used the centroid of the box to calculate biohybrid robot swimming
speeds. Distance vs time was plotted, a line was fit with average
$R^{2}=0.9894\pm 0.0141$, and the slope was used as the velocity. This program
was used to create Figure 3C-D. To determine instantaneous velocities, Blender
was used to track the location of a point chosen near the centroid of the
biohybrid robot. This tracking method better reproduces the intracycle speed
changes as the jellyfish biohybrid robot contracts and relaxes. This position
data was then imported into MATLAB where velocities were calculated as finite
differences. Horizontal velocities were assumed as proxies for out of plane
velocities and all 3 components were taken into account for all calculations.
Intracycle speed changes were found with a program in MATLAB that identifies
peaks and valley in the speed trace and calculates the difference from a peak
to the closest valley. Blender tracking data was used for Figures 2A-B, 4A,
and the overlayed data in Figures 5A-B. Average speeds for both programs were
the same. To model the plateau observed in Figure 5C, the terminal velocity
equation in terms of $C_{D}$ was fit to the experimental data with scaling
coefficients $u=0.1417+0.0225\sqrt{1/C_{D}}$ where $u$ is the observed speed.

The forebody drag coefficients $C_{D}$, used in calculating the streamlining
coefficient $S$ were found by conducting droptests of the forebodies.
Forebodies of the same form factor as those used in our experiments but at a
higher infill density were printed out of PLA such that they were negatively
buoyant. These forebodies were attached to synthetic 0.3 mm diameter line for
ease of recovery, submerged at the top of the vertical tank, and released to
sink at terminal velocity. These experiments were recorded using the same
Gopro camera setup and post-processing methods as described above. The
coefficients of drag $C_{D}$ were found using the drag equation (?) which
after accounting for buoyancy becomes

$C_{D}=2g(m_{f}-V_{f}\rho_{w})/(\rho_{w}S_{f}u^{2})$ (1)

where $g$ is the acceleration of gravity, $m_{f}$ is the mass of the forebody,
$V_{f}$ is the volume of the forebody, $\rho_{w}$ is the density of saltwater,
$S_{f}$ is the frontal surface area of the forebody impacting the flow, and
$u$ is the terminal velocity at which the forebody falls. Each forebody was
tested three times and the average velocity was used to calculate
$C_{D,1}=0.93\pm 0.05$, $C_{D,2}=0.52\pm 0.03$, and $C_{D,3}=0.17\pm 0.01$
where the error is the standard deviation of all three trials. The drag on the
jellyfish without a forebody was approximated as a flat plate with
$C_{D,0}=1.12$ (?).

### Kinematics Model

We derived a first-principles physics-based kinematics model for our biohybrid
robot by extending previous models (?, ?, ?). Starting with a force balance
from Newton’s second law

$\sum{F}=m_{R}du/dt$ (2)

where $\sum{F}$ is the sum of the forces, $m_{R}$ is the combined mass of the
biohybrid robot, $u$ is the swimming speed, and $t$ is time. Then the thrust
of the biohybrid robot $T$ is assumed to oppose drag $D$ and the unsteady
acceleration reaction $Ac$

$T-D-Ac=m_{R}du/dt$ (3)

with the full forces given as

$T=(\rho_{w}/A_{sub})(dV_{R}/dt)^{2}$ (4)

$D=1/2C_{d}\rho_{w}S_{R}u^{2}$ (5)

$Ac=\alpha\rho_{R}V_{R}du/dt$ (6)

where $\rho_{w}=1.024g/cm^{3}$ is the density of saltwater at 35 PPT and 21∘C,
$A_{sub}$ is the subumbrellar opening area of the jellyfish, $V_{R}$ is the
instantaneous volume of the biohybrid robot, $C_{d}$ is the coefficient of
drag calculated above, $S_{R}$ is the instantaneous frontal surface area of
the biohybrid robot, $\alpha=(h/r)^{1/4}$ is the added mass coefficient where
$h$ is the bell height and $r$ is the bell radius, and $\rho_{R}$ is the
density of the biohybrid robot. Then we define

$V_{R}=V_{c}+V_{0}+(dv/dt)t$ (7)

$S_{R}=S_{c}+S_{0}+(3/{2h})(dv/dt)t$ (8)

the relaxed volume of the jellyfish $V_{0}=2/3\pi hr^{2}$, the relaxed frontal
surface area of the jellyfish $S_{0}=\pi r^{2}$, and the change in volume
$dv/dt=V_{ex}V_{0}$ with a constant volume exchange during contraction
$V_{ex}$. Substituting these equations back into the force balance leads to

$(\rho_{w}/A_{R})(dV_{R}/dt)^{2}=1/2C_{d}\rho_{w}S_{R}u^{2}+(1+\alpha)m_{R}du/dt$
(9)

which was solved for $u$ in MATLAB using the ODE45 function. Jellyfish
geometry was measured to experimentally determine $h$ and $r$ for each
jellyfish, and video analysis led to approximate contraction and relaxation
times of $t_{c}=0.5$ and $t_{r}=1.5$.

scibib

## References

  1. 1.

E. Ramirez-Llodra, A. Brandt, R. Danovaro, B. De Mol, E. Escobar, C. R.
German, L. A. Levin, P. Martinez Arbizu, L. Menot, P. Buhl-Mortensen, B. E.
Narayanaswamy, C. R. Smith, D. P. Tittensor, P. A. Tyler, A. Vanreusel, and M.
Vecchione, Deep, diverse and definitely different: unique attributes of the
world’s largest ecosystem. Biogeosciences 7, 2851 (2010).

  2. 2.

L. A. Levin, N. Le Bris, The deep ocean under climate change. Science 350, 766
(2015).

  3. 3.

V. Ryabinin1, J. Barbière, P. Haugan, G. Kullenberg, N. Smith, C. McLean, A.
Troisi, A. Fischer, S. Aricò, T. Aarup, P. Pissierssens, M.Visbeck, H. O.
Enevoldsen, J. Rigaud, The UN decade of ocean science for sustainable
development. Frontiers in Marine Science 6, 470 (2019).

  4. 4.

A. C. Wölfl, H. Snaith, S. Amirebrahimi, C. W. Devey, B. Dorschel, V. Ferrini,
V. A. I. Huvenne, M. Jakobsson, J. Jencks, G. Johnston, G. Lamarche, L. Mayer,
D. Millar, T. Haga Pedersen, K. Picard, A. Reitz, T. Schmitt, M. Visbeck, P.
Weatherall, R. Wigley, Seafloor mapping–the challenge of a truly global ocean
bathymetry. Frontiers in Marine Science p. 283 (2019).

  5. 5.

H. M. Benway, L. Lorenzoni, A. E. White, B. Fiedler, N. M. Levine, D. P.
Nicholson, M. D. DeGrandpre, H. M. Sosik, M. J. Church, T. D. O’Brien, M.
Leinen, R. A. Weller, D. M. Karl, S. A. Henson1, R. M. Letelier, Ocean time
series observations of changing marine ecosystems: an era of integration,
synthesis, and societal applications. Frontiers in Marine Science 6, 393
(2019).

  6. 6.

J. Yuh, G. Marani, D. R. Blidberg, Applications of marine robotic vehicles.
Intelligent service robotics 4, 221 (2011).

  7. 7.

E. Zereik, M. Bibuli, N. Miskovi C., P. Ridao, A. Pascoal, Challenges and
future trends in marine robotics. Annual Reviews in Control 46, 350 (2018).

  8. 8.

D. Ventura, A. Bonifazi, M. F. Gravina, A. Belluscio, G. Ardizzone, Mapping
and classification of ecologically sensitive marine habitats using unmanned
aerial vehicle (UAV) imagery and object-based image analysis (OBIA). Remote
Sensing 10, 1331 (2018).

  9. 9.

N. Reul, S. Fournier, J. Boutin, O. Hernandez, C. Maes, B. Chapron, G. Alory,
Y. Quilfen, J. Tenerelli, S. Morisset, Y. Kerr, S. Mecklenburg and S. Delwart,
Sea surface salinity observations from space with the SMOS satellite: A new
means to monitor the marine branch of the water cycle. Surveys in Geophysics
35, 681 (2014).

  10. 10.

D. Parker, C. Folland, M. Jackson, Marine surface temperature: observed
variations and data requirements. Climatic change 31, 559 (1995).

  11. 11.

M. T. Kavanaugh, T. Bell, D. Catlett, M. A. Cimino, Scott C. Doney, Willem
Klajbor, Monique Messié, Enrique Montes, Frank E. Muller-Karger, Daniel Otis,
Jarrod A. Santora, Isaac D. Schroeder, Joaquin Triñanes, David A. Siegel,
Satellite remote sensing and the marine biodiversity observation network.
Oceanography 34, 62 (2021).

  12. 12.

D. L. Rudnick, Ocean research enabled by underwater gliders. Annual review of
marine science 8, 519 (2016).

  13. 13.

P. Favali, L. Beranzoli, Seafloor observatory science: A review. Annals of
geophysic (2006).

  14. 14.

A. Sahoo, S. K. Dwivedy, P. Robi, Advancements in the field of autonomous
underwater vehicle. Ocean Engineering 181, 145 (2019).

  15. 15.

M. A. Charette, W. H. Smith, The volume of Earth’s ocean. Oceanography 23, 112
(2010).

  16. 16.

R. K. Katzschmann, J. DelPreto, R. MacCurdy, D. Rus, Exploration of underwater
life with an acoustically controlled soft robotic fish. Science Robotics 3,
eaar3449 (2018).

  17. 17.

G. Picardi, M. Chellapurath, S. Iacoponi S. Stefanni, C. Laschi, M. Calisti,
Bioinspired underwater legged robot for seabed exploration with low
environmental disturbance. Science Robotics 5, eaaz1012 (2020).

  18. 18.

G. Li, X. Chen, F. Zhou, Y. Liang, Y. Xiao, X. Cao, Z. Zhang, M. Zhang, B. Wu,
S. Yin, Y. Xu, H. Fan, Z. Chen, W. Song, W. Yang, B. Pan, J. Hou, W. Zou, S.
He, X. Yang, G. Mao, Z. Jia, H. Zhou, T. Li, S. Qu, Z. Xu, Z. Huang, Y. Luo,
T. Xie, J. Gu, S. Zhu, W. Yang, Self-powered soft robot in the Mariana Trench.
Nature 591, 66 (2021).

  19. 19.

M. I. D. Carter, K. A. Bennett, C. B. Embling, P. J. Hosegood, D. J. Russell,
Navigating uncertain waters: a critical review of inferring foraging behavior
from location and dive data in pinnipeds. Movement Ecology 4, 1 (2016).

  20. 20.

E. Cuevas, F. A. Abreu-Grobois, V. Guzman-Hernandez, M. Liceaga-Correa, R. P.
Van Dam, Post-nesting migratory movements of hawksbill turtles Eretmochelys
imbricata in waters adjacent to the Yucatan Peninsula, Mexico. Endangered
Species Research 10, 123 (2008).

  21. 21.

T. A. Mooney, K. Katija, K. A. Shorter, T. Hurst, J. Fontes, P. Afonso, ITAG:
an eco-sensor for fine-scale behavioral measurements of soft-bodied marine
invertebrates. Animal Biotelemetry 3, 1 (2015).

  22. 22.

S. Fossette, K. Katija, J. A. Goldbogen, S. Bograd, W. Patry, M. J. Howard, T.
Knowles, S. H.D. Haddock, L. Bedell, E. L. Hazen, B. H. Robison, T. A. Mooney,
K. A. Shorter, T. Bastian, A. C. Gleiss, How to tag a jellyfish? A
methodological review and guidelines to successful jellyfish tagging. Journal
of Plankton Research 38, 1347 (2016).

  23. 23.

H. Chung, J. Lee, W. Y. Lee, A review: Marine bio-logging of animal behaviour
and ocean environments. Ocean Science Journal 56, 117 (2021).

  24. 24.

C. J. Sanchez, C.-W. Chiu, Y. Zhou, J. M. González, S. B. Vinson, H. Liang,
Locomotion control of hybrid cockroach robots. Journal of The Royal Society
Interface 12, 20141363 (2015).

  25. 25.

A. Bozkurt, R. F. Gilmour, A. Lal, Balloon-assisted flight of radio-controlled
insect biobots. IEEE transactions on biomedical engineering 56, 2304 (2009).

  26. 26.

S.-J. park, M. Gazzola, K. S. Park, S. Park, V. D. Santo, E. L. Blevins, J. U.
Lind, P. H. Campbell, S Dauth, A K. Capulli, F S. Pasqualini, S Ahn, A Cho, H.
Yuan, B. M. Maoz, R. Vijaykumar, J.-W. Choi, K Deisseroth, G V. Lauder, L.
Mahadevan, K. K. Parker, Phototactic guidance of a tissue-engineered soft-
robotic ray. Science 353, 158 (2016).

  27. 27.

J. C. Nawroth, H. Lee, A.W. Feinberg, C. M. Ripplinger, M. L. McCain, A.
Grosberg, J. O. Dabiri, K. K. Parker, A tissue-engineered jellyfish with
biomimetic propulsion. Nature biotechnology 30, 792 (2012).

  28. 28.

J. E. Purcell, Jellyfish and ctenophore blooms coincide with human
proliferations and environmental perturbations. Annual review of marine
science 4, 209 (2012).

  29. 29.

J. E. Purcell, R. R. Hopcroft, K. N. Kosobokova, T. E. Whitledge,
Distribution, abundance, and predation effects of epipelagic ctenophores and
jellyfish in the western Arctic Ocean. Deep Sea Research Part II: Topical
Studies in Oceanography 57, 127 (2010).

  30. 30.

N. D. Gallo, J. Cameron, K. Hardy, P. Fryer, D. H. Bartlett, L. A. Levin,
Submersible-and lander-observed community patterns in the Mariana and New
Britain trenches: influence of productivity and depth on epibenthic and
scavenging communities. Deep Sea Research Part I: Oceanographic Research
Papers 99, 119 (2015).

  31. 31.

J. E. Purcell, Environmental effects on asexual reproduction rates of the
scyphozoan Aurelia labiate. Marine Ecology Progress Series 348, 183 (2007).

  32. 32.

B. J. Gemmell, J. H. Costello, S. P. Colin, J. O. Dabiri, D Tafti, S Priya,
Passive energy recapture in jellyfish contributes to propulsive advantage over
other metazoans. Proceedings of the National Academy of Sciences 110, 17904
(2013).

  33. 33.

N. W. Xu, J. O. Dabiri, Low-power microelectronics embedded in live jellyfish
enhance propulsion. Science Advances 6, eaaz3194 (2020).

  34. 34.

N. Xu, O. Lenczewska, S. Wieten, C. Federico, J. Dabiri, Ethics of biohybrid
robotic jellyfish modification and invertebrate research. Preprint, (2020).

  35. 35.

E. S. J. Smith, G. R. Lewin, Nociceptors: a phylogenetic view. Journal of
Comparative Physiology A 195, 1089 (2009).

  36. 36.

N. W. Xu, J. P. Townsend, J. H. Costello, S. P. Colin, B. J. Gemmell, J. O.
Dabiri, Field testing of biohybrid robotic jellyfish to demonstrate enhanced
swimming speeds. Biomimetics 5, 64 (2020).

  37. 37.

A. M. Bradley, M. D. Feezor, H. Singh, F. Y. Sorrell, Power systems for
autonomous underwater vehicles. IEEE Journal of oceanic Engineering 26, 526
(2001).

  38. 38.

A. Manthiram, An outlook on lithium ion battery technology. ACS central
science 3, 1063 (2017).

  39. 39.

G. K. Batchelor, An Introduction to Fluid Dynamics (Cambridge university
press, 1967).

  40. 40.

S. F. Hoerner, Fluid-Dynamic Drag. (Bricktown, NJ, 1958).

  41. 41.

T. L. Daniel, Mechanics and energetics of medusan jet propulsion. Canadian
Journal of Zoology 61, 1406 (1983).

  42. 42.

M. J. McHenry, J. Jed, The ontogenetic scaling of hydrodynamics and swimming
performance in jellyfish (Aurelia aurita). Journal of experimental biology
206, 4125 (2003).

## Acknowledgments

The authors would like to thank Cabrillo Marine Aquarium for providing the
jellyfish used in these experiments. Funding: This work was supported by the
National Science Foundation Alan T. Waterman Award and NSF Graduate Research
Fellowship grant number DGE‐1745301. Author contributions: S.R.A. and J.O.D.
conceived of project, S.R.A. conducted experiments, S.R.A. and J.O.D. analyzed
results and wrote paper. Competing interests: The authors declare that they
have no competing financial interests. Data and materials availability:
Contact J.O.D. for data and source code.

## Supplementary materials

Video S1. Sample video of free swimming biohybrid jellyfish robot swimming
through tracer particles without attached forebody. Video is sped up at 4
times speed.
Video S2. Sample video of free swimming biohybrid jellyfish robot swimming
through tracer particles with attached $AR_{3}$ forebody. Video is sped up at
4 times speed.
Video S3. Sample timelapse video of free swimming biohybrid jellyfish robot
swimming in the 6 m saltwater facility without attached forebody. Video is
sped up at 15 times speed.
Video S4. Sample timelapse video of free swimming biohybrid jellyfish robot
swimming in the 6 m saltwater facility with attached $AR_{3}$ forebody. Video
is sped up at 15 times speed.


# Ferroelectricity in oxyfluoride Aurivillius phase Bi2TiO4F2 from first-
principles

Sarah Benomar  Laboratoire de Chimie Théorique Computationnelle et Photonique,
faculté de Chimie, Université Houari Boumedienne, BP.32, El Alia, Bab Ezzouar,
Alger, Algeria  Hania Djani  Centre de Développement des Technologies
Avancées, Cité 20 août 1956, Baba Hassen, Alger, Algeria

###### Abstract

Using first-principles calculations, we restrained the m=1 oxyfluoride
Aurivillius Bi2TiO4F2 to adopt three ordered key configurations, such as
tetragonal $I4/mmm$ phases with fluorine anions distributed either at the
octahedral apical sites or at the octahedral equatorial ones and orthorhombic
$Pma2$ phase with fluorine anions distributed over both apical and equatorial
sites. We explored the energetics of the metastable phases and their potential
for ferroelectricity and found that the least favorable configuration is
equatorial sites for fluorine and the most favorable one is the mixture
between apical and equatorial sites for fluorine. Phonon calculations
performed in the tetragonal $I4/mmm$ with fluorine at the apical sites showed
a strong c-axis polar instability $\Gamma_{3}^{-}$ that co-exists with the ab-
plane polar instability $\Gamma_{5}^{-}$, commonly found in Aurivillius oxide
phases. These two polar instabilities lead to $Pc$ phase that displays in-
plane and out-of-plane spontaneous polarizations as large as 44 $\mu$c/cm2 and
35 $\mu$c/cm2, respectively, opening new prospective design in thin films
Aurivillius FeRAM devices.

##  I introduction

Aurivillius phase structures family, with the general formula Bi2Am-1BmO3m+3
are mostly ferroelectric oxides, consisting of alternating, along the c-axis,
of [Bi2O2]+2 layers and perovskite-like (Am-1BmO3m+1)-2 blocks, with, A= Ca2+,
Ba2+, Sr2+, Bi2+ …, B= Fe3+, Ti4+, Nb5+, W6+ .., and $m$ is the number of BO6
octahedra in the perovskite-like block ($m$=1-8) and is interdependent with A
and B cations oxidation states to always provide right number of electrons to
oxygen atoms Aurivillius (1949); Perez-Mato et al. (2008); Withers et al.
(1991); Perez-Mato et al. (2004); Etxebarria et al. (2010); Hervoches et al.
(2002). Typically, oxide Aurivillius phase structures crystallize at high
temperature into centrosymmetric $I4/mmm$ tetragonal phase and adopt at room
temperature a polar orthorhombic ground state involving in-plane polar
displacement and in-plane and out-of-plane octahedral tilts about [110]t and
[001]t directions, respectively ($t$ referring to the tetragonal parent
structure). They are also known to exhibit a large spontaneous polarization Ps
$\approx$ 30-50 $\mu$c/cm2 and high Curie temperature (Tc > 350°C). Bi2W06 is
the smallest member with $m=1$. It is a one-layer strong ferroelectric with A
cation deficiency in the perovskite block, with Ps = 50 $\mu$c/cm2 and Tc =
950°C.

Oxyfluoride Aurivillius phases are also naturally occuring. The $m$=1
Bi2TiO4F2, isostructural to Bi2NbO5F and Bi2TaO5F, was first reported in 1952
Aurivillius (1953) and has been intensively studied because of its
ferroelectric and photocatalysis properties Wang et al. (2015, 2011). The
first experimental study claimed that Bi2TiO4F2 was ferroelectric with a Curie
temperature Tc $\approx$ 12°C and that the ferroelectric phase transition was
accompanied with an orthorhombic distortion of the high temperature tetragonal
paraelectric $I4/mmm$ parent phase, that lead to polar $Imm2$ phase Ismailzade
and Ravez (1978). Later Needs et $al$. Needs et al. (2005), by using the
combined refinement of the powder X-ray and neutron diffraction, refined
Bi2TiO4F2 and Bi2NbO5F into centrosymmetric tetragonal $I4/mmm$ symmetry and
the fluorine anions assigned to the equatorial sites. No evidence of an
orthorhombic distortion of the crystal structure was observed. Furthermore, E.
McCabe et $al$. McCabe et al. (2007) showed that Bi2NbO5F adopts a distorted
structure with octahedral tilts about [110]t and [001]t axis. The absence of
SHG signal led them to conclude that the structure was centrosymmetric $Pbca$
and the bond valence sum calculations indicated that fluorine anions prefer
the apical sites. However, all these efforts and numerous investigations have
not dispelled, so far, the doubts about ferroelectricity and fluorine ordering
in Bi2TiO4F2.

Here, we have explored different anions ordering configurations for Bi2TiO4F2
Aurivillius phase, with the aim to examine the effect of fluorine positions on
the structural, dielectric and dynamical properties of this compound. We have
generated three key configurations: the centrosymmetric tetragonal $I4/mmm$
phase with fluorine anions at the octahedral equatorial sites of the
perovskite-like blocks ($I4/mmm$ (Feq), see Fig. 1(a)), the centrosymmetric
tetragonal $I4/mmm$ phase with fluorine anions at the octahedral apical sites
($I4/mmm$ (Fap), see Fig. 1(b)) and the polar orthorhombic $Pma_{2}$ (Feq\+
Fap) with $fac$-arrangement, consisting of fluorine anions occupying the three
corners of one face of the octahedron (see Fig. 1(c)). Recent calculations
have recently found this latter configuration as the most stable among twenty
tested ones Morita et al. (2019).

Our study shows that Bi2TiO4F2 adopts a polar ground state of monoclinic $Pc$
symmetry with the fluorine anions distributed over the octahedral equatorial
and apical sites (Feq\+ Fap) in a $fac$-arrangement. Phonon calculations
performed in the three different configurations revealed numerous polar,
antipolar and octahedral tilts instabilities that lead to many metastable
phases and particularly, a strong polar instability, along the $c$-axis
direction was found in the $I4/mmm$ (Fap) configuration.

Figure 1:  High symmetry parent phases structures. (a) Tetragonal $I4/mmm$
(Feq), (b) Tetragonal $I4/mmm$ (Fap) and (c) Orthorhombic $Pma2$ (Feq+Fap). Bi
atoms in Purple, Ti in blue, F in light grey and O in red

##  II computational details

This work have been performed within density functional theory Kohn and Sham
(1965); Hohenberg and Kohn (1964) as implemented in ABINIT package Gonze et
al. (2002, 2009, 2005). The exchange and correlation functional was evaluated
with the generalized gradient approximation GGA-PBEsol as parameterized by
perdew, Burke and Ernzerhof Perdew et al. (1996) where Bi ($5d,6s,6p$), Ti
($3d,4s$), O and F ($2s,2p$) levels were treated as valence states. Wave
functions were expanded up to a kinetic energy cutoff (E$cut$) of 42 Hartrees.
Integrals over the Brillouin zone approximated by sums on 4 $\times$ 4
$\times$ 1 Monokhorst-Pack mesh of special $k$ points Monkhorst and Pack
(1976) was enough for good convergence of structural optimization and phonons
calculations (tests on 6 $\times$ 6 $\times$ 1 were also performed for
validation). The Broyden-Fletcher-Goldfarb-Shanno Schlegel (1982) minimization
algorithm was used to give the structural optimisation. We relaxed the
position of each individual atom until the absolute values of the forces on
these atoms were converged to less than $10^{-5}$ Ha/Bohr. Phonons, Born
effective were calculated with a variational approach to density functional
perturbation theory (DFPT) Gonze and Lee (1997); Baroni et al. (2001) and the
spontaneous polarization was calculated using the Berry phase formalism Resta
(1994).

##  III Results and Discussion

As mentioned in the introduction, we generated three high symmetry
hypothetical reference structures for Bi2TiO4F2, with different fluorine
anions ordering into the octahedral sites: $I4/mmm$ (Feq), $I4/mmm$ (Fap) and
$Pma_{2}$ (Feq+Fap). cells parameters and atomic coordinates are given in
Table LABEL:table:1. After fully relaxation, phonons calculations were
performed to identify the kind of instabilities that can lower the internal
energy of the different systems.

Table 1:  Calculated cell parameters for high symmetry tetragonal $I4/mmm$(Feq) and $I4/mmm$(Fap) and orthorhombic $Pma2$ (Feq+Fap) phases. Here, the tetragonal phases are doubled along the [110] diagonale $\sqrt{2}\times\sqrt{2}\times 1$, to allow the atomic positions to match those of the orthorhombic metastable phases. Our results are in good agreement with the experimental cell parameters for the $I4mmm$ phase obtained by Needs and $al$. Needs et al. (2005) with a=5.37$\AA$ and c=16.32$\AA$ and experimental atomic coordinates between brackets.  Structure | Atom |  | x | y | z
---|---|---|---|---|---
I4/mmm (F${}_{\textbf{eq}}$) |  |  |  |  |
a=$5.41\mathring{A}$ | Bi |  | 0.0000 | 0.0000 | 0.32654
c= $16.4\mathring{A}$ |  |  | (0.0000 | 0.0000 | 0.32829)
| Ti |  | 0.0000 | 0.0000 | 0.00000
| Feq |  | 0.2500 | 0.7500 | 0.00000
| Oap |  | 0.0000 | 0.0000 | 0.11208
|  |  | (0.0000 | 0.0000 | 0.11614)
| Ofl |  | 0.2500 | 0.7500 | 0.25000
I4/mmm (F${}_{\textbf{ap}}$) |  |  |  |  |
a=$5.42\mathring{A}$ | Bi |  | 0.0000 | 0.0000 | 0.33317
c= $15.2\mathring{A}$ |  |  | (0.0000 | 0.0000 | 0.32829)
| Ti |  | 0.0000 | 0.0000 | 0.00000
| Oeq |  | 0.2500 | 0.7500 | 0.00000
| Fap |  | 0.0000 | 0.0000 | 0.12847
|  |  | (0.0000 | 0.0000 | 0.11614)
| Ofl |  | 0.7500 | 0.5000 | 0.25000
Pma2 (F${}_{\textbf{eq}}$+F${}_{\textbf{ap}}$) |  |  |  |  |
a=$5.41\mathring{A}$ | Bi1 |  | 0.75000 | 0.22394 | 0.41923
b=$5.37\mathring{A}$ | Bi2 |  | 0.75000 | 0.25675 | 0.07642
c= $16.25\mathring{A}$ | Ti |  | 0.75000 | 0.28640 | 0.74983
| Oeq |  | 0.00000 | 0.50000 | 0.25601
| Oap |  | 0.25000 | 0.17128 | 0.35827
| O’ap |  | 0.25000 | 0.67704 | 0.64230
| Feq |  | 0.50000 | 0.00000 | 0.21796
| Fap |  | 0.75000 | 0.18477 | 0.87569
| F’ap |  | 0.75000 | 0.70640 | 0.12499
| Ofl |  | 0.00000 | 0.00000 | 0.49981

###  III.1 High symmetry I4/mmm (F${}_{\textbf{eq}}$) paraelectric phase

Phonon calculations were performed in this parent phase and numerous mode
instabilities were identified at high-symmetry points. They include $(i)$ two
[110]t in-plane polar modes ($\Gamma_{5}^{-}$): one consisting of polar
displacements of Bi against Ofl and Ti against Oap -Ti and Feq displace in the
same direction at variance of what is commonly observed in oxide Aurivillius
phases, where B cations displace against Oeq\- ($\Gamma_{5}^{-}$ at 163$i$
cm-1) and a second polar rigid layer mode (RL-mode), in which, perovskite
blocks displace against (Bi2O2) layers -as commonly observed in oxide
Aurivillius phases- ($\Gamma_{5}^{-}$ at 33$i$ cm-1). $(ii)$ [110]t in-plane
antipolar modes (M${}_{5}^{-}$, M${}_{5}^{+}$, X${}_{3}^{-}$ and
X${}_{4}^{-}$), consisting of inter-perovskite antipolar motion of Oap
(M${}_{5}^{-}$ at 137$i$ cm-1) and inter-perovskite antipolar motion of Feq
(M${}_{5}^{-}$ at 27$i$ cm-1); inter-(Bi2O2) layers antipolar motion of Bi
(M${}_{5}^{+}$); antipolar motion of Oap between in-plane linked octahedra
(X${}_{3}^{-}$ and X${}_{4}^{-}$), and $(iii)$ octahedral tilts about [110]t
direction of Oap and Feq (X${}_{3}^{+}$ and X${}_{4}^{+}$) and about [001]
direction of Feq (X${}_{2}^{+}$).

Considering the condensation of individual modes, reported in Table
LABEL:table:2, we see that the octahedral tilt about [110]t axis,
X${}_{3}^{+}$ (Cmca) drives the largest energy gain, followed by the polar
distortion $\Gamma_{5}^{-}$ ($Fmm2$), the antipolar M${}_{5}^{-}$ ($Cmcm$) and
the octahedral tilt about [100] axis X${}_{2}^{+}$, ($Pbam$). M${}_{5}^{+}$
($Cmca$) and X${}_{4}^{+}$ ($Cccm$) distortions give also a sizeable energy
gain. This hierarchy differs from what is noticed in the $m=1$, Bi2W06 in
which the polar $\Gamma_{5}^{-}$ distortion gives the largest energy gain
followed by the [001]t axis tilt, X${}_{2}^{+}$ Djani et al. (2012).

Considering the condensation of combined unstable modes, we identified many
additional metastable phases at lower energy (see Table LABEL:table:2). The
most stable phase is polar of $Pca2_{1}$ symmetry, resulting from the
combination of $\Gamma_{5}^{-}$, X${}_{3}^{+}$ and X${}_{2}^{+}$ ($M_{5}^{+}$
mode is brought by trilinear coupling $X_{2}^{+}$ $\oplus$ $X_{3}^{+}$
$\oplus$ $M_{5}^{+}$). This $Pca2_{1}$ phase is closely followed by polar $Pc$
phase, resulting from the combination of M${}_{5}^{-}$, X${}_{3}^{+}$ and
${}_{2}^{+}$ and containing insignificant contributions from M${}_{5}^{+}$ and
$\Gamma_{5}^{-}$ mode, not reported in Table LABEL:table:2.

The mode-by-mode contributions to the atomic distortion $\Delta$ of each
metastable phase with respect to $I4/mmm$ parent reference structure are also
reported in table LABEL:table:2. These contributions are evaluated with qAgate
software Bieder (2021), by expressing the distortion $\Delta$ in the basis of
phonon eigendisplacements vectors $\eta_{i}$ of the $I4/mmm$ phase (such that
$<\eta_{i}|M|\eta_{j}>=\delta_{ij}$), following the scheme explained in Ref.
Djani et al. (2012): $\Delta=A\sum_{i}\alpha_{i}\eta_{i}$, where $A$ is the
total distortion amplitude and $\alpha_{i}$ are the relative mode
contributions such that $\sum_{i}\alpha_{i}^{2}=1$. Contributions of distinct
phonon modes $i$ to the distortion $\Delta$ of a given phase correspond
therefore to the amplitudes $A\alpha_{i}$. From Table LABEL:table:2 (up part),
we notice a strong cooperation between $\Gamma_{5}^{-}$ [163$i$] and
$X_{3}^{+}$ (in $Aba2$) and $M_{5}^{-}$ [137$i$] and $X_{3}^{+}$ (in $Pc$),
revealed from the distortion amplitudes of these modes, which when combined,
became almost the double of those of single ones ( for $\Gamma_{5}^{-}$, 7.5
$\mathring{A}$ in $Fmm2$ to 11.8 $\mathring{A}$ in $Aba2$ or 13 $\mathring{A}$
in $Pca2_{1}$; for $X_{3}^{+}$, 11.1 $\mathring{A}$ in $Cmca$ to 22.1
$\mathring{A}$ in $Aba2$; for $M_{5}^{-}$, 8.2 $\mathring{A}$ in $Cmcm$ to
12.7 $\mathring{A}$ in $Pc$). These strong modes cooperations lead to the most
stable $Pca2_{1}$ and $Pc$ phases. $Pca2_{1}$ phase does not exhibit any
further phonon instabilities, so can be considered as the hypothetical ground
state for $I4/mmm$ (Feq) configuration.

###  III.2 High symmetry I4/mmm (F${}_{\textbf{ap}}$) paraelectric phase

Phonon calculations performed in the parent $I4/mmm$ (Fap) phase with fluorine
anions sitting in apical positions, showed the same pattern of unstable modes
as frequently found in oxide Aurivillius phases (see Djani et al. (2020)),
such as: a polar $\Gamma_{5}^{-}$ mode (the softest) that consists of polar
motion of Ti against both Oeq and Fap in perovskite blocks and Bi against Ofl
in (Bi2O2) layers (the Rigid-Layer mode is stable at 66 cm-1); an antipolar
$M_{5}^{-}$ mode that consists of inter-perovskite antipolar motion of Ti;
out-of-plane octahedral tilt of Oeq and Fap about [110]t axis (X${}_{3}^{+}$
and X${}_{4}^{+}$) and in-plane octahedral tilt of Oeq about [001]t axis
(X${}_{2}^{+}$). In-plane twists of Oeq ($\Gamma_{4}^{-}$ and $M_{4}^{-}$)
were also found. More importantly, an unstable polar mode of $\Gamma_{3}^{-}$
symmetry, polarized along the out-of-plane [001]t direction was also found at
a frequency of 65$i$ cm-1. This mode usually exists in Aurivillius phases, but
is never found unstable Co et al. (2019) (e.g. at 98 cm-1 in $I4/mmm$ (Feq),
at 88 cm-1 in $I4/mmm$ phase of Bi2WO6 and at 65 cm-1 in $I4/mmm$ phase of
Bi2W2O9). Its eigendisplacement involves mainly a rigid layer motion of the
(Bi2O2) layer against the perovskite block along the [001]t direction, with a
domination of Bi motion against Oeq (see Fig. 2).

Table 2: Modes contributions $A\alpha_{i}$ (A in $\mathring{A}$) (see the
text) of metastable phase, derived from condensation of individual and
combined unstable modes, wrt. paraelectric $I4/mmm$ (Feq) (up table) and
$I4/mmm$ (Fap) (down table) parent phases. Energy difference $\Delta E_{t}$
(in meV/f.u) wrt. paraelectric $I4/mmm$ (Feq) (up table) and $I4/mmm$ (Fap) is
also reported. Frequencies [$\omega$] are in cm-1 and Space group are given in
a conventional setting.

I4/mmm (F${}_{\textbf{eq}}$) |  |  | $\Gamma_{5}^{-}$ | M${}_{5}^{-}$ | X${}_{2}^{+}$ | X${}_{3}^{+}$ | X${}_{4}^{+}$ | X${}_{3}^{-}$ | M${}_{5}^{+}$ | X${}_{4}^{-}$ | $\Delta E$
---|---|---|---|---|---|---|---|---|---|---|---
phases | irreps | A |  [163$i$ |  33$i$] |  [137$i$ |  27$i$] |  [136$i$] |  [127$i$] |  [118$i$] |  [94$i$] |  [40$i$] |  [30$i$] |
$Fmm2$(42) | $\Gamma_{5}^{-}$ | 11.3 | [7.5 | 8.0] |  |  |  |  |  |  |  |  |  $-$328
$Cmcm$(63) | M${}_{5}^{-}$ | 10.7 |  |  | [8.2 | 5.5] |  |  |  |  |  |  |  $-$234
$Pbam$(55) | X${}_{2}^{+}$ | 7.2 |  |  |  |  | 5.2 |  |  |  |  |  |  $-$237
$Cmca$(64) | X${}_{3}^{+}$ | 11.6 |  |  |  |  |  | 11.1 |  |  |  |  |  $-$389
$Cccm$(66) | X${}_{4}^{+}$ | 8.1 |  |  |  |  |  |  | 7.6 |  |  |  |  $-$205
$Cmcm$(63) | X${}_{3}^{-}$ | 10.1 |  |  |  |  |  |  |  | 9.6 |  |  |  $-$132
$Cmca$(64) | M${}_{5}^{+}$ | 29.2 |  |  |  |  |  |  |  |  | 29.0 |  |  $-$212
$Cmma$(67) | X${}_{4}^{-}$ | 1.2 |  |  |  |  |  |  |  |  |  | 1.2 |  $-$0.28
$Abm2$(39) |  $\Gamma_{5}^{-}$+$X_{2}^{+}$\+ $X_{4}^{-}$ | 10.7 | [6.4 | 4.9] |  |  | 4.6 |  |  |  |  | 1.3 |  $-$423
$Aba2$(41) |  $\Gamma_{5}^{-}$+$X_{3}^{+}$ | 29.5 | [11.8 | 5.3] |  |  |  | 22.1 |  |  |  |  |  $-$624
$P$-$1$ (2) |  $M_{5}^{+}$+$X_{2}^{+}$\+ $X_{3}^{+}$ | 21.8 |  |  |  |  | 3.7 | 10.3 |  |  | 18.5 |  |  $-$467
$Pc$ (7) |  $M_{5}^{-}$+$X_{2}^{+}$\+ $X_{3}^{+}$ | 28.3 |  |  | [12.7 | 6.5] | 2.8 | 15.8 |  |  |  |  |  $-$807
$Pca2_{1}$(29) |  $\Gamma_{5}^{-}$\+ $X_{2}^{+}$\+ $X_{3}^{+}$\+ $M_{5}^{+}$ | 32.3 | [13.0 | 5.2] |  |  | 3.5 | 16.2 |  |  | 15.2 |  |  $-$822

I4/mmm (F${}_{\textbf{ap}}$) |  |  | $\Gamma_{5}^{-}$ | M${}_{5}^{-}$ | M${}_{4}^{-}$ | $\Gamma_{4}^{-}$ | $X_{3}^{+}$ | X${}_{2}^{+}$ | $\Gamma_{3}^{-}$ | X${}_{4}^{+}$ |  |  |  |  |  |  |  |  | $\Delta E$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
phases | irreps | A |  [200$i$] |  [197$i$] |  [110$i$] |  [107$i$] |  [99$i$] |  [90$i$] |  [65$i$] |  [54$i$] |  |  |  |  |  |  |  |  |
$Fmm2$(42) | $\Gamma_{5}^{-}$ | 2.2 | 1.9 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  $-$27
$Cmcm$(63) | M${}_{5}^{-}$ | 2.1 |  | 1.8 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  $-$24
$P4_{2}/nmc$(137) | M${}_{4}^{-}$ | 2.4 |  |  | 2.3 |  |  |  |  |  |  |  |  |  |  |  |  |  |  $-$15
$I$-$4m2$(119) | $\Gamma_{4}^{-}$ | 2.3 |  |  |  | 2.3 |  |  |  |  |  |  |  |  |  |  |  |  |  $-$14
$Cmca$(64) | X${}_{3}^{+}$ | 5.4 |  |  |  |  | 5.2 |  |  |  |  |  |  |  |  |  |  |  |  $-$56
$Pbam$(55) | X${}_{2}^{+}$ | 2.8 |  |  |  |  |  | 2.0 |  |  |  |  |  |  |  |  |  |  |  $-$14
$I4mm$(107) | $\Gamma_{3}^{-}$ | 8.2 |  |  |  |  |  |  | 8.1 |  |  |  |  |  |  |  |  |  |  $-$56
$Cmcm$(63) | X${}_{4}^{+}$ | 3.7 |  |  |  |  |  |  |  | 3.7 |  |  |  |  |  |  |  |  |  $-$10
$Abm2$(39) |  $\Gamma_{5}^{-}$\+ X${}_{2}^{+}$ | 3.7 | 1.9 |  |  |  |  | 2.2 |  |  |  |  |  |  |  |  |  |  |  $-$45
$Aba2$(41) |  $\Gamma_{5}^{-}$ \+ X${}_{3}^{+}$ | 6.8 | 2.2 |  |  |  | 5.9 |  |  |  |  |  |  |  |  |  |  |  |  $-$103
$Pc$(7) |  $\Gamma_{5}^{-}$+X${}_{3}^{+}$+X${}_{2}^{+}$+$\Gamma_{3}^{-}$ | 9.2 | 2.2 |  |  |  | 5.2 | 0.6 | 6.4 |  |  |  |  |  |  |  |  |  |  $-$126
$Aba2$(41) |  $\Gamma_{3}^{-}$ \+ X${}_{2}^{+}$ | 8.0 |  |  |  |  |  | 1.4 | 7.6 |  |  |  |  |  |  |  |  |  |  $-$60
$Abm2$(39) |  $\Gamma_{3}^{-}$ \+ X${}_{3}^{+}$ | 8.8 |  |  |  |  | 4.7 |  | 7.1 |  |  |  |  |  |  |  |  |  |  $-$87
$Pna2_{1}$(33) |  $\Gamma_{3}^{-}$\+ X${}_{4}^{+}$+M${}_{5}^{-}$ | 11.4 |  | 2.1 |  |  |  |  | 8.4 | 5.1 |  |  |  |  |  |  |  |  |  $-$114
$Pca2_{1}$(29) |  $\Gamma_{3}^{-}$ \+ X${}_{3}^{+}$\+ X${}_{2}^{+}$+M${}_{5}^{-}$ | 9.4 |  | 2.2 |  |  | 5.3 | 0.6 | 6.5 |  |  |  |  |  |  |  |  |  |  $-$124

As for $I4/mmm$ (Feq), we condensed individual and combined unstable modes
into the parent paraelectric phase $I4/mmm$ (Fap) and the results are also
reported in Table LABEL:table:2. It is worth noting that both total amplitude
of distortions and energy gain are very small compared to $I4/mmm$ (Feq). This
is consistent with previous works on Bi2NbO5F compound Medvedeva et al.
(1993), that showed that the energy to substitute Oap with fluorine anions is
much lower than the energy to substitute Oeq. The condensation of individual
modes gave, equally, the larger energy gain to $\Gamma_{3}^{-}$ ($I4mm$) and
X${}_{3}^{+}$ ($Cmca$), followed by polar $\Gamma_{5}^{-}$ ($Fmm2$) and
antipolar M${}_{5}^{-}$ ($Cmcm$). The condensation of combined modes showed
that the most stable phase is associated to $\Gamma_{3}^{-}$ and
$\Gamma_{5}^{-}$ polar modes in addition to X${}_{3}^{+}$ and X${}_{2}^{+}$
octahedral tilts. This phase of $Pc$ symmetry is polar along the in-plane
[110]t direction and out-of-plane [001]t direction and is considered as the
ground state for the (Fap) configuration, since further phonon calculations in
this phase do not reveal any instabilities. The condensation of polar
$\Gamma_{3}^{-}$, $M_{5}^{-}$ antipolar mode and X${}_{3}^{+}$ and
X${}_{2}^{+}$ octahedral tilts lead also to a stable phase of $Pca2_{1}$
symmetry, very close in energy to the ground state $Pc$.

In an aim to explore the origin of out-of-plane $\Gamma_{3}^{-}$ instability,
the ferroelectricity being directy related to the length and covalency of
cations-anions bonds Ravez and Simon (1988), we examined the selected bonds
lengths given in Table 3. We noticed that Fap is moving away from Ti and
getting closer to Bi (Ti-Fap=1.95$\mathring{A}$ and Bi-Fap=2.77$\mathring{A}$
in $I4/mmm$ (Fap), in contrast to Ti-Oap=1.83$\mathring{A}$ and Bi-
Oap=2.88$\mathring{A}$ in $I4/mmm$ (Feq)). Accordingly, the calculations of
the Born effective charges ($Z^{*}$), reported in Table LABEL:table:7, showed
that, on the one hand, the ionic character of the Ti-Fap bonding, marked by
$Z^{*}_{zz}$ Ti(+4.18$e$) that is close to its nominal value along the apical
direction can explain the elongation of the Ti-Fap bonds and the non
hybridization between F-$2p$ and Ti-$d$ orbitals and, on the other hand, the
strongly anomalous $Z^{*}$ Fap(-2.3$e$) and anomalous $Z^{*}$ Bi(+4.2$e$) can
be linked to an increasing hybridization between F-$2p$ and Bi-$6p$-$6s$
orbitals and hence, to a shortening of Bi-Fap bonds.

Figure 2:  Out-of-plane, [001]t direction polar mode, $\Gamma_{3}^{-}$. The
atomic displacement is showing a Rigid-Layer mode along the polar direction:
(Bi2O2) layers displace against Perovskite-like blocks, with a domination of
Oeq motion against Bi cations. Bi atoms in Purple, Ti in blue, F in light grey
and O in red.

The ionic character of the Ti-Feq is also confirmed in $I4/mmm$ (Feq) with
$Z^{*}_{xx=yy}$ Ti(+4.7$e$), that is close to its nominal value along the
equatorial octahedral plane. However, here, the Ti-Feq bond length is equal to
Ti-Oeq bond length in $I4/mmm$ (Fap), although Ti and Oeq are highly
hybridized ($Z^{*}_{xx=yy}$ Ti(+8.4$e$) and $Z^{*}_{xx=yy}$ Oeq(-4.3$e$)).
This equivalence in bond lengths, despite the clear difference in
hybridization between Ti and Feq/Oeq could be explained by the constrain on
bond lengths in the equatorial plane due to the mismatch between the
perovskite-like block and the (Bi2O2) layer.

Moreover, the calculation of the interatomic force constants (IFCs),
$C_{\alpha,\beta}(lk,l^{\prime}k^{\prime})$, through the expression
$F_{\alpha}(lk)=-C_{\alpha,\beta}(lk,l^{\prime}k^{\prime})\tau_{\beta}(l^{\prime}k^{\prime})$
Ghosez et al. (1999) which relates $\alpha$-component of the force
$F_{\alpha}(lk)$ on atom $k$ in cell $l$, to the induced displacement
$\tau_{\beta}(l^{\prime}k^{\prime})$ of atom $k^{\prime}$ in cell
$l^{\prime}$, as defined in Refs Ghosez et al. (1996), allow us to identify
which driving forces lead the system to structural instabilities Ghosez et al.
(1999); Amoroso et al. (2018). By convention when the IFC is positive it
corresponds to a destabilizing interaction. Accordingly, our results on IFC,
reported in Table 5, showed that the shortening of the Bi-Fap bonds in
$I4/mmm$ (Fap) configuration decreases the distance between Bi and Oeq atoms,
thus, increases their hybridization, as reflected by the large and positive
IFC of Bi-Oeq interaction (+0.011 Ha/Bohr2), in contrast to the IFC of Bi-Feq
that is positive but vanishingly small (+0.002 Ha/Bohr).

Table 3: Interatomic bond lengths ($\mathring{A}$) for selected atoms. structure | pairs of atoms |  | distance($\mathring{A}$)
---|---|---|---
I4/mmm (F${}_{\textbf{eq}}$) |  |  |
| Bi-Oap |  | 2.88
| Bi-Feq |  | 3.42
| Bi-Ofl |  | 2.28
| Ti-Oap |  | 1.83
| Ti-Feq |  | 1.91
I4/mmm (F${}_{\textbf{ap}}$) |  |  |
| Bi-Fap |  | 2.77
| Bi-Oeq |  | 3.17
| Bi-Ofl |  | 2.29
| Ti-Fap |  | 1.95
| Ti-Oeq |  | 1.91
Table 4: Nonzero elements of the Born effective charges tensor $Z^{*}$ (in Cartesian coordinates) for Bi2TiO4F2 in the paraelectric tetragonal $I4/mmm$ phase. The nominal charges of Bi, Ti, O and F are +3, +4, $-$2 and $-$1, respectively. structure | Atom | $Z_{xx}^{*}$ | $Z_{yy}^{*}$ | $Z_{zz}^{*}$ | $Z_{xy}^{*}$ | $Z_{yx}^{*}$ | $Z_{xz}^{*}$ | $Z_{zx}^{*}$
---|---|---|---|---|---|---|---|---
I4/mmm (F${}_{\textbf{eq}}$) |  |  |  |  |  |  |  |
| Bi | 5.2 | 5.2 | 4.0 | 0.0 | 0.0 | 0.0 | 0.0
| Ti | 4.7 | 4.7 | 9.0 | 0.0 | 0.0 | 0.0 | 0.0
| Feq |  $-$2.0 |  $-$2.0 |  $-$0.9 |  $-$1.3 |  $-$1.3 | 0.0 | 0.0
| Oap |  $-$2.4 |  $-$2.4 |  $-$4.9 | 0.0 | 0.0 | 0.0 | 0.0
| Ofl |  $-$3.0 |  $-$3.0 |  $-$2.6 | 0.0 | 0.0 | 0.0 | 0.0
I4/mmm (F${}_{\textbf{ap}}$) |  |  |  |  |  |  |  |
| Bi | 4.8 | 4.8 | 4.2 | 0.0 | 0.0 | 0.0 | 0.0
| Ti | 8.4 | 8.4 | 4.2 | 0.0 | 0.0 | 0.0 | 0.0
| Oeq |  $-$4.3 |  $-$4.3 |  $-$2.0 |  $-$2.3 |  $-$2.3 | 0.0 | 0.0
| Fap |  $-$1.7 |  $-$1.7 |  $-$2.3 | 0.0 | 0.0 | 0.0 | 0.0
| Ofl |  $-$2.9 |  $-$2.9 |  $-$2.0 | 0.0 | 0.0 | 0.0 | 0.0
Table 5: Interatomic force constants (IFC) in (Ha/$bohr^{2}$) between different selected pairs of atoms in their local coordinates system along the longitudinal ($\parallel$) direction. Dipole-Dipole (DD) and Short-Range (SR) contributions are also reported. structure |  | Total | DD | SR
---|---|---|---|---
I4/mmm (F${}_{\textbf{ap}}$) |  |  |  |
| Bi-Fap($\parallel$) | 0.011 | 0.013 |  $-$0.002
| Bi-Oeq($\parallel$) | 0.011 | 0.012 |  $-$0.001
I4/mmm (F${}_{\textbf{eq}}$) |  |  |  |
| Bi-Oap($\parallel$) | 0.015 | 0.022 |  $-$0.007
| Bi-Feq($\parallel$) | 0.002 | 0.004 |  $-$0.002

###  III.3 Spontaneous polarization

In order to estimate the spontaneous polarization Ps for the different
metastable polar phases derived from $I4/mmm$ (Feq) and $I4/mmm$ (Fap) parent
phases, we used two methods: from the Berry phase approach Resta (1994), with
the careful determination of the polarization quantum (see Djani et al.
(2012)) and from the knowledge of $Z^{*}$ and atomic displacements following
the expression.

$P_{s,\alpha}=\frac{e}{\Omega}\sum_{k,\beta}Z^{*}_{k,\alpha\beta}\delta\tau_{k,\beta}$

where $\delta\tau_{k,\beta}$ is the displacement of atom k along direction
$\beta$ from the paraelectric to the ferroelectric phase,
$Z^{*}_{k,\alpha\beta}$the Born effective charge tensor of atom k, and
$\Omega$, the unit cell volume. The results, reported in Table LABEL:table:6
showed

Table 6: Spontaneous polarization $P_{s}$ (in $\mu$c/cm2) of different metastable polar phases of Bi2TiO4F2. $P_{s}$(BP) refers to the polarization computed using the Berry-phase approach, and $P_{s}$($Z^{*}$) corresponds to the polarization deduced from the knowledge of the Born effective charges and atomic displacements. structure | phase |  | direction |  $P_{s}$($Z^{*}$) |  |  $P_{s}$(BP)
---|---|---|---|---|---|---
I4/mmm (F${}_{\textbf{eq}}$) |  |  |  |  |  |
|  $Fmm2$(42) |  | [110]t | 71 |  | 72
|  $Aba2$(41) |  | [110]t | 136 |  | 138
|  $Pca2_{1}$(29) |  | [110]t | 112 |  | 110
I4/mmm (F${}_{\textbf{ap}}$) |  |  |  |  |  |
|  $Fmm2$(42) |  | [110]t | 34 |  | 30
|  $I4mm$(107) |  | [001]t | 45 |  | 41
|  $Pc$(7) |  | [110]t | 39 |  | 44
|  |  | [001]t | 37 |  | 35

that, in the structure with (Feq) configuration, the polar phases display Ps
along the [110]t direction. Ps in $Fmm2$ phase (with only $\Gamma_{5}^{-}$
polar mode) is smaller than Ps in $Aba2$ ($\Gamma_{5}^{-}$ and $X_{3}^{+}$)
and $Pca2_{1}$ ($\Gamma_{5}^{-}$, $X_{3}^{+}$, $X_{2}^{+}$ and $M_{5}^{+}$)
phases. This is in contrast to the results obtained in Bi2WO6 Djani et al.
(2012), in which the polarization decreases from $Fmm2$ to $Aba2$ to
$Pca2_{1}$. The increasing of Ps in $Aba2$ of Bi2TiO4F2 with (Feq) is due to
the strong cooperation between $\Gamma_{5}^{-}$ and $X_{3}^{+}$ modes, as
mentioned in section A, which is weak in the case of Bi2WO6 (see Djani et al.
(2012)). The polarization can reach, in this case, a very large value of 138
$\mu$c/cm2 (see Table LABEL:table:6).In the structure with (Fap)
configuration, Ps expands along both [110]t and [001]t directions. Ps in
$I4mm$ phase (with only $\Gamma_{3}^{-}$ polar mode) is larger than Ps in
$Fmm2$ ($\Gamma_{5}^{-}$). In this (Fap) configuration, the value of Ps is
almost constant in the different polar phases, in agreement with the weak
cooperative coupling observed between the polar modes and octahedral tilts, as
mentioned in section B. Finally, the ground state $Pc$ displays Ps along two
directions:44 $\mu$c/cm2 along [110]t direction and 35 $\mu$c/cm2 along [001]t
direction. This finding is of valuable importance since ferroelectricity in
Aurivillius phases is usually in-plane [110]t oriented and not very suitable
for thin films application, such as nonvolatile memory devices, that would
prefer a polarization switching along out-of-plane [001]t direction to detect
the switching states during read and write operations Auciello et al. (1998).

###  III.4 High symmetry Pma2 (F${}_{\textbf{eq}}$+F${}_{\textbf{ap}}$)

When fluorine anions are distributed over apical and octahedral sites in
$fac$-arrangement, the highest symmetry phase is polar with orthorhombic
$Pma2$ symmetry, with always the same sort of anion (Oap or Fap) on both sides
of (Bi2O2) layer. Phonon calculations performed in this $Pma2$ phase showed
four unstable modes: three modes of $\Gamma_{4}$ symmetry at frequencies of
68$i$, 61$i$ and 33$i$ cm-1 and one mode of $\Gamma_{2}$ symmetry at 44$i$
cm-1. $\Gamma_{4}$ [68$i$] and $\Gamma_{4}$ [61$i$] are octahedral tilts about
[001] direction. Because of absence of inversion center and mirror plane
perpendicular to [001] direction, octahedral tilts within the two perovskite
blocks are decoupled: each $\Gamma_{4}$ mode defines tilts of octahedra within
one single block. $\Gamma_{4}$ [33$i$] corresponds to a translation of (Bi2O2)
layer against Oap and $\Gamma_{2}$ is an octahedral tilt about [110]
direction.

We reported in Table LABEL:table:5, the metastable phases arising from the
condensation of individual and combined modes. We noticed that the energy gain
in this (Fap+Feq) configuration is very small in comparison to (Fap) and (Feq)
configurations. The condensation of one $\Gamma_{4}$ mode brought the
contribution of the other $\Gamma_{4}$ modes and the larger energy gain (21
meV/f.u) is attributed to the phase $Pc$ involving the largest contribution of
$\Gamma_{4}$ [61$i$]. Here, we notice, that contrary to $I4/mmm$
(F${}_{\textbf{eq}}$) and $I4/mmm$ (F${}_{\textbf{ap}}$) configurations, the
condensation of [001] and [110] octahedral tilts together ($\Gamma_{4}$ \+
$\Gamma_{2}$) in $Pma2$ does not increase the system stability. In this case,
the ground state is polar of $Pc$ symmetry derived from the condensation of
unstable $\Gamma_{4}$ [001] octahedral tilt mode, but, we were unable to
calculate the spontaneous polarization, since no related centrosymmetric
parent phase was available in this configuration Spaldin (2012).

Table 7: Modes contributions $A\alpha_{i}$ (A in $\mathring{A}$) (see the text) of metastable phase, derived from condensation of individual and combined unstable modes, wrt. $Pma2$ orthorhombic parent phase. Energy difference $\Delta E_{t}$ (in meV/f.u) wrt. $Pma2$ is also reported. Frequencies [$\omega$] are in cm-1 and Space group are given in a conventional setting. |  |  |  |  | $\Gamma_{4}$ | $\Gamma_{2}$ |  | $\Delta E$
---|---|---|---|---|---|---|---|---
phase |  |  |  | A | [$i$68 $i$61 $i$33] | [$i$44] |  |
$Pc$(7) |  | $\Gamma_{4}$ |  | 7.6 | [2.7 1.9 5.5] |  |  |  $-$17
$Pc$(7) |  | $\Gamma_{4}$ |  | 7.2 | [2.7 2.4 5.6] |  |  |  $-$21
$P2$(3) |  | $\Gamma_{2}$ |  | 5.4 |  | 4.9 |  |  $-$11
$Pc$(7) |  | $\Gamma_{4}$ |  | 8.2 | [2.3 1.9 6.8] |  |  |  $-$19
$P1$(1) |  |  $\Gamma_{4}$\+ $\Gamma_{2}$ |  | 5.8 | [- 2.2 -] | 4.7 |  |  $-$16
Figure 3: Summary diagram of relevant metastable phases of Bi2TiO4F2 in
different fluorine sites ordering. In the colored boxes are reported the
symmetry of the phase (in conventional setting), the involved modes, their
irreps and frequencies (between brackets) and $\Delta E_{t}$ (in mev/f.u.)
with respect to $I4/mmm$ (Feq) taken as the zero energy.

In previous sections, we have discussed independently the key configurations
of fluorine sites ordering in Bi2TiO4F2. In Figure 3, we provided a global
view recapitulation, comparing the internal energies of the three parent
phases and the metastable phases resulting from the condensation of their
unstable modes. We found that $I4/mmm$ (Feq) is the less stable configuration,
followed by $I4/mmm$ (Fap) to finally identify the orthorhombic $Pma2$
(Fap+Feq) as the most stable configuration (consistently with the calculations
of K. Morita Morita et al. (2019)) and its derived polar monoclinic $Pc$ phase
as the most stable among the different metastables phases explored herein. Our
calculations do not fit the $I4/mmm$ experimental model of Needs et $al$.
Needs et al. (2005) with fluorine assigned to equatorial sites thanks to Bond
Valence Sum calculations. However, it is important to point out that these
authors mentioned the possibility of symmetry lowering to polar $I4mm$, via a
phase transition occuring with a very small perturbation that is not apparent
to diffraction techniques Weller et al. (2004); Knee et al. (2000); Knee and
Weller (2001); Hyatt et al. (2003). In particular, these authors exclude the
possible assignment of an equal mixture of oxide and fluorine anions over
apical and equatorial sites, although their BVS calculations showed close
values between those on mixture of fluorine and oxide anions and those on
fluorine at equatorial sites only. This excluded assignment is consistent with
our (Fap+Feq) $fac$-arrangement.

##  IV conclusion

We have shown in the $m$=1 oxyfluoride Aurivillius phase Bi2TiO4F2, that the
ground state is polar of $Pc$ symmetry and that the ordered distribution of
fluorine anions over $both$ apical and equatorial octahedral sites in the
perovskite-like blocks (F${}_{\textbf{eq}}$+F${}_{\textbf{ap}}$) is favored.
When fluorine anions are positioned at equatorial sites only
(F${}_{\textbf{eq}}$), the most favorable phase is polar of $Pca2_{1}$
symmetry, deriving from small distortions of the parent $I4/mmm$ phase,
involving in-plane polar and antipolar displacements of Bi and Ti cations and
in-plane and out-of-plane octahedral tilts. Whereas, when fluorine anions are
at apical sites only (F${}_{\textbf{ap}}$), the most favorable phase is polar
of $Pc$ symmetry, also deriving from parent $I4/mmm$ phase and involving in-
plane polar displacement of Bi and Ti cations, out-of-plane polar displacement
of Bi against Oeq and in-plane and out-of-plane octahedral tilts. In
particular, this $Pc$ phase is displaying two distinct spontaneous
polarizations: in-plane Ps= 44 $\mu$c/cm2 and out-of-plane Ps= 35 $\mu$c/cm2.
This phase is of great importance in FeRAM’s thin films applications and
appropriate strain engineering on this phase should be investigated to
increase its stability.

## References

  * Aurivillius (1949) B. Aurivillius, Arki for Kemi 1, 463 (1949).
  * Perez-Mato et al. (2008) J. M. Perez-Mato, P. Blaha, K. Schwarz, M. Aroyo, D. Orobengoa, I. Etxebarria, and A. García, Phys. Rev. B 77, 184104 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.77.184104.
  * Withers et al. (1991) R. L. Withers, J. G. Thompson, and A. D. Rae, J. Solid State Chem. 94, 404 (1991).
  * Perez-Mato et al. (2004) J. M. Perez-Mato, M. Aroyo, A. García, P. Blaha, K. Schwarz, J. Schweifer, and K. Parlinski, Phys. Rev. B 70, 214111 (2004).
  * Etxebarria et al. (2010) I. Etxebarria, J. M. Perez-Mato, and P. Boullay, Ferroelectrics 401, 17 (2010).
  * Hervoches et al. (2002) C. H. Hervoches, A. Snedden, R. Riggs, S. H. Kilcoyne, P. Manuel, and P. Lightfoot, Journal of Solid State Chemistry 164, 280 (2002), ISSN 0022-4596, URL http://www.sciencedirect.com/science/article/pii/S0022459601994733.
  * Aurivillius (1953) B. Aurivillius, Arkiv for kemi 5, 39 (1953).
  * Wang et al. (2015) W. Wang, W. Mao, Z. Ye, and J. Huang, RSC Adv. 5, 81087 (2015), URL http://dx.doi.org/10.1039/C5RA14288A.
  * Wang et al. (2011) S. Wang, B. Huang, Z. Wang, Y. Liu, W. Wei, X. Qin, X. Zhang, and Y. Dai, Dalton Trans. 40, 12670 (2011), URL http://dx.doi.org/10.1039/C1DT10889A.
  * Ismailzade and Ravez (1978) I. H. Ismailzade and J. Ravez, Ferroelectrics 21, 423 (1978), eprint https://doi.org/10.1080/00150197808237285, URL https://doi.org/10.1080/00150197808237285.
  * Needs et al. (2005) R. L. Needs, S. E. Dann, M. T. Weller, J. C. Cherryman, and R. K. Harris, J. Mater. Chem. 15, 2399 (2005), URL http://dx.doi.org/10.1039/B502499D.
  * McCabe et al. (2007) E. E. McCabe, I. P. Jones, D. Zhang, N. C. Hyatt, and C. Greaves, J. Mater. Chem. 17, 1193 (2007), URL http://dx.doi.org/10.1039/B613970A.
  * Morita et al. (2019) K. Morita, J.-S. Park, S. Kim, K. Yasuoka, and A. Walsh, The Journal of Physical Chemistry C 123, 29155 (2019), eprint https://doi.org/10.1021/acs.jpcc.9b09806, URL https://doi.org/10.1021/acs.jpcc.9b09806.
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965), URL https://link.aps.org/doi/10.1103/PhysRev.140.A1133.
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964), URL https://link.aps.org/doi/10.1103/PhysRev.136.B864.
  * Gonze et al. (2002) X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, et al., Computational Materials Science 25, 478 (2002), ISSN 0927-0256, URL http://www.sciencedirect.com/science/article/pii/S0927025602003257.
  * Gonze et al. (2009) X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Côté, et al., Computer Physics Communications 180, 2582 (2009), ISSN 0010-4655, 40 YEARS OF CPC: A celebratory issue focused on quality software for high performance, grid and novel computing architectures, URL http://www.sciencedirect.com/science/article/pii/S0010465509002276.
  * Gonze et al. (2005) X. Gonze, G. Rignanese, M. Verstraete, J. Beuken, Y. Pouillon, R. Caracas, F. Jollet, M. Torrent, G. Zerah, M. Mikami, et al., Zeitschrift für Kristallographie 220, 558 (2005), ISSN 0044-2968, © 2005 Oldenbourg Wissenschaftsverlag. Published in Zeitschrift für Kristallographie and uploaded in accordance with the publisher’s self archiving policy.
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996), URL https://link.aps.org/doi/10.1103/PhysRevLett.77.3865.
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976), URL https://link.aps.org/doi/10.1103/PhysRevB.13.5188.
  * Schlegel (1982) H. B. Schlegel, The Journal of Chemical Physics 77, 3676 (1982), eprint https://doi.org/10.1063/1.444270, URL https://doi.org/10.1063/1.444270.
  * Gonze and Lee (1997) X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997), URL https://link.aps.org/doi/10.1103/PhysRevB.55.10355.
  * Baroni et al. (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001), URL https://link.aps.org/doi/10.1103/RevModPhys.73.515.
  * Resta (1994) R. Resta, Rev. Mod. Phys. 66, 899 (1994), URL https://link.aps.org/doi/10.1103/RevModPhys.66.899.
  * Djani et al. (2012) H. Djani, E. Bousquet, A. Kellou, and P. Ghosez, Phys. Rev. B 86, 054107 (2012), URL http://link.aps.org/doi/10.1103/PhysRevB.86.054107.
  * Bieder (2021) J. Bieder, _Qt interface for agate_ (2021), URL https://doi.org/10.5281/zenodo.4606005.
  * Djani et al. (2020) H. Djani, E. E. McCabe, W. Zhang, P. S. Halasyamani, A. Feteira, J. Bieder, E. Bousquet, and P. Ghosez, Phys. Rev. B 101, 134113 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.101.134113.
  * Co et al. (2019) K. Co, F.-C. Sun, S. P. Alpay, and S. K. Nayak, Phys. Rev. B 99, 014101 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.014101.
  * Medvedeva et al. (1993) N. I. Medvedeva, S. A. Turzhevsky, V. A. Gubanov, and A. J. Freeman, Phys. Rev. B 48, 16061 (1993), URL https://link.aps.org/doi/10.1103/PhysRevB.48.16061.
  * Ravez and Simon (1988) J. Ravez and A. Simon, Ferroelectrics 81, 309 (1988), eprint https://doi.org/10.1080/00150198808008870, URL https://doi.org/10.1080/00150198808008870.
  * Ghosez et al. (1999) P. Ghosez, E. Cockayne, U. V. Waghmare, and K. M. Rabe, Phys. Rev. B 60, 836 (1999).
  * Ghosez et al. (1996) P. Ghosez, X. Gonze, and J.-P. Michenaud, Europhysics Letters (EPL) 33, 713 (1996).
  * Amoroso et al. (2018) D. Amoroso, A. Cano, and P. Ghosez, Phys. Rev. B 97, 174108 (2018).
  * Auciello et al. (1998) O. Auciello, J. F. Scott, and R. Ramesh, Physics Today 51, 22 (1998), eprint https://doi.org/10.1063/1.882324, URL https://doi.org/10.1063/1.882324.
  * Spaldin (2012) N. A. Spaldin, Journal of Solid State Chemistry 195, 2 (2012), ISSN 0022-4596, polar Inorganic Materials: Design Strategies and Functional Properties, URL https://www.sciencedirect.com/science/article/pii/S0022459612003234.
  * Weller et al. (2004) M. T. Weller, R. W. Hughes, J. Rooke, C. S. Knee, and J. Reading, Dalton Trans. pp. 3032–3041 (2004), URL http://dx.doi.org/10.1039/B401787K.
  * Knee et al. (2000) C. S. Knee, B. D. Rainford, and M. T. Weller, J. Mater. Chem. 10, 2445 (2000), URL http://dx.doi.org/10.1039/B006120O.
  * Knee and Weller (2001) C. S. Knee and M. T. Weller, J. Mater. Chem. 11, 2350 (2001), URL http://dx.doi.org/10.1039/B102092G.
  * Hyatt et al. (2003) N. C. Hyatt, J. A. Hriljac, and T. P. Comyn, Materials Research Bulletin 38, 837 (2003), ISSN 0025-5408, URL https://www.sciencedirect.com/science/article/pii/S0025540803000321.