(3 + 2)D modulation governs vacancy ordering and oxide-ion transport in γ-type BIMEVOX conductors

Abstract

Incommensurate structural modulation is a defining yet poorly understood feature of several functional solids, particularly regarding its impact on defect dynamics in fast-ion conductors. Here, using the model oxide-ion conductor Bi₂V_0.9Cu_0.1O_5.35 (BICUVOX.10), we achieve the first full determination of a (3 + 2)D incommensurately modulated structure in the well-known γ-type BIMEVOX family. Combined single-crystal and powder X-ray diffraction, neutron total scattering with reverse Monte Carlo modelling, and ab initio molecular dynamics (AIMD) reveal that the γ′-phase exhibits short-range oxygen-vacancy ordering that intrinsically causes its modulation behavior. Upon heating, this vacancy ordering transforms into the dynamically disordered γ-phase, establishing the structural origin of fast-ion conduction. The modulation waves strongly constrain oxygen motion within the vanadate layers, elevating the activation energy, while the coupled apical–equatorial positional modulations generate versatile V/Cu coordination geometries and a zig-zag oxide-ion diffusion network. In addition, Cu atoms were found to act as local vacancy traps. These findings identify vacancy-driven modulation as the mechanism governing the reversible γ′ ↔ γ transition and establish a general framework linking incommensurate structural modulation, defect organization, and ionic transport in complex solid-state conductors.

---

Introduction

The physical and chemical properties of functional solids are ultimately governed by their thermodynamically stable structures. Changes in thermodynamic stability often alter the degree of structural order, and when local ordering becomes incompatible with the underlying lattice periodicity, structural modulations arise.¹ Such modulations represent a compromise between the lattice's drive to maintain global symmetry and the local chemical tendency to minimize bonding energy. Modulations have been widely found in functional crystalline materials, including ferroelectrics,^2–4 thermoelectrics^5–7 and fast-ion conductors,^8–11 where lattice periodicity interplays with subtle electronic or ionic degrees of freedom. In particular, incommensurate structural modulations represent a special class of aperiodic order, yet their atomic-scale relationship with defect dynamics remains poorly understood. Resolving this connection is key to understanding how complex crystalline architectures enable fast ionic transport and to guiding the design of next-generation solid electrolytes and related energy materials.

Oxide ion conducting solids have important applications in solid oxide fuel cells (SOFCs), solid oxide electrolyzer cells (SOECs), oxygen separation membranes, and gas sensors.^12,13 Such materials provide an ideal platform to probe the interplay between structural modulation and defect dynamics. While most such solids (e.g., yttria stabilized zirconia and gadolinium doped ceria) show high oxide ion conductivity at high and intermediate temperatures, the BIMEVOXes have gained significant attention due to their high oxide ion conductivity at relatively low temperatures in the range of 300–600 °C.^14,15 The BIMEVOX family of solid electrolytes is derived from cation substitution of V and/or Bi in the layered Aurivillius-phase Bi₄V₂O₁₁, with many different solid solutions reported covering a large extent of the periodic table.^16–18 Of interest to most researchers is the highly conducting tetragonal γ-phase, which although only present in the parent Bi₄V₂O₁₁ above ca. 550 °C, can be preserved to room temperature through solid solution formation.¹⁹ The structure of the BIMEVOXes consists of alternating [Bi₂O₂]_nⁿ²⁺ and [ME_x^lV_1−xO_{3.5−(5−l)x/2−δ}]_nⁿ²⁻ (where l is the valency of the ME cation) layers, with a 2-dimensional diffusion pathway for oxide ions in the heavily disordered vanadate layer.¹⁶ Depending on the level of substitution and temperature, BIMEVOXes typically show three main polymorphs, namely α (monoclinic), β (orthorhombic), and γ (tetragonal), among which the γ-phase shows the highest ionic conductivity. It is convenient to relate the three phases through an orthorhombic mean cell of approximate dimensions a_m ≈ 5.53 Å, b_m ≈ 5.61 Å, and c_m ≈ 15.28 Å, with a_α = 3a_m, b_α = b_m, c_α = c_m, β » 90°; a_β = 2a_m, b_β = b_m, c_β = c_m; and a_γ = b_γ ≈ a_m/ , c_γ ≈ c_m.^16,20

The Cu-substituted system, BICUVOX (Bi₂V_1−xCu_xO_5.5−3x/2), initially reported by Abraham et al. shows extraordinarily high ionic conductivity, with the γ-phase BICUVOX.10 (x = 0.10) composition showing a conductivity at 300 °C in the order of 10⁻³ S cm⁻¹.²¹ The fully disordered γ-phase, present above ca. 500 °C, was determined to be tetragonal in space group I4/mmm. However, at lower temperatures, partial ordering occurs (γ′-phase) evidenced by the appearance of weak superlattice peaks and a change in slope in the Arrhenius plot of conductivity, with a slightly higher activation energy below ca. 500 °C.²¹ The reversibility of the γ ↔ γ′ (disorder ↔ order) transition was later confirmed using high-temperature X-ray diffraction on a single crystal sample of BICUVOX.12, and transmission electron microscopy (TEM) for BICUVOX.10.^22,23 Using neutron powder diffraction, the γ ↔ γ′ transition was also observed in other BIMEVOXes (e.g. ME = Co, Ge).^24,25 Muller et al. suggested that the γ ↔ γ′ transition in BICOVOX.15 is first order in nature.²⁶

The complex satellite reflections at low temperatures in the γ′-phase appear to be incommensurately related to the γ-phase.^27–29 In the BICUVOX system, the γ-type phase occurs in the composition range 0.07–0.135.^22,30 Pernot et al. observed weak incommensurate superlattice peaks using single crystal X-ray diffraction and they suggested that the room-temperature γ′-phase of BICUVOX.12 shows a modulation of tetragonal symmetry with a modulation vector q = 0.34(ma* + nb*), where m and n are integers and a* and b* are cell dimensions of the orthorhombic mean cell in reciprocal space.²² However, using TEM on a single crystal sample of BICUVOX.10, the modulation vector was later claimed to be along both a* and b* axes in a 3a_m × 3b_m × c_m cell of tetragonal symmetry.²³ In a (3 + 2)D modulated structure, the satellite reflections can be indexed as hklmn, where the hkl represent the main subcell reflections and m and n are integers to accompany q₁ and q₂ respectively. The atoms in superspace need to be described using five coordinates, x₁, x₂, x₃, x₄, and x₅, the latter two representing correlated modulation directions. The atomic parameters (occupancies, positional parameters, anisotropic displacement parameters, etc.) are described with the help of periodic modulation functions p(x₄, x₅) = p(x₄ + n₄, x₅ + n₅), where n₄ and n₅ are integers. Thus, any (3 + 2)D modulation properties, such as atomic displacements and interatomic distances, are also functions of x₄ and x₅. The coordinates x₄ and x₅ are also used to produce sections through a generalized density map which can also be used during structure solution as in the de Wolff method.³¹

