The Stokes–Einstein–Debye (SED) relation [1] ${D}_{r}={k}_{{\rm{B}}}T/\varsigma $ correlates the rotational diffusion constant ${D}_{r},$ rotational frictional coefficient $\varsigma ,$ Boltzmann constant k_B and temperature T. When a rigid sphere with a radius $a$ moves in a fluid with viscosity $\eta ,$ the $\varsigma $ can be described by the Stokes’ formula $\varsigma =8\pi \eta {a}^{3}.$ And the SED relation can be expressed as ${D}_{r}={k}_{{\rm{B}}}T/8\pi \eta {a}^{3}.$

Debye [1, 2] proposed that $\varsigma $ is correlated with rotational relaxation time ${\tau }_{r}$ and the SED can be expressed as the variant ${D}_{r}=1/\left[{\tau }_{rn}n\left(n+1\right)\right],$ where ${\tau }_{rn}$ is determined by the decay of the nth degree Legendre polynomials. If one assumes the effective hydrodynamic radius $a$ for a soft particle is also a constant, one gets another variant ${D}_{r}\sim T/\eta $ [3], where ‘∼’ means proportional. Furthermore one can get the third variant ${D}_{r}\sim T/{\tau }_{t}$ with the approximate relation $\eta ={G}_{\infty }{\tau }_{t}$ [4], where G_∞ is the instantaneous shear modulus presumed to be a constant, and ${\tau }_{t}$ is the structural relaxation time. Therefore, the SED relation has at least four forms, the original formula ${D}_{r}\,={k}_{{\rm{B}}}T/8\pi \eta {a}^{3}$ and the three variants ${D}_{r}\,=1/\left[{\tau }_{rn}n\left(n+1\right)\right],$ ${D}_{r}\sim T/\eta $ and ${D}_{r}\sim T/{\tau }_{t}.$ Moreover, if ${D}_{r}\,=1/\left[{\tau }_{rn}n\left(n+1\right)\right]$ is satisfied, $m\left(m+1\right){\tau }_{rm}=n\left(n+1\right){\tau }_{rn}$ should be established for different n, m. On the other hand, if the Stokes–Einstein (SE) relation ${D}_{t}={k}_{{\rm{B}}}T/6\pi \eta a$ and the assumption of constant $a$ are both satisfied, one gets more variants of SED relation by directly combination of SE and SED relation or its variants, such as the independence of the ratios ${D}_{t}{\tau }_{rn},$ ${D}_{t}/{D}_{r}$ and ${\tau }_{rn}T/\eta $ with the temperature, where ${D}_{t}$ is the translational diffusion constant.

Many studies proposed the SED relation is invalid in supercooled liquids by testing with above variants. A deviation from ${\tau }_{r1}/{\tau }_{r2}=3$ is observed in both analytical models and experiments [5–7]. De Michele and Leporini [8] found the ratio $n\left(n+1\right){\tau }_{rn}/2{\tau }_{r1}\ne 1$ for $n=2,\,3,\,4$ in a supercooled liquid consisted of rigid dumbbells interacting via a Lennard–Jones potential. They found the $n\left(n+1\right){\tau }_{rn}/2{\tau }_{r1}$ and $n\left(n+1\right){D}_{r}{\tau }_{rn}$ firstly decrease to a minimum and then increase with deceasing temperature. The ${\tau }_{r1}{D}_{t}$ in a dumbbell model interacting via a repulsive ramp like potential [9] is not a constant but temperature and density dependent. Similar phenomena are also observed in the supercooled SCP/E water [10–12]. The ratio ${D}_{r}{\tau }_{t}/T$ is almost a constant in supercooled SPC/E water [4] within 280–350 K, but increases with cooling below 280 K. However, a fractional form ${D}_{r}\sim {\left(T/{\tau }_{t}\right)}^{\xi }$ with $\xi =0.75$ is observed for the whole temperature range. The ratio ${D}_{t}/{D}_{r}$ is found not to be a constant but is decreased with decreasing temperature. Kawasaki and Kim [3] observed that ${D}_{t}/{D}_{r}$ is also decreased with cooling in TIP4P water; however, ${D}_{r}{\tau }_{rn}$ shows a reverse trend for $n=1,\,2,\,3,\,6.$ They found the ${\tau }_{r6}T/\eta $ is almost established but a fractional form ${\tau }_{r1}\sim {\left(\eta /T\right)}^{\xi }$ is observed with $\xi =0.8.$

