*Article* **Fingerprints of Topotactic Hydrogen in Nickelate Superconductors**

**Liang Si 1,2,\*, Paul Worm <sup>2</sup> and Karsten Held 2,\***


**Abstract:** Superconductivity has entered the nickel age marked by enormous experimental and theoretical efforts. Notwithstanding, synthesizing nickelate superconductors remains extremely challenging, not least due to incomplete oxygen reduction and topotactic hydrogen. Here, we present density-functional theory calculations for nickelate superconductors with additional topotactic hydrogen or oxygen, namely La1−*x*Sr*x*NiO2H*<sup>δ</sup>* and LaNiO2<sup>+</sup>*δ*. We identify a phonon mode as a possible indication for topotactic hydrogen and discuss the charge redistribution patterns around oxygen and hydrogen impurities.

**Keywords:** superconductivity; nickelates; strongly correlated electron systems

#### **1. Introduction**

Computational materials calculations have predicted superconductivity in nickelates [1] and the heterostructures thereof [2–4] since many decades, mainly based on apparent similarly to cuprate superconductors. Three years ago, superconductivity in nickelates was finally discovered in an experiment by Li, Hwang, and coworkers [5], breaking the grounds for a new age of superconductivity, the nickel age. It is marked by an enormous theoretical and experimental activity, including but not restricted to [5–31]. Superconductivity has been found by now, among others, in Nd1−*x*Sr*x*NiO2 [5,6], Pr1−*x*Sr*x*NiO2 [7], La1−*x*Ca*x*NiO2 [8], La1−*x*Sr*x*NiO2 [9], and most recently in the pentalayer nickelate Nd6Ni5O12 [10]. Figure 1 shows some of the hallmark experimental critical temperatures (*Tc*'s) for the nickelates in comparison with the preceding copper [32] and iron age [33] of unconventional superconductivity. Also shown are some other noteworthy superconductors, including the first superconductor, solid Hg, technologically relevant NbTi, and hydride superconductors [34]. The last are superconducting at room temperature [35], albeit only at a pressure of 267GPa exerted in a diamond anvil cell. All of these compounds are marked in gray in Figure 1 as they are conventional superconductors. That is, the pairing of electrons originates from the electron-phonon coupling, as described in the theory of Bardeen, Cooper, and Schrieffer (BCS) [36].

In contrast, cuprates, nickelates, and, to a lesser extent, iron pnictides are strongly correlated electron systems with a large Coulomb interaction between electrons because of their narrow transition metal orbitals. Their *Tc* is too high for BCS theory [37,38], and the origin of superconductivity in these strongly correlated systems is still hotly debated. One prospective mechanism is antiferromagnetic spin fluctuations [39–43] stemming from strong electronic correlations. Another mechanism is based on charge density wave fluctuations and received renewed interest with the discovery of charge density wave ordering in cuprates [44,45]. Dynamical vertex approximation [46–49] calculations for nickelates [27], which are unbiased with respect to charge and spin fluctuations, found that spin fluctuations dominated and successfully predicted the superconducting dome prior to experiment in Nd1−*x*Sr*x*NiO2 [6,50,51].

**Citation:** Si, L.; Worm, P.; Held, K. Fingerprints of Topotactic Hydrogen in Nickelate Superconductors. *Crystals* **2022**, *12*, 656. https:// doi.org/10.3390/cryst12050656

Academic Editor: Xiaoguan Zhang

Received: 19 April 2022 Accepted: 30 April 2022 Published: 4 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** Superconducting *Tc* vs. year of discovery for selected superconductors. The discovery of cuprates, iron pnictides and nickelates led to enormous experimental and theoretical activities. Hence, one also speaks of the copper, iron and nickel age of superconductivity.

Why did it take 20 years to synthesize superconducting nickelates that have been so seemingly predicted on a computer? To mimic the cuprate Cu 3*d*<sup>9</sup> configuration, as in NdNiO2, nickel has to be in the uncommon oxidation state Ni1<sup>+</sup>, which is rare and prone to oxidize further. Only through a complex two-step procedure, Lee, Hwang, and coworkers [52] were able to synthesize superconducting nickelates. In a first step, modern pulsed laser deposition (PLD) was used to grow a Sr*x*Nd1−*<sup>x</sup>*NiO3 film on a SrTiO2 substrate. This nickelate is still in the 3D perovskite phase—see Figure 2 (left)—with one oxygen atom too many and will thus not show superconductivity. Hence, this additional oxygen between the layers needs to be removed in a second step. The reducing agent CaH2 is used to this end within a quite narrow temperature window [52]. If all goes well, one arrives at the superconducting Sr*x*Nd1−*<sup>x</sup>*NiO2 film (top center). However, this process is prone to incomplete oxidation or to intercalate hydrogen topotactically, i.e., at the position of the removed oxygen; see Figure 2 (bottom center). Both of those unwanted outcomes are detrimental for superconductivity.

In [21,53,54], it was shown by density functional theory (DFT) calculations that NdNiO2H is indeed energetically favorable to NdNiO2 + 1/2 H. For the doped system, on the other hand, Nd0.8Sr0.2NiO2 without the hydrogen intercalated is energetically favorable. The additional H or likewise an incomplete oxidation to SrNdNiO2.5 alters the physics completely. Additional H or O0.5 will remove an electron from the Ni atoms, resulting in Ni2<sup>+</sup> instead of Ni1<sup>+</sup>. The formal electronic configuration is hence 3*d*<sup>8</sup> instead of 3*d*9, or two holes instead of one hole in the Ni *d*-shell. Dynamical mean-field theory (DMFT) calculations [21] evidence that the basic atomic configuration is the one of Figure 2 (lower right). That is, because of Hund's exchange, the two holes in NdNiO2H occupy two different orbitals, 3*dx*<sup>2</sup>*y*<sup>2</sup> and 3*d*3*z*2−*r*<sup>2</sup> , and form a spin-1. A consequence of this is that DMFT calculations predict NdNiO2H to be a Mott insulator, whereas NdNiO2 is a strongly correlated metal with a large mass enhancement of about five [21].

To the best of our knowledge, such a two-orbital, more 3D electronic structure is unfavorable for high-*Tc* superconductivity. The two-dimensionality of cuprate and nickelate superconductors helps to suppress long-range antiferromagnetic order, while at the same time retaining strong antiferromagnetic fluctuations that can act as a pairing glue for superconductivity. In experiment, we cannot expect ideal NdNiO2, NdNiO2H or NdNiO2.5 films, but most likely some H or additional O will remain in the NdNiO2 film, after the CaH2 reduction. Additional oxygen can be directly evidenced in standard X-ray diffraction

analysis after the synthesis step. However, hydrogen, being very light, evades such an X-ray analysis. It has been evidenced in nickelates only by nuclear magnetic resonance (NMR) experiments [55] which, contrary to X-ray techniques, are very sensitive to hydrogen. Ref. [56] suggested hydrogen in LaNiO2 to be confined at grain boundaries or secondary-phase precipitates. Given these difficulties, it is maybe not astonishing that it took almost one year before a second research group [6] was able to reproduce superconductivity in nickelates. Despite enormous experimental efforts, only a few groups have succeeded hitherto.

**Figure 2.** For synthesizing superconducting nickelates (1, **left**), a perovskite film of Nd(La)1−*x*Sr*x*NiO3 is grown on a SrTiO3 substrate, and (2, **center**) the O atoms between the planes are removed by reduction with CaH2. Besides the pursued nickelate Nd(La)1−*x*Sr*x*NiO2 (top center), also excess oxygen or topotactic H may remain in the film, yielding Nd(La)1−*x*Sr*x*NiO2H (**bottom center**). The excess hydrogen results in two holes instead of one hole within the topmost two Ni 3*d* orbitals (**right**). Adapted from Ref. [57].

In this paper, we present additional DFT results for topotactic hydrogen and incomplete oxygen reduction in nickelate superconductors: In Section 3, we provide technical information on the DFT calculations. In Section 3, we analyze the energy gain to topotactically intercalate hydrogen in LaNiO2 and NdNiO2. In Section 4, we analyze the phonon spectrum and identify a high-energy mode originating from the Ni-H-Ni bond as a characteristic feature of intercalated hydrogen. In Section 5, we show the changes of the charge distribution caused by topotactic hydrogen or oxygen. Finally, Section 6 provides a summary and outlook.

#### **2. Method**

*Computational details on Eb.* In both our previous theoretical study [21] and this article, the binding energy *Eb* of hydrogen atoms is computed as:

$$E\_b = E[AB\mathbf{O}\_2] + \mu[\mathbf{H}] - E[AB\mathbf{O}\_2\mathbf{H}].\tag{1}$$

Here, *E*[*AB*O2] and *E*[*AB*O2H] are the total energy of infinite-layer *AB*O2 and hydrideoxides *AB*O2H, while *μ*[H] = *E*[H2]/2 is the chemical potential of H. Note that H2 is a typical byproduct for the reduction with CaH2 and also emerges when CaH2 is in contact with H2O. Hence, it can be expected to be present in the reaction. A positive (negative) *Eb* indicates the topotactic H process is energetically favorable (unfavorable) to obtain *AB*O2H instead of *AB*O2 and H2/2.

In the present paper, we go beyond [21] that reported *Eb* of various *AB*O2 compounds by investigating *Eb* of La1−*x*Ca*x*NiO2 systems for many different doping levels. Here, the increasing Ca-doping is achieved by using the virtual crystal approximation (VCA) [58,59] from LaNiO2 (*x* = 0) to CaNiO2 (*x* = 1). For each Ca concentration, structure relaxation and static total energy calculation is carried out for La1−*x*Ca*x*NiO2 and La1−*x*Ca*x*NiO2H within the tetragonal space group *P*4/mmm. To this end, we use densityfunctional theory (DFT) [60,61] with the VASP code [62,63] and the generalized gradient approximations (GGA) of Perdew, Burke, and Ernzerhof (PBE) [64], and PBE revised for solids (PBEsol) [65]. For undoped LaNiO2, the GGA-PBEsol relaxations predict its in-plane lattice constant as 3.890 Å, which is close to that of the STO substrate: 3.905 Å. The computations for La1−*x*Ca*x*NiO2 and LaCoO2, LaCuO2, SrCoO2, and SrNiO2 are performed without spin-polarization and a DFT+*U* treatment [66], as the inclusion of Coulomb *U* and spin-polarization only slightly decrease the *Eb* by ∼5% for LaNiO2 [57]. For NdNiO2, an inevitably computational issue is the localized Nd-4 *f* orbitals. These *f*-orbitals are localized around the atomic core, leading to strong correlations. In non-spinpolarized DFT calculations, this generates flat bands near the Fermi level E*<sup>F</sup>* and leads to unsuccessful convergence. To avoid this, we employed DFT+*U* [*Uf*(Nd) = 7 eV and *Ud*(Ni) = 4.4 eV] and initialized a *G*-type anti-ferromagnetic ordering for both Nd- and Ni-sublattice in a <sup>√</sup><sup>2</sup> <sup>×</sup> <sup>√</sup><sup>2</sup> <sup>×</sup> <sup>2</sup> supercell of NdNiO2. For the Nd0.75Sr0.25NiO2 case, 25% Sr-doping is achieved by replacing one out of the four Nd atoms by Sr in a <sup>√</sup><sup>2</sup> <sup>×</sup> <sup>√</sup><sup>2</sup> <sup>×</sup> <sup>2</sup> NdNiO2 supercell.

*Computational details on phonons.* The phonon computations for LaNiO2, LaNiO2H, LaNiO2H0.125, and LaNiO2.125 are performed with the frozen phonon method using the PHONONY [67] code interfaced with VASP. Computations with density functional perturbation theory (DFPT) method [68] are also carried out for double checking. For LaNiO2 and LaNiO2H, the unit cells shown in Figure 3a,b are enlarged to a 2 × 2 × 2 supercell, while for LaNiO2H0.125 and LaNiO2.125, the phonons are directly computed with the supercell of Figure 3c,d.

*Computational details on electron density.* The electron density distributions of LaNiO2, LaNiO2H, LaNiO2H0.125, and LaNiO2.125 are computed using the WIEN2K code [69] while taking the VASP-relaxed crystal structure as input. The isosurfaces are plotted from 0.1 (yellow lines) to 2.0 (center of atoms) with spacing of 0.1 in units of e/Å2.

**Figure 3.** *Cont*.

**Figure 3.** Phonon spectra of (**a**) LaNiO2 and (**b**) LaNiO2H and in a 2 × 2 × 2 LaNiO2 supercell doped with a single (**c**) H and (**d**) O atom (i.e., LaNiO2H0.125 in (**c**) and LaNiO2.125 in (**d**)). The orange and black arrows in (**b**,**d**) represent vibrations of H and O atoms. The blue dashed oval in (**d**) labels the unstable phonon modes induced by intercalating additional O atoms in LaNiO2.

#### **3. Energetic Stability**

Figure 4 shows the results of the hydrogen binding energy *Eb* for the infinite layer nickelate superconductors Nd1−*x*Sr*x*NiO2 [5,6,50] and La1−*x*Ca*x*NiO2 [8]. To reveal the evolution of *Eb* when the *B*-site band filling deviates from their original configurations (3*d*<sup>9</sup> in LaNiO2 when *x* = 0 and 3*d*<sup>8</sup> in CaNiO2 when *x* = 1), we also show the binding energy of LaCoO2 (3*d*8), LaCuO2 (3*d*10), SrCoO2 (3*d*7), and SrNiO2 (3*d*8).

**Figure 4.** Hydrogen binding energy (*Eb*) per hydrogen in two nickelate superconductors, La1−*x*Ca*x*NiO2 and Nd1−*x*Sr*x*NiO2 vs. Sr/Ca doping concentration *x*; LaCuO2, LaCoO2, and SrCoO2 are shown for comparison. Slightly above 10% (Sr,Ca)-doping infinite layer nickelates are energetically more stable. Note that the doping changes the filling of the *B*-3*d* orbital. To study the relationship between *Eb* and the types of *B*-site elements, *Eb* of several other *AB*O2 compounds is computed: LaCoO2, LaCuO2, SrCoO2 and SrNiO2. Note that this changes the filling of *B*-3*d* orbital within a large range: e.g., 3*d*<sup>8</sup> for LaCoO2 and 3*d*<sup>9</sup> for LaNiO2.

Let us start with the case of La1−*x*Ca*x*NiO2 [8]. Here, the unoccupied La-4 *f* orbitals make the computation possible even without spin-polarization and Coulomb *U* for La-4 *f* , whereas for NdNiO2, this is not practicable due to Nd-4 *f* flat bands near E*F*. Positive (negative) *Eb* above (below) the horizontal line in Figure 4 indicates topotactic H is energetically favorable (unfavorable). When *x* = 0, i.e., for bulk LaNiO2, the system tends to confine H atoms, resulting in oxide-hydride *AB*O2H with *Eb* = 157 meV/H. As the concentration of Ca increases, *Eb* monotonously decreases, reaching −248 meV for the end member of the doping series, CaNiO2. The turning point between favorable and unfavorable topotactic H inclusion is around 10% to 15% Ca-doping. Let us note that *Eb* = 0 roughly agrees with the onset of superconductivity, which for Ca-doped LaNiO2 emerges for *x* > 15% Ca-doping [8].

To obtain *Eb* in NdNiO2 a much higher computational effort is required: firstly, the Nd-4 *f* orbitals must be computed with either treating them as core-states or including spin-splitting. Secondly, for the spin-polarized DFT(+*U*) calculations, an appropriate (anti-)ferromagnetic ordering has to be arranged for both Ni- and Nd-sublattices. In oxide-hydride *AB*O2H compounds, the *δ*-type bond between Ni and H stabilizes a *G*type anti-ferromagnetic order by driving the system from a quasi two-dimensional (2D) system to a three dimensional (3D) one [21]. Given the large computational costs of *Eb* for Nd1−*x*Sr*x*NiO2 by using anti-ferromagnetic DFT+*U* calculations for both Nd-4 *f* (*U* ∼7 eV) and Ni-3*d* (*U* = 4.4 eV) orbitals, we merely show here the results of NdNiO2 (*x* = 0), Nd0.75Sr0.25NiO2 (*x* = 0.25), and SrNiO2 (*x* = 1), which are adopted from [21]. With 25% Sr-doping, the *Eb* of NdNiO2 is reduced from 134 meV to −113 meV. Please note that *Eb* of (Nd,Sr)NiO2 is slightly smaller than in (La,Ca)NiO2, at least in the low doping range. This can be explained by shorter lattice constants in NdNiO2, in agreement with the finding [21] that compressive strain plays an important role in reducing *Eb*.

One can speculate that this suppression of topotactic hydrogen may also play a role when comparing the recently synthesized (Nd,Sr)NiO2 films on a (LaAlO3)0.3(Sr2TaAlO6)0.7 (LSAT) substrate [51] with the previously employed SrTiO3 (STO) substrate [50]. Lee et al. [51] reported cleaner films without defects and also a higher superconducting transition temperature *Tc* ∼20 K for the LSAT film as compared to *Tc* = 15 K and plenty of stacking fault defects for the STO substrate [50]. As for (La,Ca)NiO2, *Eb* = 0 falls in the region of the onset of the superconductivity for (Sr,Nd)NiO2, which is *x* ∼10% Sr-doping in LSAT-strained defect-free films [51] and *x* ∼12.5% at SrTiO3-substrate states [50]. Topotactic hydrogen might play a role in suppressing superconductivity in this doping region.

In Figure 4, we further show additional infinite layer compounds LaCoO2, LaCuO2, SrCoO2, and SrNiO2 for comparison. Their *Eb* is predicted to be 367, −42, 69, and −134 meV, respectively. Combining the results of LaNiO2 and CaNiO2, we summarize several tendencies on how to predict *Eb* of *AB*O2: (1) the strongest effect on *Eb* is changing the *B*-site element. However, this seems unpractical for nickelate superconductors as the band filling is strictly restricted to be 3*d*9−*<sup>x</sup>* (*<sup>x</sup>* ∼0.2). For both trivalent (La, Nd) and bivalent (Sr, Ca) cations, *Eb* decreases when the *B*-site cation goes from early to late transition metal elements, e.g., from LaCoO2 (3*d*8) to LaNiO2 (3*d*9) to LaCuO2 (3*d*10). (2) Compressive strains induced by either substrate or external pressure can effectively reduce *Eb*, and we believe that this might be used for growing defect-free films. (3) According to our theoretical calculations, *Eb* mainly depends on lattice parameters and band filling of the *B*-site 3*d*-orbitals, but much less on magnetic ordering and Coulomb interaction *U*.

#### **4. Phonon Dispersion**

As revealed by previous DFT phonon spectra calculations [16], NdNiO2 is dynamically stable. One of the very fundamental questions would be whether topotactic H from overreacted reduction and/or O from unaccomplished reductive reactions affect the lattice stability. To investigate this point, we perform DFT phonon calculations and analyze the lattice vibration induced by H/O intercalation, as shown in Figure 3.

The phonon spectrum of LaNiO2 (Figure 3a) is essentially the same as in Ref. [16]; all the phonon frequencies are positive, indicating it is dynamically stable. Its upmost optical phonon at around 14 to 16 THz can be identified with the recent experimental resonant inelastic x-ray scattering (RIXS) data [70] showing a weakly dispersing optical phonon at ∼60 meV ≈ 15 THz. In Figure 3b, the oxides-hydride LaNiO2H is also predicted to be dynamically stable. Please note that the phonon dispersions between 0 and 20 THz are basically the same as those in LaNiO2 (Figure 3a; note the different scale of the *y*axis). However, one can see new, additional vibration modes from the light H-atoms at frequencies of ∼27 THz and ∼43 THz. Among these vibrations, the double degenerate mode at lower frequency is generated by an in-plane (*xy*-plane) vibration of the topotactic H atom. There are two such in-plane vibrations of H atoms, either along the (100) or (110) direction (and symmetrically related directions), as indicated by the orange arrows in Figure 3b. The mode located at the higher frequency ∼43 THz is, on the other hand, formed by an out-of-plane (*z*-direction) vibration and is singly degenerate.

We explain these phonon modes in detail by computing the bonding strength between H-1*s*–Ni-*dz*<sup>2</sup> and H-1*s*–La-*dxy* orbitals. Our tight-binding calculations yield an electron hopping term of −1.604 eV between H-1*s* and Ni-*dz*<sup>2</sup> , while it is −1.052 eV from La-*dxy* to H-1*s*. That is, the larger H-1*s*–Ni-*dz*<sup>2</sup> overlap leads to a stronger *δ*-type bonding and, together with the shorter *c*-lattice constant, to a higher phonon energy. Additionally, the shorter *c*-lattice in LaNiO2 should also play a role at forming a stronger H-1*s*–Ni-*dz*<sup>2</sup> bond.

In our previous analysis of the band character for LaNiO2H [21], the H-1*s* bands were mainly located at two energy regions: a very flat band that is mostly from the H-1*s* itself at ∼−7 to −6 eV, and a hybridized band between H-1*s* and Ni-*dz*<sup>2</sup> at ∼−2 eV. Together with the higher phonon energy, this indicates that the topotactic H atoms are mainly confined by a Ni sub-lattice via bonding and anti-bonding states formed by H-1*s* and Ni-*dz*<sup>2</sup> orbitals, instead of the La(Nd) sub-lattice.

The complete (full) topotactic inclusion of H, where all vacancies induced by removing oxygen are filled by H, is an ideal limiting case. Under varying experimental conditions, such as chemical reagent, substrate, temperature, and strain, the H-topotactic inclusion may be incomplete, and thus *AB*O2H*<sup>δ</sup>* (*δ* < 1) may be energetically favored. Hence, we also compute the phonon spectrum at a rather low H-topotactic density: LaNiO2H0.125, achieved by including a single H into 2 × 2 × 2 LaNiO2 supercells as shown in Figure 3c. Moreover, such a local H defect, as revealed by the positive frequency at all *q*-vectors in the lower panel of Figure 3c, does not destroy the dynamical stability of the LaNiO2 crystal. In fact, the only remarkable qualitative difference between the complete and 12.5% topotactic H case is the number of phonon bands at 0 THz to 20 THz. This is just a consequence of the larger 2 × 2 × 2 LaNiO2 supercell, with eight times more phonons. Some quantitative differences can be observed with respect to the energy of the phonon mode: The out-ofplane vibration energy is enhanced from ∼43 THz in LaNiO2H (Figure 3b) to ∼47 THz in LaNiO2H0.125 (Figure 3b), and the in-plane vibration mode frequency is reduced from ∼27 THz in LaNiO2H (Figure 3b) to ∼21 THz in LaNiO2H0.125 (Figure 3c). This is because the H-intercalation shrinks the local *c*-lattice, i.e., the distance between two Ni atoms separated by topotactic H, from 3.383 Å in (LaNiO2H: Figure 3b) to 3.327 Å (LaNiO2H0.125: Figure 3c). The bond length between H and La is, on the other hand, slightly increased from 2.767 Å in (LaNiO2H: Figure 3b) to 2.277 Å (LaNiO2H0.125: Figure 3c). This lattice compression (enlargement) explains the enhancement (reduction) for the out-of-plane (in-plane) phonon frequencies (energies).

These results pave a new way to detect the formation of topotactic H in infinite nickelate superconductors: by measuring the phonon modes. The existence of localized phonon modes with little dispersion at ∼25 THz and ∼45 THz indicates the presence of topotactic hydrogen, which otherwise would be extremely hard to detect. These frequencies correspond to energies of 103 meV and 186 meV, respectively, beyond the range <80 meV measured for La1−*x*Sr*x*NiO2 in [71].

Lastly, we further study the case representing an incompleted reduction process: LaNiO2.125, achieved by intercalating a single O into a 2 × 2 × 2 LaNiO2 supercell (LaNiO2.125: Figure 3d). As the same consequence of employing a supercell in phonon computation, the number of phonon bands is multiplied by a factor of 8 in the frequency region between 0 THz to 20 THz. One obvious difference between undoped LaNiO2 (Figure 3a) and LaNiO2.125 (Figure 3d) is that the additional O leads to an unstable phonon mode near *q* = *X*(*π*,0,0) (blue region in Figure 3d). This phonon mode is formed by an effective vibration of the additional O along the *xy* plane in the (001) or (110) direction (and symmetrically related directions depending on the exact *q*-vector) of locally cubic coordinate. Such a mode is related to the structural transition from cubic *P*m-3m to a *R*-3c rhombohedral phase as in bulk LaNiO3, with the Ni-O-Ni bond along the *z*-direction deviating from 180◦. Our simulations for other concentrations of additional O atoms (not shown) also indicate that incomplete oxygen reduction reactions generally result in local instabilities of LaNiO2<sup>+</sup>*<sup>δ</sup>* with *δ* > 0.

#### **5. Charge Distribution**

In this section, we perform electron density calculations for LaNiO2, LaNiO2H, LaNiO2H0.125, and LaNiO2.125 compounds to investigate the bond types resulting from intercalated H and O atoms. Figure 5a,b show the electron density of LaNiO2 at the NiOplane and La-plane (light green planes of the top panels). In Figure 5a, a strong Ni-O bond is observed, while the low electron density between each Ni-O layer reveals a very weak inter-layer coupling, indicating the strong quasi-2D nature of the infinite layer nickelates. In Figure 5b, no bonds are formed between the La (Nd) atoms. The *A*-site rare-earth elements merely play the role of electron donors.

**Figure 5.** DFT calculated valence charge density of (**a**,**b**) LaNiO2, (**c**,**d**) LaNiO2H, and a LaNiO2 supercell doped with a single (**e**,**f**) H and (**g**,**h**) O atom. For each compound, the charge density of (020) and (001) planes are shown in panels (**a**,**c**,**e**,**g**) and (**b**,**d**,**f**,**h**), respectively. The La, Ni, O, and H atoms are labeled by blue, green, red, and black circles, respectively.

Figure 5c,d present the electron density of LaNiO2H along the same planes. In the NiOH-plane of Figure 5c, the comparison to Figure 5c shows that intercalated H boosts a 3D picture with an additional *δ*-type bond formed by Ni-*dz*<sup>2</sup> and H-1*s* orbitals (black circle). Along the LaH-plane (Figure 5d), *δ*-type bonds are formed by the orbital overlap between La-*dxy* and H-1*s* orbitals. For LaNiO2 with partial topotactic H (LaNiO2H0.125 in Figure 5e,f), the additional H atoms play similar roles at the Ni-H and the La-H bonds as in LaNiO2H. The Ni and La atoms without H in between are similar to those in Figure 5a,b, and those with H are akin to Figure 5c,d. This indicates that the effects induced by topotactic H are indeed very local, i.e., they only affect the the nearest Ni and La atoms.

In Figure 5g,h, for LaNiO2.125, the additional O increases the local *c*-lattice (Ni-Ni bond length via the additional O) from the LaNiO2 value of 3.338 Å to 4.018 Å which is even larger than the DFT-relaxed value of LaNiO3: 3.80 Å. This lattice expansion can be clearly seen in Figure 5g. The large electron density between Ni and O along the *z*-direction indicates the strength of this Ni-O bond in the *z*-direction is comparable with the ones along *x*/*y* directions. From Figure 5h, we conclude that similar La-O bonds are formed after intercalating additional O atoms, and the La-La distance is shrunken by the additional O atom from 3.889 Å (LaNiO2) to 3.746 Å between the La atoms pointing to the additional O. However, from the electron density plot, the La-O bond strength seems not stronger than the La-H bonding in Figure 5c,e. This can be explained by the fact that both O-*px* and -*py* orbitals do not point to orbital lobes of La-*dxy*, leading to a comparable bond strength as the La-H bond in LaNiO2H*δ*.

#### **6. Conclusions and Outlook**

Our theoretical study demonstrates that the parent compounds of infinite-layer nickelate superconductors, LaNiO2 and NdNiO2, are energetically unstable with respect to topotactic H in the reductive process from perovskite La(Nd)NiO3 to La(Nd)NiO2. The presence of H, which reshapes the systems from *AB*O2 to the hydride-oxide *AB*O2H, triggers a transition from a quasi-2D strongly correlated single-band (*dx*<sup>2</sup>−*y*<sup>2</sup> ) metal, to a two-band (*dx*<sup>2</sup>−*y*2+*dz*<sup>2</sup> ) anti-ferromagnetic 3D Mott insulator. Our predictions [21] have been reproduced by other groups using DFT+*U* calculations for other similar *AB*O2 systems [53,54]. The recent experimental observation [70] of Ni2<sup>+</sup> (3*d*8) in nickelates indicates the existence of topotactic H, as do NMR experiments [55]. The presence of H and its consequence of a 3D Mott-insulator is unfavorable for the emergence of superconductivity in nickelates. However, it is difficult to detect topotactic H in experiment. Three factors contribute to this difficulty: (1) the small radius of H makes it hard to be detected by commonly employed experimental techniques such as X-ray diffraction and scanning transmission electron microscopy (STEM). (2) As revealed by our phonon calculations, the dynamical stability of La(Nd)NiO2 does not rely on the concentration of intercalated H atoms. Hence, the same infinite-layer structures should be detected by STEM even in the presence of H. (3) As revealed by electron density distributions, the topotactic H does not break the local crystal structure either (e.g., bond length and angle); the H atoms merely affect the most nearby Ni atoms via a Ni-*dz*<sup>2</sup> -H-1*s δ*-bond. This is different if we have additional O atoms instead of H: O atoms do not only induce a dynamical instability but also obviously change the local crystal by enlarging the Ni-Ni bond length and angle visibly. Oxygen impurities also lead to unstable phonon modes in LaNiO2+*<sup>δ</sup>* and thus a major lattice reconstruction.

The ways to avoid topotactic H revealed by our calculations are: in-plane compressive strains and bivalent cation doping with Sr or Ca. This draws our attention to the recently synthesized (Nd,Sr)NiO2 films [51], which were grown on a (LaAlO3)0.3(Sr2TaAlO6)0.7 (LSAT) instead of a SrTiO3 (STO) substrate, inducing an additional 0.9% compressive strain. These new films were shown to be defect-free and with a considerably larger superconducting dome from 10% to 30% Sr-doping and a higher maximal *Tc* ∼20 K [51], compared to 12.5%–25% Sr-doping and *Tc* ∼15 K for nickelate films grown on STO which show many stacking faults [5,6,50]. The compressive strain induced by replacing the STO

substrate (*a* = 3.905 Å) by LSAT (*a* = 3.868 Å) may turn the positive *Eb* to negative, thus contributing to suppressing defects and recovering a single *dx*<sup>2</sup>−*y*<sup>2</sup> -band picture.

Besides avoiding topotactic H, compressive strain is also predicted as an effective way to enhance *Tc*. Previous dynamical vertex approximation calculations [27,57] reveal that the key to enhance *Tc* in nickelates is to enhance the bandwidth *W* and reduce the ratio of Coulomb interaction *U* to *W*. Based on this prediction, we have proposed [27,57] three experimental ways to enhance *Tc* in nickelates: (1) In-plane compressive strain, which can indeed be achieved by using other substrates having a smaller lattice than STO, such as LSAT (3.868 Å), LaAlO3 (3.80 Å), or SrLaAlO4 (3.75 Å). The smaller in-plane lattice shrinks the distance between Ni atoms and thus increases their orbital overlap, leading to a larger *W* and a smaller *U*/*W*. Recent experimental reports have confirmed the validity of this approach by growing (Nd,Sr)NiO2 on LSAT [51] and Pr0.8Sr0.2NiO2 on LSAT [72]. (2) Applying external pressure on the films plays the same role as in-plane strain for the, essentially 2D, nickelates. This has been experimentally realized in [73]: under 12.1 GPa pressure, *Tc* can be enhanced monotonously to 31 K without yet showing a saturation. (3) Replacing 3*d* Ni by 4*d* Pd. In infinite-layer palladates such as NdPdO2 or LaPdO2 and similar compounds with 2D PdO2 layers and separating layers between them, the more extended 4*d* orbitals of Pd are expected to reduce *U*/*W* from *U*/*W* ∼7 for nickelates to *U*/*W* ∼6 for palladates. Further experimental and theoretical research on the electronic and magnetic structure and the superconductive properties of palladates are thus worth performing.

**Author Contributions:** L.S. and K.H. conceptionalized the study; L.S. performed the DFT calculations; L.S., K.H. and P.W. contributed to the writing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been supported by the Austrian Science Funds (FWF) through project P 32044.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data will be made available upon reasonable request.

**Acknowledgments:** We thank M. Kitatani, J. Tomczak, O. Janson and Z. Zhong for valuable discussions and the Austrian Science Funds (FWF) for funding through project P 32044. L.S. also thanks the starting funds from Northwest University. Calculations have been done on the Vienna Scientific Clusters (VSC).

**Conflicts of Interest:** The authors declare no competing interests.

#### **References**


#### *Review* **Superconductivity and Charge Ordering in BEDT-TTF Based Organic Conductors with** *β--***-Type Molecular Arrangement**

**Yoshihiko Ihara 1,\*,† and Shusaku Imajo 2,†**


**\*** Correspondence: yihara@phys.sci.hokudai.ac.jp

† These authors contributed equally to this work.

**Abstract:** Exotic superconductivity that appears near the charge ordering instability has attracted significant interest since the beginning of superconducting study. The discovery of possible coexistence of charge ordering and superconductivity in cuprates and kagome metals has further fascinated researchers in recent years. In this review, we focus on the BEDT-TTF-based organic superconductor with *β*---type molecular packing sequence, which shows the charge ordering transition in the very vicinity of superconducting transition, and summarize the experimental results reported up to the present. At the charge ordering temperature, ultrasonic measurement detects the softening of the crystal lattice, and 13C-NMR measurement shows an increase in nuclear spin-lattice relaxation rate divided by temperature 1/*T*1*T*. These results suggest that low-energy dynamics are activated near the charge ordering transition, leading us to invoke the charge-fluctuation mediated superconducting pairing mechanism.

**Keywords:** superconductivity; charge order; molecular conductor

**1. Introduction**

Metallic conductivity in a material is introduced by doping carriers to the conduction band. When the conduction band is empty, a material simply shows an insulating behavior and is referred to as the band insulator. Even with the carriers in the conduction band when the doping level is 1/2, the onsite electron-electron correlations disturb the itinerancy of carriers to cause another insulating state, known as the Mott–Hubbard insulator [1,2]. Prominent many-body effects near the half filling result in fascinating physics, such as metalinsulator transition and unconventional superconductivity. Moreover, at the doping level of 1/4, long-range Coulomb interactions enforce the carriers to localize at every second site, leading again to an insulating state with charge ordering [3,4]. The electrons near this charge ordered state would also host exotic ground states through the long-range electron-electron correlation effect. The effect of charge ordering on superconductivity has been intensively studied since the discovery of charge ordering in copper-oxide superconductors [5,6]. Recent discovery of superconductivity in the charge ordered state of kagome metal CsV3Sb5 further promotes the discussion on the novel superconducting mechanism mediated by charge fluctuations [7–9]. The increased charge fluctuations near the valence criticality in CeCu2Si2 were also proposed to support superconductivity [10]. Above all, organic conductors that show both superconductivity and charge ordering have been the most intensively studied for decades [11–19]. Although coexistence between superconductivity and charge instability has been found in a variety of materials with different electronic orbitals, no universal understanding on its mechanism has been given. In this review article, we focus on a BEDT-TTF-based organic conductor, which shows the superconducting and charge ordering transitions at almost the same temperature, and summarize the experimental studies reported thus far.

**Citation:** Ihara, Y.; Imajo, S. Superconductivity and Charge Ordering in BEDT-TTF Based Organic Conductors with *β*---Type Molecular Arrangement. *Crystals* **2022**, *12*, 711. https://doi.org/ 10.3390/cryst12050711

Academic Editor: Andrej Pustogow

Received: 19 April 2022 Accepted: 11 May 2022 Published: 17 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The BEDT-TTF based organic conductors with the chemical formula (BEDT-TTF)2*X* (*X* = monovalent anion) ideally possess 1/4 carriers per single BEDT-TTF site. In the celebrated *κ*-type salts, however, the conduction band is constructed from effective (BEDT-TTF)2 molecular orbital formed by the strongly dimerized two BEDT-TTF molecules, and thus the half filling per (BEDT-TTF)2 dimer results in Mott physics [20]. The 1/4 filling per BEDT-TTF site is restored by arranging the BEDT-TTF molecules in a way that prevents dimerization, for instance, *θ*-type and *β*-type molecular arrangements [18,19,21,22]. The charge ordered state is found in these series of organic conductors, and, most importantly, superconductivity appears near the charge criticality. In this respect, our interest is directed to the *β*---type molecular arrangement, because the charge ordering temperature *T*CO 8.5 K was found in the proximity to the superconducting transition temperature *T*<sup>c</sup> 7 K in *β*---(BEDT-TTF)4[(H3O)Ga(C2O4)3]·PhNO2 [23]. Here, we summarize the electronic properties of *β*---(BEDT-TTF)4[(H3O)*M*(C2O4)3]·*G* salts with several metallic ions *M* and guest molecules *G* and review the detailed experiments performed for a salt with the specific combination of *M* = Ga and *G* = PhNO2.

This review article is constructed in the following way. In Section 2, we summarize the crystal structure and superconducting transition temperatures reported for a series of *β*---(BEDT-TTF)4[(H3O)*M*(C2O4)3]·*G* salts. (Hereafter, these chemical formulae are abbreviated as *β*---*M*/*G*.) Section 3 addresses the chemical pressure effects introduced by the guest molecules. The quantum oscillation, elastic constants, 13C-NMR, and EPR measurements conducted to investigate the charge ordered state in *β*---Ga/PhNO2 are summarized in Sections 4–7, respectively. In Section 8, we present the superconducting properties measured by the temperature, field, and field-orientation dependence of heat capacity. Finally, in Section 9 a summary of our current understanding on the basis of reported results and perspective will be given.

