Control of superexchange interactions with DC electric fields

We discuss DC electric-field controls of superexchange interactions. We first present generic results about antiferromagnetic and ferromagnetic superexchange interactions valid in a broad class of Mott insulators, where we also estimate typical field strength to observe DC electric-field effects: $\sim 1~\mathrm{MV/cm}$ for inorganic Mott insulators such as transition-metal oxides and $\sim 0.1~\mathrm{MV/cm}$ for organic ones. Next, we apply these results to geometrically frustrated quantum spin systems. Our theory widely applies to (quasi-)two-dimensional and thin-film systems and one-dimensional quantum spin systems on various lattices such as square, honeycomb, triangular, and kagome ones. In this paper, we give our attention to those on the square lattice and on the chain. For the square lattice, we show that DC electric fields can control a ratio of the nearest-neighbor and next-nearest-neighbor exchange interactions. In some realistic cases, DC electric fields make the two next-nearest-neighbor interactions nonequivalent and eventually turns the square-lattice quantum spin system into a deformed triangular-lattice one. For the chain, DC electric fields can induce singlet-dimer and Haldane-dimer orders. We show that the DC electric-field-induced spin gap $\propto |\boldsymbol E|^{2/3}$ in the Heisenberg antiferromagnetic chain will reach $\sim 10~\%$ of the dominant superexchange interaction in the case of a spin-chain compound $\mathrm{KCuMoO_4(OH)}$ when the DC electric field of $\sim 1~\mathrm{MV/cm}$ is applied.


I. INTRODUCTION
Controlling quantum states of matter has been a longstanding subject of condensed-matter physics and other related fields. Historically, the condensed-matter-physics community has made much effort for the challenging task of searching for novel quantum states of matter such as, in quantum magnetism, spin liquids [1][2][3][4][5][6] and spinnematic states [7][8][9][10][11][12]. This task includes a search for experimental realizations of theoretical models that can host such quantum states. However, the synthesis of a compound that faithfully realizes the theoretical model does not mean realizing the quantum phase of our interest, which depends on the parameters of the compound. In general, it is even more challenging to obtain a model compound with a parameter set suitable for realizing the desired phase.
Microscopic controls of quantum states by external forces then become important. External forces can induce a quantum phase transition into the phase of our interest and can bring us the microscopic information of the quantum state through the response to the external forces. An external DC (i.e., static) magnetic field is a typical example. The DC magnetic field adds the Zeeman interaction to the Hamiltonian and induces quantum phase transitions such as Bose-Einstein condensation of magnetic excitations [13]. The pressure is another interesting example. The pressure can change the microscopic parameters of the Hamiltonian without introducing additional interactions [14][15][16].
Despite these recent developments, microscopic DC electric-field controls, which should be theoretically simpler than AC ones, are less considered in quantum spin systems partly because these systems are usually realized in Mott insulators where the charge degree of freedom is frozen. The DC electric field actually affects the Hamiltonian of quantum spin systems, as shown in Fig. 1, because the exchange interaction has an electronic origin. The superexchange interaction of spins comes from hoppings of electrons carrying the spin. It also motivates us to study DC electric-field effects that the DC field is free from heating effects that the AC one inevitably induces.
Previously, some of the authors investigated DC electric-field controls of "direct superexchange" interactions [41], originating from direct hoppings between magnetic ions, by starting from a fundamental electron model, the Hubbard model [42]. Reference [42] dealt with the DC electric potential as a site-dependent on-site potential. This treatment of the DC electric field applies to general site-dependent potentials and is convenient to discuss several related phenomena on an equal footing, as shown in Fig. 2. Note that strong DC electric fields are often required to induce sizable effects on magnetism, where the surface geometry [ Fig. 2 (b)] [43,44] and the single-cycle terahertz (THz) laser pulse [ Fig. 2 (c)] [45][46][47] are helpful for that purpose. DC electric field generated by these methods can be deemed spatially uniform. Note that we can also apply a local DC electric field to the material, for example, by using a needle-like device [Fig. 2 (d)] [39,40] [48] As of this writing, DC electric-field controls of superex- The antiferromagnetic superexchange interaction (3.8) between two magnetic ions (the orange balls) is plotted as a function of E/U d , where is the distance between the j = 0 (j = 1) and j = 1/2 sites, and U d is the on-site Coulomb repulsion at the d-orbital sites. We assume that the two magnetic-ion sites and the ligand (the blue ball) site are linearly aligned along the unit vector e0,1 that connects the j = 0, 1 sites. The DC electric field is applied either parallel (E ) or perpendicular (E ⊥ ) to a unit vector e0,1. (b) When the ligand site is located to form the right angle as shown in the left panel, the superexchange interaction between the magnetic ions can become ferromagnetic. The ferromagnetic interaction (3.11) is plotted in the right panel. When the DC electric field is applied parallel (perpendicular) to e0,1, we take = ( = ⊥ , respectively). For (a) and (b), we used typical values, U d = 5 eV, Up = 1 eV, and ∆ dp = 2 eV. Also, we took JH = Up = 1 eV for (b). The hoppings t0 and t1 can be arbitrary when we consider ratios JA(E)/JA(0) and JF(E)/|JF(0)| [see Eqs. (3.8) and (3.11)]. change interactions mediated by nonmagnetic ions are largely unexplored despite their importance to a broad class of Mott-insulating materials such as transitionmetal oxides. Reference [42] discussed DC electric-field controls of exchange interactions, which is focused on controlling the spatial anisotropy by the DC electric field. It was not yet discussed how to control exchange interactions by keeping the spatial anisotropy intact. This paper develops a theory of DC electric-field controls of superexchange interactions in geometrically frustrated quantum spin systems, starting from simple electron models. We are mainly focused on controlling microscopic Hamiltonian without affecting the dimensionality of the sample, in contrast to our previous work [42]. Our theory is widely applicable to (quasi-)two-dimensional thin-film and one-dimensional quantum spin systems on various lattices such as square, honeycomb, triangular, and kagome ones. We discuss basic exemplary applications of our theory to frustrated quantum spin systems on the square lattice and those on the chain. This paper is organized as follows. In Sec. II, we define two models that give a firm foundation for DC electric- FIG. 2. Three exemplary methods to impose a site-dependent potential to a sample. (a) We can employ a field-effect transistor to generate the DC electric field indicated by an orange arrow [37,38]. (b) A surface potential at a boundary of two different materials can be used as a source of the gradient of the potential V (x) due to the explicit inversion-symmetry breaking at the surface. (c) A single-cycle terahertz laser pulse can also be a DC electric-field source within a time scale sufficiently shorter than the pulse width (see Sec. VI D.). (d) A needle-like electrode device such as scanning tunneling microscopes can yield a spatially local DC electric field [39,40]. The dashed curves depict electric force lines.
field controls of superexchange interactions in specific cases. In Sec. III, we perform fourth-order perturbation expansions on those two models and derive superexchange interactions in generic forms. Also, we estimate typical values of DC electric-field strength from our results. We apply generic results of Sec. III to specific situations in Secs. IV and V. Section IV is devoted to geometrically frustrated quantum spin systems on the square lattice, where the nearest-neighbor J 1 and nextnearest-neighbor J 2 exchange interactions compete with each other (see Sec. IV A). First, we deal with a simple toy model in Sec. IV B to demonstrate controlling a ratio of J 1 /J 2 by a DC electric field without affecting the spatial anisotropy. Next, we investigate a more realistic situation corresponding to a compound BaCdVO(PO 4 ) 2 [49] (Sec. IV C). To get insight into the essential effects of DC electric fields, we employ an electron model that significantly simplifies the structures of those compounds. The simplified model that emulates the experimental reality is studied in Secs. IV D and IV E. Section V discusses another application of our generic results to frustrated quantum spin chains formed on a CuO 2 chain. We show that the DC electric field induces an alternation of nearest-neighbor exchange interactions, called bond alternation. We discuss experimentally observable consequences of the DC electric-field-induced bond alternation, which turns out to differ from the DC magneticfield-induced one. Section VI discusses other major DC electric-field effects not incorporated in our analysis. We summarize the paper in Sec. VII. 3. (a) An electron configuration of an unperturbed ground state of the model (2.2) is schematically drawn. Bars represent the unperturbed energy levels of the j = 0, 1 2 , and 1 site from left to right, whose energy differences are depicted as those among heights. Each bar corresponds to a single orbital. Namely, each bar can have two electrons with the opposite spins. Since the p orbital at the j = 1/2 site is fully occupied, the eigenenergy is raised by the repulsive potential energy Up. (b) Example of the d and p orbitals for the model (2.2), where the j = 0 and 1 sites have the d x 2 −y 2 orbital and the j = 1/2 site has the px orbital. The hopping amplitudes satisfy −t0 = t1 = t.