Stillinger and co-workers [13] found the scaled ratio ${D}_{t}/{D}_{r}$ is almost equal to 1.0 in ortho-terphenyl (OTP) within 260–346 K, but it deviates from 1.0 below 260 K and the deviation gets larger with cooling. Moreover, a controversial result is observed that the directly simulated ${D}_{t}/{D}_{r}$ shows an opposite trend with the data deduced from ${\tau }_{r2}$ below 260 K. Sillescu and co-workers [14] observed the ${D}_{r}$ of OTP molecule is still proportional to ${\eta }^{-1}$ down to the glass transition temperature ${T}_{g},$ but the ${D}_{t}$ is only proportional to ${\eta }^{-1}$ for $T\gt 1.2{T}_{g}$ and ${D}_{t}\sim {\eta }^{-0.75}$ for $T\lt 1.2{T}_{g}$. Cicerone and Ediger [15] exploited four probes in OTP. They found the ${D}_{t}\sim T/\eta $ and ${D}_{t}{\tau }_{r}$ are size dependent, the SE variant ${D}_{t}\sim T/\eta $ is broken down with cooling and the SED variant ${D}_{t}{\tau }_{r}$ for the small probe increases almost two orders of magnitude as ${T}_{g}$ is approached. The different behavior observed in SED and SE relation is usually called the decoupling of the transitional and rotational motion. There exist several different explanations, such as the mechanism changes for the transitional diffusion arise around $1.2{T}_{g}$ [14, 16, 17], spatially heterogeneous dynamics [15], various sized structured domain [18] or fluidized domain [19] with different dynamics forms in supercooled liquids.

Although there are many studies that suggest the breakdown of SED relation in supercooled liquids, no study directly tests the original SED relation ${D}_{r}={k}_{{\rm{B}}}T/8\pi \eta {a}^{3}.$ This is questionable because the equivalence of the variants to the original SED formula is on the basis of various assumptions, while there is some evidence showing that those assumptions may not be valid all the time. For instance, the $a$ for organic molecules varies with the volume fraction in their diluted solutions [20], and the behavior of ions in aqueous solutions is observed to deviate from the SE relation by taking $a$ as a constant but the original SE relation actually holds if a is allowed to change [21–25]. Moreover, there exist simulations [26, 27] indicate the Einstein relation ${D}_{t}={k}_{{\rm{B}}}T/\alpha $ is valid for several supercooled liquids, and the $a$ should be considered to be varied with temperatures on the basis of the validity of the Stokes’ formula $\alpha =6\pi \eta a.$ The relation $\eta ={G}_{\infty }{\tau }_{t}$ is approximately established only when the memory effect is exponential [2]; however, it is found that the structural relaxation follows a non-exponential decay in supercooled liquids due to the dynamic heterogeneity [28–30]. In addition, there exist many studies that show the breakdown of SE relation variants [4, 31–35], it is questionable to test the SED relation by using the combination of SED and SE relation or its variants. So the validity of the SED relation is still elusive and it is necessary to consider the original formula. In this work, we explore the SED relation from its original form and variants to verify their validity by performing molecular dynamics (MD) simulations with the Lewis–Wahnstrom model of OTP [34, 36].

