Effective Cubic Coupling Between Modes?

[irrotational]

I would like to compute a quantity that can measure the (cubic) coupling between modes (I'm interested in $\Gamma$-point modes only), having already obtained the third-order force constants (FC3) and phonon eigenvectors. Something like:

$$ A^{(3)}_{\nu \lambda \sigma} = \sum_{\alpha \beta \gamma} \sum_{i j k} \Phi^{\alpha \beta \gamma}_{i j k} M^{\alpha \beta \gamma}_{i j k} \ \eta^{\nu}_{\alpha i} \ \eta^{\lambda}_{\beta j} \ \eta^{\sigma}_{\gamma k} $$

Where:

$$ \Phi^{\alpha \beta \gamma}_{i j k} $$

is the FC3 tensor (from .fcs file), and:

$$ M^{\alpha \beta \gamma}_{i j k} $$

is the FC3 multiplicity tensor (also from .fcs file), and:

$$ \eta^{\nu}_{\alpha i} $$

is a (normalized) phonon eigenvector at the $\Gamma$-point (from .evec file). $(\nu,\lambda,\sigma)$ index phonon branches, $(\alpha,\beta,\gamma)$ index atoms and $(i,j,k)$ index Cartesian directions.

Or perhaps it is more appropriate to use this, which accounts for the frequency/mass factors in the mode length scale:

$$ A^{(3)}_{\nu \lambda \sigma} = \frac{1}{\sqrt{\omega_{\nu} \ \omega_{\lambda} \ \omega_{\sigma}}} \sum_{\alpha \beta \gamma} \sum_{i j k} \frac{ \Phi^{\alpha \beta \gamma}_{i j k} M^{\alpha \beta \gamma}_{i j k} \ \eta^{\nu}_{\alpha i} \ \eta^{\lambda}_{\beta j} \ \eta^{\sigma}_{\gamma k} } {\sqrt{m_{\alpha} \ m_{\beta} \ m_{\gamma}}} $$

Both approaches are essentially projecting the (Cartesian) FC3 tensor onto the normal modes. The hope is that modes with 'strong cubic coupling' will have a large A value, whereas 'weakly cubic coupled' modes will not. Of course, this can be extended to higher orders (FC4, etc).

I am certain that this will have been discussed in the literature, but I cannot find a rigorous definition, and am uncertain about the mass/frequency factors. The idea is to investigate which modes could potentially contribute to third-harmonic generation.

[ttadano]

Options to printing out the anharmonic coupling coefficients are already implemented in ALAMODE, albeit not supported officially. Please add the following options with MODE=RTA and run anphon.

&kpoint
2
1 1 1 # Just gamma point is fine for the current purpose.
/

&analysis
 PRINTV3 = 2 # 0: Not computed, 1: Print V3 (with frequency, mass factors), 2: Print \Phi_3 (w/o frequency, mass factors)
 PRINTV4 = 2 # The same for the fourth-order couplings
 QUARTIC = 1 # Necessary for computing the fourth-order
 KS_INPUT = mode.in # a file defining target modes to compute the couplings
 PROJECTION_AXES = 1 0 0, 0 1 0 # Sometimes helpful to project degenerate modes along specific directions
/

In mode.in, please list the phonon momentum and branch index of the mode $q_1$ you want to compute the couplings $\Phi(q_1,q_2,q_3)$, etc. as follows:

15
0 0 0 1
0 0 0 2
0 0 0 3
0 0 0 4
0 0 0 5
0 0 0 6
0 0 0 7
0 0 0 8
0 0 0 9
0 0 0 10
0 0 0 11
0 0 0 12
0 0 0 13
0 0 0 14
0 0 0 15

The first row defines the number of target modes, and the following 15 lines specify (qx, qy, qz) and the branch index (1-base).

After the anphon calculation finishes, the coupling terms are saved in PREFIX.V4.XXX or PREFIX.Phi4.XXX where XXX represents the index of the target mode. For example, below is the example output of Phi4 computed for lowest-lying optical mode of cubic SrTiO3:

# xk =               0              0              0
# mode = 4
# Frequency = 127.391
## Matrix elements Phi4 for given mode
## q1, j1, omega(q1j1) (cm^-1), q2, j2, omega(q2j2) (cm^-1), q3, j3, omega(q3j3) (cm^-1), Phi4(qj,q1j1,q2j2,q3j3) (Ry/(u^{1/2}Bohr)^{4}), multiplicity
    1    1    1.09737e-10    1    1    1.09737e-10    1    1    1.09737e-10    4.42723e-21   -1.12829e-24    1
    1    1    1.09737e-10    1    1    1.09737e-10    1    2    1.09737e-10   -1.77451e-22   -9.90776e-26    1
    1    1    1.09737e-10    1    1    1.09737e-10    1    3    1.09737e-10    9.49443e-20   -3.33297e-23    1
    1    1    1.09737e-10    1    1    1.09737e-10    1    4        127.391    1.57593e-19   -3.03119e-23    1
    1    1    1.09737e-10    1    1    1.09737e-10    1    5        127.391     2.6777e-22   -2.97621e-25    1
    1    1    1.09737e-10    1    1    1.09737e-10    1    6        127.391   -5.33185e-21    9.05408e-25    1
    1    1    1.09737e-10    1    1    1.09737e-10    1    7        183.098    1.71485e-19   -3.00777e-23    1
    1    1    1.09737e-10    1    1    1.09737e-10    1    8        183.098    -6.7424e-23    6.54282e-27    1
    1    1    1.09737e-10    1    1    1.09737e-10    1    9        183.098   -4.64226e-23    5.17693e-27    1
    1    1    1.09737e-10    1    1    1.09737e-10    1   10        269.071    7.12564e-21   -2.93637e-24    1
    1    1    1.09737e-10    1    1    1.09737e-10    1   11        269.071   -8.91432e-22   -1.84891e-25    1
    1    1    1.09737e-10    1    1    1.09737e-10    1   12        269.071    6.95612e-22   -1.21882e-26    1
    1    1    1.09737e-10    1    1    1.09737e-10    1   13        548.093   -1.05744e-19    2.10425e-23    1
    1    1    1.09737e-10    1    1    1.09737e-10    1   14        548.093   -1.51701e-21    8.23223e-26    1
    1    1    1.09737e-10    1    1    1.09737e-10    1   15        548.093   -6.26322e-21    3.35537e-24    1

    1    1    1.09737e-10    1    2    1.09737e-10    1    1    1.09737e-10   -7.63865e-22    7.98099e-27    1
    1    1    1.09737e-10    1    2    1.09737e-10    1    2    1.09737e-10    1.04311e-21    1.87353e-26    1
    1    1    1.09737e-10    1    2    1.09737e-10    1    3    1.09737e-10   -1.99424e-22   -3.23542e-24    1
    1    1    1.09737e-10    1    2    1.09737e-10    1    4        127.391    9.94129e-22   -7.64668e-24    1
    1    1    1.09737e-10    1    2    1.09737e-10    1    5        127.391   -2.61138e-21     3.1633e-24    1
    1    1    1.09737e-10    1    2    1.09737e-10    1    6        127.391     7.7752e-22    3.22376e-26    1
    1    1    1.09737e-10    1    2    1.09737e-10    1    7        183.098   -1.15378e-21   -5.99167e-24    1
    1    1    1.09737e-10    1    2    1.09737e-10    1    8        183.098    4.05252e-22   -5.59434e-26    1
    1    1    1.09737e-10    1    2    1.09737e-10    1    9        183.098   -2.05387e-22    6.04247e-27    1
    1    1    1.09737e-10    1    2    1.09737e-10    1   10        269.071    -6.0208e-22    2.67008e-26    1
    1    1    1.09737e-10    1    2    1.09737e-10    1   11        269.071     2.0398e-21    4.01009e-25    1
    1    1    1.09737e-10    1    2    1.09737e-10    1   12        269.071    2.79639e-21    5.83056e-25    1
    1    1    1.09737e-10    1    2    1.09737e-10    1   13        548.093   -9.83531e-22    5.22683e-24    1
    1    1    1.09737e-10    1    2    1.09737e-10    1   14        548.093   -1.89714e-21    3.83145e-26    1
    1    1    1.09737e-10    1    2    1.09737e-10    1   15        548.093    7.02582e-22   -1.32904e-25    1

The last two float values are real and imaginary parts of the coupling coefficient.