#### **2. Crystal Structure and Fermi Surface**

*β*---(BEDT-TTF)4[(H3O)Ga(C2O4)3]PhNO2, which we focus on in this review, belongs to the 4:1 *β*---type BEDT-TTF-tris(oxalato)metallate salts known as the Day series [24,25]. This 4:1 *β*---type family, having the formula (BEDT-TTF)4*AM*(C2O4)3*G*, is synthesized by means of electrochemical oxidation of the organic BEDT-TTF donor with 18-crown-6 ether and counter molecule *AM*(C2O4)3 in solvents *G* under galvanostatic current. Martin and co-workers successfully obtained pseudo-*κ* [26–28], *α*-*β*-- [26,29], and 2:1-*β*- with 18-crown-6 [26,30–32] phases in addition to the 4:1 *β*---phase by optimizing the synthesis conditions for each phase. For the 4:1 *β*---phase, independent studies by Akutsu, Coronado, Martin, and Prokhorova and their co-workers report an extensive variety of *M*(=Fe3<sup>+</sup>, Cr3<sup>+</sup>, Ga3<sup>+</sup>, Rh3<sup>+</sup>, Ir3<sup>+</sup>, Ru3<sup>+</sup>, Al3<sup>+</sup>, Co3<sup>+</sup>) and *G*(=C5H5N, DMF, PhCN, PhNO2, PhBr, PhCl, PhF, PhI, CH2Cl2) [23,25–28,33–38]. When *M* = Ga and *G* = PhNO2, the 4:1 *β*---salt is often obtained as the minor product with the major product of the semiconducting pseudo-*κ* phase crystals. *β*---Ga/PhNO2 salt crystallizes as black distorted hexagonal rods, whereas the pseudo-*κ* phase crystallizes as dark-brown diamond-shaped plates.