The present work is based on our previous work [27], the OTP model adopted [34, 36] and simulation details are the same. The frictional coefficient $\alpha $ and viscosity $\eta $ are directly taken from [27], which are also plotted in figures 1(a) and (b) for convenience. The three sites of the OTP molecule are named after A, B and A; their charges are ${q}_{A}=0$ and ${q}_{B}=0,$ respectively. To explore the possible influences of torque on the SED variants, other two charged systems are simulated with ${q}_{A}\,=0.02,\,{q}_{B}=-0.04$ and ${q}_{A}=0.04,\,{q}_{B}=-0.08$ in unit e, respectively. To improve statistics, seven independent trajectories have been simulated to determine the structural relaxation time ${\tau }_{t},$ rotational diffusion constant ${D}_{r},$ rotational correlation time ${\tau }_{rn}$ and rotational non-Gaussian parameter ${\alpha }_{2}\left(t\right).$

The structural relaxation is described by the self-intermediate scattering function [37]where N is the number of molecules, wavevector $k=14.5\,{\mathrm{nm}}^{-1}$ corresponding to the first maximum of the static structure factor, ${\mathop{r}\limits^{\longrightarrow}}_{j}\left(t\right)$ is the position of center of mass for the jth molecule, $\langle \rangle $ denotes time average, and ${\tau }_{t}$ is determined by ${F}_{s}\left(k,{\tau }_{t}\right)={e}^{-1}.$

The rotational diffusion constant ${D}_{r}$ is calculated via its asymptotic relation with the rotational mean square displacement (RMSD) [3, 4]where $\left\langle {\mathop{\varphi }\limits^{\longrightarrow}}^{2}\left({\rm{\Delta }}t\right)\right\rangle =\left(1/N\right){\displaystyle {\sum }_{i=1}^{N}\left|{\mathop{\varphi }\limits^{\longrightarrow}}_{i}\left(t+{\rm{\Delta }}t\right)-{\mathop{\varphi }\limits^{\longrightarrow}}_{i}\left(t\right)\right|}^{2}$ is the RMSD for the displacement ${\mathop{\varphi }\limits^{\longrightarrow}}_{i}\left({\rm{\Delta }}t\right).$ The ${\vec{\varphi }}_{i}\left({\rm{\Delta }}t\right)\,={\int }_{t}^{t+{\rm{\Delta }}t}{\vec{\omega }}_{i}\left(t\right){\rm{d}}t,$ where the angular velocity ${\mathop{\omega }\limits^{\longrightarrow}}_{i}\left(t\right)$ of ith molecule is in magnitude $\left|{\mathop{\omega }\limits^{\longrightarrow}}_{i}\left(t\right)\right|={\cos }^{-1}\left[{\mathop{e}\limits^{\longrightarrow}}_{i}\left(t\right)\cdot {\mathop{e}\limits^{\longrightarrow}}_{i}\left(t+{\rm{\Delta }}t\right)\right]$ and direction ${\mathop{e}\limits^{\longrightarrow}}_{{\omega }_{i}}={\mathop{e}\limits^{\longrightarrow}}_{i}\left(t\right)\times {\mathop{e}\limits^{\longrightarrow}}_{i}\left(t+{\rm{\Delta }}t\right)/{\rm{\Delta }}t,$ ${\mathop{e}\limits^{\longrightarrow}}_{{\omega }_{i}}$ is the unit vector of the bisector of ith angle ∠ABA. A time interval ${\rm{\Delta }}t=0.01$ ps is adopted to calculate the RMSD.

The rotational correlation time ${\tau }_{rn}$ is calculated via the rotational correlation function [1, 3]where ${P}_{n}\left(x\right)$ is the nth order Legendre polynomial, and ${\tau }_{rn}$ is determined by ${C}_{n}\left({\tau }_{rn}\right)={e}^{-1}.$

The rotational dynamics are heterogeneous and are characterized by the rotational non-Gaussian parameter [38]

To examine the SED relation and its variants, the viscosity $\eta ,$ frictional coefficient $\alpha ,$ rotational relaxation time ${\tau }_{rn}$ for n = 1, 2, 6, structural relaxation time ${\tau }_{t}$ and rotational diffusion constant ${D}_{r}$ at different temperature T are calculated and plotted in figure 1.

