Nonlinear effects in many-body van der Waals interactions

Van der Waals interactions are ubiquitous and they play an important role for the stability of materials. Current understanding of this type of coupling is based on linear response theory, while optical nonlinearities are rarely considered in this context. Many materials, however, exhibit strong optical nonlinear response, which prompts further evaluation of dispersive forces beyond linear response. Here we present a $\textit{Discrete Coupled Nonlinear Dipole}$ approach that takes into account linear and nonlinear properties of all dipolar nanoparticles in a given system. This method is based on a Hamiltonian for nonlinear dipoles, which we apply in different systems uncovering a complex interplay of distance, anisotropy, polarizibilities, and hyperpolarizabilities in the vdW energy. This investigation broadens our basic understanding of dispersive interactions, especially in the context of nonlinear materials.

Introduction -Van der Waals (vdW) interactions play an essential role for the stability of materials with chemically inert components [1,2].The vdW interaction is important for the stacking patterns of layered materials [3,4], and it can even lead to significant changes in the electronic structure as is the case for the directto-indirect band gap transition when comparing a monolayer MoS 2 and its bulk counterpart [5].First principles computational methods for calculating vdW interactions within linear response theory have also been developed [6,7] with implementations in codes, such as VASP [8,9] and Quantum Espresso [10].In these packages, however, vdW calculations rely on linear response theory, which may be problematic for materials with optical nonlinearities.
Recently, several noncentrosymmetric transition metal dichalcogenides [11][12][13][14] or Weyl semimetals [15][16][17][18] have been found to have much enhanced second-order nonlinear hyperpolarizability, which is important in novel applications for the coherent control of spin and valley polarized currents.Many of these materials can also exhibit large third-order optical nonlinearities, which is beneficial for ultrafast optics and plasmonics [19,20].Given that optical properties play an essential role in the collective nature of vdW interactions, one expects that hyperpolarizbilities will also affect the coupling between materials.Indeed, recent electrodynamical studies have shown that even at large separations the interaction between macroscopic bodies experiences magnitude and/or scaling law modulations depending on their nonlinear properties [21][22][23].VdW quadrupolar contributions have also been examined with the random phase approximation showing that these beyond-dipolar effects cannot be neglected in molecular systems [24].
Main-stream computational schemes currently do not take into account optical nonlinearities [8][9][10], which may be a hurdle for a more accurate calculations for electronic structure properties and phenomena directly related to vdW interactions.It is known that even within linear response theory, the vdW coupling energy can have vastly different magnitude, scaling law, and/or sign than the typical attractive two-body London-type 1/R 6 behavior [25][26][27].A striking example are carbon materials, where the vdW forces can be quite different for carbon nanotubes, graphene layers, or graphene nanoribbons [28][29][30].This is a consequence of the complex interplay between many factors: separation, atomic distributions, optical response, and many-body effects.
Here we present the Discrete Coupled Nonlinear Dipole method as a generalized framework for computing dispersive interactions between nanoparticles.In this manybody approach the linear and nonlinear optical properties for all interacting nanoparticles are taken into account explicitly.For this purpose, we obtain the quantum mechanical Hamiltonian Ĥ for interacting nonlinear dipoles fluctuating around their equilibrium positions, which was previously unavailble.Our scheme relies on a diagnonalization procedure coupled with perturbation theory to obtain the eigenmodes of Ĥ whose sum is the vdW interaction energy.The method allows a miscroscopic dissection of dispersive coupling in terms of anisotropy, tensorial nature of the optical response, inversion symmetry, and distance separations for interacting materials.
Theoretical Model -We consider a system composed of N nanoparticles with no external fields applied.Each nanoparticle has a dipolar moment pi (i = 1, 2, ..., N ) generated by the induced electric field from the other nanoparticles.The Hamiltonian of this system is Ĥ = pi T ij pj , where the displacement matrix is determined by the separa- Since previously the quantum mechanical Hamiltonian for an individual nanoparticle Ĥj is available only arXiv:2307.13607v1[cond-mat.mtrl-sci]23 Jun 2023 from linear response, we begin with the derivation of Ĥj for a single nonlinear dipole.The starting point is a classical dipole for particle i in an electric field E: p i = α α α i E + β β β i EE + γ γ γ i EEE + ..., where α α α i is the polarizability tensor of rank 2 from linear response, while β β β i , γ γ γ i , ... are second-order, third-order,... hyperpolarizabilities (rank 3, 4,... tensors) [31][32][33][34].From the classi-cal Hamilton and the equations of motion for p i and its canonical momentum π i taken under equilibrium conditions, a self-consistent nonlinear equation for p i is found.Further, assuming small dipole moments oscillating with a characteristic frequency ω 0 , the solutions to the Hamilton equations combined with a standard quantization procedure are used to transform the classical Hamiltonian into its quantum mechanical equivalent, Ĥ= 1 2 where m, n, q, p = 1, 2, . . ., 3N denote the degrees of freedom associated with x, y, z directional property components of the N nanoparticles.Also, m+2 A detailed derivation of Ĥj and Ĥ is given in the Supplementary Information [35].The first term in Eq. ( 1) contains the kinetic energy of the dipoles described by a 1 × 3N column matrix for the canonical momenta Π with components Πm .The second term corresponds to the linear response properties of the fluctuating dipolar moments arranged in a 1 × 3N column matrix P with components Pm [36,37].The third term contains the β β β properties only, while the last term includes both β β β and γ γ γ nonlinearities.
The Hamiltonian in Eq. ( 1) sets the stage for the vdW energy calculations that takes into account the linear and nonlinear response of each nanoparticle.The schematic flowchart in Fig. 1 gives the main steps of the method discussed in what follows.After identifying the optical properties of each nanoparticle, we represent Pm = P 0,m + Qm , where the equilibrium position P 0,m is found from ∂ Ĥ ∂ Pm = 0 and Qm is the fluctuation of each dipole around P 0,m .This enables transforming the Hamiltonian into Ĥ = E 0 + Ĥh + Ĥanh where E 0 = Ĥ(P 0,m ) is the minimum free energy at equilibrium, Ĥh contains Πm Πn and Qm Qn terms, while Ĥanh consists of Qn Qm Qq and Qn Qm Qq Qp terms [38].From the diagonalization of Ĥh , we find a system of coupled equations for the dipolar fluctuations, [A mn + 2B mnq P 0,q + 3G mnqp P 0,q P 0,p ] F nl u l = ω 2 u l (5) where the eigenvectors u l satisfy the normalization condition F mn u mq u np = δ qp .The eigenvalues from Eq. ( 5) constitute the zero-point modes which give the zeroth order ground state energy E (0) 2 ω n .Further, considering that linear response is always stronger than the optical nonlinearities, Ĥanh is treated as a perturbation finding that where ∆E (2) 2nd is the second-order correction coming from Qn Qm Qq terms in Ĥanh (the first order correction is zero) and ∆E (1) 3rd is the first-order correction from Qn Qm Qq Qp terms in Ĥanh (details in the Supplementary Information [35]).The approach in Fig. 1 constitutes the Discrete Coupled Nonlinear Dipole method, which can now be applied to any system composed of discrete nonlinear dipoles.Note that when B mnq = G mnqp = 0, one recovers the linear Discrete Coupled Dipole method [39], which was already successfully used in first principles many-body schemes [40][41][42][43][44] for calculating vdW energy within linear dipolar interactions.The generalized method described here provides means to account explicitly for linear and nonlinear properties in many-body vdW interactions.
In addition to R 0 and g γ , we introduce . The vdW dependence on the optical axis relative orientation is tracked by the angle ψ (inserts in Fig. 2d-f ).The results in Fig. 2d-f show that the vdW energy de-pends collectively on the anisotropy (parameter ), the optical nonlinearities (g β and g γ parameters), and the optical axis orientation (angle ψ).For particles separated along the x-axis, the nonlinear response leads to a reduced attraction as compared to the London formula.This reduction is strongest for g β and g γ with opposite signs (Fig. 2d).For particles separated along the y-axis and ψ = π 2 , g β and g γ lead to stronger vdW attraction, while particles on the z-axis and ψ = π 2 experience reduced attraction when compared with U L (Fig. 2f).It is also possible to obtain a repulsive vdW energy when the particles are on the y-axis but their optical axis have a relative enagle ψ = π (Fig. 2e).The numerical results given in Fig. 2d-f are consistent with the analytical expressions found in the limit of R R 0 (see in Section S-V and Fig. S2 on Supplementary Information [35] for detailed derivations and comparisons), The above expression reflects the complexity of the vdW interaction for anisotropic particles.There is always a London-like term U L (R), but rescaled by a factor that depends nontrivially on , g β , g γ , and ψ.The vdW energy also has a term entirely due to the second order nonlinearity, which has a much longer outreach range due its 1/R 3 dependence.This term can be positive or negative (dictated mainly by ψ) and it can be a decisive factor for the energy behavior especially at larger R. Eq. ( 8) further shows that by tailoring the properties of the nanoparticles, the London contribution can be suppressed in which case the interaction is dominated by the long-ranged 1/R 3 dependence as determined by the second order nonlinearity.
The Discrete Coupled Nonlinear Dipole model is also applied to realistic systems taking azulene (C 10 H 8 ) and 4-dimethylamino-4'-nitrostilbene (C 16 H 16 N 2 O 2 ) as representatives for nonlinear dipolar nanoparticles (Fig. 2gl).In addition to different anisotropy, each molecule has hyperpolarizabilities with different tensor components as given in [46,47].Our numerical results show that the interaction has similar features as in the simplified model used earlier.The optical nonlinearity strengthens the vdW attraction in most cases, although repulsion is pos-sible for two azulenes providing they are on the x− axis with an opposite direction of their optical axis (Fig. 2g).In addition to the importance of the different optical properties, the vdW interaction is also a many-body phenomenon.The collective nature of dispersive interaction has been examined by several computational methods at the linear response level [6,7], however within this approach one can also account for the effect of optical nonlinearities.Using the scheme from Fig. 1, we calculate the vdW energy between two chains each with 10 nanoparticles separated by a distance a along the x-axis.In the case of isotropic nanoparticles, Fig. 3a,b shows that U vdW is attractive but g γ < 0 reduces its strength, while g γ > 0 enhances it.In the case of azulene chains whose relative angles are ψ = 0 and π, the results given in Fig. 3d,e show that the interaction is repulsive.This effect is dominated by the second order nonlinearity, since neglecting g β in the simulations turns U vdW into attraction.Note that this is unlike the two-particle interaction (Fig. 2g) for which the U vdW was found to always be attractive.