The 4:1 *β*---salts crystallize in the monoclinic space group *C*2/*c*. As shown in Figure 1, this series has a two-dimensional layered structure composed of the BEDT-TTF layers and counter layers. The two crystallographically independent BEDT-TTF molecules, A and B, are stacked in the *β*---type arrangement with weak tetramerization without strong dimerization in the *a*-*b* plane. In the counter anion layers, *M*3<sup>+</sup> and *A*<sup>+</sup> form a hexagonal frame bridged by (C2O4) <sup>2</sup>−, and *G* occupies the cavity. The weak tetramerization of the quarter-filled BEDT-TTF system leads to the semimetallic band structure with the compensated Fermi surface of the electron and hole pockets, as calculated by the extended-Hu¨ ckel tight-binding method [39] with a unit cell transformation from the C-centered lattice to the primitive lattice [40,41], shown in Figure 1d,e. The small Fermi pockets indicate a small number of itinerant carriers, which weakens the Coulomb screening effect, and therefore, the long-range Coulomb repulsion effectively induces the instability of charge distribution.

**Figure 1.** (**a**) Interlayer packing structure of *β*---(BEDT-TTF)4*AM*(C2O4)3*G*. The red rhombi indicate tetramers of the BEDT-TTF molecules, and the blue dashed lines represent the counter anion layers. (**b**) Arrangement of the BEDT-TTF molecules in the *β*---type packing motif. In this 4:1 *β*---phase, there are the two crystallographically independent BEDT-TTF molecules, A and B. The arrows signify transfer integrals between them along the diagonal (red), stacking (green), and horizontal (blue) directions. (**c**) Honeycomb cavity composed of the monovalent cation *A*<sup>+</sup> and tris(oxalato)metallate *M*(C2O4)3 2− in the counter anion layer. The guest molecule *G* occupies this hexagonal vacancy. (**d**,**e**) Fermi surface (**d**) and band structure (**e**) of *β*---(BEDT-TTF)4[(H3O)Ga(C2O4)3]PhNO2 (Ga/PhNO2) derived from the band calculation using the extended-Huckel tight-binding method [ ¨ 39–41].

#### **3. Chemical Pressure Effect**

Because electronic states depend on the transfer integral between the BEDT-TTF molecules, each *M*/*G* should exhibit different low-temperature behavior through the chemical pressure effect. Figure 2a shows the temperature dependence of the out-of-plane resistance *R*(*T*) reduced by *R*(300 K) for several *M*/*G* [41]. In the inset, the low-temperature region is enlarged to show the superconducting transition clearly. Whereas *R*(*T*)/*R*(300 K) above 100 K is moderate in all salts, the low-temperature *R*(*T*) behavior strongly depends on the chemical pressures. For the higher-*T*<sup>c</sup> salts, such as Fe/PhCN and Ga/PhNO2, an abrupt increase in the resistance is observed at approximately 10 K, whereas the lower-*T*c salts and non-superconducting salts do not show such semiconducting behavior. A simple deduction is that the factor producing this semiconducting nature raises *T*c. Namely, the origin of the change in the electronic state around 10 K should be connected to the origin of superconductivity. To organize these results, the chemical substitution effect is summarized in view of the effective pressure applied to the electronic state. Assuming that the size of the hexagonal cavity in the counter layer governs the arrangement of the BEDT-TTF molecules and the distances between them, the size of *G* occupying the vacancy must be important. Because the guest molecules *G* orient toward the *b* axis as shown in Figure 1, the *b*-axis length should dominantly control the chemical pressure. Indeed, the longer *G*, such as PhCN and PhNO2, gives the longer *b*-axis length compared to the shorter *G*, such as C5H5N and CH2Cl2 (Figure 2c). To ascertain the *b*-axis length dependence of the electronic state, we plotted *T*c as a function of *b*-axis length in Figure 2c. The values of *T*c for the series of salts reported in Refs. [23,28,34–36,38,41–43] were determined by the resistivity data presented in these studies. The large positive correlation coefficient *R*∼0.84 indicates that the electronic state is certainly under the influence of the *b*-axis length.

To understand the low-temperature electronic state in more detail, the magnetoresistance measured at 1.5 K [41] is displayed in Figure 3a. For some salts, superconductivity is observed and immediately suppressed below a few teslas. At higher fields (>10 T), Shubnikov–de Haas (SdH) oscillations are observed. Note that Fe/PhCN and Ga/PhNO2 do not exhibit large SdH signals as compared to other salts; however, small oscillations

certainly appear in the higher-field region as enlarged in the inset. The Fourier transform spectra of the SdH oscillations are shown in Figure 3b. As expected for a compensated metal, Fe/PhCN and Ga/PhNO2 show only one signal at the SdH frequencies of 195 T and 220 T, respectively. These results are consistent with the band calculation [40,41], suggesting that the cross-sectional area of the Fermi pockets *A*FS is approximately 10% of that of the Brillouin zone *A*BZ. In the case of the other salts, some other additional peaks are observed in the low-frequency region (∼50 T), which means that a split of the Fermi surface occurs and changes the topology of the Fermi surface. Indeed, the band calculations for some salts [41] indicate that an additional tiny hole pocket appears around the *Y* point depending on transfer integrals. To study the chemical pressure effect on the SdH oscillations simply, we here focus on the signal of the electron pocket at the *M* point (marked with each circle). Referring to the reported data by Bangura, Coldea, Prokhorova, Uji, and ourselves [37,40–44], we show the ratio *A*FS/*A*BZ (c) and the effective electron mass *m*∗/*m*<sup>e</sup> (d) as a function of *b* in Figure 3. These results also clarify the *b*-axis length dependence of the electronic state. As the elongation of the *b*-axis length with the longer *G* should reduce the amplitude of the transfer integrals in the molecular stacking, it is reasonable that increasing *b* leads to the diminution in *A*FS/*A*BZ and augmentation of *m*∗/*m*e. This perspective is consistent with the variation of band structure calculated as a function of long-range electron correlations *V* normalized by the band width *W* [41]. That is to say, the electron correlations leading to the charge disproportion facilitate superconductivity.

**Figure 2.** (**a**) Temperature dependence of the reduced electrical resistance *R*(*T*)/*R*(300 K) for some salts [41]. The inset is an expanded view below 20 K. The arrows indicate the superconducting transition temperatures *T*c determined by the onset of the drop of the resistance. (**b**) Schematic comparison of the *b*-axis length of some salts depending on the length of the guest molecule *G*. (**c**) *T*c vs. *b* of the reported *M*/*G* salts. The *T*c of each salt was determined from the resistivity data reported in Refs. [23,28,34–36,38,41–43]. This plot gives a large positive correlation coefficient *R*∼0.84. The bars indicate the onset and zero-resistivity temperatures.