The variant ${D}_{r}\sim T/\eta $ behaves as a fractional form ${D}_{r}\sim {\left(T/\eta \right)}^{{\xi }_{1}}$ with ${\xi }_{1}=0.61$ as shown in figure 2(a), which deviates 0.39 from the exact result and indicates the breakdown of ${D}_{r}\sim T/\eta .$ The exponent ${\xi }_{1}=0.61$ is smaller than the $\chi \approx 0.9$ in ${D}_{t}\sim {\left(T/\eta \right)}^{\chi }$ in [27], which implies ${D}_{t}/{D}_{r}\sim {a}^{2}\sim {\left(T/\eta \right)}^{\chi -{\xi }_{1}}$ should decrease with decreasing temperature. The decreasing of ${D}_{t}/{D}_{r}$ with cooling is also observed in the simulation of the TIP4P water [3] and OTP [13]. Similar breakdown is observed in ${D}_{r}\sim T/{\tau }_{t}$ but the fractional form ${D}_{r}\sim {\left(T/{\tau }_{t}\right)}^{{\xi }_{2}}$ has an exponent ${\xi }_{2}=0.49$ as plotted in figure 2(b). The ${\xi }_{1}=0.61$ and ${\xi }_{2}=0.49$ show the relation $\eta ={G}_{\infty }{\tau }_{t}$ is not exact, the similar phenomena are observed in testing SE relation in KSCN aqueous solutions [23], TIP5P water [26] and OTP [34]. The $\eta ={G}_{\infty }{\tau }_{t}$ is an approximate relation and is only valid when the system gets an exponential decay. However, the supercooled liquids usually follow a non-exponential decay for the dynamic heterogeneity [28, 30].

The ${D}_{r}=1/\left[{\tau }_{rn}n\left(n+1\right)\right]$ is tested by ${D}_{r}\sim {\tau }_{rn}^{-1}$ for n = 1, 2, 6. The three are all in fractional forms as ${D}_{rn}\sim {\tau }_{rn}^{-{\xi }_{3}}$ with ${\xi }_{3}\approx 0.9$ as plotted in figure 2(c). The breakdown is small and ${D}_{r}\sim {\tau }_{rn}^{-1}$ is approximately valid. It is similar as that observed by Cicerone and Ediger [15] with four probes in OTP. A little different from that observed by Sillescu and co-workers [14], who found the variant ${D}_{t}/{D}_{r}$ is valid even close to the ${T}_{g}$ with ${D}_{r}$ deduced from ${\tau }_{r}.$ The differences may be resulted from the different methods adopted to get ${\tau }_{r}$ in simulation and experiment. Different from the TIP4P water [3], the ${D}_{r}\sim {\tau }_{r1}^{-{\xi }_{3}}$ is breakdown with ${\xi }_{3}\approx 0.8$ and ${D}_{r}\sim {\tau }_{r6}^{-{\xi }_{3}}$ is almost established with ${\xi }_{3}\approx 1.0.$ Comparing ${D}_{t}\sim {\left(T/\eta \right)}^{\chi }$ and ${D}_{r}\sim {\tau }_{r2}^{-{\xi }_{3}},$ the two exponents are almost equal. However, the increases of $\eta $ are much faster than ${\tau }_{r2}$ with decreasing temperature as shown in figure 1, and one can explain why an opposite trend is observed in ${D}_{t}/{D}_{r}\,{\rm{vs}}\,T$ for the simulated result and data deduced from ${\tau }_{r2}$ in [13].

By assuming establishment of the Stokes’ formula $\alpha =6\pi \eta a,$ the $a$ is deduced by $a\sim \alpha /\eta $ and the original form ${D}_{r}={k}_{{\rm{B}}}T/8\pi \eta {a}^{3}$ can tested by ${D}_{r}{\alpha }^{3}\sim {\left(T{\eta }^{2}\right)}^{{\xi }_{4}}.$ The directly fitted exponent ${\xi }_{4}$ is 1.06, which is so close to the exact result marked by the red dashed with ${\xi }_{4}=1.0$ as plotted in figure 2(d). The result indicates ${D}_{r}={k}_{{\rm{B}}}T/8\pi \eta {a}^{3}$ is valid. Comparing with ${D}_{r}\sim T/\eta ,$ one can see the adopted assumption of constant $a$ is not appropriate in testing the SED relation. The result is consistent with the [27], the correlation among molecules plays an important role in $a$ and makes it vary with temperature. The SE relation for OTP is valid after taking the changes of $a$ into account, although its variants are all broken down.