Conclusion -
The Discrete Coupled Nonlinear Dipole method presented here establishes a computational framework of vdW interactions that takes into account linear and nonlinear optical properties from a microscopic perspective.This approach relies on a newly derived quantum mechanical Hamiltonian of a nonlinear dipole combined with a perturbation theory treating linear properties as dominant when compared to second and third order optical nonlinearities.Given that vdW interactions are truly collective phenomena, many factors need to be considered for calculating interaction energies in a given system, as shown by advanced computational methods based entirely on linear response theory.Our method gives the means to add another layer of complex-ity concerning the role of hyperpolarizability in a system with many interacting components.Indeed, the second order nonlinear hyperpolarizbility plays an especially prominent in the dispersive coupling.In many cases, it can lead to a significant reduction of the vdW attraction, while in other situations it can result in repulsion.This is a consequence of the complex and un-intuitive interplay between anisotropy, relative position of the nanoparticles, and the particular tensor structure of the optical nonlinearity.Perhaps, the most striking signature of β is its much long-ranged outreach (1/R 3 ) as compared to the characteristic two-body London interaction.It is even possible to tailor the materials properties to completely suppress the 1/R 6 behavior and have a situation when the energy is completely dominated by the second order nonlinearity.The many-body effects bring additional complexity to the interaction, as found for the azulenes, for example: an attraction between two molecules with ψ = 0 can turn into a repulsion between chains of the same molecules.
The newly found quantum mechanical Hamiltonian for nonlinear dipoles and the novel computational schemes for interactions between nonlinear dipoles indeed open up new directions for vdW phenomena and their close relations to materials.In addition to the specific molecules examined here with a diverse structure of their β and γ tensors, other particles, such as alkali-halide clusters, have giant polarizabilities resulting in R 0 being 3-4 times bigger than the physical size of the cluster [48].For such materials, taking into account the optical nonlinearities may be crucial for the overall stability of the interacting system.This is also relevant to noncentrosymmetric materials with giant optical nonlinearities.In light of its inter-dependence on α, β, γ, the vdW interaction and its effects on stacking patterns, interlayer distances, bond lengths, etc... must be re-examined for a more accurate description.Our Discrete Coupled Nonlinear Dipole method can easily be applied to such situations and give further qualitative and quantitative understanding of the collective competition of linear and nonlinear many-body effects in vdW interactions in direct relation to equilibrium separations and binding and adhesion energies in materials.This approach can serve as a stepping stone for future design of main-stream computational schemes for more accurate simulations, which might be necessary for nonlinear materials with chemically inert components.

3 being the integer part of m+2 3 tracks
the nanoparticle and its m th dipolar component.A generalized Kronecker delta notation δ abcd... = 1 when a = b = c = d = ... 0 otherwise and the Einstein summation rule are also used in Eqs.(1)-(4).
FIG. 1. Schematic flowchart of the Discrete Coupled Nonlinear Dipole method.