II. GENERIC MODELS
Our argument is based on a degenerate perturbation theory of single-band or multi-band Hubbard models with a site-dependent potential, whose Hamiltonian H is given by the following generic form: (2.1) Here, H U represents potential terms such as the on-site Coulomb repulsion and site-dependent potentials, and H t represents hopping terms of electrons. Note that we include the DC electric field in the site-dependent potential. Throughout the paper, H t is regarded as a perturbation to H U . We take two simple and generic models. Both models are made of three sites, two of which have d orbitals, and the other site in the middle has p orbitals. The number of sites is minimized but large enough to yield the superexchange interaction. To discuss the superexchange interaction between d electrons, we assume that each d orbital is half occupied, and the p orbitals are fully occupied in the subspace of the degenerate unperturbed ground states. The difference of two models lies in a degeneracy of p orbitals at the middle site (Figs. 3 and 4).
The first model (Fig. 3) is the three-site single-band Hubbard model with the following Hamiltonian:  (2.5). Since the p orbitals are half occupied, the Coulomb exchange interaction acts on the p 1 2 ,µ orbitals (µ = 0, 1) and lift their degeneracy by reconstructing them. The reconstructed p orbitals at the j = 1/2 site are p 1 2 ,± = p 1 2 ,0 ± p 1 2 ,1 . The eigenenergies of the two p orbitals differ by 2JH − Up when they are half occupied, that is, when the lower-energy level p 1 2 ,+ is fully occupied and the higher-energy one p 1 2 ,− is empty. (c) Example of the d and p orbitals. Here, we took d0 = d1 = d 3x 2 −r 2 , p 1 2 ,0 = px, and p 1 2 ,1 = py. There are two hoppings, t0 = −t between the d 3x 2 −r 2 and px orbitals and t1 = 1 2 t between the d 3x 2 −r 2 and py orbitals.

2)
where the jth site has a d orbital for j = 0, 1 and a p orbital for j = 1/2. d j,σ and p j,σ are annihilation operators of d and p electrons with the spin σ =↑, ↓.
Hopping amplitudes, t 0 , t 1 ∈ C are supposed to be complex to incorporate U(1) fluxes, if necessary. V j represents the site-dependent and spin-independent potential of electrons, into which the DC electric potential enters. Generally, p and d orbitals have different eigenenergies. These atomic orbital eigenenergies are incorporated into our model through the on-site potential term V j . The last term of Eq. (2.3) represents the Zeeman energy, where h is the DC magnetic field. Though h = gµ B B with B = µ 0 H is the Zeeman splitting rather than the DC magnetic field H, we roughly identify h and H, and call them simply the DC magnetic field in this paper. Both the external magnetic field is spatially uniform. If nonuniform, the exchange interactions would depend on the magnetic field. We assume that the on-site repulsions U d > U p > 0 are much larger than |V j | and |t j | to validate a degenerate perturbation expansion about |t j |/U d and |t j |/U p . The condition of small |V j | is required to ensure that the low-energy physics involves magnetic excitations only. The same parameter conditions apply to the other model given below.
The second model (Fig. 4) has two p orbtals at the j = 1/2 site:

5)
The j = 1/2 site has two p orbitals labeled by the index µ = 0, 1. The operator s µ = 2 s,s =↑,↓ p † 1 2 ,s,µ σ ss p 1 2 ,s ,µ denotes the spin-1/2 hosted by the p orbital µ = 0, 1, where σ = (σ x , σ y , σ z ) is a set of the Pauli matrices. The operator s µ makes sense only when the p orbitals µ = 0, 1 are half occupied. We include the interaction, −J H s 0 · s 1 , only when both the two p orbitals are half occupied and otherwise may drop it from Eq. (2.6). The coupling −J H < 0, which is the ferromagnetic direct exchange between the p orbitals and related to the Hund's rule, affects the p orbitals only when they are half occupied. Though the inter-band Coulomb potentials are omitted in this model for simplicity, it is straightforward to take them into account. For simplicity, we focus on a situation where the eigenenergy, E p , of the empty p orbital is equal to or lower than the eigenenergy, E d , of the empty d orbital in both models, that is, E p ≤ E d . These eigenenergies are encoded in the on-site potential V j . Here, we introduce a parameter ∆ j (E) for j = 0, 1, E denotes the DC electric field, and ∆ j (0) = E d − E p represents the eigenenergy difference of the d orbital at the j site from the p orbital aside from the Coulomb repulsion energy. For example, when we apply E = E in Fig. 1 (a), we obtain ∆ 0 (E ) = E d − E p − |e| E and ∆ 1 (E ) = E d − E p + |e| E , where −e < 0 is the electron charge.
Though we use parameters that meet the condition, ∆ j (0) ≥ 0 (j = 0, 1), (2.9) for simplicity throughout this paper, nothing forbids real materials from violating the inequality (2.9). We emphasize that our results also hold for ∆ j (0) < 0. The eigenenergy differences can be shifted by the amount of the Coulomb repulsive energy when the orbitals are partially or fully occupied [Figs. 3 (a) and 4 (a,b)]. We conclude this section by mentioning the orbital degeneracy. We assumed above the single d orbital at each magnetic-ion site, namely, the absence of the d-orbital degeneracy. This assumption can be easily relaxed. Let us consider situations where only one d orbital at each magnetic-ion site can participate in hoppings between the intermediate nonmagnetic-ion site because of symmetries. For example, in the setup of Fig. 4 (c), the p x orbital at the j = 0 has a nonzero hopping amplitude to the d 3x 2 −r 2 orbital at the j = 1/2 site and has vanishing hopping amplitude to the other four d orbitals, d xy , d yz , d zx , and d y 2 −z 2 . The low-energy physics in these multi-d-orbital cases are also effectively described by the models (2.2) and (2.5) unless symmetry-breaking structural distortion occurs.

III. SUPEREXCHANGE INTERACTIONS
This section describes the degenerate perturbation theory of the two models (2.2) and (2.5) and shows that the former (the latter) model gives rise to a Heisenberg exchange interaction with the antiferromagnetic (ferromagnetic, respectively) coupling.