The above results indicate the variants ${D}_{r}\sim T/\eta $ and ${D}_{r}\sim T/{\tau }_{t}$ are definitely invalid, and whether ${D}_{r}\sim {{\tau }_{rn}}^{-1}$ is really breakdown. If ${\mathop{\varphi }\limits^{\longrightarrow}}_{i}\left({\rm{\Delta }}t\right)$ follows a Gaussian distribution, the ${D}_{r}\sim {\tau }_{r1}^{-1}$ is an exact result due to $\left\langle {P}_{1}\left(\cos \,\varphi \right)\right\rangle \,=\left\langle {{\rm{e}}}^{{\rm{i}}\varphi }\right\rangle ={{\rm{e}}}^{-\left\langle {\varphi }^{2}\right\rangle /2}.$ However, the rotational dynamics are also observed to deviate from Gaussian and show a heterogeneous dynamics as the translational dynamics [8, 38–40]. The non-Gaussian parameter ${\alpha }_{2}\left(t\right)$ plotted in figure 3(a) is similar as that observed in SPC/E water [38], which deviates from zero and the maximum increases with cooling. The result indicates the system gets more heterogeneous dynamics at a lower temperature. Therefore the ${D}_{r}\sim {\tau }_{r1}^{-1}$ should not be exactly established.

On the other hand, if a molecule rotates without net external torque, the probability distribution of the chosen unit vector $\mathop{e}\limits^{\longrightarrow}$ is $\rho \left[\mathop{e}\limits^{\longrightarrow}\left(t\right)\right]=\displaystyle \displaystyle {\sum }_{n,m}{{\rm{e}}}^{-n\left(n+1\right){D}_{r}t}{{\rm{Y}}}_{n}^{m}\left[\mathop{e}\limits^{\longrightarrow}\left(t\right)\right]{{\rm{Y}}}_{n}^{-m}\left[\mathop{e}\limits^{\longrightarrow}\left(0\right)\right]$ [2], and ${D}_{r}=1/\left[{\tau }_{rn}n\left(n+1\right)\right]$ is an exact result, where ${{\rm{Y}}}_{n}^{m}$ is the spherical harmonic function. Because of the dynamic heterogeneity, the system is heterogeneous both in dynamics and structure, and the more mobile particles are likely to from a string structure [30, 38]. So the net torque applied on a molecule may not be zero due to the interaction among molecules even without an external applied torque. To explore the possible influence introduced by the interaction among molecules, we simulate two other systems consisting of polar molecules. The three sites of OTP in the two systems are taken small electric charges with ${q}_{A}=0.02,\,{q}_{B}=-0.04$ and ${q}_{A}=0.04,{q}_{B}=-0.08,$ respectively. The three systems are labeled as $q=0,$ $q=0.02$ and $q=0.04,$ respectively, for simplicity.

Figure 3(b) shows the ${\alpha }_{2}\left(t\right)$ for $q=0.02$ and $q=0.04$ also deviate from Gaussian but not much differences are observed compared with the $q=0.$ However, the fractional form ${D}_{rn}\sim {\tau }_{rn}^{-{\xi }_{3}}$ for n = 1, 2, 6 has a smaller exponent ${\xi }_{3}$ for a larger q as identified by the figures 3(c) and (d). For instance, the ${\xi }_{3}$ for n = 1 is 0.913 for $q=0,$ 0.82 for $q=0.02,$ and 0.78 for $q=0.04.$ The results indicate the interaction among molecules play an important role in the breakdown of ${D}_{r}\sim {\tau }_{rn}^{-1}$ other than the heterogeneous dynamics, which may lead the probability distribution of the chosen unit vector $\mathop{e}\limits^{\longrightarrow}$ deviate from $\rho \left[\mathop{e}\limits^{\longrightarrow}\left(t\right)\right]\,=\displaystyle {\sum }_{n,m}{{\rm{e}}}^{-n\left(n+1\right){D}_{r}t}{{\rm{Y}}}_{n}^{m}\left[\mathop{e}\limits^{\longrightarrow}\left(t\right)\right]{{\rm{Y}}}_{n}^{-m}\left[\mathop{e}\limits^{\longrightarrow}\left(0\right)\right].$