**Figure 3.** (**a**) Magnetoresistance of some salts at 1.5 K [41]. Each dataset has an offset for clarity. The inset is an enlarged plot for Fe/PhCN to highlight the quantum oscillations. (**b**) Fourier spectra of the SdH oscillations shown in (**a**). The peaks marked with the circles originate from the electron Fermi pocket around the M point shown in the inset. (**c**,**d**) Ratio of the size of the electron Fermi pocket *A*FS and the Brillouin zone *A*BZ (**c**) and the effective mass (**d**) as a function of the *b*-axis length. These results were reported in Refs. [37,40–44].

#### **4. Charge Ordered State Viewed from High-Field Quantum Oscillations**

In the last section, we elucidated that the long-range electron correlations are important for superconductivity. Considering that electron correlations of *π* electrons in the quarter-filled *β*---type salts induce charge instability in itinerant carriers according to *V*/*W*, the charge degrees of freedom must be discussed. As introduced above, Ga/PhNO2 is a good target for this discussion due to the coexistence of the charge ordering and superconductivity, and thus, we hereafter focus on Ga/PhNO2.

Figure 4a shows the SdH oscillatory component in the electrical transport Δ*R*osc/*R* up to 60 T at various temperatures [45]. At low temperatures, the amplitude of the observed SdH oscillations dwindles above 40 T and cannot be explained by a single-component SdH oscillation. This behavior can be understood by assuming that the two oscillations interfere with each other and render the oscillation beating. Indeed, using the two-component Lifshitz–Kosevich (LK) formula, the field dependence can be reproduced, as described by the dotted curves. In Figure 4b, we present the thermal variation in the frequency *F* of one of the SdH signals. At higher temperatures, the *F* is approximately 220 T, which agrees with the value discussed above [40,41,43]. Below ∼8 K, where the charge disproportion develops, the value of *F* decreases. This means that unbalance between the hole and electron Fermi pockets manifests in the charge ordered state, and the different cross sections of Fermi surface yield the beating of the SdH oscillations. Considering the proportional relation between *F* and *A*FS, the result indicates that the charge ordering reduces the number of itinerant carriers. Naturally, the charge ordering affects the effective electron mass. Figure 4c displays the amplitude *A* at 42 T as a function of temperature. The dotted curve represents a fit to the LK formula with *m*∗/*m*<sup>e</sup> = 1.6. Above 8 K, the obtained data have larger values than the fitting curve. To make this deviation clearer, the ratio of the data and fit, *A*/*A*(*m*∗/*m*<sup>e</sup> = 1.6), is shown in the inset. As a smaller *m*∗/*m*<sup>e</sup> causes a larger *A*, the behavior indicates that the charge ordering enhances *m*∗/*m*e. As the charge ordering should connect with the localization of the electrons in the strong electron correlations, these results seem to lend support to the present scenario.

**Figure 4.** (**a**) Oscillatory component of the high-field magnetoresistance Δ*R*osc/*R* [45]. The dotted curves are fits to the two-component LK formula. The dashed envelope is a visual guide for the amplitude of the oscillations. (**b**,**c**) Temperature variation in the frequency *F* (**b**) and the amplitude *A* (**c**) of one of the SdH signals. The dotted curve in (**c**) is a calculation obtained by the LK formula with *m*∗/*m*e = 1.6. The inset displays the ratio of *A* of the data and the simulation when *m*∗/*m*e = 1.6.

#### **5. Elastic Response to Charge Instability**

As charge ordering transition can be detected by elastic properties through the coupling between an electric quadrupole moment and strain of lattice [46], the ultrasonic properties are measured by our group and the results reported in Ref. [45] are presented in Figure 5. The relative change in the longitudinal ultrasonic attenuation Δ*α* shows the broad

maximum around 9 K. The relative change in the elastic constant Δ*C*L/*C*<sup>L</sup> also suggests that the lattice is softened in this temperature region. The absence of the magnetic field dependence in this temperature range indicates that this anomaly is not related to superconductivity and magnetic degrees of freedom. Therefore, the lattice softening is induced by the development of the pure charge fluctuations that result in the static charge ordering at 8.5 K. This viewpoint shows a good agreement with the results of the 13C-NMR [47,48] as well as the high-field SdH [45] studies discussed in Sections 4 and 6. Note that the small softening observed at 6 K (inset) is attributable to the superconducting transition due to the suppression in a magnetic field.

**Figure 5.** (**a**,**b**) Low-temperature elastic properties, (**a**) the relative change in the ultrasonic attenuation Δ*α*, and (**b**) the elastic constant Δ*C*L/*C*L, as a function of temperature [45]. The red and blue curves are the data at 0 and 5 T, respectively. The arrow in (**a**) indicates the transition temperature of the charge ordered state *T*CO. The inset in (**b**) is the enlarged plot around *T*c. To emphasize the superconducting transition, a linear extrapolation estimated from the behavior of the normal state is shown as a dashed line in the inset.

#### **6. Charge Ordering Observed by 13C-NMR Spectroscopy**

NMR spectroscopy is a powerful tool to measure the disproportionate local charges at BEDT-TTF molecules. In this section, we explore the static properties of the charge-ordered state by reviewing the highly resolved NMR spectra taken at 15 T and precisely quantify the site-charge modulation that appears near the superconducting transition [47,48].

When the external field is applied along the *b* axis, a single 13C-NMR peak was observed at 20 K, as shown in the inset of Figure 6a. The single-peak spectrum was broadened at low temperatures, and a clearly resolved two-peak structure was observed at 1.6 K. We determined the NMR shift *δ* from the peak positions and display the temperature dependence of the NMR shift for each peak in Figure 6a. The peak splitting was observed below 8.5 K both at 15 T and 8 T. This peak splitting cannot be explained by the effects of a superconducting transition, because *T*c is suppressed to 3 K at 15 T. Only a barely resolved kink in the NMR shift was observed at *T*c. It is clear that the temperature variation of *δ* and the peak splitting in units of parts per million (ppm) are identical between 8 and 15 T. The field-independent peak separation confirms that electron spins are in the paramagnetic state and have polarization proportional to the external magnetic field. This indicates the absence of a spontaneous internal field, which could have been generated by some field-independent magnetic ordering. We identified that the transition at 8.5 K is the order in the charge degrees of freedom, in which the molecular site charges deviate from the formal value of 0.5*e*. We assigned the broad and the sharp peak as the NMR signal from the charge-rich (R) and charge poor (P) sites, respectively, according to the discussion given later. The peak separation is related to the charge imbalance between the R and P sites and is proportional to the order parameter of the charge-ordered state.

**Figure 6.** (**a**) Temperature dependence of the NMR shift at 15 T and8T[48]. The vertical axis is the peak shift from tetramethylsilane (TMS). The field-independent peak separation at the lowest temperature indicates the absence of any fixed internal field that could be associated with a magnetic transition. Inset shows 13C-NMR spectra at various temperatures. The peak splitting starts below 8 K, and a clearly split two-peak structure was observed at 1.6 K. The ratio of the integrated area between the sharp (cyan) and broad (orange) peaks was evaluated as 0.4:1. (**b**) Temperature dependences of 1/*T*1*T* measured at 3.6 and7T[47]. The peaks at *T*CO were found in both fields, suggesting unconventional electronic state affected by long-range electron-electron correlations.

The dynamical properties of electrons associated with the charge ordering transition are investigated by measuring the nuclear spin-lattice relaxation rate 1/*T*<sup>1</sup> at 3.6 and 7 T [47]. As shown in Figure 6b, 1/*T*1*T* increases with decreasing temperature, forming a peak at *T*CO. In general, 1/*T*1*T* is proportional to the square of the density of states in the Fermi liquid state and is temperature independent. The temperature dependence of 1/*T*1*T* is introduced by the enhanced magnetic fluctuations in the vicinity of magnetic criticality. For charge ordering *β*---Ga/PhNO2, however, weak magnetic fluctuations cannot induce a strong temperature dependence in 1/*T*1*T*. Charge fluctuations should be enhanced around *T*CO, but they are not directly coupled with 13C nuclear spins, because 13C nuclei with a nuclear spin *I* = 1/2 do not carry an electric quadrupole moment, which can interact with charge fluctuations and cause relaxation of nuclear magnetization. The coupling between charge and magnetic fluctuations is required to increase 1/*T*1*T* at *T*CO. One possible interpretation is that the fluctuations in local spin density, which are generated by charge density fluctuations, create fluctuating hyperfine fields at the 13C site. Direct observation of charge fluctuations by quadrupolar nuclear spins is desired to reveal the mechanisms of spin-charge coupling.

Below *T*CO, the Fermi liquid behavior in 1/*T*1*T* is absent until the superconducting state emerges. The non-Fermi liquid behavior is consistent with the semiconducting resistivity at the corresponding temperature range just above *T*c. Because 1/*T*1*T* decreases below *T*CO following a power-law close to *T*2, a clear anomaly associated with a superconducting transition was not observed at *T*c. The temperature dependence of quasi particle density of states is more clearly observed from the heat capacity measurement (Section 8) [45].

The NMR intensity is proportional to the number of 13C nuclei on the molecular sites. In the charge-ordered state, if the number of R sites were equal to that of P sites (PR pattern), the NMR intensity of the narrow peak would be comparable to that of the broad peak. However, as shown in the inset of Figure 6a, the intensity ratio between the peaks was 0.4:1, which indicates that there are at least twice as many R sites as P sites in the charge-ordered state (PRR pattern) [48]. The crystal structure of the *β*---Ga/PhNO2 salt consists of two crystallographically independent BEDT-TTF molecules, each of which forms a pair with another molecule connected by inversion symmetry. The simplest PR charge pattern can be attributed to the inversion symmetry breaking between a pair of BEDT-TTF molecules. The experimentally suggested PRR pattern addresses a more complicated charge pattern, possibly originating from the competition between the long and short range Coulomb interactions, which will be discussed later.

The large charge modulation was confirmed by the field angle dependence of NMR spectra at the lowest temperature [48]. The angle dependence of the NMR shift originates from the anisotropic dipolar Knight shift from *π* electrons and from the anisotropy of the chemical shift. The Knight (chemical) shift becomes maximum (minimum) when the external field is along the long axis of the *π* orbital, which is perpendicular to the BEDT-TTF molecular plane. To obtain the pure Knight shift contribution, the chemical shift contribution should be subtracted from the total NMR shift using a chemical shift tensor for BEDT-TTF molecules [49,50]. In Figure 7a, the Knight shift is plotted as a function of the angle between the field and the normal to the molecular plane. The amplitude of this angle dependence is proportional to the spin density in the *π* orbital. The contrasting behavior for the two independent sites confirms that the site charge is strongly modified from the formal value (0.5*e*). In particular, a very weak angle dependence for the P sites indicates that they possess a small site charge.

**Figure 7.** (**a**) Field direction dependence of Knight shift [48]. The field direction is defined with respect to the molecular long axis. The large amplitude for R sites indicates rich charge filling to the BEDT-TTF *π* orbital. (**b**,**c**) Schematic images of charge patterns in the conducting plane. The yellow balls on BEDT-TTF molecules indicate charge-rich (R) sites. The vertical P-R-R pattern can be aligned horizontally to make stripes (**b**), or shifted by one site to form a honeycomb structure (**c**).

Now we discuss a possible charge pattern in the ordered state. The NMR intensity ratio between the R and P sites suggests a PRR structure. In the two-dimensional conducting plane, we can assume PRR stripe and honeycomb structures, as shown in Figure 7b,c. The possibility of interlayer charge ordering was excluded because of the large energy cost for a charge-rich plane. To get insight into the energetic stability of the structure, we calculated the average Coulomb potential (per site) between neighboring BEDT-TTF sites, by placing point charges on molecular sites,

$$E(q) = \frac{e^2}{2\pi\epsilon\_0} \langle r^{-1} \rangle \left( 1 - Zq^2 \right),\tag{1}$$

where *<sup>r</sup>*−<sup>1</sup> = 1/8 ∑ *r*−<sup>1</sup> *<sup>i</sup>* is the average of the inverse intermolecular distance over eight neighboring BEDT-TTF sites, and *qe* is the charge deviation of the P site from 0.5*e* [48]. The coefficient *Z*, which depends on the charge pattern, is calculated to be 0.21 and 0.41 for the PRR stripe and honeycomb structures, respectively. The larger *Z* value for the honeycomb structure suggests that the honeycomb structure is energetically more stable, and is thus expected to be realized. In terms of the nearest-neighbor Coulomb potential, however, the conventional PR stripe structure with equal P and R sites should be the most

stable charge pattern, whereas this possibility is excluded by our NMR data. It is thus clear that other interactions should be taken into account. From theoretical studies on other BEDT-TTF salts with *θ*-type molecular packing, we learn that the PRR charge pattern can be stabilized in a wide parameter range in the static limit of the extended Hubbard model [51]. The exotic charge pattern appears in a *θ*-type structure, because the in-plane structure is closer to the triangular than to the square lattice. Considering the in-plane structure of *β*---Ga/PhNO2 salt to be a squeezed triangular lattice, the exotic charge pattern can be induced by the competition between the off-site and the on-site Coulomb interactions. As in *β*---Ga/PhNO2 salt, the charge ordering transition occurs in a metallic state, the effect of the transfer integral should certainly be taken into account [52] to explain the charge pattern and understand the superconducting pairing mechanism inside the charge ordered state.

#### **7. EPR Measurement to Detect SC/CO Coexistence**

NMR spectroscopy is one of the most powerful techniques to investigate the electronic state from a microscopic viewpoint. However, because of an insufficient spectrum resolution [inset of Figure 6a], we were not able to exclude the possibility of phase segregation. As an alternative probe, the electron paramagnetic resonance (EPR) experiment was performed by our group [53].

The EPR signal in *β*---Ga/PhNO2 salt originates from the *π* electrons in the highest occupied molecular orbital of the BEDT-TTF molecules. We succeeded in detecting the charge anomaly on the clearly resolved EPR spectrum by taking advantage of the anisotropy of *g* factors and the bilayer crystal structure of *β*---Ga/PhNO2 salt. Thus, the EPR experiment allows us to observe the phase segregation, if any, as an additional component of the EPR spectrum. The present results, which are explained by a single EPR contribution at any temperature, clearly evidence a uniform coexistence between the superconducting and charge ordering states.

We measured the temperature dependence of the EPR spectrum in the field applied parallel to the *b* axis (*b*-axis field) and to the direction 45◦ rotated from the *b* axis (45◦ field). In the 45◦ field, *T*<sup>c</sup> is suppressed below 3 K by a field of approximately 300 mT. At the lowest temperature in the 45◦ field, a two-peak EPR spectrum was observed, as shown in Figure 8b. With increasing temperature, the peak separation becomes small, and a single peak was observed at temperatures higher than 8.5 K. We found a trace of the two-peak structure at 8 K as the wiggle around the center of the spectrum at 344.4 mT. Therefore, the peak positions were determined by the two-peak Lorentzian fit for the spectra below 8 K [solid symbols in Figure 8c] and by the single-peak Lorentzian fit above 8.5 K [open symbols in Figure 8c]. The abrupt increase in the peak separation below 8 K clearly evidences a phase transition. This anomaly agrees with the charge ordering transition at *T*CO = 8.5 K determined from the ultrasonic [45] (Section 5) and 13C-NMR measurements [47,48] (Section 6) . The origin of the EPR peak splitting associated with charge ordering is discussed in Ref. [53].

As the charge ordering anomaly is successfully observed in the EPR spectrum, we compare the EPR spectra in the *b*-axis and 45◦ fields to unveil the relationship between the charge ordering and superconducting states. Typical EPR spectra for *b*-axis fields are presented in Figure 8a. In the *b*-axis field of 345 mT, *T*<sup>c</sup> does not change because of the extremely high upper critical field (*Bc*<sup>2</sup> > 30 T). The effect of the superconducting transition was observed in the EPR spectrum as a reduction of the integrated intensity below 7 K [Figure 8d]. Such a decrease in intensity was not observed in the 45◦ field, because *T*<sup>c</sup> is suppressed below 3 K. This result confirms that the electronic spins that would show superconductivity in zero field contribute to the EPR intensity when superconductivity is suppressed by the 45◦ field. Thus, if the superconducting part of the sample did not show the charge ordering transition, which is the case for the macroscopic phase segregation, an additional EPR peak originating from the electrons in a normal metallic state should be observed at the center of the two-peak spectrum. However, such an extra contribution was not observed at the lowest temperature of 3.6 K, as shown at the bottom of Figure 8b. The clear two-peak spectrum in the 45◦ field allows us to conclude that the superconducting state coexists uniformly with the charge ordered state. We note that the EPR intensity decreases gradually below *T*c in the *b*-axis field, and finite intensity remains even at 3.6 K. This behavior contrasts with the conventional behavior expected for a homogeneous superconducting state, in which EPR signal should disappear. The EPR in the superconducting state may originate from the nearly localized electrons in the charge ordered state, for instance, the pin site in the pinball liquid state [54–56].

**Figure 8.** Temperature dependence of EPR spectra in *b*-axis field (**a**) and 45◦ field (**b**) [53]. No spectrum splitting was observed in *b*-axis field, whereas clear spectrum splitting was observed in 45◦ field below *T*CO = 8.5 K. The EPR intensity decreases in *b*-axis field below *T*<sup>c</sup> 7 K, but finite intensity was observed even at 4 K. (**c**) EPR peak position for 45◦ field determined by Lorentzian fit. (**d**) Temperature dependence of EPR intensity normalized at 12 K.

#### **8. Gap Symmetry of the Charge-Fluctuation-Mediated Superconductivity**

In the previous sections, we understand the properties of the charge ordering, and that superconductivity emerges from the charge ordered state. Proximity of these transition temperatures implies that superconductivity is mediated by the charge fluctuations; nevertheless, we need to know the exact nature of the superconducting state itself for the direct determination of the pairing mechanism of superconductivity. The superconducting gap function is one of the most important pieces of information for discussing the pairing mechanism. Heat capacity measurement, sensitive to low-energy excitations, is a powerful method that can scrutinize this through temperature, field, and field-angle dependence [57–60].

Figure 9a shows the low-temperature heat capacity in various fields measured by our group [45]. The zero-field data indicate that the finite intercept, which corresponds to the residual electronic heat capacity, exists. This means that the electronic state is inhomogeneous even at 0 T. Considering the cleanness of the electronic state confirmed by the SdH oscillation, the coexistence with the normal state should be an intrinsic characteristic, which is consistent with the results of the previous NMR [47,48] and EPR studies [53] discussed in Sections 6 and 7. With increasing field, the electronic heat capacity is recovered because of the suppression of superconductivity. As this field-dependent contribution relates to superconductivity, we here plot the superconducting electronic heat capacity *C*sc, subtracting the lattice heat capacity and residual electronic component as a function of reduced temperature *T*/*T*<sup>c</sup> in Figure 9b. The low-temperature *C*sc exhibits neither *T*<sup>2</sup> nor *T*<sup>3</sup> dependence (dashed lines) observed in nodal superconductors, such as the *κ*-type organic superconductors, but is well reproduced by the solid curve for the anisotropic full-gapped model using <sup>Δ</sup>0[1 + *<sup>A</sup>*cos(4*φ*)]1/2, where <sup>Δ</sup><sup>0</sup> = 2.5*k*B*T*<sup>c</sup> [61] and *<sup>A</sup>* = 0.64 (Δmax/Δmin∼2.2). Magnetic-field dependence of *Cp* [45] agrees with the full-gapped model because of the observation of *H*-linear behavior at low fields. Thus, the superconducting state of Ga/PhNO2 should be fully gapped.

**Figure 9.** (**a**) Low-temperature heat capacity in various magnetic fields plotted as *Cp*/*T* vs. *T*<sup>2</sup> [45]. (**b**) Electronic heat capacity related to the superconductivity *C*sc at 0 T. The dashed line represents *T*-squared relation, and the solid curve is a calculation for the anisotropic *s*-wave superconductivity using the gap function Δ0[1 + *A*cos(4*φ*)]1/2. (**c**) In-plane angle-resolved heat capacity in a field of 0.5 T rotated from the *a* axis (0◦) to the *b* axis (90◦). (**d**) Temperature dependence of the fourfold term *C*4/*T*. The inset is a schematic illustration of the anisotropic *s*-wave gap.

To deepen the understanding of the gap symmetry in more detail, in Figure 9c we present the in-plane field-angle dependence of heat capacity in a field of 0.5 T [45]. The angular dependence is simply described by the equation, *Cp*/*T* = [*C*<sup>0</sup> + *C*2cos(2*φ*) + *C*4cos(4*φ*)]/*T* (solid curve). Taking account of the twofold rotational symmetry of the crystal structure [23], it is natural that the anisotropy of the Fermi velocity yields the twofold component *C*2cos(2*φ*). As the crystal structure does not have the fourfold symmetry, the fourfold term *C*4cos(4*φ*) should reflect the anisotropic part of the gap function. Hence, to shed light on the fourfold component, we present the temperature dependence of *C*4/*T* in Figure 9d. The negative values mean that the minima of the fourfold term are located in the directions of the crystal axes. Such a fourfold periodic term is often observed in *d*-wave superconductivity according to the anisotropy of the gap function with four gap nodes. The correspondence between the angle-resolved heat capacity and the anisotropy of the gap function has been theoretically and experimentally verified [57–59]. Based on the zero-energy Doppler effect, the gap-node positions can be determined by the sign of *C*4/*T* in a low-energy limit. However, in the present case, the fourfold component disappears at low temperatures. This disappearance is the feature of the anisotropic *s*-wave gap because the finite gap minima prohibit the quasiparticle excitations by the zero-energy Doppler effect [57]. This result also suggests that the present superconductivity is fully gapped, albeit anisotropically.

The full-gap superconductivity is distinct from the well-known dimer-Mott type organic superconductivity exhibiting *d*-wave symmetry with line nodes [60]. For the dimer-Mott system, the superconductivity is mediated by antiferromagnetic spin fluctuations, which originate from the strong on-site Coulomb repulsion *U*. Naturally, the on-site pairing is disadvantageous when utilizing *U* as an attractive force for the Cooper pairing, and the gap function must have nodes that render the sign of the gap function reverse. Similar to the dimer-Mott system, even for the quarter-filled system including the *β*---phase, *d*-wave symmetry is favored by spin fluctuations when *U* is sufficiently stronger than *V* [62]. However, with increasing *V*, the pairing around (*π*,*π*) in the momentum space is suppressed [62], and the on-site pairing interaction becomes more advantageous to reduce

the Coulomb repulsion between the nearest neighbor sites and the second nearest neighbor sites. Note the emergent gap symmetry strongly depends on the lattice geometry [62–64]; however, the charge fluctuations should facilitate the on-site pairing. Although there is no specific theoretical prediction for Ga/PhNO2, the anisotropic *s*-wave gap function does not contradict the scenario that superconductivity in Ga/PhNO2 arises from the charge fluctuations. To verify this picture, further detailed theoretical research is desirable in the future.

#### **9. Summary**

We have reviewed the experimental studies to explore the electronic properties of the organic superconductor with *β*---type molecular arrangement. Among various combinations between metallic ions *M* and guest molecules *G*, Ga/PhNO2 was found to show a remarkable electronic state, in which high superconducting transition sets in at a temperature very close to the charge ordering transition (Section 2). At the charge ordering transition, 13C-NMR experiments revealed that the local site charge at each BEDT-TTF molecule deviates from their average value of 0.5*e* (Section 6). Associated with this charge order, the Fermi surface is deformed, as evidenced from the magnetoresistance experiment (Section 4). This is the smoking gun evidence for the charge ordering on the metallic background, which can maintain high conductivity even with a partial carrier localization. It is noteworthy that the EPR measurement suggests a uniform electronic state, excluding a possibility of phase segregation between the charge ordered insulator and metallic parts of the sample (Section 7).

At the charge ordering temperature, softening of crystal lattice was observed from the ultrasonic experiment (Section 5). 13C-NMR measurement also detected an increase in the low-energy fluctuations in the temperature dependence of 1/*T*1*T* (Section 6). As the peak in the temperature dependence of ultrasonic attenuation and 1/*T*1*T* appears at a temperature very close to the superconducting transition temperature, one would naturally speculate that the low-energy fluctuations introduced near the charge criticality assists the superconducting pairing interaction.

The superconducting state was investigated by heat capacity measurements (Section 8). The temperature dependence and field direction dependence of *Cp* suggest an anisotropic *s*-wave superconducting state. This superconducting gap symmetry is preferable for the charge-fluctuation-induced superconductivity. However, to identify the novel superconducting mechanism in *β*---Ga/PhNO2 and to uncover the effect of charge fluctuations on superconductivity, further studies from both theory and experiments are required.

**Author Contributions:** Y.I. and S.I. have contributed equally to this article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by the Suhara Memorial Foundation and by the Japan Society for the Promotion of Science KAKENHI Grant Nos. 23740249, 19H01832, 20K14406.

**Data Availability Statement:** The authors declare that the data supporting the findings of this study are available within the paper. Further information can be provided by Y.I. or S.I.

**Acknowledgments:** We would like to acknowledge H. Akutsu, A. Akutsu-Sato, A. L. Morritt, L. Martin, Y. Nakazawa, R. Kurihara, T. Yajima, Y. Kohama, M. Tokunaga, K. Kindo, H. Seki, A. Kawamoto, M. Jeong, H. Mayaffre, C. Berthier, M. Horvatic, K. Moribe, and S. Fukuoka for collaborations. ´

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

