Rectification of Vibrational Energy Transfer in Driven Chiral Molecules

TL;DR

This paper demonstrates that combining molecular chirality with phase-controlled driving can rectify vibrational energy transfer, with effects depending on frequency, polarization, and temperature, including at room temperature.

Contribution

It introduces a novel method to control vibrational energy flow in chiral molecules using phase-controlled driving, supported by classical models.

Findings

Rectification depends on driving frequency and phase.

Effect persists at room temperature.

Directionality can be controlled by polarization.

Abstract

We show that the combination of molecular chirality and phase-controlled driving can lead to rectification of vibrational energy transfer. We demonstrate this effect using classical models of (1) a single helical chain and (2) a more realistic model of polyethylene double helix. We examine the effect of the driving frequency, polarization, and temperature on this phenomenon. Notably, we find that the direction and magnitude of the observed directionality preference depend on the driving frequency and phase, and that the effect persists at room temperature.

Figures

a and 0.34 for Fig. 2c). On the other hand, linear driving along a certain direction (in this case along $x=y$) does not have a directional preference, as the orange lines are roughly symmetric about the center of the chain. For the angular momentum observable (Fig. 2(b) and 2(d)), we observe the opposite angular momentum polarization on the opposite side of the chain, a result that persists regardless of the driving phase. This observation is in agreement with our recent study of the angular momentum observable in response to white-noise driving [[23](https://arxiv.org/html/2509.00324v2#bib.b

Equations

$$H_1$$

$$V_{\nu }\left(\{{\vec{x}}_j\}\right)=\sum _{j=1}^{N-\nu }\frac{1}{2}K{\left[|({\vec{x}}_j+{\vec{r}}_j)-({\vec{x}}_{j+\nu }+{\vec{r}}_{j+\nu })|-a_{\nu }\right]}^2$$

$${\vec{r}}_j=\left(r_{\perp }\cos (2j\pi /3),r_{\perp }\sin (2j\pi /3),j{r}_{\parallel }\right)$$

$$r_{\perp }=\sqrt{3}{a}_1/6$$

$$r_{\parallel }=\sqrt{3}{a}_1/2.$$

$$a_2=\sqrt{13}{a}_1/2$$

$$a_3=3\sqrt{3}{a}_1/2.$$

$$V_{{\nu }_h}\left(\{{\vec{x}}_j\}\right)=\sum _{j=1}^{N-\nu }\frac{1}{2}({\vec{x}}_j-{\vec{x}}_{j+\nu }{)}^T{\text{K}}_{j,j+\nu }({\vec{x}}_j-{\vec{x}}_{j+\nu })$$

$${\text{K}}_{j,j+1}={\mathbf{R}}_z^{-1}(2\pi j/3){\mathbf{R}}_y^{-1}(\pi /3){\text{K}}_0{\mathbf{R}}_y(\pi /3){\mathbf{R}}_z(2\pi j/3),$$

$${\text{K}}_{j,j+2}=$$

$${\text{K}}_{j,j+3}={\mathbf{R}}_y^{-1}(\pi /2){\text{K}}_0{\mathbf{R}}_y(\pi /2),$$

$${\text{K}}_0=\left(\begin{array}{ccc}K& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right).$$

$$\begin{array}{rl}{H}_2& =\sum _{j=1}^{2N}\frac{{\vec{p}}_j^2}{2m}+\sum _{j=1}^{2N-2}{k}_b(l_j-l_0{)}^2+\sum _{j=1}^{2N-4}{k}_{\theta }({\theta }_j-{\theta }_0{)}^2+\\ & \sum _{j=1}^{2N-6}\sum _{n_j}^4\frac{C_n}{2}\left[1+(-1{)}^{n_{j-1}}\cos (n_j{\varphi }_j)\right]+\\ & \sum _{j=1}^{2N}\sum _{i\ne j}4{ϵ}_{ij}\left[(\frac{\sigma }{r_{ij}}{)}^{12}-(\frac{\sigma }{r_{ij}}{)}^6\right],\end{array}$$

$$\begin{array}{r}\vec{F}(t)=F_0(\text{cos}\omega t,\text{cos}(\omega t+\varphi ),0),\end{array}$$

$$J_L$$

$$J_R$$

$$m{\ddot{\vec{x}}}_1=-{\text{K}}_{1,2}({\vec{x}}_1-{\vec{x}}_2)-{\text{K}}_{1,3}({\vec{x}}_1-{\vec{x}}_3)-{\text{K}}_{1,4}({\vec{x}}_1-{\vec{x}}_4)$$

$$+{\vec{F}}_1-m{\gamma }_1{\dot{\vec{x}}}_1,$$

$$m{\ddot{\vec{x}}}_j=-\sum _{\nu =1,2,3}\left[{\text{K}}_{j-\nu ,j}({\vec{x}}_j-{\vec{x}}_{j-\nu })+{\text{K}}_{j,j+\nu }({\vec{x}}_j-{\vec{x}}_{j+\nu })\right]$$

$$+{\vec{F}}_j-m{\gamma }_j{\dot{\vec{x}}}_j,$$

$$⋮$$

$$m{\ddot{\vec{x}}}_N=-{\text{K}}_{N-1,N}({\vec{x}}_N-{\vec{x}}_{N-1})-{\text{K}}_{N-2,N}({\vec{x}}_N-{\vec{x}}_{N-2})$$

$$-{\text{K}}_{N-3,N}({\vec{x}}_N-{\vec{x}}_{N-3})+{\vec{F}}_N-m{\gamma }_N{\dot{\vec{x}}}_N,$$

$$m\ddot{\mathbb{X}}=\mathbb{K}\cdot \mathbb{X}-m\Gamma \cdot \dot{\mathbb{X}}+\mathbb{F},$$

$$\mathbb{X}={\left({\vec{x}}_1,{\vec{x}}_2,{\vec{x}}_3,...,{\vec{x}}_N\right)}^T,$$

$$\Gamma =\mathrm{Diag}\left[{\gamma }_1,{\gamma }_1,{\gamma }_1,...,{\gamma }_j,{\gamma }_j,{\gamma }_j,...,{\gamma }_N\right],$$

$$\mathbb{F}={\left[{\vec{F}}_1,{\vec{F}}_2,{\vec{F}}_3,...,{\vec{F}}_N\right]}^T.$$

$$-m{\omega }^2\tilde{\mathbb{X}}=\mathbb{K}\cdot \tilde{\mathbb{X}}-im\omega \Gamma \cdot \tilde{\mathbb{X}}+\tilde{\mathbb{F}},$$

$$\tilde{\mathbb{X}}(\omega )=\tilde{\mathbb{G}}(\omega )\tilde{\mathbb{F}}(\omega )$$

$$\tilde{\mathbb{G}}(\omega )=[-\mathbb{K}+im\omega \Gamma -m{\omega }^{2}1{]}^{-1}.$$

$$x_j(t)$$

$$y_j(t)$$

$$z_j(t)$$

$$J_L$$

$$J_R$$

Full text

††preprint: AIP/123-QED

I Introduction

The correlation between linear and angular momenta of charge or energy carriers in chiral media, often referred to as linear-angular momentum locking [1, 2, 3], has been recognized as an important factor governing electronic and spin transport properties. More recently, attention has also shifted to such correlations for nuclear vibrations with angular momentum, often referred to as chiral phonons, which have been observed and analyzed in a variety of systems [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 11, 4, 12, 5]. While most studies have focused on such motions in crystalline systems, our earlier studies have revealed that angular–linear momentum locking can arise not only in crystalline solids but also in isolated chiral molecules [14, 15]. Interest in chiral phonons has been further fueled by proposals implicating them mechanistically in the chiral-induced spin selectivity (CISS) effect [16, 10, 17, 18, 19]. It has also been suggested [20] that this correlation between linear and angular momentum can lead to preferred directionality of energy transfer in such systems but, since phonons with positive and negative angular (and corresponding linear) momenta appear symmetrically a crystal spectral distribution, the actual realization of this phenomenon is not assured.

In this paper we describe analytical and numerical studies of directional energy transport in chiral molecules. Our primary finding is that the combination of molecular chirality and phase-controlled driving can lead to highly directionally asymmetric energy flow without any asymmetry in the system-bath couplings.

The relationship between heat transport and chirality has previously been analyzed from several perspectives; for example, theoretical analyses based on the Boltzmann transport equation suggest that thermal gradients in chiral solids can serve as a source of nuclear angular momentum whose directionality depends on the handedness of the chiral solid [21, 22]. Our recent work has predicted a similar effect, showing that heat transport across chiral molecular chains induces nuclear angular momentum, and that the effect exhibits robust properties, persisting in the steady state in both harmonic and anharmonic model systems [23]. Directional energy transport under mechanical or optical driving is another intriguing manifestation of dynamical chirality that provides a novel route towards directional control of energy flow. We exhibit this effect in two model systems described in previous work [23, 24, 14, 20] (which the single chain model is enhanced with including next-nearest-neighbour and next-next-nearest-neighbour coupling), and we examine its dependence on various parameters. Notably, we find that the effect persists at room-temperature.

II Molecular Models

The models chosen for this study are (1) a single (generally anharmonic) helical chain [23, 20] that we study also in the harmonic approximation, and (2) an anharmonic polyethylene double helix [24, 14, 23].

The first model (henceforth referred to as Model 1) is a chain of atomic sites defined by the Hamiltonian [25, 23]

[TABLE]

[TABLE]

where $N$ is the number of atoms in the chain, $m$ is the atomic mass (taken to be equal across all atoms), $\vec{r_{j}}$ are the equilibrium positions of the atoms and $\vec{x_{j}}$ are the displacements of the atoms (note that this Hamiltonian, as well as the Hamiltonian of Model 2, Eq. 13, are defined for the bare system without any external forces and without coupling to any dissipative bath). The model defines a single-chain helix obtained as follows: Taking the $x$-axis as the direction of a given bond, the next bond is generated by rotating the given bond around the $y$-axis by $60^{\circ}$ followed by rotation around the $z$-axis by $120^{\circ}$. The chirality is enforced by the structural dependent equilibrium bond lengths of nearest-neighbor $a_{1}$, next nearest-neighbors $a_{2}$ and $a_{3}$. The position of each atom can be expressed as

[TABLE]

where

[TABLE]

[TABLE]

Then we have

[TABLE]

[TABLE]

Making the harmonic approximation [23] leads to

[TABLE]

where,

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

Here, $\mathbf{R}_{y}$, $\mathbf{R}_{z}$ denote rotations about the $y$ and $z$-axes, respectively. Note that the harmonic model used in Ref. [23] is a simplified version of this Hamiltonian that includes only the nearest-neighbor interactions (Eq. 8, $\nu=1$). In this harmonic approximation the energy transport dynamics can be evaluated analytically.

The second model (referred to below as Model 2) we consider is a polyethylene double-helix [24, 14] with more realistic bond, angle, dihedral, and inter-chain potentials, and thus can be analyzed only using numerical molecular dynamics only. The Hamiltonian is given by

[TABLE]

where $N$ is the polymerization of a single chain of the double-helix, $\dot{\vec{x}}_{j}$ is the velocity of atom $j$, $l_{j}$,$\theta_{j}$, and $\phi_{j}$ are the values of the $j$th bond-length, angle, and dihedral respectively, and $r_{ij}$ are the inter-atom distances. The parameters $k_{b}$ and $k_{\theta}$ are the spring constants for the bonds and angles, and $C_{n}$ are constants that define the dihedral potentials. The $\sigma$ and $\epsilon_{ij}$ are Lennard-Jones parameters to account for non-bonded interactions (such as between the strands in opposite chains). As in Ref. [26] the above parameters were chosen to fit observed physical properties [27, 28].

Figure 1 shows the construction used in this study, where the chain constructed either according to model (1) or (2) above is oriented along the $z$-direction. A section in the middle of the chain (a single central atom in the case of (1) and the two central atoms in each strand in the case of (2)) is subjected to a coherent driving force in the plane perpendicular to the chain length

[TABLE]

where $F_{0}$ is the amplitude of the driving force, $\omega$ is its frequency, and $\phi$ is a phase lag between the $x$ and $y$-components.

We defined the site-resolved kinetic temperature $T_{j}$ of atom $j$ by equipartition as $T_{j}=\frac{2}{3k_{B}}\langle E_{k,j}\rangle=\frac{m}{3k_{B}}\langle\dot{\vec{x}}_{j}^{2}\rangle$. In most of the calculations described below, the temperature at the two chain ends was taken as $T_{cold}=0~\mathrm{K}$. This choice is easily implemented by subjecting the segments at the chain ends to damping forces of the form $-\gamma m\dot{\vec{x}}_{j}$ where $\gamma$ is a damping constant. Once the system reaches steady-state, the heat currents flowing in either direction are then obtained from

[TABLE]

where $S_{L}$ and $S_{R}$ are groups of atomic sites at the left and right chain ends respectively. These groups were taken large enough (10-20 sites) to insure that the temperature fell to zero by the ends of the chain (see Fig. 2). Interestingly, in one numerical test reported below, we found that the heat rectifying behavior still persists when $T_{\text{cold}}=300$ K.

This heat flux observable defined in Eq. (15) can be obtained numerically via standard velocity Verlet integration of classical trajectories. In the case of model (1) this is achieved via an in-house code using integration timestep $\Delta t=0.25~\mathrm{fs}$, while for (2) it is achieved using the LAMMPS code (the version corrected for heat flux [29, 30]) using integration timestep $\Delta t=1.25$ fs.

In the case of the harmonic model (1), the heat flux can be evaluated analytically as well using the methods discussed in Ref. [23]. In particular, the equations of motion are derived from Eq. (1), supplemented by an external driving force $\vec{F}_{j}$ on each atom j and the dissipation $\vec{f}_{j}=-m\gamma_{j}\dot{\vec{x}}_{j}$

[TABLE]

can be cast in the compact form

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The expression for the strength tensor $\mathbb{K}$ can be found in the appendix section I. Note that in our construction, the elements of $\mathbb{\Gamma}$ are nonzero only for the indices $j\in S_{L}\cup S_{R}$, and the elements of $\mathbb{F}$ are only nonzero for the central driven atom(s), where the respective values have been specified above.

The steady state reached while driving with frequency $\omega$ is described by the Fourier transform of Eq. (17)

[TABLE]

where the tildes denote Fourier conjugates in frequency space, e.g $\tilde{\mathbb{X}}=\tilde{\mathbb{X}}(\omega).$ The solution to Eq. 17 is

[TABLE]

where $\tilde{\mathbb{G}}$ is the Green’s function

[TABLE]

Thus, using Eq. (14) and Eq. (23) in Eq. (22) and then inverting the Fourier transform, the motion of each atom $j$ is obtained

[TABLE]

Using Eq. (24) in Eq. (15) and averaging over time, we obtain

[TABLE]

which is the desired analytical result.

The dimensionless ratio $J_{R}/J_{L}$, computed using either of the methods outlined above, is a measure of thermal rectification; in what follows we report this quantity for both models as a function of the parameters that characterize the driving force in Eq. (14).

III Results

As outlined in the introduction, heat transport in chiral environment generates excitations of atomic motions with angular momentum (“chiral phonons”) whose sign depends on the direction of heat current and the molecular handedness. This behavior, together with the thermal distribution (defined as $\mathrm{(2/3k_{B})\langle E_{K}}\rangle$ where $\mathrm{E_{K}}$ is the atomic kinetic energy) along the chain of Fig. 1 is demonstrated in Fig. 2 for (a-b) Model 1 and (c-d) Model 2. On the one hand, focusing on the temperature profile (Fig. 2(a) and 2(c)), we clearly observe that when the driving force is circular (in this case clockwise), the temperature on the left-hand side is significantly higher than the right-hand side. The same effect holds for the heat flux. The higher temperature to the left of the driving point at Z=0 indicates a larger energy flux going leftward (the heat flux ratio was calculated to be ($J_{L}/J_{R}$=0.021 for Fig. 2a and 0.34 for Fig. 2c). On the other hand, linear driving along a certain direction (in this case along $x=y$) does not have a directional preference, as the orange lines are roughly symmetric about the center of the chain. For the angular momentum observable (Fig. 2(b) and 2(d)), we observe the opposite angular momentum polarization on the opposite side of the chain, a result that persists regardless of the driving phase. This observation is in agreement with our recent study of the angular momentum observable in response to white-noise driving [23], and likewise can be explained as a linear-angular momentum locking effect.

The directional energy transport characterizing these systems is more clearly seen in Figures 3 – 6. Fig. 3 and 4 show the ratio $J_{\text{R}}/J_{\text{L}}$ as a function of the driving frequency for models 1 and 2, respectively. The strong frequency dependence of the rectification phenomenon follows directly from the solution of Eq. 22, which yields the density of phonon states as well as the phonon polarization in the phonon band structure. A detailed plot of the density of states as well as the phonon band structure for the single chain is shown in Ref. [23]. The currents are larger in one direction than in the other when the polarization of the phonon is strongly coupled with the group velocity (so that excitation with a certain polarization and frequency of the system will lead to the phonon generation with certain direction of the group velocity). These results were calculated for the driving phase $\phi=\pi/2$ and amplitude $F_{0}=0.5~\mathrm{kcal}~\mathrm{mol}^{-1}$Å*-1*). Note that Fig. 3 presents numerical results obtained for the full Model 1 as well as results obtained from the analytical solution of the harmonic approximation. The excellent agreement seen indicates the validity of the harmonic approximation for the chosen small amplitude driving.

Remarkably, we find that at various driving frequencies, both ratios $J_{\text{R}}/J_{\text{L}}\gg 1$ and $J_{\text{R}}/J_{\text{L}}\ll 1$ are observed, and such values can differ from the unity by a factor of $10^{2}$–$10^{3}$. We note the spectral quality of this rectifying behavior: the same structure will favor heat transport in different directions depending on the frequency. Furthermore, we find that changing either the helical enantiomer (results not shown) or the polarization of the driving (see discussion in Fig. 4 below) will reverse the directional preference at a particular frequency (i.e. $J_{\text{R}}/J_{\text{L}}\rightarrow J_{\text{L}}/J_{\text{R}}$).

Interestingly, for the case of the double helical Model 2 with chain length $N=98$, a significant rectification is observed (Fig. 4(a-b)) only in the highly twisted (15 twists) structure, while for a mildly twisted (3-twists) structure $J_{\text{R}}/J_{\text{L}}$ is close to 1 at all frequencies.

Fig. 5 show similar results for $J_{R}/J_{L}$, now plotted against the driving phase $\phi$ (see Eq. 14) at fixed frequency. Results are shown for the two models in panels (a) and (b) respectively. We observe the greatest difference from unity when the driving phase is equal to $\pm\pi/2$, and that reversing the phase ($\phi\rightarrow-\phi$) results in the reciprocal rectification ratio ($J_{\text{R}}/J_{\text{L}}\rightarrow J_{\text{L}}/J_{\text{R}}$). The peaks at $\pm\pi/2$ driving phase reflect the nature of circular motion of chiral phonons.

Using the analytical results obtained for the harmonic approximation to Model 1 we can also examine the combined effect of the driving frequency ($\omega$) and phase ($\phi$). Fig. 6 shows a heat map showing the rectification displayed against both $\omega$ and $\phi$. As for Fig. 3, the phase dependence of the rectification effect is related to the structure of chiral phonon modes and the frequency dependence results from both the polarization of the phonon band structure and the phonon density of states.

Finally, to address the extent to which the findings outlined above can be expected to influence room-temperature experiments, we conducted sample Langevin molecular dynamics simulations with $T_{\text{cold}}=300$ K according to the methods described in Ref. [23]. We used the $N=98$ polyethylene double-chain model containing 15 twists described above, and the driving parameters were set to $\phi=\pi/2$, $\omega=79.5$ cm*-1* (2.38 THz), and $F_{0}=4.0~\mathrm{kcal}/\mathrm{mol}^{-1}$Å${}^{-1}.$ We performed five independent simulations of length 10 ns, and we obtain a ratio of $J_{\text{R}}/J_{\text{L}}=32.5\pm 3.1.$ This result strongly suggests that the phenomenon observed in this study is relevant to room-temperature experiment.

IV Conclusion

In this work, we analytically and computationally study the heat rectification phenomenon in a chiral system with phase-controlled driving. We examined two structures: a single helical chain (studied without and with the harmonic approximation) and a polyethylene double helix. We found significant asymmetric heat flux under circularly polarized driving at the center of the system. The strong frequency dependence of this phenomenon allows us to modify the rectification effect by tuning the driving frequency, sometimes even switching the direction of the heat transfer, which is also achievable by inverting the polarization of the circular driving. The effect is highly phase dependent, as it decreases when the driving becomes more linearly polarized and is maximal when the driving is circularly polarized in the plane perpendicular to the chain direction. The effect is also correlated with the structural chirality, such that the effect is much weaker (orders of magnitude) in less chiral systems (a polyethylene double helix with 3 twists) than in more chiral systems (a polyethylene double helix with 15 twist). The agreement between analytical results and molecular dynamics in the harmonic model ensures that the effect does not hinge on numerical artifacts. The effect is robust enough that in a double helix system with realistic interaction potentials and room-temperature baths, we still get a remarkable heat flux ratio (e.g. $J_{R}/J_{L}\sim 30$ at $\omega=79.5~\mathrm{cm}^{-1}$), which is a large effect and indicates the feasibility of experimental realization, perhaps by local driving using circularly polarized light.

Future work will focus on quantum-regime extensions with electron-phonon and electron-spin coupling. Although in our simulations the drive is applied at the chain center as a computational convenience, local circularly polarized driving force can be applied to one end of the molecular chain by introducing near-field circularly polarized light (CPL) and with the form

[TABLE]

which is effectively equivalent to the external force in our computational setup above. Changing the enantiomer will also lead to a similar reversal as changing the polarization of the driving. This chiral rectifier/diode effect gives rise to the possibility of building controllable molecular energy transport networks, and possibly CPL catalysis in chiral systems.

ACKNOWLEDGMENTS

This paper is submitted in honor of Prof. YiJing Yan and his seminal contribution to the research of molecular dynamical processes. The research of A.N. is supported by the National Science Foundation award number 2451953. The research of J.E.S. is supported by the Air Force Office of Scientific Research under award number FA9550-23-1-0368. E.A acknowledges the support of the Massachusetts Institute of Technology Department of Chemistry and Department of Chemical Engineering. We thank Claudia Climent for many helpful discussions.

AUTHOR DECLARATION

Conflict of Interest

The authors have no conflicts to disclose.

DATA AVAILABILITY

The data that support the findings of this study are available from corresponding authors upon reasonable request. Structures corresponding the double-chain model examined above are available at https://github.com/eabes23/polymer_twist/.

APPENDIX: MODEL PARAMETERS

In Model 1, the strength constant for the coupling term is $K=500~\mathrm{N/m}$, the bond length for the nearest coupling is $a_{1}=1.5~$Å, $a_{2}=2.7~$Å, $a_{3}=4.0~$Å. The mass of the atom is $m=12~\mathrm{AMU}$.

The strength tensor can be expressed as

[TABLE]

In Model 2, the force field parameters used are $k_{b}=450$ kcal mol*-1Å-2*, $l_{0}=1.54$ Å, $k_{\theta}=62.1$ kcal mol*-1rad-2*, $\theta_{0}=114^{\circ}$, $C_{1}=1.4110$ kcal mol*-1*, $C_{2}=-0.2708$ kcal mol*-1*, $C_{3}=3.143$ kcal mol*-1*, $C_{4}=0.0$ kcal mol*-1*, $\sigma=3.95$ Å, and $\epsilon_{ij}=0.0912$ kcal mol*-1* (for $|\vec{r}_{i}-\vec{r}_{j}|<14\text{\AA },=0$ otherwise). The particle mass is taken $m=14$ AMU, modeling each CH2 bead of the polyethylene chain as a single interaction site.