To further verify the validity of ${D}_{r}\sim {\tau }_{rn}^{-1},$ the relation $n\left(n+1\right){\tau }_{rn}=2{\gamma }_{n}{\tau }_{r1}$ are adopted. If ${D}_{r}\sim {\tau }_{rn}^{-1}$ is valid, ${\gamma }_{n}=1.0$ and otherwise ${\gamma }_{n}\ne 1.0.$ Figure 4 shows the exponents ${\gamma }_{2}$ ≈ 1.194 for $q=0,$ 1.288 for $q=0.02,$ and 1.295 for $q=0.04.$ The ${\gamma }_{6}$ ≈ 0.495 for $q=0$，0.584 for $q=0.02,$ and 0.561 for $q=0.04.$ Both ${\gamma }_{2}$ and ${\gamma }_{6}$ deviate significantly from the exact result ${\gamma }_{n}=1.0$ and imply the breakdown of ${D}_{r}\sim {\tau }_{rn}^{-1}.$ Moreover, the systems with $q\ne 0$ have a larger value of ${\gamma }_{2}$ and ${\gamma }_{6}$ than q = 0, and the larger q has a larger ${\gamma }_{2},$ which also signifies the importance of interaction for the breakdown of ${D}_{r}\sim {\tau }_{rn}^{-1}.$ The result is consistent with the data plotted in figure 3. Combining the results given by figures 2, 3 and 4, we conclude the variant ${D}_{r}\sim {\tau }_{rn}^{-1}$ is explicitly broken down.

In summary, we have examined the rationality of the SED relation and its variants in OTP liquids by performing MD simulations. Our results indicate ${D}_{r}\sim T/\eta $，${D}_{r}\sim T/{\tau }_{t}$ and ${D}_{r}\sim {{\tau }_{rn}}^{-1}$ with n = 1，2，6 are all broken down and in fractional forms. The breakdown of the variant ${D}_{r}\sim T/\eta $ is due to the assumption of constant a is not satisfied. And the variant ${D}_{r}\sim T/{\tau }_{t}$ further adopts a not usually established approximation relation $\eta ={G}_{\infty }\tau .$ The ${D}_{r}\sim {{\tau }_{r1}}^{-1}$ is only approximately established and its breakdown is resulted from the rotational heterogeneous dynamics. The ${D}_{r}\sim {{\tau }_{rn}}^{-1}$ with n = 1，2，6 for $q\ne 0$ gets a stronger breakdown than the $q=0.$ So the interaction among molecules also plays an important role in the breakdown of ${D}_{r}\sim {{\tau }_{rn}}^{-1}$ other than the heterogeneous dynamics. Although the three variants are all broken down, the original SED relation ${D}_{r}={k}_{{\rm{B}}}T/8\pi \eta {a}^{3}$ tested by ${D}_{r}{\alpha }^{3}\sim T{\eta }^{2}$ is established after considering the changes of $a.$ The result is consistent with our previous work for the SE relation in OTP [27], which shows the SE relation is valid by taking the changes of $a$ into account. So no decoupling of translation and rotational motion is observed for OTP within 260–400 K. Our simulations indicate that the $a$ is such an important parameter that is closely connected with conclusions drawn on the validity of SED relation, and the assumption of constant $a$ should be carefully evaluated when testing the SED relation as well as the approximation relations adopted in the variants of SED relation. To explore the applicability of the main idea proposed in this work, our future work could be considering more supercooled liquids, such as the supercooled water.