Over the years, although many BIMEVOX compositions have been developed and reported, the details of the modulated structure of the γ′-phase are still unclear owing to the difficulty of structure determination based on the weak satellite reflections. The fact that the modulation peaks are more easily observed in neutron diffraction compared to X-ray diffraction indicates that the modulation in the γ′-phase likely involves ordering in the oxide-ion sublattice.^24,32 Crystallographic analysis based on neutron diffraction has shown that many BIMEVOX systems (e.g., ME = Co, Ni, Mg) have oxide-ion vacancies concentrated in the equatorial plane of the V/ME-O_x polyhedra.^32–34 Oxygen/vacancy ordering has been difficult to study because of the diffuse character of the oxide ion distribution in the vanadate layers, for example, as shown in Fourier maps for γ-Bi₂Mg_0.15V_0.85O_5.35 based on neutron diffraction data.²⁶

In previous work, we have shown the complex distribution of oxide ions and oxide ion vacancies in the vanadate layer using reverse Monte Carlo (RMC) analysis of total neutron scattering data in the Ge and Sn substituted BIMEVOX systems BIGEVOX and BISNVOX.^20,25 Two limiting models, namely the equatorial vacancy (EV) model and the apical vacancy (AV) model, were proposed to predict the solid solution limits of BIMEVOXes depending on the preferred geometry of the dopant cations.¹⁶ Based on our initial neutron diffraction work^32,33 and experimental studies confirming predominantly 2-dimensional ionic conduction in BIMEVOXes,³⁵ we have previously proposed a conduction mechanism involving reorientation and coordination number changes of vanadate polyhedra resulting in diffusion of equatorial oxide ion vacancies within the vanadate layer of BIMEVOX systems.¹⁶ However, our recent studies have revealed a combination of both equatorial and apical vacancies around the metal cations in the vanadate layer, suggesting a more complex mechanism.^20,25 Indeed, in a molecular dynamics study of the parent material, Bi₄V₂O₁₁, a diffusion pathway involving both equatorial and apical oxide ion sites was identified.¹⁹

These findings have prompted a re-examination of BICUVOX.10, arguably one of the most important BIMEVOX compositions. Many questions regarding the structure of BICUVOX.10 remain. For example, the superlattice structure remains unresolved, and it is yet to be proven that the observed modulation in the γ′-phase is truly related to oxygen atom/vacancy ordering. It is also unclear as to whether BICUVOX.10 shows similar differences in local structure between γ- and γ′-phases to those seen in the BIGEVOX system. Additionally, while it is likely that the conduction mechanism in BICUVOX.10 resembles that in Bi₄V₂O₁₁, it is unclear whether this involves both apical and equatorial vacancies and what role the dopant cation plays in this mechanism. We answer these complex questions in the present work, using a combination of powerful techniques including single crystal/powder X-ray diffraction, high-resolution neutron powder diffraction, reverse Monte Carlo (RMC) analysis of neutron total scattering, and A.C. impedance spectroscopy. These have revealed details of how the complex modulation behavior governs ionic dynamics in the γ′-phase, the temperature-dependent γ ↔ γ′ phase evolution, and short-range vacancy ordering in this important oxide ion conductor. The concept of modulation–vacancy coupling opens a new space for designing functional materials where modulation structure is not an obstacle but a controllable higher-dimensional ordering state of matter.

Methods

Materials and synthesis

Bi₂V_0.9Cu_0.1O_5.35 was synthesized by solid-state reaction. Stoichiometric quantities of Bi₂O₃ (Aldrich, 99.9%), V₂O₅ (Aldrich, 98%), and CuO (Aldrich, 99.9%) were ground together in an agate mortar for 15 min using ethanol as a dispersant. After drying at ambient temperature, the mixture was transferred to a platinum crucible and heated in air at 650 °C for 12 h. After cooling, the powder was then reground and reheated to 850 °C for 24 h before slow cooling to room temperature. To prepare high-quality single crystals, the powder obtained after heating at 650 °C was reground and reheated to 910 °C for 2 h in an O₂-rich (95%) atmosphere, followed by slow cooling to 850 °C at a cooling rate of −1 °C h⁻¹ and then to room temperature at a rate of −3 °C h⁻¹.

Characterization

Single-crystal X-ray diffraction (SCXRD) was performed on a Rigaku XtaLAB Rapid-S X-ray diffractometer at 100 K, with a 4-circle kappa goniometer, an IP detector, and Mo-Κα radiation (λ = 0.7107 Å). Data reduction was performed in the CrysAlisPro software, and the modulated structure was refined using the JANA2020 software.³⁶

Powder X-ray diffraction (PXRD) patterns were collected on a PANalytical X'Pert Pro diffractometer, using Ni-filtered Cu-Kα radiation (λ = 1.5418 Å) with an X'Celerator detector over the angular range (2θ) 5° to 120° with a step width of 0.03342° and an equivalent scan time of 200 s per step. High temperature PXRD experiments were performed on the same instrument using an Anton Paar HTK-16 furnace over the temperature range from 100 to 700 °C with a 50 °C interval on heating and cooling and a dwell time of 90 min at each temperature. Repeat experiments showed no significant differences in transition temperatures.

Neutron total scattering experiments were performed on the Polaris time-of-flight powder diffractometer at the ISIS Facility, Rutherford Appleton Laboratory, UK. Five detector banks were used for data collection, and the detector parameters are given elsewhere.²⁵ The powdered samples were measured in an 8 mm diameter vanadium can in an evacuated furnace. Data were collected at room temperature and from 100 to 700 °C in intervals of 50 °C. To allow for analysis of total neutron scattering, longer data collections were carried out at room temperature and 700 °C, equivalent to proton beam charges of ca. 1000 µA h each. At other temperatures, shorter data collections of ca. 30 µA h were performed to monitor the phase evolution. Data for the empty instrument (3500 µA h), an 8 mm diameter V-rod (3500 µA h), an empty vanadium can (ca. 800 µA h) and the empty furnace (ca. 400 µA h) were collected at room temperature for normalization and data correction purposes.