A. Definition of the effective Hamiltonian
The unperturbed ground state is 2 N -fold degenerate for the two models (2.2) and (2.5), where N = 2 is the number of d-orbital sites. Let us define a projection operator P onto the subspace of the Hilbert space spanned by the degenerate ground state. Q = 1 − P is a projection operator onto its complementary space spanned by the unperturbed excited states.
We carry out the perturbation expansion as follows. First, we perform the Schrieffer-Wolf canonical transformation [50] on the full Hamiltonian (2.1), H → e η He −η , with an anti-Hermitian operator η. The effective Hamiltonian H eff that describes the low-energy physics is then defined as H eff = P e η He −η P. (3.1) The unitary operator e η keeps the excitation spectrum of the model unchanged but can simplify the effective Hamiltonian (3.1). η is determined so that the e η He −η is commutative with P [51]. Generic forms of the perturbation expansion of H eff up to the sixth order are listed in the appendix of Ref. [51].
In our case, a relation, P H t P = 0, simplifies the perturbation significantly. Up to the fourth-order perturbation expansion, the effective Hamiltonian (3.1) is expanded as where E g is the unperturbed ground-state energy and ( 1 The first-and third-order terms vanish trivially in our models. The zeroth-order term (3.3) is mostly constant but contains one important term, the Zeeman energy: S j = 2 s,s =↑,↓ d † j,s σ ss d j,s is the S = 1/2 spin operator at the j = 0, 1 site. Note that the spins of the p orbitals do not appear in Eq. (3.6) since they are fully occupied in the unperturbed ground-state subspace. The DC magnetic field h appears only in the zeroth-order term (3.6) because it is spatially uniform and our model Hamiltonians do not have spin-orbit couplings. The secondorder term (3.4) is a constant that has no impact on low-energy physics and thus discarded hereafter.

B. Antiferromagnetic superexchange
Performing the fourth-order perturbation expansion on the model (2.2) (see the Appendix Sec. A 1), we obtain the following effective Hamiltonian: where J A is the antiferromagnetic superexchange interaction, where the DC electric field E enters into Eq. (3.8) through the energy difference ∆ j (E). When E = 0, the Heisenberg superexchange coupling (3.8) is antiferromagnetic. In fact, is positive when ∆ j (0) ≥ 0 and U d > U p > 0. The latter inequalities are the case. The former condition (2.9) is likely to be the case but can be violated, as we mentioned below Eq. (2.8). Small DC electric field modifies the strength of the antiferromagnetic exchange coupling (3.8) with keeping its sign.

C. Ferromagnetic superexchange
When applied to the model (2.5), the fourth-order perturbation expansion (appendix A 2) leads to the following effective Hamiltonian of the model (2.5): The sign of J F (E) is determined by that of J H . The coupling J H must be positive since it represents the di-rect Coulomb exchange interaction [52]. The additional also required to guarantee that the right hand side of Eq. (3.11) is negative. The inequality is usually satis- . When the inequality (2.9) is violated, the superexchange coupling (3.11) can become antiferromagnetic. Therefore, the superexchange interaction J F (E) is ferromagnetic.

D. Estimates in typical situations
Here, we estimate the DC electric-field dependence of the antiferromagnetic (3.8) and ferromagnetic (3.11) superexchange interaction parameters for typical situations. As Fig. 1 shows, the DC electric field affects the antiferromagnetic (ferromagnetic) superexchange interaction drastically when the DC electric field is applied parallel to e (e ⊥ ), where the unit vector e (e ⊥ ) is parallel (perpendicular, respectively) to the vector that connects the j = 0 site and the j = 1 site.
Hence, to estimate the typical field strength, we give our attention to the following situations: for the antiferromagnetic case ( Fig. 3) and for the ferromagnetic case (Fig. 4). Note that the unit vector e ⊥ , perpendicular to e , is on the plane where j = 0, 1 2 , and 1 sites are put. J A (E) of Eq. (3.8) and J F (E) of Eq. (3.11) are plotted for E/U d in Fig. 1, where E is either E ⊥ or E , and is a length scale between the d-orbital and the p-orbital sites. As shown in Fig. 1 Supposing transition-metal oxides, we employ their typical values U d = 5 eV, U p = 1 eV, ∆ dp = 2 eV, J H = 1 eV, and = 5Å (e.g., Ref. [53]). Hopping amplitudes t 0 and t 1 can be arbitrary, though they should be perturbative.
J A (E) and J F (E) depend on the direction of the DC electric field. J A (E) is sensitive to E but independent of E ⊥ , which is obvious because the latter case only shifts V j uniformly for j = 0, 1 2 , 1. By contrast, J F (E) is more sensitive to E ⊥ than E . To increase J A (E) and J F (E) by 1 % of their original values at E = 0, we need |E |/U d ≈ 0.05 for the former and ⊥ |E ⊥ |/U d ≈ 0.003 for the latter case, namely, for the antiferromagnetic superexchange and for the ferromagnetic superexchange. These values are large but feasible with current experimental techniques.
As summarized in Ref. [37], the DC electric field in the conventional field-effect transistor setup [ Fig. 2 (a)] reaches the order of 1 MV/cm. Also, recently developing techniques, such as the electric double layer transistors, are capable of inducing the electric field stronger than 10 MV/cm [37,38]. Despite these technical developments, realization of the DC electric field with O(1) MV/cm is still challenging. It is thus important to look into other approaches that facilitate the achievement of the large DC electric field. There are basically two options. One is to use an alternative sample with more suitable parameters. The other is to use an alternative method to induce the DC electric field.
The simplest way in the first option is to reduce the thickness of the sample. If we can use a thin-film sample, we will be able to use larger DC electric field though the electric field is then applied only in a direction perpendicular to the thin film [ Fig. 2 (a)]. Reducing the sample thickness meets our purpose of controlling microscopic parameters without changing the composition of the compound. However, for comparison, it is now worth considering the use of a completely different sample to enhance the DC electric-field effects. Recall that the superexchange interactions (3.8) and (3.11) are the functions of E/U d . Instead of increasing |E|, we may increase /U d . Organic Mott insulators will be suitable to this purpose for their long and weak U d [54]. Let us assume, for example, U d = 1 eV, U p = 0.5 eV, ∆ dp = 0.8 eV, = 10Å (e.g., Refs. [55,56]). Then, for the antiferromagnetic superexchange and for the ferromagnetic superexchange. The last option we consider here is to employ an alternative DC electric-field source, the THz laser pulse (Fig. 2), instead of using modifying the sample [57]. The THz laser can induce a larger electric field than the DC one [58]. We can regard the THz laser pulse effectively as a DC electric field in a shorter time than the temporal width of the pulse. To adopt the THz laser pulse for the DC electric-field control of quantum magnetic states, we need to use a quantum magnet with fast enough spin dynamics (see Sec. VI D).

IV. FRUSTRATED FERROMAGNETS ON SQUARE LATTICE
Sections IV and V are devoted to applications of the generic results of Sec. III to simple frustrated quantum spin systems.

A. Model
This section deals with a frustrated quantum spin system on the square lattice. We give our attention to a frustrated spin-1/2 J 1 -J 2 model [5,8,[60][61][62][63][64][65][66][67][68][69][70][71][72][73] with the following Hamiltonian, where J 1 and J 2 represent the nearest-neighbor and next-nearest-neighbor interactions, respectively. r, r 1 ( r, r 2 ) denotes a pair of a nearest-neighbor (nextnearest-neighbor, respectively) bond connecting two sites r = (r x , r y ) and r = (r x , r y ) of the square lattice We can apply the DC magnetic field h to the model (4.1) through the zeroth-order term (3.6), but in what follows, we focus on the h = 0 case. The model (4.1) shows a rich ground-state phase diagram [8,59] depending on the ratio of J 1 and J 2 [ Fig. 5 (b)], which contains the spin-liquid phase and the spin-nematic phase. Hereafter, we employ the unit system = e = a 0 = 1 for simplicity, where −e < 0 is the electron charge and a 0 is the lattice spacing.

B. Toy model
Before addressing an experimentally feasible case, we first consider a simple toy model on a lattice of Fig. 6 (a) to get insight into the DC electric-field effects on J 1 /J 2 . The lattice contains edge-sharing pyramids, where the empty balls form a square lattice. Every crossing of the square lattice has one magnetic-ion site with a single d orbital that hosts an S = 1/2 spin. The top of the pyramid has one ligand site, which is depicted as a filled ball in Fig. 6 (a). We assume that this ligand site hosts two p orbitals. Similarly to the model (2.5), these p orbitals are fully occupied and degenerate in the unperturbed ground-state subspace. When the system is excited to a state with the half-filled p orbitals, the Coulomb exchange J H lifts the orbital degeneracy. For simplicity, we consider only two hoppings: t on the square edge and t that climbs up the pyramid to the top [ Fig. 6 (b)].
As we did in Sec. III, we minimize the number of sites down to 5 [ Fig. 6 (b)] and discuss nearest-neighbor and next-nearest-neighbor exchange interactions of spins.
The five-site model on the single pyramid has the following Hamiltonian: where d † j,σ for j = 1, 2, 3, and 4 and p † 0,σ,µ for µ = 0, 1 are creation operators of the d-orbital electron and the p-orbital one, respectively. The index µ = 0, 1 distinguishes two p orbitals. Note that we employed the periodic boundary condition in Eq. (4.2), d † j+4,σ = d † j,σ . Here, we apply the DC electric field to this system so that the electric field points perpendicular to the square lattice. Namely, we assume the following relations among on-site potentials, V j : where ∆ dp = ∆ 0 (0) = E d − E p ≥ 0 is the eigenenergy difference of the d orbital and the p orbital and > 0 denotes the height of the p-orbital site from the squarelattice plane. The electric field E is applied so that where e z is perpendicular to the basal square lattice of the pyramids (Fig. 6). The situation (4.7) is in contrast to Ref. [42] where the DC electric field is applied within the square-lattice plane. We designed the Hamiltonian (4.2) so that the superexchange interaction along the path 1 → 0 → 2 is ferromagnetic and that along the path 1 → 0 → 3 is antiferromagnetic.
with coupling constants, within the fourth-order perturbation expansion about the hoppings. Note that we implicitly assumed that the second-order direct superexchange and the fourthorder superexchange interactions are comparable with each other. In other words, we assumed a relation O(|t|) = O(|t | 2 ) so that the leading terms of the "direct superexchange" and superexchange interactions are of the same order. The first term of J 1 (E) is the direct superexchange interaction, and the second term is the ferromagnetic superexchange interaction. J 2 (E) is the antiferromagnetic superexchange interaction. We can strengthen (E > 0) or weaken (E < 0) the exchange interactions thanks to the explicit breaking of the z → −z inversion symmetry. Since the electric field (4.7) is antisymmetric in this inversion, the superexchange interactions (4.9) and (4.10) can also contain antisymmetric terms under the z → −z inversion, which is in contrast to the results of Ref. [42]. Figure 7 shows how E controls the ratio J 1 (E)/J 2 (E).

C. Motivation from experiments
We demonstrated that the DC electric field indeed controls the ratio of the nearest-neighbor and next-nearestneighbor interactions of the toy model (4.2). Since both J 1 and J 2 are antiferromagnetic in this model, the DC electric field does not drive the system into the spinnematic phase [ Fig. 5 Many experiments were reported about J 1 −J 2 squarelattice quantum magnets with J 1 < 0 < J 2 , for example, BaCdVO(PO 4 ) 2 [49,84,85] and AMoOPO 4 Cl (A = K, Rb) [83]. The former compound has (J 1 , J 2 ) ≈ (−3.6 K, 3.2 K) [49], which is in the vicinity of the phase boundary between the canted antiferromagnetic phase and the spin-nematic phase. The latter compound has (J 1 , J 2 ) ≈ (−2 K, 19 K) for A = K and (0 K, 29 K) for A = Rb [83], both of which are supposed to be deep in the canted antiferromagnetic phase.
As we show below, these compounds have characteristic crystal structures that allow for the DC electric-field control of J 1 /J 2 in a different mechanism from that for the pyramid model (4.2). In what follows, we describe the essential characteristics of the crystal structure relevant to our purpose and build a model that simplifies the crystal structure without interfering with the essence. Figure 8 (a) shows a schematic picture of the crystal structure of BaCdVO(PO 4 ) 2 viewed along the c axis. The single layer that hosts the J 1 -J 2 model is projected onto the ab plane in this figure. Large squares and small gray squares represent VO 4 pyramids and PO 4 tetrahedra, respectively [ Fig. 8   D in Fig. 8 (a).
The single layer of the J 1 -J 2 Heisenberg model (4.1) is composed of U and D pyramids connected by the PO 4 tetrahedra [49]. It is crucial in the following argument that the c coordinate of the magnetic ion depends on whether it belongs to the pyramid U or D [ Fig. 8 (b)].
To change the ratio J 1 /J 2 , we apply a uniform DC electric field, to the system, where e c is the unit vector along the c axis. If we define the origin of the electric potential at the height of the P ions, the V ion of the upward pyramid feels the electric potential −E and that of the downward pyramid feels +E with 2 > 0 being the height of the pyramid. Though the external field is uniform, the magnetic ion at r = (r x , r y ) feels a staggered electric potential (−1) rx+ry E due to the staggered structure of the pyramids [ Fig. 8 (a)].  1 and MII,2). There is a direct hopping tM−M between two magnetic ions that belong to the same kind of pyramid. To hop between two magnetic ions on the two different kinds of pyramids, the electron first hops to the ligand ion (the empty ball) at the center of the tetrahedron and next to the other magnetic ion. The hopping amplitude between the magnetic ion and the ligand is given by tM−L.

D. Simplified electron model
To investigate DC electric-field effects on BaCdVO(PO 4 ) 2 , we develop a model that emulates BaCdVO(PO 4 ) 2 with a much simpler structure. The compound BaCdVO(PO 4 ) 2 has complicated exchange paths between neighboring magnetic V ions. Starting from a V ion site, the path goes through an O ion, a P ion, and another O ion before reaching the other V ion. The miscellaneous O, Ba, and Cd atoms are omitted in Fig. 8. In our simplified model, the miscellaneous atoms of BaCdVO(PO 4 ) 2 are removed. Namely, we consider more direct hoppings between the magnetic ions along paths such as V → V or V → P → V as shown in Fig. 9 instead of V → O → P → O → V. The U and D pyramids are replaced by the magnetic ions with the same labels. Figure  Our simplified model inherits the essential structure of BaCdVO(PO 4 ) 2 . First, we keep the ligand ion at the center of a tetrahedron formed by the each of four magnetic-ion sites to yield superexchange interactions (Fig. 9). Second, our model has the alternating structure that the magnetic-ion site r = (r x , r y ) is located Δ dp = 0 Δ dp = 1 In panels (a) and (b), we take ∆ dp = 0 and ∆ dp = 1, respectively. The other parameters are taken in common as U d = 1, Up = 0.5, JH = 0.5, t = 0.01, and t = 0.12. The ratio δJ2(E)/J2(E) grows with an increase of E/U d for the two cases, where the square lattice becomes the anisotropic triangular lattice. In particular, panel (b) shows that the DC electric field changes only δJ1(E)/J2(E). above (below) the ligand site for r x + r y is even (odd) [ Fig. 8 (a)].
Following the generic argument of Sec. III, we minimize the number of sites and consider a five-site electron model with a Hamiltonian, (4.14) The model has a d orbital at four magnetic-ion sites at the vertices of the tetrahedron and two orbitals at the ligand site at the center. d † j,α,σ represents the annihilation operator of electron at the magnetic-ion site j with the spin σ. The index α = U, D indicates the upward and downward magnetic ions, respectively. Accordingly, the on-site potential is given by V U =V − E and V D =V + E with an E-independent constantV . The p L,α,σ operator annihilates the p electron at the ligand site. As the hopping term (4.14) shows, the p orbitals labeled by the index α = U, D have a hopping term from/to the d orbital with the same α index. −J H < 0 is the Coulomb exchange interaction between s α = 1 2 s,s p † L,α,s σ s,s p L,α,s for α = U, D at the ligand site.

E. Effective spin model
The low-energy effective Hamiltonian of the model (4.12) is given by Let us embed the effective model (4.15) on the tetrahedron into the square lattice of Fig. 5 (a). J 1 (E) and J 2 (E) are the nearest-neighbor and the next-nearestneighbor interactions, respectively. The perturbation expansion gives a side effect that the effective Hamiltonian (4.15) contains the third interaction of δJ 2 (E), which makes the two diagonal bonds nonequivalent (the middle panel of Fig. 10). The appearance of nonzero δJ 2 (E) is due to the explicit violation of an inversion symmetry by E. When E = 0, the bond connecting M U,1 and M U,2 is symmetrically equivalent to that connecting M D,1 and M D,2 . The DC electric field along the c axis makes these bonds nonequivalent, that is, δJ 2 (E) = 0. Figure 11 shows the DC electric-field dependence of the exchange interaction for some parameter sets.
The nearest-neighbor interaction is a superexchange interaction, and the next-nearest-neighbor interaction is mainly a direct superexchange interaction though its subleading terms are superexchange ones as shown below. Up to the fourth order of the perturbation expansion, J 1 (E), J 2 (E), and δJ 2 (E) are given by (4.17) (4.18) to make the leading terms of the exchange and superexchange interactions comparable. The side effect of nonzero δJ 2 (E) is an interesting phenomenon by itself. In the extreme situation of δJ 2 = J 2 , the spin model turns into a deformed triangular-lattice Heisenberg antiferromagnet (Fig. 10) [86][87][88]. Each triangle unit has one J 2 bond and two J 1 bonds. We can thus expect that the DC electric field changes the squarelattice antiferromagnet eventually into the triangularlattice one. The deformed triangular-lattice antiferromagnet also exhibits a rich phase diagram under magnetic fields [87]. It will be interesting to apply the DC electric and magnetic fields simultaneously to the model (4.12), searching for E-induced quantum phase transitions.
We emphasize that the argument in Secs. IV D and IV E also applies to AMoOPO 4 Cl. There are minor differences between BaCdVO(PO 4 ) 2 and AMoOPO 4 Cl from the DC electric-field viewpoint. AMoOPO 4 Cl has the Mo ion whose 4f orbital contributes to quantum magnetism, however, which can be dealt with on equal footing with the model (4.12). The orbital hosting the S = 1/2 spin does not need to have the d-orbital symmetry though we used the notation of the d orbital symbolically. Our model here as well as the generic ones of Sec. II assume that the orbital at the magnetic-ion site is nondegenerate and has a large on-site repulsion, U d > U p . The 4f orbital fulfills this assumption.
BaCdVO(PO 4 ) 2 and AMoOPO 4 Cl have the same characteristics in response to the DC electric field for the following reasons. First, AMoOPO 4 Cl also has two kinds of magnetic-ion sites with two different electric potential energies. Instead of the upward and downward pyramids, AMoOPO 4 Cl has two kinds of octahedra shifted in the opposite direction along the c axis that results in the same staggered electric-field potential (−1) rx+ry E as BaCdVO(PO 4 ) 2 . Second, the same simplification of the hopping paths to Fig. 9 applies to AMoOPO 4 Cl as well as BaCdVO(PO 4 ) 2 . We can expect that the simplification works even better in AMoOPO 4 Cl thanks to the outspread probability distribution of the 4f orbital. We thus conclude that the results (4.16), (4.17), and (4.18) will thus hold also for AMoOPO 4 Cl.

A. Experimental realizations
This section discusses another important frustrated quantum spin system of an S = 1/2 J 1 -J 2 spin chain described by the Hamiltonian The dimer is a nonmagnetic singlet type for J 1 > 0 and a magnetic triplet type for J 1 < 0 [91]. Some materials with CuO 2 chains [103-109] are known to be described by the J 1 − J 2 chain with J 1 < 0. When −2.7 < J 1 /J 2 < 0, this spin chain has the spin nematic phase in the vicinity of the fully polarized phase forced by the strong magnetic field h [ Fig. 12 (b)]. When the magnetic field h is in the vicinity of the saturation field, the phase diagram becomes rich. The phase diagram contains quasi-long-range multipolar order phases of spins, the quadrupolar, octapolar, and hexadecapolar phases for J1 < 0 [10,90]. In particular, the quadrupolar phase is known as the (quasi-long-range) spin-nematic phase [10,11,73,[92][93][94][95][96][97][98]. There are one-and two-component TL-liquid phases for J1/J2 > 4 and 0 < J1/J2 < 4, respectively [99,100].

B. Goodenough-Kanamori rule
The model (5.1) can experimentally be realized, for example, in a CuO 2 chain of Fig. 13 (a), where the Cu ion carries the S = 1/2 spin. Figure 13 (a) shows a generic case that two O ions, O I and O II , are nonequivalent. Two O ions can be equivalent but, in what follows, are supposed to be nonequivalent to derive the antiferromagnetic J 1 and ferromagnetic J 1 from the single model.
The sign J 1 is determined by the angle ∠Cu − O − Cu. The two O ions mediate the nearest-neighbor superexchange interaction. According to the Goodenough-Kanamori rule [110][111][112], the superexchange interaction is ferromagnetic when an angle ∠Cu − O − Cu is 90 • and antiferromagnetic when the angle is 180 • . We can easily understand this difference by looking into the d and p orbitals of Fig. 13 (b). When the angle is 90 • , there are no hoppings between the d x 2 −y 2 orbital of the Cu ion and the p x orbital of the O ion. Then, the superexchange interaction is well described by the model (2.5) with the two degenerate p orbitals at the oxygen site. On the other hand, when the angle is 180 • , the single p orbital allows for hoppings of electrons from the O ion to the two Cu ions on both sides. Then, the superexchange interaction is modeled by Eq. (2.2) with the nondegenerate p orbital at the oxygen site [see also Fig. 3 (b) and Fig. 4 (c)]. Generally, the angles are intermediate, when the superexchange process also has an intermediate character of the two models (2.2) and (2.5). For simplicity, we consider an ideal situation where the superexchange interaction via the O I ion is described by the model (2.5) with t 0 = t 1 = t and that via the O II ion is described by the model (2.2) with t 0 = t 1 = t [ Fig. 13 (a)]. It follows from the above modeling that We assume that the two oxygens are exactly in the middle of the two nearest-neighbor Cu ions along the chain direction, e .

C. Twofold screw symmetry
When the CuO 2 plane has two nonequivalent oxygen sites, the CuO 2 chain has a twofold screw symmetry along the chain direction, a combination of the discrete translation symmetry in the chain direction and the π spatial rotation around that direction. Let us denote the discrete translation as T 1 and the π rotation as R π . The CuO 2 chain of Fig. 13 (a) has neither T 1 nor R π symmetries though it has the T 1 R π symmetry. Nevertheless, the low-energy spin-chain model (5.1) in the absence of the electric field has both the T 1 and R π symmetries instead when the electric potentials at the two oxygen sites are exactly balanced.
The DC electric field can violate this balance, replacing the emergent high symmetry of the spin chain by the original T 1 R π one of the CuO 2 chain. Let us apply the following DC electric field to the system, along the direction perpendicular to the CuO 2 chain and on the CuO 2 plane [ Fig. 13 (a)]. We denote the electric potentials at the Cu, O I , and O II sites as 0, ± I E, and ∓ II E, respectively. The sign of the latter two potentials is positive if the O ion is below the Cu ion in the e ⊥ direction and negative if the O ion is above the Cu ion. We can assume I > II to be consistent with the assumption The effective spin-chain model under the DC electric field (5.2) has the following Hamiltonian, where the last term, called the bond alternation, is a direct consequence of the imbalance of the electric potentials of O I and O II . The nearest-neighbor interactions, J 1 (E) and δJ 1 (E), are given by We emphasize that Eqs.  (Figs. 14 and 15). This E dependence is consistent with the twofold screw symmetry of the CuO 2 chain. For instance, since R π maps E = Ee ⊥ to −E, the T 1 R π operator maps the bond alternation Hence, the Einduced bond alternation has the twofold screw symmetry though it breaks the T 1 symmetry.
Note that the right-hand side of Eq. (5.5) does not vanish for I = II , though it should from the symmetry viewpoint. This inconsistency comes from the initial assumption that we made. We assumed that the superexchange process via the O I ion and that via the O II ion are supposed to be nonequivalent from the beginning. Accordingly, Eqs. (5.4) and (5.5) hold when the difference I − II is large enough to justify the nonequivalent superexchange processes. This condition I,II is convenient for our purpose of the DC electric-field control. J 2 (E) depends on E in a complex manner. However, its DC electric-field dependence is less important than the induction of the bond alternation because the latter is much more relevant in the renormalization-group sense [113,114].
Before describing the E-induced bond alternation effects, we comment on the field direction. If we apply the DC electric field along the chain direction (e of Fig. 13), the DC electric field does not yield the bond alternation because it keeps the balance of the electric potentials at the two oxygen sites. The DC electric field E e purely changes the ratio J 1 (E)/J 2 (E), keeping δJ 1 (E) = 0, which will be relevant to studies of the quasi-long-range spin-nematic phase of J 1 − J 2 spin chains.
Δ dp = 0 Δ dp = 0.5 FIG. 14. The E-induced bond alternation δJ1(E) is plotted for the DC electric potential IE for the J1, J2 > 0 case with parameters, U d = 1, Up = 0.5, JH = 0.5, t = 0.04, t = 0.1, and II/ I = 0.60. The ratio II/ I is determined from the data for KCuMoO4(OH) [109]. We took ∆ dp = 0 for the upper panel and ∆ dp = 0.5 for the lower panel. When both J 1 and J 2 are positive, the ground-state phase diagram for E = h = 0 contains two phases: a Tomonaga-Luttinger (TL) liquid phase [113,114] and a spontaneously dimerized phase. A quantum critical point, J 2 /J 1 = α c ≈ 0.2411 [115,116], separates these two phases: the TL-liquid phase for J 2 /J 1 ≤ α c and the spontaneously dimerized phase for α c < J 2 /J 1 . The former is gapless, and the latter is gapped. Turning on the DC electric field (5.2), we can introduce the bond alternation to the J 1 -J 2 spin chain. Trivially, the bond alternation turns the spontaneously dimerized phase into an induced dimerized phase, which we call an E-induced dimerized phase, by lifting the ground-state degeneracy of the spontaneously dimerized phase [ Fig. 13 (c2)]. More interestingly, the bond alternation drives the TL liquid into the E-induced dimerized phase.
The dimer order parameter D(E) characterizes the Einduced dimer phase. When the ground state belongs to the TL-liquid phase for E = 0, the dimer order parameter is obviously zero: D(0) = 0. As soon as we apply Δ dp = 0 Δ dp = 0.5 FIG. 15. The E-induced bond alternation δJ1(E) is plotted for the DC electric potential IE for the J1 < 0 < J2 case with parameters, U d = 1, Up = 0.5, JH = 0.5, t = 0.1, t = 0.01, and II/ I = 0.81. The ratio II/ I is determined from the data for NaCuMoO4(OH) [109]. We took ∆ dp = 0 for the upper panel and ∆ dp = 0.5 for the lower panel.
the DC electric field to the TL liquid, we can open the excitation gap. The gap opening due to the bond alternation is well described by the sine-Gordon theory with the following Hamiltonian [113].
where v is the spinon velocity and g(E) ∝ δJ 1 (E) ∝ |E| is the coupling constant that leads to the spin gap. φ and θ are related to the spin operator S j as [113,114,117] with nonuniversal constants a 1 , b 0 , and b 1 . It is exactly known that the sine-Gordon theory (5.6) gives the following the dimer order parameter D(E) [118]: where the parameter d z relates the spin operator and the φ boson via (−1) j S p j S p j+1 = d z sin( √ 2πφ) + · · · for p = x, y, z in the absence of the DC electric field [119,120] and Γ(z) is the gamma function. The parameter ∆(E) represents the lowest-energy excitation gap induced by E. Note that v and d z are calculated in the model with E = 0 since the E dependence of these quantities merely leads to hardly observable corrections to the right hand side of Eq. (5.9), which we will discuss later. The fieldtheoretical result (5.9) becomes accurate when the ratio J 1 (0)/J 2 (0) is close to the critical value α c . When J 1 (0)/J 2 (0) ≈ α c , the parameters v and d z are given by v = 1.174J 1 (0) [121] and d z = 0.182 [119].
One can obtain the exact excitation gap, ∆(E), of the sine-Gordon theory (5.6) [118,122,123]: The gap (5.10) and the dimer order (5.9) are plotted in Fig. 16. We can see the power-law behaviors (5.11) and (5.12) near E = 0. The power laws (5.11) and (5.12) hold near the origin but break down at some field strength, say, for | I E/U d | ≈ 0.1 when ∆ dp = 0 (Fig. 16), because of a nonlinear E dependence of J 1 (E) and δJ 1 (E) [Eqs. (5.4) and (5.5)]. Since U d is typically ∼ 5 eV for Cu [124][125][126], we can expect this nonlinear effect for |E| > ∼ 34 MV/cm, which is extremely strong. The nonlinear effect will be hardly observable. At the same time, this estimation of |E| justifies our perturbative treatment of the DC electric field in the field-theoretical argument such as Eq. (5.9).
This E-induced dimerization differs significantly from a magnetic-field-induced dimerization that one of the authors found recently [127]. Though a fourfold screw symmetry is essential to cause the magnetic-field-induced dimerization, such a complication is not required in the E-induced dimerization. (the dashed lines). We used the same parameters as those in Fig. 14. 2. J1 < 0: Haldane dimers When J 1 < 0 < J 2 , the ground-state phase diagram at E = h = 0 contains the ferromagnetically ordered phase, a vector-chiral phase, and the Haldane-dimer phase [91]. When 0 < J 1 /J 2 1, the ground state is expected to belong to the Haldane-dimer phase, where the triplet dimer is formed on the nearest-neighbor bond. Reference [91] also discusses the effects of the bond alternation on the ground-state phase diagram of the J 1 -J 2 chain. The bond alternation weakens the geometrical frustration between the nearest-neighbor and the nextnearest-neighbor exchange interactions. While the J 1 -J 2 spin chain with J 1 < 0 < J 2 has the spontaneous Haldane-dimer phase with the twofold ground-state degeneracy, the J 1 -J 2 -δJ 1 chain has the induced Haldanedimer phase with the unique and gapped ground state. In the latter phase, the triplet dimer is formed on the bond with the stronger nearest-neighbor exchange interaction, namely, on the J 1 + δJ 1 bond for E > 0 and on the J 1 − δJ 1 bond for E < 0 [ Fig. 13 (c3)].
Differently from the J 1 > 0 case, the DC electric field does not open the spin gap for large |J 1 |. When J 1 < 0 and |J 1 | is large enough, the ground state has the spontaneous ferromagnetic order accompanied by a gapless nonrelativistic Nambu-Goldstone mode, which is robust against the small bond alternation. When |J 1 | is much smaller than J 2 , the DC electric field opens the gap exactly in the same way for both the J 1 < 0 and J 1 > 0 cases.
When max{|J 1 |, |δJ 1 |} J 2 , we can regard the J 1 -J 2 spin chain as weakly coupled spin chains. Bosonizing each spin chain, we obtain an excitation gap [128], instead of Eq. (5.10). Accordingly, the dimer order parameter follows Therefore, we find for weak E Though we obtained different power laws from Eqs. (5.11) and (5.12), it is not attributed to the sign of J 1 . In fact, even when J 1 > 0, we will find the power laws (5.15) and (5.16) if J 1 and |δJ 1 | are much smaller than J 2 . In this weak |J 1 | region, we can find the difference due to the sign of J 1 only in the nature of the dimer order whether it is the Haldane dimer or the singlet one.

E. Proposals for experiments
The following are our proposals for experiments. The CuO 2 chain of Fig. 13 is realized, for example, in ACuMoO 4 (OH) (A = Na, K) [109,129,130]. (J 1 , J 2 ) are given by (−51 K, 36 K) for A = Na and (238 K, 0) for A = K. We can expect the E-induced dimerization for A = K and the E-induced Haldane-dimer phase for A = Na For A = K, it will be intriguing to observe the 2/3power-law dependence (5.11) of the excitation gap ∆(E). We expect that the DC electric field with 1 MV/cm opens the gap ∼ 22 K in KCuMoO 4 (OH), which reaches 8.6 % of J 1 (0). Spectroscopic methods such as the electron spin resonance (ESR) [131,132], the inelastic neutron scattering [133], and the Raman scattering [134,135] will be suitable candidates for observing those characteristics.
There is another interesting phenomenon specific to KCuMoO 4 (OH). The compound KCuMoO 4 (OH) has a staggered Dzyaloshinskii-Moriya (DM) interaction, j (−1) j D · S j × S j+1 , in addition to the Hamiltonian (5.3). The DM interaction originates from the spin-orbit coupling. Though we did not take the spin-orbit coupling into account thus far in this paper, we can approximately apply our theory to this compound since the DM interaction is weak. As far as the leading effects caused by the perturbative DM interaction and the perturbative electric field are concerned, we can still adapt our theory to those systems with the weak DM interaction. If one wishes to discuss DC electric-field controls of DM interactions, one has to take into account the spin-orbit coupling in the electric models such as Eq. (2.2) and (2.5), which goes beyond the scope of this paper.
The staggered DM interaction yields a staggered magnetic-field interaction (say, −h s j (−1) j S x j ) when the external DC magnetic field h is applied in a direction perpendicular to D [136,137].
The resultant staggered magnetic field is proportional to h and |D|: h s ∝ |D|h/J. As well as the bond alternation δJ 1 (E), the staggered magnetic field h s also generates the excitation gap. The staggered magnetic field induces the Néel order in a direction perpendicular to the external DC magnetic field [137]. When the DC electric and magnetic fields are applied to the spin chain simultaneously and are tuned, a crossover will occur between the E-induced dimerized phase and the h sinduced Néel phase. The entanglement of the singlet dimer is protected by one of the time-reversal, the D 2 spin-rotation, and the bond-centered inversion symmetries [138][139][140]. All these symmetries are explicitly broken in the simultaneous presence of the DC electric and magnetic fields. Namely, nothing forbids the smooth deformation of the E-induced dimer-ordered ground state into the H-induced Néel ordered ground state. It will be interesting to look into whether the crossover actually occurs on the two-dimensional parameter space shown in Fig. 17.
Due to the staggered DM interaction, the field dependence of the excitation gap will be tractable in ESR measurements [141][142][143]. One significant difference of these two phases lies in a localized excited state at the chain end, called the boundary bound state [143], exists only in 17. Possible E-H ground-state phase diagram of KCuMoO4(OH). When E = 0 and hs = 0 (= h), the Einduced singlet-dimer phase is realized. By contrast, when E = 0 and hs = 0, the hs-induced Néel phase is realized. These ground states will be smoothly connected to each other because the symmetries protecting these gapped quantum phases are explicitly broken.
the h s -induced Néel phase [143]. With an increase of |E|, the boundary bound state will be eventually lost though it will survive for a while in the vicinity of the horizontal axis of Fig. 17.
For A = Na, the nearest-neighbor interaction is ferromagnetic. The ratio J 1 (0)/J 2 (0) ≈ −1.4 implies that the ground state at zero magnetic field belongs to the Haldane-dimer phase (Fig. 1 of Ref. [91]). Experimentally, DC electric-field effects on the Haldane-dimer phase can be observed as an increase of the excitation gap similarly to the J 1 > 0 case. We can also observe an S = 1/2 edge state. The spin-1 Haldane phase is a topological phase protected by the Z 2 × Z 2 spin-rotation symmetry [144,145]. Since the DC electric fields keep the Z 2 × Z 2 symmetry, the E-induced Haldane-dimer phase has doubly degenerate edge states on each chain end, that is, an S = 1/2 edge spin. The edge-spin degrees of freedom can be observed by, for example, the ESR spectroscopy [146,147]. An increase in the excitation gap immediately means a decrease in the correlation length. The decrease of the correlation length would affect the ESR spectrum of the spin chain with a finite chain length [147].

VI. OTHER DC ELECTRIC-FIELD EFFECTS
This section is devoted to discussions on other major DC electric-field effects that we have not dealt with in this paper. We can incorporate some effects into our model with slight modifications and some with substantial changes. The renormalization of the dielectric constant and the structural distortion fall into the former. The spin-orbit coupling falls into the latter and requires the substantial change of the model. In what follows, we briefly discuss these three effects and another important effect of the THz laser pulse.

A. Dielectric constant
We implicitly assumed that the electrons feel the external DC electric field itself. In real materials, the electron is surrounded by various charges that can screen or enhance the external electric field. Generally, the dielectric function represents how the external DC electric field is screened or enhanced. In particular, the dielectric constant in the material gets shifted from its vacuum value. We can incorporate this effect into our model by regarding E(r) in our model as the actual DC electric field that the electrons in materials actually feel.

B. Structural distortion
The strong DC electric field can possibly distort the crystal structure. Still, the E-induced change of the onsite potential V j turns out to be dominant, as shown below. The structural distortion will make the hopping amplitude t j depend on E by modifying the lattice spacing. We did not include this effect in our analysis thus far but can immediately include it without any problem. Namely, we just replace the constant t j by an E-dependent function t j (E). If one wishes to predict the precise E dependence of the hopping amplitude, one needs to model how the DC electric field distorts the lattice.
Besides, the structural distortion can lower the crystalline symmetry. The symmetry lowering leads to an important effect of switching on hoppings between d and p orbitals that were forbidden by symmetries in the absence of the DC electric field. Then, we need to include d-or porbital degeneracy explicitly into the model Hamiltonian. However, since such a newly introduced hopping amplitude is proportional to the DC electric-field strength, the inclusion of the degenerate d or p orbitals will become important only under extremely strong DC electric fields. Let us denote the additional hopping amplitude as t j (E) for j = 0, 1. Note that t j (0) = 0 by definition. It is straightforward to include t into our results. We can obtain the correction by t j (E) to Eqs. (3.8) and (3.11) by replacing t j in these relations to t j totally or partially. The fourth-order perturbation expansion shows that this correction to the superexchange coupling is an even order of |E|/U d because an electron in a p orbital hopped to a d orbital must come back to the same p orbital from the same d orbital in the Mott-insulating phase. As we saw throughout the paper, the strong DC electric field of O(1) MV/cm is already required to change the superexchange interaction by a few percent. Accord-ingly, we need much stronger DC electric field to observe the nonlinear change of the superexchange. Under such situations, the second-or fourth-order effect will be relevant only for O(10) MV/cm for inorganic materials and for O(1) MV/cm for organic ones. There will be a chance in organic materials to observe the nonlinear field effects including that by the symmetry-lowering structural distortion. However, it will be difficult to distinguish the structural distortion effect from other nonlinear terms in Eqs. (3.8) and (3.11). Therefore, the inclusion of the structural distortion keeps our result intact unless the subleading nonlinear field effects are concerned.
Finally, we comment on a possible spontaneous formation of the electric polarization due to the E-driven structural phase transition. Such a spontaneous polarization can become large even if |E|/U d 1. This large polarization can potentially lead to a large correction to the superexchange coupling. It will be interesting in the future to investigate such E-driven phase-transition effects.

C. Spin-orbit coupling
The spin-orbit coupling can dramatically change the spin Hamiltonian by adding to the spin Hamiltonian anisotropic interactions such as the DM interaction. It needs straightforward but lengthy calculations to discuss the effects of the spin-orbit coupling on electric-field controls of magnetism [148]. We will discuss in a subsequent paper [148] the combination effect of the spin-orbit coupling and the DC electric field.

D. THz laser pulse
The single-cycle THz laser pulse can be deemed the effective DC electric field when the relevant time scale of the spin system is fast enough. We need this assumption of the fast spin dynamics to guarantee that the quantum spin system quickly reaches the equilibrium before the pulse laser disappears. The typical time scale ranges from ∼ 0.1 ps to ∼ 1 ns [149][150][151][152][153][154][155][156][157]. For example, the spin dynamics with the time scale ∼ 0.1 − 1 ps is much faster than the single cycle (∼ 10 ps) of the THz laser pulse with the 0.1 THz frequency [45][46][47]. Then, we can regard the single-cycle THz laser pulse as the effective DC electric field. On the other hand, if we apply the laser pulse with the 10 THz frequency to the spin system with the time scale > ∼ 10 ps, the pulse (∼ 0.1 ps) disappears much faster than the equilibration of the spin system. Then, the single-cycle THz laser pulse can be approximated as a delta-function electric field, E(t) = E 0 δ(t). Even when the spin dynamics is much slower than the THz laser pulse width, the electron dynamics can be much faster than the pulse width. Indeed, the hopping amplitude t ∼ O(1) eV gives the time scale O(1) fs = O(10 −3 ) ps. Therefore, we can incorporate the THz laser pulse as an effective DC electric field into our model even when the pulse width is too short to equilibrate the spin system. We then notice an interesting possibility. If we derive the effective spin model, the THz laser pulse turns into the delta-function potential, δJ(E 0 )δ(t), that instantaneously modifies the exchange interaction at the time t = 0. It is an interesting direction of future studies to investigate dynamical effects caused by the instantaneous potential.
To conclude, we emphasize two points. The THz laser pulse can be deemed the effective DC electric field. The THz laser pulse can equilibrate the spin system and otherwise induces the instantaneous modification of the exchange interaction, which depends on materials.

VII. SUMMARY
This paper discussed DC electric-field controls of superexchange interactions in Mott insulators. We first presented the generic results (3.8) and (3.11) about the superexchange interaction by the fourth-order degenerate perturbation expansion of the two basic electron models. We considered the antiferromagnetic and ferromagnetic superexchange interactions and obtained their E dependence.
The other part of the main text was devoted to applications of the generic results to basic geometrically frustrated quantum spin systems. While the results of Sec. III applies to various quantum spin systems regardless of the lattice and the dimensionality, we give our attention to low-dimensional quantum spin systems in this application part to investigate how DC electric fields applied perpendicular to the system control the superexchange interactions without interfering with the spatial anisotropy. In contrast to this paper, the previous paper [42] discussed how to control the "direct superexchange" interaction [41] and introduce the spatial anistoropy to the system through the DC electric field parallel to the system. These two theories complement each other.
As the first application, we considered J 1 -J 2 frustrated quantum spin systems on the square lattice in Sec. IV. We first demonstrated in the toy model how the DC electric field adjusts the parameter J 1 (E)/J 2 (E) that determines the fate of the ground state. Next, we applied the generic results of Sec. III to the more realistic model that emulates two compounds BaCdVO(PO 4 ) 2 [49] and AMoOPO 4 Cl [83]. Combined to their crystal structures, the DC electric field breaks the inversion symmetry that introduces the nonequivalence of the next-nearestneighbor interactions, J 2 (E)±δJ 2 (E) though the original model exactly has δJ 2 (0) = 0 at E = 0. An increase of the nonequivalence by δJ 2 (E) turns the frustrated square-lattice quantum antiferromagnet eventually into the deformed triangular-lattice antiferromagnet, which will offer a unique experimental method to change the lattice structure effectively.
We also discussed applications to geometrically frustrated one-dimensional quantum spin systems, the J 1 -J 2 spin chains. We assumed that there are two oxygen sites between the two nearest-neighbor magnetic ion sites, such as the CuO 2 chain. When two oxygen sites feel different electric potentials, the DC electric field breaks the one-site translation symmetry down to the two-site one. In other words, the number of spins per unit cell is doubled. The DC electric field then yields the bond alternation δJ 1 (E) j (−1) j S j · S j+1 while δJ 1 (0) = 0. The appearance of the bond alternation drives the quantum critical phase of the spin chain into a unique gapped quantum phase with a dimerization, which is the first proposal of the DC electric-field-induced dimerization.
The E-induced dimer phase is the singlet-dimer one for J 1 > 0 or the Haldane-dimer (triplet-dimer) one for J 1 < 0. The E dependence of the excitation gap will be experimentally visible in, for example, the ESR, the inelastic neutron scattering, and the Raman scattering experiments, which give evidence of the growth of the E-induced dimer orders.
We also briefly investigated other major DC electricfield effects that were not incorporated in our analyses. We hope that this paper stimulates further theoretical and experimental studies on electric-field controls of quantum magnetism.
In addition, some products of operators at d orbitals are rewritten in terms of the S = 1/2 spin operators S j = 1 2 σ,σ d † j,σ σ σσ d j,σ . The contribution (A2) of the process (a) is thus reduced to a simple form, We can obtain H eff;b 4 similarly: Combining Eqs. (A4) and (A5), we reach the final result of the antiferromagnetic superexchange coupling (3.8), which is consistent with the special case of ∆ 0 = ∆ 1 and t 0 = t 1 ∈ R [52].

Ferromagnetic superexchange interaction
The ferromagnetic superexchange coupling (3.11) is similarly derived. Two processes contribute to the fourth-order term (3.5) as shown in Fig. 19. The process (a) of Fig. 19 exchanges spins at j = 0 and j = 1 sites. On the other hand, the process (b) of Fig. 19 does not. Nevertheless, the latter must be taken into account because it gives rise to a term S z 0 S z 1 , as we will see soon. A major difference of the present model (2.5) from the previously dealt one (2.2) comes from the Coulomb exchange (the ferromagnetic direct exchange), J H , that reconstructs the unperturbed eigenstates of the degenerate p orbitals only when they are half occupied. We assumed the presence of two degenerate p orbitals p x and p y . Without the hopping terms, the eigenstates are given by a product state |φ d |φ p , where |φ µ denotes the eigenstate at the µ = d, p orbital. |φ d is further split into a product of local states at two dorbital sites. The same applies to |φ p except for the case mentioned above. When the p orbitals are half occupied, their eigenstates are reconstructed as the singlet |φ p = (|↑ x ↓ y − |↓ x ↑ y )/ √ 2 and triplets, |φ p = |↑ x ↓ y , (|↑ x ↓ y + |↓ x ↑ y )/ √ 2, and |↓ x ↓ y . Here, |σ x σ y denotes the eigenstate of the p orbitals. Note that the triplets have the eigenenergy lower than the singlet by 2J H (> 0).  . 19. Two typical fourth-order processes, (a) and (b), of the perturbative expansion for the model (2.5), similarly drawn to Fig. 18. The red rounded square represents a reconstruction of the p-orbital eigenstates by the Coulomb exchange JH [cf. Fig. 4 (b)]. The p-orbital states surrounded by this curve is either symmetrized or antisymmetrized. The former corresponds to one of the triplet states and the latter to the singlet state. This reconstruction permits hoppings, for instance, from the j = 0 site to the y orbital at the j = 1/2 site, that were forbidden for JH = 0.
The last two terms represent spin-dependent hoppings that appear as a direct consequence of nonzero J H . They vanish in the J H → 0 limit. Discarding the p operators and translating the d operators into spin operators with the aid of the projection P , we obtain