. They can be obtained through the second-order response of the density matrix with the spin-orbit mean-field operator as the perturbation. The two parameters κ m,ξ and λ mn,ξ are intra-and interstate coupling constants, respectively, that can be determined numerically with the use of wave function overlaps between states at geometries with different displacements ±δQ ξ as follows