To check the reproducibility of neutron scattering data, neutron powder diffraction data and total scattering data were also collected independently on the SuperHRPD and NOVA diffractometers at the Materials and Life Science Experimental Facility (MLF), J-PARC. For total scattering measurements, the sample was contained in a V–Ni alloy cell of 6 mm diameter. Scans were carried out at room temperature for ca. 6 h and the data from the 90° bank were corrected by subtracting the intensities from the empty furnace and V–Ni alloy cell.

For electrical measurements, powder samples were pelletized at 150 MPa in a 13 mm diameter cylindrical die and then sintered at 800 °C for 5 h, followed by slow cooling. The obtained pellets were cut and then polished into blocks of approximate dimensions 2 mm × 3 mm × 5 mm and coated with platinum electrodes by cathodic discharge. Impedance was measured on a Novocontrol Alpha analyzer with a ZG4 extension interface over the frequency range of 1.0 Hz to 1.0 MHz at temperatures from 100 °C up to 820 °C at intervals of 30 °C, with 1.0 h stabilization time at each temperature. The ionic transference numbers were determined by a modified EMF method on ceramic membranes in the concentration cell: O₂|Pt|pellet|Pt|air. The automated measurement system consisted of a Solartron 1260 analyzer, a Keithley 2400LV current–voltage source, and a Keithley 2181A voltage meter. RAPIDOX sensors were used to measure oxygen partial pressure.

Scattering data analysis

For average structure analysis, the diffraction data were analyzed using Rietveld refinement with the GSAS suite of programs and the Z-Rietveld program.^37–39 Where possible joint refinements using both X-ray and neutron diffraction data were performed. In the neutron experiment, a slight reduction of the sample occurred at 700 °C, as a result of the long data collection time under vacuum. This was evident as a small increase in volume compared to the corresponding X-ray data making the combined X-ray and neutron refinement at this temperature inaccurate. Hence, at 700 °C structure refinement proceeded using only the neutron data. Unless otherwise specified, the results of the data analysis on neutron diffraction in the following are all from the Polaris diffractometer. The γ- and γ′-phase structures were refined using the tetragonal model in space group I4/mmm, with a = 3.9274 Å and c (stacking direction) = 15.4274 Å.

The neutron total scattering structure functions, S(Q), and the total pair distribution functions, G(r), were produced using the software GudrunN.⁴⁰ The local structure was modelled using the reverse Monte Carlo (RMC) method with the RMCprofile software.⁴¹ The S(Q) and G(r) functions were fitted with the neutron Bragg data used as a long-range order constraint. The initial model was constructed based on a 16 × 16 × 4 supercell of the ideal tetragonal model. Ten parallel configurations were prepared each with V atoms randomly replaced by Cu and oxygen vacancies randomly located on the oxygen equatorial positions in the vanadate layer to achieve the correct stoichiometry. A soft bond valence summation (BVS) constraint was used along with a series of bond-stretching pseudopotential constraints for metal–oxygen pairs to avoid unrealistically short bonds. Cation swapping was tested but found to have no significant influence on the fits; therefore, only translational movements of atoms were permitted. To constrain the movement of the vanadium atoms in the RMC calculations, a V–O partial pair correlation obtained from an ab initio simulation based on the refined tetragonal model was applied as an additional dataset to fit.

DFT calculations

Density functional theory (DFT) calculations were carried out using the Vienna ab initio Simulation Package (VASP),⁴² employing the projector augmented-wave (PAW) method with a plane-wave energy cutoff of 450 eV. The exchange–correlation interactions were treated using the Perdew–Burke–Ernzerhof (PBE)⁴³ functional within the generalized gradient approximation (GGA).⁴⁴ Brillouin zone sampling was performed at the Γ-point, appropriate for the supercell sizes used. Ab initio molecular dynamics (AIMD) simulations were conducted in the canonical ensemble using the Nosé–Hoover thermostat⁴⁵ at a constant temperature of 1100 K, with a time step of 1 fs and a Nosé mass corresponding to 40 fs. Three independent AIMD simulations were performed on 2 × 2 × 1 supercells, each comprising 100 ps of trajectory data. Each supercell contained 150 atoms and was chosen such that each cation had a first cationic coordination shell composed of distinct ions; in particular, the nearest cation neighbors on opposite sides of a given ion were not periodic images of one another. Because the diffusion barriers under investigation are governed by relatively local interactions and because image–image correlations were thus minimized by construction, the supercell size used is considered to be reasonably accurate in describing the properties discussed in this study as in previous studies on related systems.¹⁹ The initial 10 ps of each run were excluded from the analysis to allow for equilibration.

Both DFT- and RMC-based models were used to calculate potential energy landscapes of the studied compounds. To this end, a 3D histogram of the local ionic density G_x,y,z was computed for each model. Only oxide ions coordinated to V or Cu were considered. The sampling region was defined as a box extending 2 Å in the a- and b-directions, and 3 Å in the c-direction from each V/Cu site. A grid of 30 × 30 × 40 sampling points was used. The resulting density histograms, G_x,y,z, were smoothed using Gaussian smearing with a width of 0.1 Å. The corresponding probability-based energy distribution, E_x,y,z, was obtained via the Boltzmann relation E_x,y,z = −k_BT ln(G_x,y,z), where k_B is the Boltzmann constant and T is the temperature. The energy landscape was normalized such that the minimum value of E_x,y,z was set to zero. Diffusion pathways between saddle points (energy barriers separating oxide ion sites) and oxide ion sites were determined using the steepest descent algorithm.

Results and discussion

Modulated structure of the γ′-phase

