Supplementary Material

The optimized lattice parameters, the different continuous switching paths with different selected coordinate origins for tetragonal BaTiO3 and PbTiO3, the eigenvectors of the three Γ15 modes for cubic BaTiO3 and PbTiO3, the anharmonic coupling coefficients under η = 0.0 ~ 0.022 for cubic BaTiO3, the eigenvalue κ of the ξ in the range of η = 0.0 ~ 0.022 for cubic BaTiO3, and the Wannier-interpolated electronic band structures projected on O-2pz orbitals of BaTiO3 at the energy minimum of u for η = 0.01 and η = 0.02. Additional calculations within GGA-PBEsol functionals and evaluations of the in-plane shrinkage when applying uniaxial tensile strains along z-axis are demonstrated as well.