Single crystal diffraction images of BICUVOX.10 at 100 K are shown in Fig. 1. While the main reflections of the γ′-phase can be indexed on an average tetragonal cell (a ≈ 3.92 Å and c ≈ 15.40 Å) in space group I4/mmm, this only accounted for 18% of the total observations. In order to accurately describe the structure of the γ′-phase, modulation vectors need to be considered. Indexing of the satellite reflections within the tetragonal subcell yielded 2D incommensurate modulation vectors of q₁ = (0.364, 0.364, 0) and q₂ = (−0.364, 0.364, 0), as illustrated in Fig. 1b. These findings establish a (3 + 2)-D superlattice model with superspace group I4/mmm(a,a,1)0s0s(−a,a,1)0s00.

Fig. 1 Single-crystal diffraction patterns. (a) Reconstruction of the hk0 reciprocal space layer for BICUVOX.10. (b) Projected diffraction spots in the tetragonal unit cell with satellites indexed using (3 + 2)D modulation vectors.

The initial average structure was solved in a routine way with the SUPERFLIP method integrated in JANA2020.³⁶ The positions of Bi and O(1) were readily located within the tetragonal framework. The derived difference Fourier map (Fig. S1a) exhibited diffuse electron density in the vanadate layer, indicative of significant site disorder. Three distinct oxygen sites were identified in the average (I4/mmm) structure (Fig. 2a) as seen in other BIMEVOX systems.³² Bi atoms coordinate to four O(1) atoms with shorter bonds (ca. 2.3 Å) and four O(4) atoms with longer bonds (ca. 2.6 Å), while V/Cu atoms coordinate to equatorial O(3) and apical O(2) and O(4). O(4) can be interpreted as a disordered site of O(2), where the short interatomic distance precludes simultaneous occupancy of these two sites. While this model adequately described the average structure, substantial residual diffuse electron density remained in the vanadate layer (Fig. S1b), highlighting strong positional modulations. To simplify the modulation model, O(4) was removed, and its occupancy was reassigned to O(2) to maintain charge balance. Positional modulation waves were applied to all atoms, while occupational modulation waves were refined for V/Cu, O(2), and O(3). Crystallographic and refinement parameters for the modulated structure are summarized in Table S1, with basic atomic parameters in Table S2 and a complete overview of the refined modulation functions in Table S3.

Fig. 2 Average and modulated structures of BICUVOX.10 at 25 °C. (a) An average tetragonal unit cell highlighting the complex coordination geometries of Bi and V/Cu. (b) Projection of the incommensurately modulated structure model along the [010] direction which covers 11 subcells along a, showing a nearly complete period of the modulated structure (the strong modulation region is represented by the highlighted red strip). (c) Modulated bond length of atom pairs plotted along the 1st modulation axis, t, with the 2nd modulation axis u = 0. (d) Displacement magnitude of atoms with respect to their original position in the basic structure plotted along t with u = 0. (e) Bond valence sums for V and Cu as a function of t at different positions in u. For clarity the modulation behaviors along the t-axis are expanded to show two periods.

Fig. 2b shows the refined modulated structure covering a nearly complete wavelength of the modulation vector along the a-axis (across 11 subcells). The corresponding movie of the modulation behavior is shown in Video S1. Pronounced site disorder is seen in the vanadate layer for O(2) and O(3), while the bismuthate layer remains more ordered. The disorder is reflected in the modulation of the M–O(2) and M–O(3) bond lengths (M = V/Cu), compared to the less modulated Bi–O(1) bond length along t, as shown in Fig. 2c (at u = 0). Here, t and u represent the two independent modulation axes in superspace structure, associated with q₁ and q₂, respectively. The corresponding movies for modulated bond lengths along t at various values of u are shown in Videos S2 and S3. It is also noted that Bi–O(2) modulates between a bonding interaction at ∼2.4 Å and a longer non-bonding interaction at ∼2.9 Å. Fig. 2d shows the atomic displacement for each atom in the modulated structure with respect to their average positions in the subcell. It is clearly seen that the displacement magnitudes of O(2) and O(3) are more significant than other atoms (Video S4). Fig. 2e shows the bond valence sums for V/Cu along t at different u values. The fluctuation of effective charges indicates that the coordination environment of V/Cu changes along t. The de Wolff sections (Fig. S2) show continuous electron densities for Bi, V and O(1), whereas O(3) shows discontinuous electron density and O(2) exhibits a zig-zag electron density pattern, confirming strong positional modulation. The overall modulation is highly intricate: while the bismuthate layer remains relatively rigid, severe modulations are evident in the vanadate layer, particularly for oxygen atoms, as highlighted in Video S1. This is consistent with the high mobility of oxygen in the direction parallel to the vanadate layer, which underpins the material's high oxide-ion conductivity, further supported by a.c. impedance spectroscopy on single-crystal samples.³⁵

Temperature-dependent phase evolution

With increasing temperature, the main reflections for BICUVOX.10 gradually shift towards larger d-values due to thermal expansion, as observed in a series of neutron powder diffraction and X-ray patterns from 25 to 700 °C (Fig. S3a and b). In addition to the main reflections, a set of superlattice peaks are also present in the patterns (indicated by arrows) at temperatures below 500 °C but are absent at and above this temperature (Fig. S3a). The high-resolution neutron diffraction data (Fig. S3c) also confirm this behavior. Similar behavior has been reported in other γ-phase BIMEVOXes (ME = Mg, Ge), attributed to the transition from the ordered γ′-phase to the disordered γ-phase.^25,46 To assess structural changes during phase evolution, all datasets were fitted using the average tetragonal model, although this model represents a commensurate subcell model that ignores the superlattice peaks. Fig. 3a shows the evolution of the tetragonal lattice parameters and unit cell volume. The γ′ ↔ γ phase transition is marked by a subtle increase in the a-parameter and a slight decrease in the c-parameter near 500 °C, the latter indicating reduced interlayer spacing between the bismuthate and vanadate layers upon disordering. A more direct measure of interlayer proximity is provided by the Bi–M (V/Cu) distance along the c-axis (d_Bi–M, Fig. S4a), which exhibits a sharp drop at the transition (Fig. S4b). This indicates enhanced covalency between the bismuthate and vanadate layers in the γ-phase, where Bi atoms adopt a more cubic-like coordination involving neighboring apical oxygen atoms from the vanadate layer.

Fig. 3 Structural differences between the γ- and γ′-phases. (a) Variation of lattice parameters and unit cell volume upon heating (converted to the orthorhombic mean cell; error bars smaller than marker size). (b) Partial pair correlation functions (M–O) for BICUVOX.10 in the γ′- and γ-phases. (c) Atomic density maps of the supercell configuration folded onto a single subcell for each element type at 25 °C. (d) Coordination numbers of atoms within their first coordination shell. (e) Radial variation of equatorial vacancy (EV) shell content, expressed as a ratio relative to the EV content of the first shell at 25 and 700 °C (y = 0 indicates complete randomness).

PXRD data collected during cooling confirm the reversibility of the phase transition and the absence of hysteresis (Fig. S3d). The refined neutron diffraction patterns at 25 °C (Fig. S5a) and 700 °C (Fig. S5b), along with the PXRD pattern at 25 °C (Fig. S5c), further validate the adequacy of the tetragonal average model. Refinement and structural parameters are listed in Tables S4 and S5. During refinement, the total oxygen content in the vanadate layer was initially constrained to 3.35 atoms per V/Cu atom using a soft chemical restraint on O(2), O(3), and O(4). The results indicate that O(2) and O(4) together account for approximately 2 oxygen atoms per V/Cu, implying that vacancies are localized predominantly at the O(3) site. Although only average crystal structural models (space group I4/mmm) were used at this stage for both phases, it is evident that the occupancy of apical-position O(2) is higher in the γ-phase than that in the γ′-phase, while the trend for the disordered apical position O(4) is reversed, suggesting changes in the local structure between the γ- and γ′-phases. Additionally, the isotropic thermal displacements of atoms in the γ-phase are generally larger than in the γ′-phase, in most cases more than would be expected from increased thermal vibration at the higher temperature, particularly for O(3), reflecting increased structural disorder. The mean contact distances of Bi–O (2.3683(7) Å) and M–O (1.854(3) Å) at 25 °C increase significantly at 700 °C (Table S6), consistent with thermal expansion.

Short-range structure

The fitted neutron total scattering data (S(Q) and G(r)) for BICUVOX.10 at 25 °C and 700 °C are shown in Fig. S6a–d, with the derived partial pair distribution functions (PDFs) of M–O correlations displayed in Fig. 3b. In the γ′-phase (25 °C), the Bi–O correlation peak is broad and centered at ca. 2.24 Å, consistent with the short Bi–O(1) bond in the average model (ca. 2.33 Å, Table S6). In the γ-phase (700 °C), the Bi–O peak broadens further and becomes asymmetric, developing a shoulder at ca. 2.5 Å, due to the closer approach of O(4) to Bi. The mean and mode contact distances for M–O and O–O pairs are summarized in Table S7. As expected, the mode distance for Bi–O slightly decreases at 700 °C, while the mean distance increases from 2.318 Å at 25 °C to 2.334 Å at 700 °C, aligning with the average structural changes (Table S6). At 25 °C, Cu–O bond lengths are slightly longer than those of V–O, but both contract at 700 °C, indicating the development of a more tetrahedral geometry within the disordered γ-phase. Representative RMC configurations at 25 °C and 700 °C (Fig. S7) show significantly greater disorder of vanadate layers compared to the bismuthate layers, corroborated by the oxygen atom density distributions in the vanadate layer within the folded tetragonal cell (i.e., the RMC supercell folded back onto the crystallographic subcell) shown in Fig. 3c (25 °C). The O(2) and O(3) atoms display diffuse distributions, forming continuous nuclear density between the 4e and 8g sites in agreement with the strong positional modulations observed in the modulated structure.

Fig. 3d shows the coordination numbers (CNs) of atoms within their first coordination shell at 25 °C and 700 °C. The average CN of Bi slightly increases at 700 °C due to the closer contact of O(4) from the vanadate layer, while the average CNs of V, Cu and O decrease. The mean CN of V/Cu can be calculated as 0.9 × 4.97 + 0.1 × 4.57 = 4.93 at 25 °C and 0.9 × 4.7 + 0.1 × 4.35 = 4.67 at 700 °C. The method for theoretical CN calculation is detailed in our previous work.²⁵ For BICUVOX.10, the theoretical CN of V/Cu is 4.7. The higher mean CN at low temperature exceeds the theoretical value because, in the average structure, apical oxygen atoms are non-bridging, whereas in the local structure, these atoms become bridging due to severe displacement. At high temperature, the mean CN closely matches the theoretical value, reflecting the fully disordered γ-phase structure. Table S8 summarizes the occupancies of apical and equatorial oxygen atoms and the corresponding number of vacancies in the RMC models at 25 and 700 °C. The RMC models suggest that while at room temperature there are some apical vacancies (8.6% of apical sites are vacant), the number of apical vacancies decreases significantly at 700 °C (3.9%) while the equatorial site vacancy concentration shows the opposite trend. This suggests that, at elevated temperatures, oxide ions move from equatorial to apical positions, reducing the mean CN. The observed apical-to-equatorial (AV–EV) vacancy migration and enhanced equatorial vacancy population at high temperatures are crucial for facilitating continuous oxygen ion pathways, directly correlating with the improved ionic conductivity.

The structural modulation in the γ′-phase is closely linked to oxygen atom displacements with higher-dimensional periodicity. This phase is thus partially ordered, driven by specific oxygen vacancy distributions. In such an incommensurately modulated structure, the positional modulation functions describe the probability distribution of atomic positions within a higher-dimensional periodic potential landscape. In the γ′-phase, the significant positional modulation of oxide ions reflects spatially correlated displacements associated with short-range vacancy ordering. Oxygen diffusion occurs dynamically within this modulated energy landscape and the vacancy ordering confines oxygen motion, so increasing the effective diffusion barrier. Vacancy ordering was analyzed via the radial distribution of equatorial vacancies around M in the vanadate layer. Vacancy counts in each coordination shell were normalized to the first shell and compared against random distributions, yielding relative ratios (Fig. 3e). The y = 0 line indicates complete randomness. For reference, Fig. 3e also includes the initial random distributions at 25 and 700 °C. As expected, simulations which started from random distributions evolved distinct vacancy orderings at both temperatures. Notably, the γ′-phase exhibits stronger vacancy ordering than the γ-phase, underscoring the link between vacancy order and the partially ordered γ′-phase. Both phases show a negative deviation near 4 Å, indicating a depletion of vacancies close to M atoms. At this distance, two types of vacancy pairs are observed (Fig. S8): one across the MO_x polyhedron and another across adjacent voids. Analysis reveals that >90% of these pairs occur across void spaces rather than within a single polyhedron, suggesting that vacancy ordering preferentially aligns along 〈100〉 or 〈010〉 directions.

Conductivity measurements

The total conductivity for BICUVOX.10 during the first heating and cooling cycle was evaluated by a.c. impedance spectroscopy and the results (Fig. S9a) show a step-like increase in conductivity is observed around 500 °C upon heating, likely due to insufficient preheating of the platinum electrodes before measurement. During cooling, no abrupt change is detected; instead, conductivity decreases monotonically with decreasing temperature. A distinct difference in activation energy is observed between the high-temperature HT (450–800 °C) and low-temperatures LT (100–450 °C) regimes, with ΔE_HT = 0.397(6) eV and ΔE_LT = 0.625(5) eV, respectively. This difference is attributable to the γ ↔ γ′ transition, which, as shown by structural analyses, occurs near 450–500 °C. The ordering of oxide ions in the γ′-phase leads to a reduction in activation energy at lower temperatures. Ionic transference number measurements (Fig. S9b) confirm the nearly purely ionic nature of the conductivity in BICUVOX.10, particularly below 650 °C. At higher temperatures, a minor contribution from electronic conductivity becomes apparent but remains below 10% of the total conductivity.

Diffusion pathways

Ab initio molecular dynamics (AIMD) simulations were used to characterize the oxide-ion diffusion pathways. Fig. 4a shows the collapsed oxide-ion trajectories from a representative AIMD run. For clarity, a low-pass filter with a cutoff frequency of 2 THz was applied to suppress thermal vibrations and highlight long-range ion motion. The trajectories indicate that the vast majority of diffusion events correspond to jumps from apical to equatorial sites (A–E), as well as between neighboring equatorial sites (E–E). In contrast, diffusion involving O(1) sites is rare, with only a few instances observed across a total of 300 ps of AIMD data. A more detailed analysis, in which A–E and E–E jumps were distinguished explicitly, shows that A–E events occur approximately eight times more frequently than E–E events. Unfiltered trajectories were used to compute the diffusion coefficient, D = 2.8 × 10⁻¹⁰ m² s⁻¹, and the corresponding bulk conductivity, σ = 0.046 S cm⁻¹ (calculated using the Nernst–Einstein equation). The conductivity value compares well with an experimental value of 0.045 S cm⁻¹ for total conductivity at 800 °C.

Fig. 4 Oxide-ion diffusion pathway and mechanism. (a) Collapsed oxide-ion trajectories from a representative AIMD run, processed using a low-pass filter with a 2 THz cutoff to enhance visual clarity. Red lines indicate oxide-ion migration pathways, while purple and blue spheres represent Bi and V/Cu sites, respectively. Iso-surfaces of potential energy derived from (b) AIMD and (c) RMC models at 700 °C, respectively, showing the diffusion pathways for oxide ions involving the apical (A) and equatorial (E) sites (specified by local energy minima), with (d) and (e) the calculated energy barriers for oxide ion migration along the two pathways indicated by blue (from apical to equatorial sites) and green lines (between equatorial sites). (f) and (g) cross-sections of the potential energy distribution around Cu and V ions in (f) AIMD and (g) RMC models, taken along both [100] and [110] crystallographic directions. The color scale indicates the local potential energy level. (h) Proposed oxide-ion conduction mechanism in BICUVOX.10.

The ionic conductivity mechanism in BICUVOX.10 was investigated by analyzing the potential energy surfaces derived from atomistic configurations obtained using both RMC and DFT methods. Fig. 4b shows the oxide ion potential energy iso-surface calculated from AIMD trajectories, with that derived from the RMC model in Fig. 4c. Both apical and equatorial sites are clearly visible, with the latter further split into two subsites, which is consistent with the site splitting of O(3) in the crystallographic model (Fig. 2a). Energy barriers for ion migration from apical to equatorial sites, as well as between adjacent equatorial sites, were calculated and are presented in Fig. 4d and e for the AIMD and RMC models, respectively. The results indicate that the energy barrier for jumps between apical and equatorial sites is lower, by 0.2 eV in the DFT model and 0.05 eV in the RMC model, compared to that between neighboring equatorial sites. The discrepancies in the energy profiles are mainly due to differences in the energy for barrier 1, with reasonably close agreement for barrier 2. The differences are attributable to the differences in approaches used in the two methods. The AIMD method basically yields a collection of snapshots of ion positions in a small supercell as a function of time and reflects the ion dynamics, while the RMC method gives a time averaged picture of a large collection of ions reflecting variation in the local structure. The data sets used in the two methods are thus very different, i.e. energetic information in the case of AIMD and diffraction data in the case of RMC. Thus, it is not surprising that there are differences in the absolute values of energies, but the fact that the same basic trends were observed in both methods is remarkable. These findings suggest that ionic transport is dominated by jumps between apical and equatorial positions. This is in agreement with studies on the parent compound Bi₄V₂O₁₁.¹⁸ Additional calculations focused on vanadium- and copper-centered coordination environments reveal that the energy barriers for oxide ion migration are significantly lower in the vicinity of vanadium atoms compared to copper. This trend is most pronounced in the DFT-based data but is also evident in the results derived from RMC calculations. Representative cross-sections of the corresponding potential energy surfaces are shown in Fig. 4f and g. More importantly, the results suggest that while introduction of copper into bismuth vanadate introduces charge carriers in the form of oxide ion vacancies, diffusion of vacancies around Cu is more limited with these atoms acting as vacancy traps, which is reflected in lower oxide ion coordination numbers in Fig. 3d. Thus, a conduction mechanism in BICUVOX.10 is proposed in Fig. 4g, which is more complex than the original model,³⁵ involving a more zig-zag like pathway between equatorial and apical oxide ion sites that facilitates long range diffusion. In addition, although site splitting above and below the equatorial plane is visible in the RMC model at room temperature (Fig. 4h), it is much less so at 700 °C, reflecting the slightly reduced value of the z-coordinate and increased isotropic thermal parameter for O(3) at the higher temperature in the crystallographic model (Table S5).

The current results collectively indicate that in the γ′-phase, in which short-range vacancy ordering is present, the available oxygen sites for diffusion are decreased in number, which in turn reduces long-range diffusion pathway connectivity, thereby suppressing diffusion. Upon heating across the γ′ ↔ γ transition, the structural modulation weakens and the vacancy order melts, leading to a smoother and more percolative energy landscape that enables extended zig-zag diffusion pathways. Both γ- and γ′-phases show zig-zag pathways between equatorial and apical oxide ion sites, with oxide ions preferentially diffusing around V atoms rather than Cu, because of the vacancy trapping effect of the latter.

Conclusions

In this work, we have successfully clarified the complex structural characteristics and phase evolution of Bi₂V_0.9Cu_0.1O_5.35 (BICUVOX.10), with a particular focus on the transition between the ordered γ′-phase and the disordered γ-phase. For the first time, the γ′-phase superlattice has been fully resolved as a (3 + 2)D modulated structure in superspace group I4/mmm(a,a,1)0s0s(−a,a,1)0s00 with 2D incommensurate modulation vectors. Comprehensive neutron/X-ray diffraction and total scattering analyses revealed that positional and occupational modulations of oxygen atoms within the vanadate layers are key to the formation of diverse local coordination geometries and drive the disorder–order phase transition. Temperature-dependent studies demonstrated a reversible γ ↔ γ′ transition near 500 °C accompanied by measurable changes in lattice parameters, Bi–M interlayer distances, and local vacancy distributions. Notably, short-range ordering of oxygen vacancies was clearly observed in the γ′-phase but diminished in the γ-phase, highlighting its critical role in modulating ionic conductivity. The AV-EV transformation enhances continuous oxide ion pathways and therefore facilitates oxide ion conduction. Conductivity measurements confirmed that high oxide-ion conductivity is dominated by ionic transport below 650 °C, with a clear reduction in activation energy in the disordered γ-phase.

The AIMD studies reveal that oxide ion diffusion is predominantly 2-dimensional but involves both equatorial and apical oxide ion sites in a zig-zag like pathway within the vanadate/cuprate layer, with Cu atoms acting as defect traps. Bearing in mind that a similar pathway has been reported in the unsubstituted compound, it is likely that this mechanism is common to most BIMEVOX systems. Two main differences are seen in the present study on the Cu substituted system compared to previous work on the unsubstituted parent compound; firstly, that equatorial to equatorial jumps in the vanadate layer appear to be more common in the Cu doped system and secondly that the Cu dopant causes a narrowing of the diffusion pathways leading to reduced diffusion around the dopant.

The present work has demonstrated, using a comprehensive range of computational and experimental techniques, that a detailed atomistic understanding of the interplay between long- and short-range order can be achieved even in a complex incommensurately modulated system such as a BIMEVOX. The modulated structure can be understood as the average structural representation of an underlying dynamic process. The modulation encodes the spatial organization of oxide ion/vacancy distributions, which directly governs the pathways and energetics of oxide ion diffusion. The work underscores the impact of both local order and structure modulation on ion diffusion in solids and provides a methodology for interpreting these features and assessing their significance. By establishing a direct link between lattice modulation, defect distribution, and ionic transport, this work promotes a shift in the design concept of ionic conductors from tuning chemical composition toward deliberate control of higher dimensional dynamic structural order.

Author contributions

Y. Yue, P. Miao, V. Petricek, and I. Abrahams conceptualized the idea. F. Qu, A. Dzięgielewska and A. Koedtruad performed investigations. Y. Yue, F. Qu, M. Krynski, and V. Petricek carried out formal analysis. Y. Yue, V. Petricek, M. Krynski, P. Miao and I. Abrahams contributed to funding acquisition. Y. Yue and F. Qu finished the original draft and the manuscript was revised by M. Malys, S. Hull, Q. F. Zheng, F. Krok, V. Petricek, I. Abrahams and P. Miao. All authors have discussed the results and approved the submission.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data collected at the ISIS facility are available at https://doi.org/10.5286/ISIS.E.24089876. Crystal structure data are available at the Bilbao Incommensurate Crystal Structure Database (deposition code: Yz0X5FUw02Z). Other data relating to this work are available from the authors.

Supplementary information (SI): figures, tables and videos. See DOI: https://doi.org/10.1039/d6ta00458j.

Acknowledgements

This work was financially sponsored by the National Natural Science Foundation of China (No. 12305343), the Guangdong Basic and Applied Basic Research Foundation (2022A1515110210, 2022B1515120014, 2023B0303000003, 2023B1515120060) and the Guangdong Innovative Entrepreneurial Research Team Program (2021ZT09C539). Crystallographic analysis at the Institute of Physics of the Czech Academy of Sciences was supported by MGML (mgml.eu) as part of the Czech Research Infrastructures program (No. LM2023065). The authors greatly appreciate the Science and Technology Facilities Council (STFC) for beamtime awarded on Polaris at ISIS (RB1110197), and the Japan Proton Accelerator Research Complex (J-PARC) for beamtime on NOVA (No. 2022B0049) and superHRPD (No. 2023A0143). This work was supported by the LAB-TECH of Excellence grant under the Excellence Initiative – Research University at Warsaw University of Technology.

References

1. T. Janssen, G. Chapuis and M. d. Boissieu, Aperiodic Crystals from Modulated Phases to Quasicrystals, Oxford University Press Inc., New York, United States, 2007.
2. A. Y. Borisevich, E. Eliseev, A. Morozovska, C.-J. Cheng, J.-Y. Lin, Y.-H. Chu, D. Kan, I. Takeuchi, V. Nagarajan and S. V. Kalinin, Nat. Commun., 2012, 3, 775.
3. Z. Sun, J. Li, C. Ji, J. Sun, M. Hong and J. Luo, J. Am. Chem. Soc., 2017, 139, 15900–15906.
4. D. Rusu, J. J. Peters, T. P. Hase, J. A. Gott, G. A. Nisbet, J. Strempfer, D. Haskel, S. D. Seddon, R. Beanland and A. M. Sanchez, Nature, 2022, 602, 240–244.
5. J. Liang, J. Liu, P. Qiu, C. Ming, Z. Zhou, Z. Gao, K. Zhao, L. Chen and X. Shi, Nat. Commun., 2023, 14, 8442.
6. J.-W. G. Bos, H. Zandbergen, M.-H. Lee, N. Ong and R. Cava, Phys. Rev. B: Condens. Matter, 2007, 75, 195203.
7. M. Roger, D. Morris, D. Tennant, M. Gutmann, J. Goff, J.-U. Hoffmann, R. Feyerherm, E. Dudzik, D. Prabhakaran and A. Boothroyd, Nature, 2007, 445, 631–634.
8. J. Li, F. Pan, S. Geng, C. Lin, L. Palatinus, M. Allix, X. Kuang, J. Lin and J. Sun, Nat. Commun., 2020, 11, 4751.
9. C. Li, S. S. Pramana, R. D. Bayliss, C. P. Grey, F. Blanc and S. J. Skinner, Chem. Mater., 2020, 32, 2292–2303.
10. J. Wind, J. Polt, Z. Zhang, D. A. Blom, T. Vogt, R. L. Withers and C. D. Ling, Chem. Mater., 2017, 29, 9171–9181.
11. F. S. Hempel, W. A. Sławiński, B. Arstad and H. Fjellvåg, Chem. Mater., 2024, 36, 11084–11094.
12. J. B. Goodenough, Nature, 2000, 404, 821–823.
13. M. Yashima, T. Tsujiguchi, Y. Sakuda, Y. Yasui, Y. Zhou, K. Fujii, S. Torii, T. Kamiyama and S. J. Skinner, Nat. Commun., 2021, 12, 556.
14. A. Dziegielewska, M. Malys, W. Wróbel, S. Hull, Y. Yue, F. Krok and I. Abrahams, Solid State Ionics, 2021, 360, 115543.
15. Y. Yue, A. Dzięgielewska, F. Krok, R. M. Whiteley, H. Toms, M. Malys, H. Yan and I. Abrahams, J. Phys. Chem. C, 2022, 126, 2108–2120.
16. I. Abrahams and F. Krok, J. Mater. Chem., 2002, 12, 3351–3362.
17. B. Singh, S. Ghosh, S. Aich and B. Roy, J. Power Sources, 2017, 339, 103–135.
18. A. Agnaou, W. Mhaira, R. Essalim, A. Turino, A. Ammar, M. Zamama, M. Alga and F. Mauvy, J. Solid State Chem., 2023, 318, 123730.
19. H. J. Stroud, C. E. Mohn, J. A. Hernandez and N. L. Allan, Philos. Trans.: Math., Phys. Eng. Sci., 2021, 379, 20200430.
20. Y. Yue, A. Dzięgielewska, M. Zhang, S. Hull, F. Krok, R. M. Whiteley, H. Toms, M. Malys, X. Huang and M. Krynski, Chem. Mater., 2022, 35, 189–206.
21. F. Abraham, J. Boivin, G. Mairesse and G. Nowogrocki, Solid State Ionics, 1990, 40, 934–937.
22. E. Pernot, M. Anne, M. Bacmann, P. Strobel, J. Fouletier, R. Vannier, G. Mairesse, F. Abraham and G. Nowogrocki, Solid State Ionics, 1994, 70, 259–263.
23. P. Kurek, P. Pongratz and M. Breiter, Solid State Ionics, 1998, 113, 615–621.
24. F. Krok, I. Abrahams, D. Bangobango, W. Bogusz and J. Nelstrop, Solid State Ionics, 1996, 86, 261–266.
25. Y. Yue, A. Dzięgielewska, S. Hull, F. Krok, R. M. Whiteley, H. Toms, M. Malys, M. Zhang, H. Yan and I. Abrahams, J. Mater. Chem. A, 2022, 10, 3793–3807.
26. C. Muller, M. Anne and M. Bacmann, Solid State Ionics, 1998, 111, 27–36.
27. R. Vannier, G. Mairesse, G. Nowogrocki, F. Abraham and J. Boivin, Solid State Ionics, 1992, 53, 713–722.
28. R. Vannier, G. Mairesse, F. Abraham and G. Nowogrocki, J. Solid State Chem., 1993, 103, 441–446.
29. S. Lazure, R.-N. Vannier, G. Nowogrocki, G. Mairesse, C. Muller, M. Anne and P. Strobel, J. Mater. Chem., 1995, 5, 1395–1403.
30. C. K. Lee, G. S. Lim and A. R. West, J. Mater. Chem., 1994, 4, 1441.
31. P. De Wolff, Acta Crystallogr., Sect. A:Found. Crystallogr., 1974, 30, 777–785.
32. I. Abrahams, F. Krok and J. Nelstrop, Solid State Ionics, 1996, 90, 57–65.
33. I. Abrahams, J. Nelstrop, F. Krok and W. Bogusz, Solid State Ionics, 1998, 110, 95–101.
34. I. Abrahams, F. Krok, M. Malys and A. Bush, J. Mater. Sci., 2001, 36, 1099–1104.
35. P. Kurek, J. Dygas and M. Breiter, J. Electroanal. Chem., 1994, 378, 77–83.
36. V. Petříček, L. Platinus, J. Plášil and M. Dušek, Z. Kristallogr. - Cryst. Mater., 2023, 238, 271–282.
37. A. C. Larson and R. B. Von Dreele, GSAS Generalised Structure Analysis System, Report Los Alamos National Laboratory Report LAUR-86-748, Los Alamos National Laboratory, 1987.
38. B. H. Toby, J. Appl. Crystallogr., 2001, 34, 210–213.
39. M. Y. R. Oishi, Y. Nishimaki, S. Torii, A. Hoshikawa, T. Ishigaki, T. Morishima, K. Mori and T. Kamiyama, Nucl. Instrum. Methods Phys. Res., Sect. A, 2009, 600, 94–96.
40. A. K. Soper, GudrunN and GudrunX: programs for correcting raw neutron and X-ray diffraction data to differential scattering cross section, Report Rutherford Appleton Laboratory Technical Report: RAL-TR-2011- 013, Science & Technology Facilities Council: Swindon, U.K., 2011.
41. M. G. Tucker, D. A. Keen, M. T. Dove, A. L. Goodwin and Q. Hui, J. Phys.: Condens. Matter, 2007, 19, 335218–335234.
42. J. H. G. Kresse, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558.
43. J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys., 1996, 105, 9982–9985.
44. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
45. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
46. I. Abrahams, F. Krok, M. Malys and A. J. Bush, J. Mater. Sci., 2001, 36, 1099–1104.

---