<figure><img src="image_1.png" /></figure>

HAL
open science

Nonlinear effects in manybody van der Waals
interactions

Dai-Nam Le, Pablo Rodriguez-Lopez, Lilia Woods

► To cite this version:

Dai-Nam Le, Pablo Rodriguez-Lopez, Lilia Woods. Nonlinear effects in manybody van der Waals interactions. Physical Review Research, 2024, 6 (1), pp.013289. 10.1103/PhysRevResearch.6.013289. hal-04753576

HAL Id: hal-04753576

https://hal.science/hal-04753576v1

Submitted on 25 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Nonlinear effects in manybody van der Waals interactions

Dai-Nam Le,¹,* Pablo Rodriguez-Lopez,²,³ and Lilia M. Woods ¹,†

¹Department of Physics, University of South Florida, Tampa, Florida 33620, USA

²Área de Electromagnetismo and Grupo Interdisciplinar de Sistemas Complejos (GISC),
Universidad Rey Juan Carlos, 28933 Móstoles, Madrid, Spain

³Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-University of Montpellier, F-34095 Montpellier, France

<figure><img src="image_3.png" /></figure>

(Received 1 June 2023; accepted 26 January 2024; published 15 March 2024)

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 *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, polarizabilities, and hyperpolarizabilities in the vdW energy. This investigation broadens our basic understanding of dispersive interactions, especially in the context of nonlinear materials.

DOI: 10.1103/PhysRevResearch.6.013289

I. 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 direct-to-indirect band-gap transition when comparing a monolayer MoS₂ and its bulk counterpart [5]. Several first-principles 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].

Recently, several noncentrosymmetric transition-metal dichalcogenides [11–14] or Weyl semimetals [15–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 hyperpolarizabilities 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–23]. VdW quadrupolar contributions have also been examined based on 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–10], which may be a hurdle for more accurate calculations of 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–27]. A striking example are carbon materials, where the vdW force can be quite different for carbon nanotubes, graphene layers, or graphene nanoribbons [28–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 many-body 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 $\hat{H}$ for interacting *nonlinear dipoles* fluctuating around their equilibrium positions, which was previously unavailable. Our scheme relies on a diagonalization procedure coupled with perturbation theory to obtain the eigenmodes of $\hat{H}$ whose sum is the vdW interaction energy. The method allows a microscopic dissection of dispersive coupling in terms of anisotropy, tensorial nature of the optical response, and distance separations for interacting materials.

<footnote>*dainamle@usf.edu</footnote>

<footnote>†lmwoods@usf.edu</footnote>
## II. THEORETICAL MODEL

We consider a system composed of N nanoparticles with no external fields applied. Each nanoparticle i has a dipolar moment $\hat{\mathbf{p}}^i$ ($i = 1, 2, \dots, N$) generated by the induced electric field from the other nanoparticles. The Hamiltonian of this system is $\hat{H} = \sum_{i=1}^{N} \hat{H}^i - \frac{1}{2} \sum_{i,j=1}^{N} \hat{\mathbf{p}}^i \mathbf{T}^{ij} \hat{\mathbf{p}}^j$, where the displacement matrix $\mathbf{T}^{ij} = \frac{1}{4\pi\epsilon_0} \frac{3\mathbf{R}^i/\mathbf{R}^j - (R^{ij})^2}{(R^{ij})^3}$ is determined by the separation vector $\mathbf{R}^{ij} = \mathbf{R}^i - \mathbf{R}^j$.

Since previously the quantum mechanical Hamiltonian for an individual nanoparticle $\hat{H}^i$ is available only from linear response, we begin with the derivation of $\hat{H}^i$ for a single nonlinear dipole. The starting point is a classical dipole for

particle i in an electric field $\mathbf{E}$, such that $\mathbf{p}^i = \alpha^i \mathbf{E} + \beta^i \mathbf{EE} + \gamma^i \mathbf{EEE} + \dots$, where $\alpha^i$ is the polarizability tensor of rank 2 from linear response, while $\beta^i, \gamma^i, \dots$ are second-order, third-order, ... hyperpolarizabilities (rank 3, 4, ... tensors) [31–33]. From the classical Hamilton and the equations of motion for $\mathbf{p}^i$ and its canonical momentum $\pi^i$ taken under equilibrium conditions, a self-consistent nonlinear equation for $\mathbf{p}^i$ is found. Further, assuming small dipole moments oscillating with a characteristic frequency $\omega^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:

$$ \hat{H} = \frac{1}{2} \mathbb{F}_{mn} \hat{\Pi}_m \hat{\Pi}_n + \frac{1}{2} A_{mn} \hat{P}_m \hat{P}_n + \frac{1}{3} \mathbb{B}_{mnq} \hat{P}_m \hat{P}_n \hat{P}_q + \frac{1}{4} \mathbb{G}_{mnqp} \hat{P}_m \hat{P}_n \hat{P}_q \hat{P}_p, \quad (1) $$

$$ \mathbb{F}_{mn} = [(\omega^0)^2 \boldsymbol{\alpha}] \delta_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor} - \mathbf{T}_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor} (1 - \delta_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor}), \quad (2) $$

$$ \mathbb{B}_{mnq} = -[\alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1}] \delta_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor \lfloor \frac{q+2}{3} \rfloor}, \quad (3) $$

$$ \mathbb{G}_{mnqp} = [2\alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1} - \alpha^{-1} \otimes \gamma \otimes \alpha^{-1} \otimes \alpha^{-1} \otimes \alpha^{-1}] \delta_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor \lfloor \frac{q+2}{3} \rfloor \lfloor \frac{p+2}{3} \rfloor}, \quad (4) $$

where $m, n, q, p = 1, 2, \dots, 3N$ denote the degrees of freedom associated with $x, y, z$ directional property components of the N nanoparticles. A similar index notation is implied for the tensors $\boldsymbol{\alpha}$ with components $\alpha_{mn}, \boldsymbol{\beta}$ with components $\beta_{mnp}$, and $\boldsymbol{\gamma}$ with components $\gamma_{mnpq}$. Also, $\lfloor \frac{m+2}{3} \rfloor$ being the integer part of $\frac{m+2}{3}$ tracks the nanoparticle and its mth dipolar component. A generalized Kronecker delta notation $\delta_{abcd...} = \begin{cases} 1 & \text{when } a=b=c=d=\dots \text{ and the Einstein summation rule are} \\ 0 & \text{otherwise} \end{cases}$ and the Einstein summation rule are also used in Eqs. (1)–(4). A detailed derivation of $\hat{H}^i$ and $\hat{H}$ is given in Appendix A. The first term in Eq. (1) contains the kinetic energy of the dipoles described by a $1 \times 3N$ column matrix for the canonical momenta $\hat{\Pi}$ with components $\hat{\Pi}_m$. The second term corresponds to the linear response properties of the fluctuating dipolar moments arranged in a $1 \times 3N$ column matrix $\hat{P}$ with components $\hat{P}_m$ [34,35]. The third term contains only second order nonlinearities, while the last term includes both second and third order 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 $\hat{P}_m = P_{0,m} + \hat{Q}_m$, where the equilibrium position $P_{0,m}$ is found from $\frac{\partial \hat{H}}{\partial \hat{P}_m} = 0$ and $\hat{Q}_m$ is the fluctuation of each dipole around $P_{0,m}$. This enables transforming the Hamiltonian into $\hat{H} = E_0 + \hat{H}_h + \hat{H}_{anh}$ where $E_0 = \tilde{H}(P_{0,m})$ is the minimum free energy at equilibrium, $\hat{H}_h$ contains $\hat{\Pi}_m \hat{\Pi}_n$ and $\hat{Q}_m \hat{Q}_n$ terms, while $\hat{H}_{anh}$ consists of $\hat{Q}_n \hat{Q}_m \hat{Q}_q$ and $\hat{Q}_n \hat{Q}_m \hat{Q}_q \hat{Q}_p$ terms [36]. From the diagonalization of $\hat{H}_h$, we find a system of coupled equations for the dipolar fluctuations

$$ [A_{mn} + 2\mathbb{B}_{mnq}P_{0,q} + 3\mathbb{G}_{mnqp}P_{0,q}P_{0,p}]\mathbb{F}_{nl}\mathbf{u}_l = \omega^2\mathbf{u}_l, \quad (5) $$

where the eigenvectors $\mathbf{u}_l$ satisfy the normalization condition $\mathbb{F}_{mn}\mathbf{u}_{mq}\mathbf{u}_{np} = \delta_{qp}$. The eigenvalues from Eq. (5) constitute the

zero-point modes which give the 0th order ground state energy $E_{g,s}^{(0)} = \sum_{n=1}^{3N} \frac{1}{2} \hbar \omega_n$. Further, considering that linear response is always stronger than the optical nonlinearities, $\hat{H}_{anh}$ is treated perturbatively:

$$ E_{g,s} \approx E_0 + E_{g,s}^{(0)} + \Delta E_{2nd}^{(2)} + \Delta E_{3rd}^{(1)}, \quad (6) $$

where $\Delta E_{2nd}^{(2)}$ is the second-order correction coming from $\hat{Q}_n \hat{Q}_m \hat{Q}_q$ terms in $\hat{H}_{anh}$ (the first order correction is zero) and $\Delta E_{3rd}^{(1)}$ is the first-order correction from $\hat{Q}_n \hat{Q}_m \hat{Q}_q \hat{Q}_p$ terms in $\hat{H}_{anh}$ (details in Appendix B).

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

<figure><img src="image_4.png" /><figcaption>FIG. 1. Schematic flowchart of the discrete coupled nonlinear dipole method.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 2. VdW energy rescaled by $\hbar\omega^0$ vs nanoparticle separation rescaled by $R_0 = (\alpha/2\pi\epsilon_0)^{1/3}$ is shown in all panels. Two isotropic nanoparticles with (a) $g_\gamma = 0.1$; (b) $g_\gamma = -0.1$; (c) $g_\gamma = -0.15$ with numerical (num), second relation in Eq. (7) (approx.), and London formula $U_L = -(3/16)\hbar\omega^0(R_0/R)^6$ results. Two anisotropic nanoparticles with $\epsilon = 0.3$, $g_\beta = 0.3$, and $g_\gamma = 0, \pm 0.1$ with a relative angle between their optical axis (d) $\psi = 0$ (solid), $\psi = \pi$ (dashed) and separation along x axis; (e) $\psi = 0$ (solid), $\psi = \pi$ (dashed) and separation along y axis; (f) $\psi = \pi/2$ with separation along x axis (solid) and $\psi = \pi/2$ with separation along z axis (dashed). Two azulenes ($C_{10}H_8$) with (g) $\psi = 0$ (blue) and $\psi = \pi$ (red) and separation along x axis; (h) $\psi = 0$ (blue) and $\psi = \pi$ (red) and separation along y axis; (i) $\psi = \pi/2$ (blue) and $\psi = -\pi/2$ (red) and separation along y axis. Two 4-dimethylamino-4'-nitrostilbenes ($C_{16}H_{16}N_2O_2$) with (j) $\psi = 0$ (blue) and $\psi = \pi$ (red) and separation along x axis; (k) $\psi = 0$ (blue) and $\psi = \pi$ (red) and separation along y axis; (l) $\psi = \pi/2$ (blue) and $\psi = -\pi/2$ (red) and separation along y axis. Numerical calculations with all molecular linear and nonlinear properties (full lines), neglecting all nonlinear properties (harmonic $\epsilon \neq 0$), neglecting only second hyperpolarizability (anharmonic $\beta = 0$) are shown. Here, $g_\gamma = \gamma\hbar\omega^0/\alpha^2$ and $g_\beta = \beta(\hbar\omega^0)^{1/2}/(\alpha(1-\epsilon))^{3/2}$. Corresponding vdW forces are also given in Fig. 6 in Appendix C.</figcaption></figure>

when $B_{mnq} = G_{mnqp} = 0$, we recover the linear discrete coupled dipole method obtaining the results for vdW energies of carbon nanostructures and finite nanoparticle chains reported in Refs. [37–39]. Note that the discrete dipole representation within linear response has been successfully used in first-principles vdW many-body schemes [40–43]. The generalized method described here provides means to account explicitly for linear and nonlinear properties in many-body vdW interactions.

### III. ISOTROPIC PARTICLES

The method is now applied to two identical isotropic nanoparticles placed at (0,0,0) and (R, 0, 0) [insert in Fig. 2(a)] whose linear polarizabilities are $\alpha_{\mu\nu}^1 = \alpha_{\mu\nu}^2 = \alpha\delta_{\mu\nu}$. For the sake of simplicity, we further take that $\gamma_{\mu\nu\lambda\rho}^1 = \gamma_{\mu\nu\lambda\rho}^2 = \frac{\gamma}{3}(\delta_{\mu\nu}\delta_{\lambda\rho} + \delta_{\mu\lambda}\delta_{\nu\rho} + \delta_{\mu\rho}\delta_{\nu\lambda})$. In this case $\beta^1 = \beta^2 = 0$ due to the presence of inversion symmetry in isotropic materials [33]. Results from the numerical calculations for distances larger than the size of the dipole $R_0 = \sqrt[3]{\frac{\alpha}{2\pi\epsilon_0}}$ are shown in Figs. 2(a)–2(c), where $g_\gamma = \frac{\gamma\hbar\omega^0}{\alpha^2}$ “measures” the relative nonlinear-to-linear response strength (for materials, $|g_\gamma| < 1$ [33]). The numerical calculations can also be put in the perspective of the obtained analytical expressions at the $R \to R_0^+$ and $R \gg R_0$ asymptotic limits

$$U_{\text{vdW}}(R) \approx \begin{cases} g_{\gamma} \frac{R_0^6/R^6}{1-R_0^6/R^6} U_L(R), & R \to R_0^+ \\ (1+\frac{15}{2}g_{\gamma}) U_L(R), & R \gg R_0 \end{cases}. \quad (7)$$

where $U_L(R) = -\frac{3}{16}\hbar\omega^0 \frac{R_0^6}{R^6}$ is the standard London result for two atoms [44] (details in Appendix C). In the case of $\gamma > 0$, there is a stronger attraction [Fig. 2(a)] for all R when compared to $U_L$. For $\gamma < 0$, however, the energy is repulsive for $R \to R_0^+$. On the other hand, the long-ranged $U_{vdW}$ is

attractive for $g_\gamma > -2/15$ and repulsive otherwise. In the first case, the energy acquires a shallow bound state [Fig. 2(b)], while in the second case, the interaction is repulsive in the entire distance regime [Fig. 2(c)]. Thus, even though $U_{vdW}$ at $R \gg R_0$ has a rescaled $1/R^6$ London-like behavior, the interplay between the sign and magnitude of the nonlinearity and the short-range energy behavior can have profound effects on the vdW energy.

### IV. ANISOTROPIC PARTICLES

We also consider anisotropic nanoparticles taken to have two distinct optical axes with $\alpha_{yy} = \alpha_{zz} = \alpha_{\parallel} = \alpha$ and $\alpha_{xx} = \alpha_{\perp} = \alpha(1-\epsilon)$ where $0 \le \epsilon \le 1$. Many anisotropic materials, such as noncentrosymmetric transition-metal dichalcogenides, polar metals, and Weyl semimetals [11–18], have broken inversion symmetry and they exhibit strong second- and third-order nonlinear response. Although much interest has been generated due to their potential applications for optical rectification and second-harmonic generation applications, for example, here we show that the effect of nonlinearity is “felt” at a fundamental level concerning dispersive interactions.

To illustrate these nonlinear effects in the vdW energy, we first focus on two identical nanoparticles with $\beta_{xxx} = \beta$ (all other $\beta$ components are zero) and $\gamma_{yyyy} = \gamma_{zzzz} = 3\gamma_{yyzz} = \gamma$ (all other $\gamma$ components are zero). In addition to $R_0$ and $g_\gamma$, we introduce $g_\beta = \frac{\beta(\hbar\omega^0)^{1/2}}{(\alpha(1-\epsilon))^{3/2}}$. The vdW dependence on the optical axis relative orientation is tracked by the angle $\psi$ [inserts in Figs. 2(d)–2(f), also see Fig. 4]. The results in Figs. 2(d)–2(f) show that the energy depends collectively on the anisotropy (parameter $\epsilon$), the optical nonlinearities ($g_\beta$ and $g_\gamma$ parameters), and the optical axis orientation (angle $\psi$). 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. 2(d)]. For particles separated along the y axis and $ψ = \frac{\pi}{2}$, $g_β$ and $g_γ$ lead to stronger vdW attraction, while particles on the z axis and $ψ = \frac{\pi}{2}$ experience reduced attraction when compared with $U_L$ [Fig. 2(f)]. 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 angle $ψ = \pi$ [Fig. 2(e)]. The complexity of the vdW interaction for two anisotropic nanoparticles can be assessed by considering the asymptotic $R \gg R_0$ limit (see Appendix C for detailed derivations and comparisons):

$$ U_{vdW}(\mathbf{R}) \approx \left\{ \begin{array}{ll} \left(1 - \frac{2\epsilon(2-\epsilon)}{3} + 2g_{\gamma} - \frac{653g_{\beta}^2(1-2\epsilon)}{162}\right) U_L(R) & \\ \mp \frac{g_{\beta}^2 \hbar \omega^0}{16 R_0^3}, & \mathbf{R} = \mathbf{Re}_x, \psi = 0(-), \pi(+) \\ \\ \left(1 - \frac{\epsilon(2-\epsilon)}{6} + 5g_{\gamma} - \frac{653g_{\beta}^2(1-2\epsilon)}{648}\right) U_L(R) & \\ \pm \frac{g_{\beta}^2 \hbar \omega^0}{32 R_0^3}, & \mathbf{R} = \mathbf{Re}_y, \psi = 0(+), \pi(-). \quad (8) \\ \\ \left(1 - \frac{5\epsilon}{6} + \frac{(7-5\epsilon)g_{\gamma}}{2} - \frac{5153g_{\beta}^2(1-\epsilon)}{2592}\right) U_L(R) & \\ \mathbf{R} = \mathbf{Re}_x, \psi = \frac{\pi}{2} \\ \\ \left(1 - \frac{\epsilon}{3} + (5-\epsilon)g_{\gamma} - \frac{653g_{\beta}^2(1-\epsilon)}{648}\right) U_L(R) & \\ \mathbf{R} = \mathbf{Re}_y, \psi = \frac{\pi}{2} \end{array} \right. $$

One finds that the London-like term $U_L(R)$ is rescaled by factors that depend nontrivially on $\epsilon$, $g_β$, $g_γ$, and $ψ$ for the various orientations. However, the vdW energy acquires a term entirely due to the second-order nonlinearity, which has a much longer outreach range due to 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$. Equation (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_2O_2$) as representatives for nonlinear dipolar nanoparticles [Figs. 2(g)-2(l)]. The molecular anisotropy is accompanied by hyperpolarizabilities with various tensor components as given in Refs. [45,46]. 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 possible for two azulenes providing they are on the x axis with an opposite direction of their optical axis [Fig. 2(g)].

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 based on the discrete dipole model at the linear response level [6,7], however within our approach one can also account for the effect of optical nonlinearities. Using the scheme from Fig. 1,

<figure><img src="image_6.png" /><figcaption>FIG. 3. VdW energy rescaled by $\hbar\omega^0$ vs nanoparticle separation rescaled by $R_0 = (\alpha/2\pi\epsilon_0)^{1/3}$ is shown in all panels. Two chains where each one has 10 isotropic molecules along x axis and separated by (a) $a = 1.5R_0$; (b) $a = 1.75R_0$ for $g_γ = 0$ (black), $g_γ = 0.1$ (red), and $g_γ = -0.1$ (blue). Two chains where each one has 10 azulenes along x axis and separated by (c) $a = 1.5R_0$; (d) $a = 1.75R_0$ with $ψ = 0$ (red) and $ψ = 0$ (blue). Results for $g_β = g_γ = 0$ (black) and $g_β = 0, g_γ \neq 0$ (green) are also shown.</figcaption></figure>

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, Figs. 3(a) and 3(b) show 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 Figs. 3(d) and 3(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. 2(g)] for which the $U_{vdW}$ was found to always be attractive.

## V. 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. The systems studied here helped us uncover that optical nonlinearities can rescale the London attraction and their signs can even cause vdW repulsion. The $1/R^3$ distance behavior due to second order hyperpolarizability indicates a much longer-ranged outreach when compared to the typical $1/R^6$ distance dependence. Nonlinearities are also expected to matter for vdW interactions in other nanoparticle systems, including alkali-halide clusters whose giant polarizabilities result in $R_0$ being 3–4 times bigger than the physical size of the cluster [47]. For such materials, considering the optical nonlinearities may be crucial for the overall stability of the interacting system.

The discrete coupled nonlinear dipole method can also accommodate the application of external fields, which can
result in significant modifications of dispersive interactions [48,49]. For example, magnetoelectric hyperpolarizabilities of individual nanoparticles $\alpha^i(\mathbf{B}_0)$, $\beta^i(\mathbf{B}_0)$, and $\gamma^i(\mathbf{B}_0)$ [50,51] can be used to study the effects of a static magnetic field $\mathbf{B}_0$. Adding a term $\hat{H}_0 = -\sum_{i=1}^N \mathbf{E}_0 \cdot \hat{\mathbf{p}}^i$ [52,53] to the Hamiltonian takes into account a static electric field $\mathbf{E}_0$, which ultimately affects the equilibrium dipole moments $P_0$, but not the overall computational scheme.

The impact of the discrete coupled nonlinear dipole model can be even broader when considering future developments of first-principles methods for vdW calculations. The discrete dipole model within linear response is currently employed in many-body dispersion energy calculations with several options available in state of the art codes. Such schemes have been applied successfully to calculate vdW energies in a variety of materials [6,7,41,42]. In light of this study, the vdW interaction and its effects on stacking patterns, interlayer distances, bond lengths, etc., however, must be re-examined

for a more accurate description. Our theoretical framework can serve as a stepping stone for future design of mainstream computational schemes for more accurate simulations, which might be necessary for nonlinear materials with chemically inert components.

ACKNOWLEDGMENTS

We acknowledge support from the U.S. Department of Energy under Grant No. DE-FG02-06ER46297. P.R.-L. acknowledges support from Ministerio de Ciencia e Innovación (Spain), Agencia Estatal de Investigación, under project NAUTILUS (PID2022-139524NB-I00), from AYUDA PUENTE, URJC, from QuantUM program of the University of Montpellier and the hospitality of the Theory of Light-Matter and Quantum Phenomena group at the Laboratoire Charles Coulomb, University of Montpellier, where part of this work was done.

APPENDIX A: HAMILTONIAN OF NONLINEAR DIPOLE SYSTEMS

1. Hamiltonian for a single nonlinear dipole

Our objective here is to derive the quantum mechanical Hamiltonian for a single nonlinear dipole responding to an electric field. We begin with the classical Hamiltonian $H$ written in the following form [31,32]:

$$H = \frac{1}{2} F_{\mu\nu} \pi_{\mu} \pi_{\nu} + \frac{1}{2} A_{\mu\nu} p_{\mu} p_{\nu} + \frac{1}{3} B_{\mu\nu\lambda} p_{\mu} p_{\nu} p_{\lambda} + \frac{1}{4} \Gamma_{\mu\nu\lambda\rho} p_{\mu} p_{\nu} p_{\lambda} p_{\rho} + \dots, \quad (A1)$$

where $p_\mu$ is the dipolar moment and $\pi_\mu$ is its corresponding canonical momentum. To determine the tensors with components $A_{\mu\nu}$, $B_{\mu\nu\lambda}$, and $\Gamma_{\mu\nu\lambda\rho}$, we examine the dipolar response due to an electric field $\mathbf{E}$

$$p_{\mu} = \alpha_{\mu\nu} E_{\nu} + \beta_{\mu\nu\lambda} E_{\nu} E_{\lambda} + \gamma_{\mu\nu\lambda\rho} E_{\nu} E_{\lambda} E_{\rho} + \dots, \quad (A2)$$

where $\alpha_{\mu\nu}$ are the components of the polarizability tensor of rank 2 from linear response, while $\beta_{\mu\nu\lambda}$, $\gamma_{\mu\nu\lambda\rho}$, ... are the components of the second-order, third-order, ... hyperpolarizabilities (rank 3, 4, ... tensors). Considering the static free energy $\mathcal{H} = H - p_\mu E_\mu$, one can write the Hamilton's equations of motion,

$$\left\{ \begin{aligned} \frac{d p_{\mu}}{dt} &= \frac{\partial \mathcal{H}}{\partial \pi_{\mu}} \\ \frac{d \pi_{\mu}}{dt} &= -\frac{\partial \mathcal{H}}{\partial p_{\mu}} \end{aligned} \right. \Rightarrow \left\{ \begin{aligned} \frac{d p_{\mu}}{dt} &= F_{\mu\nu} \pi_{\nu} \\ \frac{d \pi_{\mu}}{dt} &= -A_{\mu\nu} p_{\nu} - B_{\mu\nu\lambda} p_{\nu} p_{\lambda} - \Gamma_{\mu\nu\lambda\rho} p_{\nu} p_{\lambda} p_{\rho} - \dots + E_{\mu}. \end{aligned} \right. \quad (A3)$$

For a stationary dipole $\frac{d\pi_\mu}{dt} = 0$, thus we find that $E_\mu = A_{\mu\nu}p_\nu + B_{\mu\nu\lambda}p_\nu p_\lambda + \Gamma_{\mu\nu\lambda\rho}p_\nu p_\lambda p_\rho + \dots$ whose substitution in Eq. (A2) yields

$$\begin{aligned} p_{\mu} &= \alpha_{\mu\nu}(A_{\nu\lambda}p_{\lambda}) + \alpha_{\mu\nu}(B_{\nu\rho\sigma}p_{\rho}p_{\sigma}) + \beta_{\mu\nu\lambda}(A_{\nu\rho}p_{\rho})(A_{\lambda\tau}p_{\tau}) + \alpha_{\mu\nu}(\Gamma_{\nu\sigma\tau\xi}p_{\sigma}p_{\tau}p_{\xi}) \\ &\quad + \beta_{\mu\nu\lambda}(A_{\nu\sigma}p_{\sigma})(B_{\lambda\tau\xi}p_{\tau}p_{\xi}) + \beta_{\mu\nu\lambda}(B_{\nu\tau\xi}p_{\tau}p_{\xi})(A_{\lambda\sigma}p_{\sigma}) + \gamma_{\mu\nu\lambda\rho}(A_{\nu\sigma}p_{\sigma})(A_{\lambda\tau}p_{\tau})(A_{\rho\xi}p_{\xi}) + \dots \\ &= \alpha_{\mu\nu}A_{\nu\lambda}p_{\lambda} + [\alpha_{\mu\nu}B_{\nu\rho\tau} + \beta_{\mu\nu\lambda}A_{\nu\rho}A_{\lambda\tau}]p_{\rho}p_{\tau} + [\alpha_{\mu\nu}\Gamma_{\nu\sigma\tau\xi} + \beta_{\mu\nu\lambda}A_{\nu\sigma}B_{\lambda\tau\xi} + \beta_{\mu\nu\lambda}B_{\nu\tau\xi}A_{\lambda\sigma} \\ &\quad + \gamma_{\mu\nu\lambda\rho}A_{\nu\sigma}A_{\lambda\tau}A_{\rho\xi}]p_{\sigma}p_{\tau}p_{\xi} + \dots. \end{aligned} \quad (A4)$$

From here, we obtain

$$\Rightarrow \begin{cases} A_{\mu\nu} = (\alpha)^{-1}_{\mu\nu} \\ B_{\mu\nu\lambda} = -(\alpha)^{-1}_{\mu\sigma} \beta_{\sigma\tau\xi} (\alpha)^{-1}_{\tau\nu} (\alpha)^{-1}_{\xi\lambda} \\ \Gamma_{\mu\nu\lambda\rho} = -(\alpha)^{-1}_{\mu\sigma} \gamma_{\sigma\tau\xi\eta} (\alpha)^{-1}_{\tau\nu} (\alpha)^{-1}_{\xi\lambda} (\alpha)^{-1}_{\eta\rho} + (\alpha)^{-1}_{\mu\sigma_1} \beta_{\sigma_1\tau_1\xi_1} (\alpha)^{-1}_{\tau_1\nu} (\alpha)^{-1}_{\xi_1\sigma_2} \beta_{\sigma_2\tau_2\xi_2} (\alpha)^{-1}_{\tau_2\lambda} (\alpha)^{-1}_{\xi_2\rho} \\ \qquad + (\alpha)^{-1}_{\mu\sigma_1} \beta_{\sigma_1\tau_1\xi_1} (\alpha)^{-1}_{\tau_1\sigma_2} (\alpha)^{-1}_{\xi_1\nu} \beta_{\sigma_2\tau_2\xi_2} (\alpha)^{-1}_{\tau_2\lambda} (\alpha)^{-1}_{\xi_2\rho} \end{cases} \quad (A5)$$

Additionally, since $\frac{d\pi_\mu}{dt} = 0$, we find the following second-order differential equation from Eq. (A3) assuming now small fluctuations under no external fields

$$\frac{d^2 p_{\mu}}{dt^2} = -F_{\mu\nu}A_{\nu\lambda}p_{\lambda} + \mathcal{O}(p^2), \quad (A6)$$
where only terms linear in $p_μ$ are retained. Further, taking into consideration that the fluctuating dipole is essentially a harmonic oscillator with frequency $ω^0$, this differential equation is solved:

$$ \frac{d^2 p_μ}{dt^2} = -(ω^0)^2 p_μ \Rightarrow F_{μν}A_{νλ} = (ω^0)^2 δ_{μλ} \Rightarrow F_{μν} = (ω^0)^2(A)^{-1}_{μν} = (ω^0)^2 α_{μν}. \quad (A7) $$

Now, replacing the classical dipole and canonical momenta with their quantum mechanical operator equivalents $p \to \hat{p}$ and $\pi \to \hat{\pi}$ and taking into account that $[\hat{p}_μ, \hat{\pi}_ν] = i\hbar δ_{μν}$, one arrives at the quantized Hamiltonian for single dipolar nanoparticle

$$ \hat{H} = H(\hat{\pi}, \hat{p}) = \frac{1}{2}(ω^0)^2 α_{μν} \hat{\pi}_μ \hat{\pi}_ν + \frac{1}{2}α_{μν}^{-1} \hat{p}_μ \hat{p}_ν + \frac{1}{3}B_{μνλ} \hat{p}_μ \hat{p}_ν \hat{p}_λ + \frac{1}{4}Γ_{μνλρ} \hat{p}_μ \hat{p}_ν \hat{p}_λ \hat{p}_ρ + \dots. \quad (A8) $$

## 2. Hamiltonian for a system of N nonlinear dipoles

The dipole-dipole interaction between two nanoparticles (labeled as i, j) separated by $R^{ij} = R^i - R^j$ is $\hat{U}^{ij} = U(\hat{p}^i, \hat{p}^j, R^{ij}) = -T_{μν}(R^{ij})\hat{p}_μ^i \hat{p}_ν^j$, where $T_{μν}(R) = \frac{3R_μR_ν - R^2δ_{μν}}{4πε_0 R^5}$ is the usual displacement tensor [54]. Thus, the Hamiltonian for a system of N dipolar nanoparticles includes the Hamiltonian for each nanoparticle from Eq. (A8) and their dipole-dipole coupling

$$ \hat{H} = \sum_{i=1}^{N} \hat{H}^{i} + \sum_{i=1}^{N} \sum_{j=i+1}^{N} \hat{U}^{ij} = \frac{1}{2} \mathbb{F}_{mn} \hat{\Pi}_{m} \hat{\Pi}_{n} + \frac{1}{2} \mathbb{A}_{mn} \hat{P}_{m} \hat{P}_{n} + \frac{1}{3} \mathbb{B}_{mnq} \hat{P}_{m} \hat{P}_{n} \hat{P}_{q} + \frac{1}{4} \mathbb{G}_{mnqp} \hat{P}_{m} \hat{P}_{n} \hat{P}_{q} \hat{P}_{p}, $$

$$ \mathbb{F}_{mn} = (ω^{0,j})^2 α_{μν}^j δ_{ij}, \quad \text{where } m = 3i + μ, n = 3j + ν, \quad (A9) $$

$$ \mathbb{A}_{mn} = (α^j)^{-1}_{μν} δ_{ij} - T^{ij}_{μν} (1 - δ_{ij}), \quad \text{where } m = 3i + μ, n = 3j + ν, \quad (A10) $$

$$ \mathbb{B}_{mnq} = B^j_{μνλ} δ_{ijk} \quad \text{where } m = 3i + μ, n = 3j + ν, q = 3k + λ, \quad (A11) $$

$$ \mathbb{G}_{mnqp} = \Gamma^j_{μνλρ} δ_{ijkl} \quad \text{where } m = 3i + μ, n = 3j + ν, q = 3k + λ, p = 3l + ρ. \quad (A12) $$

Note that in the above $\hat{P} = [\hat{p}^1, \hat{p}^2, \dots, \hat{p}^N]^T = [\hat{P}_1, \hat{P}_2, \dots, \hat{P}_{3N}]^T$ and $\hat{\Pi} = [\hat{\pi}^1, \hat{\pi}^2, \dots, \hat{\pi}^N]^T = [\hat{\Pi}_1, \hat{\Pi}_2, \dots, \hat{\Pi}_{3N}]^T$. Einstein notation for repeating indices is also implied in the final expression for $\hat{H}$.

Using a tensor product notation for $B_{μνλ}$ and $Γ_{μνλρ}$ in Eq. (A5), we find

$$ \begin{aligned} B &= -\alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1} \\ \Gamma &= 2\alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1} - \alpha^{-1} \otimes \gamma \otimes \alpha^{-1} \otimes \alpha^{-1} \otimes \alpha^{-1}. \end{aligned} \quad (A13) $$

Taking advantage of a generalized $δ$-function defined as $δ_{abcd...} = \begin{cases} 1 & \text{when } a=b=c=d=\dots, \\ 0 & \text{otherwise} \end{cases}$, Eqs. (A9)–(A12) can be represented as

$$ \hat{H} = \frac{1}{2} \mathbb{F}_{mn} \hat{\Pi}_m \hat{\Pi}_n + \frac{1}{2} \mathbb{A}_{mn} \hat{P}_m \hat{P}_n + \frac{1}{3} \mathbb{B}_{mnq} \hat{P}_m \hat{P}_n \hat{P}_q + \frac{1}{4} \mathbb{G}_{mnqp} \hat{P}_m \hat{P}_n \hat{P}_q \hat{P}_p, \quad (A14) $$

$$ \mathbb{F}_{mn} = [(ω^0)^2 α] δ_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor}, \quad \mathbb{A}_{mn} = [α^{-1}] δ_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor} - \mathbb{T}^{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor} (1 - δ_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor}), \quad (A15) $$

$$ \mathbb{B}_{mnq} = -[\alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1}] δ_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor \lfloor \frac{q+2}{3} \rfloor}, \quad (A16) $$

$$ \mathbb{G}_{mnqp} = [2\alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1} \otimes \beta \otimes \alpha^{-1} \otimes \alpha^{-1} - \alpha^{-1} \otimes \gamma \otimes \alpha^{-1} \otimes \alpha^{-1} \otimes \alpha^{-1}] δ_{\lfloor \frac{m+2}{3} \rfloor \lfloor \frac{n+2}{3} \rfloor \lfloor \frac{q+2}{3} \rfloor \lfloor \frac{p+2}{3} \rfloor}, \quad (A17) $$

which is Eq. (1) in the main text.

For dipoles fluctuating around an equilibrium $P_0$, one can further express $\hat{P}_m = P_{0,m} + \hat{Q}_m$ with $\hat{Q}_m$ being the fluctuation operator. The $P_0 = [P_{0,1}, P_{0,2}, \dots, P_{0,3N}]^T$ is found from

$$ \frac{\partial \hat{H}}{\partial \hat{P}_m}(P_{0,m}) = 0 \Rightarrow \mathbb{A}_{mn} P_{0,n} + \mathbb{B}_{mnq} P_{0,n} P_{0,q} + \mathbb{G}_{mnqp} P_{0,n} P_{0,q} P_{0,p} = 0. \quad (A18) $$

The above equation shows that in general $P_{0,m} \neq 0$. Such an effect is a consequence of the $\mathbb{B}$ and $\mathbb{G}$ terms and it does not occur for linear dipoles. At equilibrium, the free energy of the system becomes $E_0 = \hat{H}(P_{0,m})$, which also contains information about the nonlinear dipolar properties. Consequently, the Hamiltonian for the fluctuating dipoles transforms into

$$ \hat{H} = E_0 + \hat{H}_h + \hat{H}_{anh}, \quad (A19) $$

$$ E_0 = \frac{1}{2} \mathbb{A}_{mn} P_{0,m} P_{0,n} + \frac{1}{3} \mathbb{B}_{mnq} P_{0,m} P_{0,q} + \frac{1}{4} \mathbb{G}_{mnqp} P_{0,m} P_{0,n} P_{0,q} P_{0,p}, \quad (A20) $$

$$ \hat{H}_h = \frac{1}{2} \mathbb{F}_{mn} \hat{\Pi}_m \hat{\Pi}_n + \frac{1}{2} [\mathbb{A}_{mn} + 2 \mathbb{B}_{mnq} P_{0,q} + 3 \mathbb{G}_{mnqp} P_{0,q} P_{0,p}] \hat{Q}_m \hat{Q}_n, \quad (A21) $$

$$ \hat{H}_{anh} = \frac{1}{3} [\mathbb{B}_{mnq} + 3 \mathbb{G}_{mnqp} P_{0,p}] \hat{Q}_m \hat{Q}_n \hat{Q}_q + \frac{1}{4} \mathbb{G}_{mnqp} \hat{Q}_m \hat{Q}_n \hat{Q}_q \hat{Q}_p. \quad (A22) $$
APPENDIX B: EIGENMODES AND VDW ENERGY OF NONLINEAR DIPOLES

The Hamiltonian in Eq. (A19) sets the stage for a many-body nonlinear vdW energy dispersion based on zero-point eigenmode summation supplemented by a perturbation theory for $\hat{H}_{\text{anh}}$, which was conveniently separated from the linear effects in $\hat{H}_h$ assumed to be dominant.

We first consider $\hat{H}_h$ by applying an orthogonal transformation $U$, such that

$$ \hat{H}_h = \frac{1}{2} \hat{\Pi}_m \hat{\Pi}_m + \frac{1}{2} \tilde{A}_{mn} \hat{Q}_m \hat{Q}_n, \quad (B1) $$

$$ \hat{Q} = M\hat{Q} = F^{1/2}U\hat{Q}, \quad \hat{\Pi} = F^{-1/2}U\hat{\Pi} = (M^{-1})^T\hat{\Pi}, \quad (B2) $$

$$ \begin{aligned} \tilde{A}_{mn} &= (U^T F^{1/2})_{mk} [A_{kl} + 2B_{klq} P_{0,q} + 3G_{klqp} P_{0,q} P_{0,p}] (F^{1/2} U)_{ln} \\ &= (F^{-1/2} U)^{-1}_{mk} \{ [A_{kl} + 2B_{klq} P_{0,q} + 3G_{klqp} P_{0,q} P_{0,p}] F_{lj} \} (F^{-1/2} U)_{jn}. \end{aligned} \quad (B3) $$

The transformation for which $\tilde{A}$ is diagonal $\tilde{A} = \begin{bmatrix} \omega_1^2 & 0 & \cdots & 0 \\ 0 & \omega_2^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & \cdots & \omega_{3N}^2 \end{bmatrix}$ is associated with the secular equation

$$ [A_{mn} + 2B_{mnq}P_{0,q} + 3G_{mnqp}P_{0,q}P_{0,p}]F_{nl}u_l = \omega^2 u_l, \quad (B4) $$

where $\omega^2$ are the eigenvalues and $u_l$ are the corresponding eigenvectors. The Hamiltonian in Eq. (B1) can further be written in terms of creation and annihilation operators $\hat{a}_n^\dagger$, $\hat{a}_n$, such that $\hat{Q}_n = \sqrt{\frac{\hbar}{2\omega_n}}(\hat{a}_n^\dagger + \hat{a}_n)$, $\hat{\Pi}_n = \sqrt{\frac{\hbar\omega_n}{2}}(\hat{a}_n^\dagger - \hat{a}_n)$,

$$ \hat{H}_h = \frac{1}{2} \hat{\Pi}_n \hat{\Pi}_n + \frac{1}{2} \omega_n^2 \hat{Q}_n \hat{Q}_n = \hbar \omega_n (\hat{a}_n^\dagger \hat{a}_n + \frac{1}{2}). \quad (B5) $$

Thus the ground-state energy is composed of the zero-point modes $E_{g.s}^{(0)} = \langle 000 \dots 0 | \hat{H}_h | 000 \dots 0 \rangle = \sum_{n=1}^{3N} \frac{1}{2} \hbar \omega_n$. The diagonalization procedure we have described above constitutes the discrete coupled dipole method utilized in many-body dispersion schemes for ab initio simulations [40–43]. Here, however, the dipolar nonlinearities are also taken into account in the eigenmodes and eigenstates of the secular equation Eq. (B4).

The contributions from $\hat{Q}_n \hat{Q}_m \hat{Q}_q$ terms in Eq. (A22) and denoted as $\hat{H}_{\text{anh}}^{(2)}$ are now treated within perturbation theory. Using the creation and annihilation operators, we write

$$ \begin{aligned} \hat{H}_{\text{anh}}^{(2)} &= \frac{1}{3} \sqrt{\frac{\hbar}{2\omega_n}} \frac{\hbar}{2\omega_m} \frac{\hbar}{2\omega_q} \tilde{B}_{nmq} (\hat{a}_n^\dagger + \hat{a}_n)(\hat{a}_m^\dagger + \hat{a}_m)(\hat{a}_q^\dagger + \hat{a}_q), \\ \tilde{B}_{nmq} &= [\tilde{B}_{mnq} + 3G_{mnqp}P_{0,q}] M_{n'n} M_{m'm} M_{q'q}. \end{aligned} \quad (B6) $$

Given that the first-order correction to the ground-state energy is zero, the second-order perturbation yields

$$ \begin{aligned} \Delta E_{2nd}^{(2)} &= \sum_{(k_1, k_2, \dots, k_{3N}) \neq (0, 0, \dots, 0)} \frac{|\langle k_1, k_2, \dots, k_{3N} | \hat{H}_{2nd} | 000 \dots 0 \rangle|^2}{E_{g.s}^{(0)} - E_{k_1, k_2, \dots, k_{3N}}^{(0)}} \\ &= -\frac{1}{36} \sum_{n=1}^{3N} \frac{\hbar^2}{\omega_n^4} \left[ (\tilde{B}_{nnn})^2 + \frac{9}{2} \left( \tilde{B}_{nnn} + \sum_{l \neq n} \sqrt{\frac{\omega_n}{\omega_l}} \tilde{B}_{nnl} \right)^2 + \sum_{l \neq n}^{3N} \frac{9\omega_n^2}{(2\omega_n + \omega_l)\omega_l} (\tilde{B}_{nnl})^2 \right] \\ &\quad -\frac{1}{72} \sum_{\substack{n,m,q=1 \\ \text{pairwise different}}}^{3N} \frac{\hbar^2 (\tilde{B}_{nmq})^2}{\omega_n \omega_m \omega_q (\omega_n + \omega_m + \omega_q)}, \end{aligned} \quad (B7) $$

where the excited state is $|k_1, k_2, \dots, k_{3N}\rangle = \frac{1}{\sqrt{k_1! k_2! \dots k_{3N}!}} (a_1^\dagger)^{k_1} (a_2^\dagger)^{k_2} \dots (a_{3N}^\dagger)^{k_{3N}} |000 \dots 0\rangle$. Calculations for matrix elements in creation-annihilation can be found in textbooks such as Ref. [55].

The contributions from $\hat{Q}_n \hat{Q}_m \hat{Q}_q \hat{Q}_p$ terms in Eq. (A22) and denoted as $\hat{H}_{\text{anh}}^{(3)}$ are also written in terms of creation and annihilation operators

$$ \begin{aligned} \hat{H}_{\text{anh}}^{(3)} &= \frac{1}{4} \sqrt{\frac{\hbar}{2\omega_n}} \frac{\hbar}{2\omega_m} \frac{\hbar}{2\omega_q} \frac{\hbar}{2\omega_p} \tilde{G}_{nmqp} (\hat{a}_n^\dagger + \hat{a}_n)(\hat{a}_m^\dagger + \hat{a}_m)(\hat{a}_q^\dagger + \hat{a}_q)(\hat{a}_p^\dagger + \hat{a}_p), \\ \tilde{G}_{nmqp} &= G_{n'n'q'p'} M_{n'n} M_{m'm} M_{q'q} M_{p'p}. \end{aligned} \quad (B8) $$
From first-order perturbation theory, we find

$$ \Delta E_{3rd}^{(1)} = \langle 000 \dots 0 | \hat{H}_{3rd} | 000 \dots 0 \rangle = \frac{3}{8} \sum_{n,m=1}^{3N} \frac{\hbar^2}{\omega_n \omega_m} \tilde{\mathcal{G}}_{nmnm} = \frac{3}{8} \sum_{n=1}^{3N} \frac{\hbar^2}{\omega_n^2} \left[ \tilde{\mathcal{G}}_{nnnn} + \sum_{l \neq n} \frac{\omega_n}{\omega_l} \tilde{\mathcal{G}}_{nnll} \right]. \quad (B9) $$

Thus, the total ground-state energy becomes

$$ \begin{aligned} E_{g,s} & \approx E_0 + E_{g,s}^{(0)} + \Delta E_{2nd}^{(2)} + \Delta E_{3rd}^{(1)} + \dots \\ & = \sum_{n=1}^{3N} \frac{1}{2} \hbar \omega_n - \frac{1}{36} \sum_{n=1}^{3N} \frac{\hbar^2}{\omega_n^4} \left[ (\tilde{\mathbb{B}}_{nnn})^2 + \frac{9}{2} \left( \tilde{\mathbb{B}}_{nnn} + \sum_{l \neq n} \sqrt{\frac{\omega_n}{\omega_l}} \tilde{\mathbb{B}}_{nnl} \right)^2 + \sum_{l \neq n} \frac{9 \omega_n^2}{(2\omega_n + \omega_l)\omega_l} (\tilde{\mathbb{B}}_{nnl})^2 \right] \\ & \quad - \frac{1}{72} \sum_{\substack{n,m,q=1 \\ \text{pairwise different}}}^{3N} \frac{\hbar^2 (\tilde{\mathbb{B}}_{nmq})^2}{\omega_n \omega_m \omega_q (\omega_n + \omega_m + \omega_q)} + \frac{3}{8} \sum_{n=1}^{3N} \frac{\hbar^2}{\omega_n^2} \left[ \tilde{\mathcal{G}}_{nnnn} + \sum_{l \neq n} \frac{\omega_n}{\omega_l} \tilde{\mathcal{G}}_{nnll} \right] + \dots \end{aligned} \quad (B10) $$

## APPENDIX C: VDW ENERGY AND FORCE OF IDENTICAL NANOPARTICLES

### 1. vdW energy of two identical isotropic nanoparticles with a third-order nonlinearity

Here we give the detailed application of the discrete coupled nonlinear dipole method to two identical isotropic nanoparticles placed at (0,0,0) and (R, 0, 0). The response properties are taken as $\alpha_{\mu\nu}^1 = \alpha_{\mu\nu}^2 = \alpha\delta_{\mu\nu}$ and $\gamma_{\mu\nu\lambda\rho}^1 = \gamma_{\mu\nu\lambda\rho}^2 = \frac{\gamma}{3}(\delta_{\mu\nu}\delta_{\lambda\rho} + \delta_{\mu\lambda}\delta_{\nu\rho} + \delta_{\mu\rho}\delta_{\nu\lambda})$. One notes that $\beta = 0$ due to the presence of inversion symmetry in isotropic materials. By defining $R_0 = \sqrt[3]{\frac{\alpha}{2\pi\epsilon_0}}$ and $g_\gamma = \frac{\gamma\hbar\omega^0}{\alpha^2}$, we find

$$ \mathbb{F}_{nm} = (2\pi\epsilon_0 R_0^3)(\omega^0)^2 \delta_{nm}, \quad (C1) $$

$$ \tilde{\mathcal{G}}_{nmqp} = \begin{cases} -\frac{g_\gamma}{3(\hbar\omega)(4\pi^2\epsilon_0^2 R_0^6)} (\delta_{\mu\nu}\delta_{\lambda\rho} + \delta_{\mu\lambda}\delta_{\nu\rho} + \delta_{\mu\rho}\delta_{\nu\lambda}) & \text{when } (n, m, q, p) = 3k + (\mu, \nu, \lambda, \rho) \\ 0 & \text{otherwise} \end{cases}. \quad (C2) $$

The equilibrium conditions from Eqs. (A18), (A20) are obtained as

$$ P_{0,1} = P_{0,4} = \sqrt{\frac{3(2\pi\epsilon R_0^3)(\hbar\omega)}{-g_\gamma}} s(R), \quad P_{0,2} = P_{0,3} = P_{0,5} = P_{0,6} = 0, \quad (C3) $$

$$ E_0 = -\frac{\hbar\omega}{4(-g_\gamma)}s^2(R), \quad (C4) $$

$$ s(R) = \begin{cases} \frac{R_0^3}{R^3} - 1 & \text{when } R < R_0 \text{ and } g_\gamma < 0 \\ 0 & \text{otherwise} \end{cases}. \quad (C5) $$

The explicit matrix representation of the discrete coupled nonlinear dipole method Eq. (B4) is found as

$$ \left[ \begin{matrix} 1 + 3s(R) & 0 & 0 & -\frac{R_0^3}{R^3} & 0 & 0 \\ 0 & 1 + s(R) & 0 & 0 & \frac{R_0^3}{2R^3} & 0 \\ 0 & 0 & 1 + s(R) & 0 & 0 & \frac{R_0^3}{2R^3} \\ -\frac{R_0^3}{R^3} & 0 & 0 & 1 + 3s(R) & 0 & 0 \\ 0 & \frac{R_0^3}{2R^3} & 0 & 0 & 1 + s(R) & 0 \\ 0 & 0 & \frac{R_0^3}{2R^3} & 0 & 0 & 1 + s(R) \end{matrix} \right] \left[ \begin{matrix} u_{1,i} \\ u_{2,i} \\ u_{3,i} \\ u_{4,i} \\ u_{5,i} \\ u_{6,i} \end{matrix} \right] = \left[ \frac{\omega_i^2}{(\omega^0)^2} \right] \left[ \begin{matrix} u_{1,i} \\ u_{2,i} \\ u_{3,i} \\ u_{4,i} \\ u_{5,i} \\ u_{6,i} \end{matrix} \right], \quad (C6) $$

whose eigenvalues are

$$ \begin{align} \omega_1 &= \omega^0 \sqrt{1 - \frac{R_0^3}{R^3} + 3s(R)}, & \omega_2 &= \omega^0 \sqrt{1 - \frac{R_0^3}{2R^3} + s(R)}, & \omega_3 &= \omega^0 \sqrt{1 - \frac{R_0^3}{2R^3} + s(R)}, \nonumber \\ \omega_4 &= \omega^0 \sqrt{1 + \frac{R_0^3}{2R^3} + s(R)}, & \omega_5 &= \omega^0 \sqrt{1 + \frac{R_0^3}{2R^3} + s(R)}, & \omega_6 &= \omega^0 \sqrt{1 + \frac{R_0^3}{R^3} + 3s(R)}. \tag{C7} \end{align} $$
From Eqs. (B7) and (B9), the vdW energy is obtained as

$$ 
\begin{aligned}
U_{\text{vdW}}(R) &= E_{\text{g,s}}(R) - E_{\text{g,s}}(+\infty) \\
&= \frac{1}{2}\hbar\omega^0 \left\{ (f_0(R) - 6) - g_{\gamma} \left(f_1(R) - \frac{15}{2}\right) + \frac{1}{2g_{\gamma}}s^2(R) + g_{\gamma}s(R)f_2(R) \right\},
\end{aligned} 
\quad (C8) $$

where $E_{g,s}$ is the ground-state energy in Eq. (B10) and we have defined the following dimensionless functions

$$ f_0(R) = \sqrt{1 - \frac{R_0^3}{R^3} + 3s(R)} + \sqrt{1 + \frac{R_0^3}{R^3} + 3s(R)} + 2\sqrt{1 - \frac{R_0^3}{2R^3} + s(R)} + 2\sqrt{1 + \frac{R_0^3}{2R^3} + s(R)}, \quad (C9) $$

$$ f_1(R) = \left( \frac{1}{\sqrt{1 - \frac{R_0^3}{2R^3} + s(R)}} + \frac{1}{\sqrt{1 + \frac{R_0^3}{2R^3} + s(R)}} \right)^2 + \frac{3}{8} \left( \frac{1}{\sqrt{1 - \frac{R_0^3}{R^3} + 3s(R)}} + \frac{1}{\sqrt{1 + \frac{R_0^3}{R^3} + 3s(R)}} \right)^2 \\ + \frac{1}{2} \left( \frac{1}{\sqrt{1 - \frac{R_0^3}{2R^3} + s(R)}} + \frac{1}{\sqrt{1 + \frac{R_0^3}{2R^3} + s(R)}} \right) \left( \frac{1}{\sqrt{1 - \frac{R_0^3}{R^3} + 3s(R)}} + \frac{1}{\sqrt{1 + \frac{R_0^3}{R^3} + 3s(R)}} \right), \quad (C10) $$

$$ f_2(R) = \frac{1}{8(1 - \frac{R_0^3}{R^3} + 3s(R))^2} \left[ 33 + \frac{4\sqrt{\frac{1 - \frac{R_0^3}{R^3} + 3s(R)}{1 - \frac{R_0^3}{2R^3} + s(R)}}\sqrt{\frac{1 - \frac{R_0^3}{R^3} + 3s(R)}{1 + \frac{R_0^3}{2R^3} + s(R)}}\sqrt{\frac{1 - \frac{R_0^3}{R^3} + 3s(R)}{1 + \frac{R_0^3}{R^3} + 3s(R)}}}{\sqrt{1 - \frac{R_0^3}{2R^3} + s(R)}} + \sqrt{\frac{1 + \frac{R_0^3}{2R^3} + s(R)}{1 - \frac{R_0^3}{R^3} + 3s(R)}} + \sqrt{\frac{1 + \frac{R_0^3}{R^3} + 3s(R)}{1 - \frac{R_0^3}{R^3} + 3s(R)}} \right. \\ \left. + \frac{6\left(1 + 4\sqrt{\frac{1 - \frac{R_0^3}{R^3} + s(R)}{1 - \frac{R_0^3}{2R^3} + 3s(R)}}\right)}{\left(1 + \frac{R_0^3}{2R^3} + s(R)\right)^{3/2}\left(1 + 2\sqrt{\frac{1 - \frac{R_0^3}{2R^3} + s(R)}{1 - \frac{R_0^3}{R^3} + 3s(R)}}\right)} + \frac{27\left(1 + 4\sqrt{\frac{1 + \frac{R_0^3}{R^3} + 3s(R)}{1 - \frac{R_0^3}{R^3} + 3s(R)}}\right)}{\left(1 + \frac{R_0^3}{2R^3} + s(R)\right)^{3/2}\left(1 + 2\sqrt{\frac{1 + \frac{R_0^3}{R^3} + 3s(R)}{1 - \frac{R_0^3}{R^3} + 3s(R)}}\right)} \right]. \quad (C11) $$

For large distances $R \gg R_0$, $s(R) = 0$, the vdW energy becomes

$$ U_{\text{vdW}}(R) \approx -\frac{3\hbar\omega R_0^6}{16 R^6} \left(1 + \frac{15}{2}g_{\gamma}\right) = \left(1 + \frac{15}{2}g_{\gamma}\right) U_L(R), \quad (C12) $$

where $U_L$ is London's result [44]. On the other hand, if $R \to R_0^+$, $s(R) = 0$ and the vdW energy becomes

$$ U_{\text{vdW}}(R) \approx -\frac{3g_{\gamma}\hbar\omega}{16} \frac{1}{1 - \frac{R_0^3}{R^3}}. \quad (C13) $$

## 2. vdW energy of two identical anisotropic nanoparticles

We also consider the application of the discrete coupled nonlinear dipole method to two identical anisotropic nanoparticles placed at (0,0,0) and $(R_x, R_y, R_z)$. The linear response is taken as $\alpha_{yy} = \alpha_{zz} = \alpha_{\parallel} = \alpha$ and $\alpha_{xx} = \alpha_{\perp} = \alpha(1 - \epsilon)$ where $0 \le \epsilon \le 1$ while their nonlinearity is captured by $\beta_{xxx} = \beta$ (all other $\beta$ components are zero) and $\gamma_{yyyy} = \gamma_{zzzz} = 3\gamma_{\gamma yzz} = \gamma$ (all other $\gamma$ components are zero). In addition to $R_0$ and $g_{\gamma}$, another dimensionless constant is also used: $g_{\beta} = \frac{\beta(\hbar\omega^0)^{1/2}}{(\alpha(1-\epsilon))^{3/2}}$. We examine three different geometrical cases as given in Fig. 4.

For the particle at $\mathbf{R}$ whose optical axis are rotated by an angle $\psi$ with respect to the particle in the origin, we find that

$$ \mathbf{T}^{12}(\mathbf{R}) = \frac{1}{4\pi\epsilon_0 R^5} \begin{bmatrix} 2R_x^2 - R_y^2 - R_z^2 & 3R_x R_y & 3R_x R_z \\ 3R_y R_x & 2R_y^2 - R_x^2 - R_z^2 & 3R_y R_z \\ 3R_z R_x & 3R_z R_y & 2R_z^2 - R_x^2 - R_y^2 \end{bmatrix}, \quad (C14) $$

and

$$ \mathbb{B}_{nmq} = -\frac{g_{\beta}}{\sqrt{(1 - \epsilon)^3 (\hbar\omega^0) (2\pi\epsilon_0 R_0^3)^3}} \begin{cases} \delta_{mnq1} & \text{when } 1 \le n, m, q \le 3 \\ 0 & \text{when } 4 \le n, m, q \le 6, \\ 0 & \text{otherwise} \end{cases}, \quad (C15) $$
<figure><img src="image_7.png" /><figcaption>FIG. 4. Schematics of two identical uniaxial anisotropic particles. The x axis of second nanoparticle is specified by the relative angle $\psi$, while the position of this molecule is defined by ($R_x, R_y, R_z$). The considered geometrical cases are (a) $\mathbf{R} = (R, 0, 0)$, (b) $\mathbf{R} = (0, R, 0)$, (c) $\mathbf{R} = (0, 0, R)$.</figcaption></figure>

$$G_{nmqp} = G_{nmqp}^{\beta} + G_{nmqp}^{\gamma}$$

$$G_{nmqp}^{\beta} = \frac{2g_{\beta}^2}{(\hbar\omega^0)big(2\pi\epsilon_0 R_0^3)^2(1-\epsilon)^2} \times \begin{cases} \delta_{nmqp1} & \text{when } 1 \le n, m, q, p \le 3 \\ 0 & \text{when } 4 \le n, m, q, p \le 6 \\ 0 & \text{otherwise} \end{cases}$$

$$G_{nmqp}^{\gamma} = -\frac{g_{\gamma}}{3(\hbar\omega^0)(2\pi\epsilon_0 R_0^3)^2} \sum_{a=2,3} \sum_{b=2,3} \sum_{c=2,3} \sum_{d=2,3} (\delta_{ab}\delta_{cd} + \delta_{ac}\delta_{bd} + \delta_{ad}\delta_{bc}) \times \begin{cases} \delta_{na}\delta_{mb}\delta_{qc}\delta_{pd} & \text{when } 1 \le n, m, q, p \le 3 \\ 0 & \text{when } 4 \le n, m, q, p \le 6. \\ 0 & \text{otherwise} \end{cases} \quad (C16)$$

where the rotation matrix is

$$\mathcal{R}_z(\psi) = \begin{bmatrix} \cos \psi & -\sin \psi & 0 \\ \sin \psi & \cos \psi & 0 \\ 0 & 0 & 1 \end{bmatrix}. \quad (C17)$$

<figure><img src="image_8.png" /><figcaption>FIG. 5. VdW energy rescaled by $\hbar\omega^0$ vs nanoparticle separation rescaled by $R_0 = (\alpha/2\pi\epsilon_0)^{1/3}$ is shown in all panels. Two anisotropic nanoparticles with $\epsilon = 0.3$, $g_{\beta} = 0.3$, $g_{\gamma} = 0$, $\pm 0.1$ with relative angle between their optical axes (a),(c) $\psi = 0$, (b),(d) $\psi = \pi$, and (e),(f) $\psi = \pi/2$ and separation along (a),(b),(e) x axis, (c),(d) y axis or (f) z axis. Discrete data corresponds to numerical calculations, while the solid lines are from the asymptotic expressions in Eq. (8). Note that $\beta_{xxx} = \beta$ and $\gamma_{yyyy} = \gamma_{zzzz} = 3\gamma_{yyzz} = \gamma$, and all other hyperpolarizability components are zero.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 6. VdW force rescaled by $\hbar\omega^0 R_0^{-1}$ vs nanoparticle separation rescaled by $R_0 = (\alpha/2\pi\epsilon_0)^{1/3}$ is shown in all panels. Two isotropic nanoparticles with (a) $g_\gamma = 0.1$; (b) $g_\gamma = -0.1$; (c) $g_\gamma = -0.15$ with numerical (num), Eq. (C12) (approx.), and London formula $U_L = -(3/16)\hbar\omega^0(R_0/R)^6$ results. Two anisotropic nanoparticles with $\epsilon = 0.3$, $g_\beta = 0.3$, and $g_\gamma = 0, \pm 0.1$ with a relative angle between their optical axis (d) $\psi = 0$ (solid), $\psi = \pi$ (dashed) and separation along x axis; (e) $\psi = 0$ (solid), $\psi = \pi$ (dashed) and separation along y axis; (f) $\psi = \pi/2$ with separation along x axis (solid) and $\psi = \pi/2$ with separation along z axis (dashed). Two azulenes ($C_{10}H_8$) with (g) $\psi = 0$ (blue) and $\psi = \pi$ (red) and separation along x axis; (h) $\psi = 0$ (blue) and $\psi = \pi$ (red) and separation along y axis; (i) $\psi = \pi/2$ (blue) and $\psi = -\pi/2$ (red) and separation along y axis. Two 4-dimethylamino-4'-nitrostilbenes ($C_{16}H_{16}N_2O_2$) with (j) $\psi = 0$ (blue) and $\psi = \pi$ (red) and separation along x axis; (k) $\psi = 0$ (blue) and $\psi = \pi$ (red) and separation along y axis; (l) $\psi = \pi/2$ (blue) and $\psi = -\pi/2$ (red) and separation along y axis. Numerical calculations with all molecular linear and nonlinear properties (full lines), neglecting all nonlinear properties (harmonic $\epsilon \neq 0$), neglecting only second hyperpolarizability (anharmonic $\beta = 0$) are shown. Here, $g_\gamma = \gamma\hbar\omega^0/\alpha^2$ and $g_\beta = \beta(\hbar\omega^0)^{1/2}/(\alpha(1-\epsilon))^{3/2}$.</figcaption></figure>

The diagonalization procedure of $\hat{H}_h$ in Eq. (B4) followed by Taylor expansions up to the second order of $\epsilon$ yields $\omega_1, \omega_2, \omega_3, \omega_4, \omega_5, \omega_6$. For the considered positions given in Fig. 4, we obtain

$$ \mathbf{R} = (R, 0, 0) = R\mathbf{e}_x \\ \omega_{1,4} = \omega^0 \sqrt{1 \pm \frac{1}{2} \frac{R_0^3}{R^3}} \qquad (C18) $$

$$ \omega_{2,5} = \omega^0 \sqrt{1 \pm \frac{1}{2} \frac{R_0^3}{R^3} \left(1 - \frac{\epsilon(4+\epsilon)}{8} \sin^2 \psi\right)} \qquad (C19) $$

$$ \omega_{3,6} = \omega^0 \sqrt{1 \mp \frac{R_0^3}{R^3} \left(1 - \epsilon + \frac{\epsilon(4-\epsilon)}{8} \sin^2 \psi\right)}, \qquad (C20) $$

$$ \mathbf{R} = (0, R, 0) = R\mathbf{e}_y \\ \omega_{1,4} = \omega^0 \sqrt{1 \pm \frac{1}{2} \frac{R_0^3}{R^3}} \qquad (C21) $$

$$ \omega_{2,5} = \omega^0 \sqrt{1 \pm \frac{1}{2} \frac{R_0^3}{R^3} \left(1 - \epsilon + \frac{\epsilon(4-\epsilon)}{8} \sin^2 \psi\right)} \qquad (C22) $$

$$ \omega_{3,6} = \omega^0 \sqrt{1 \mp \frac{R_0^3}{R^3} \left(1 - \frac{\epsilon(4+\epsilon)}{8} \sin^2 \psi\right)}, \qquad (C23) $$

$$ \mathbf{R} = (0, 0, R) = R\mathbf{e}_z \\ \omega_{1,4} = \omega^0 \sqrt{1 \pm \frac{1}{2} \frac{R_0^3}{R^3} \left(1 - \epsilon + \frac{\epsilon(4-\epsilon)}{8} \sin^2 \psi\right)} \qquad (C24) $$

$$ \omega_{2,5} = \omega^0 \sqrt{1 \pm \frac{1}{2} \frac{R_0^3}{R^3} \left(1 - \frac{\epsilon(4+\epsilon)}{8} \sin^2 \psi\right)} \qquad (C25) $$

$$ \omega_{3,6} = \omega^0 \sqrt{1 \mp \frac{R_0^3}{R^3}}. \qquad (C26) $$
Also, using Eqs. (B7), (B9), we arrive at

$$U_{\text{vdW}}(\mathbf{R}, \psi) = E_{g,s}(\mathbf{R}, \psi) - E_{g,s}(+\infty, \psi), \qquad (C27)$$

where the ground-state energy Eq. (B10) is given as

$$E_{g,s}(\mathbf{R}, \psi) \approx \frac{\hbar\omega^0}{2} \{f_0(\mathbf{R}, \psi) - g_\gamma f_1(\mathbf{R}, \psi) + g_\beta^2 f_2(\mathbf{R}, \psi)\}. \qquad (C28)$$

The explicit forms of $f_0(\mathbf{R}, \psi)$, $f_1(\mathbf{R}, \psi)$, and $f_2(\mathbf{R}, \psi)$ can be obtained numerically, but are not given here due their extremely long expressions. For the asymptotic limit $R \gg R_0$ with small $\epsilon$, $g_\gamma$, $g_\beta \ll 1$, we obtain the expressions in Eq. (8).

Results for the interaction energies of two anisotropic particles with relative separations along the x, y, and z-axis (Fig. 4) are given in Fig. 5. The shown computations for the chosen parameters depict numerical and asymptotic calculations for a complementary comparison with the case of isotropic particles in Figs. 2(a)–2(f).

### 3. Nonlinear effects in the vdW force

In order to analyze the attractive and repulsive features of vdW interaction, we also calculate vdW force from obtained vdW energy for our examples given in main text. For fixed orientation of nanoparticles/molecules, the radial vdW force is defined as the radial gradient of vdW potential energy

$$F_{\text{vdW}}(R) = -\frac{d}{dR}U_{\text{vdW}}(R). \qquad (C29)$$

In Fig. 6 results for the calculated force for all cases in Fig. 2 in the main text are given.

[1] M. N. Gjerding, A. Taghizadeh, A. Rasmussen, S. Ali, F. Bertoldo, T. Deilmann, N. R. Knøsgaard, M. Kruse, A. H. Larsen, S. Manti, T. G. Pedersen, U. Petralanda, T. Skovhus, M. K. Svendsen, J. J. Mortensen, T. Olsen, and K. S. Thygesen, Recent progress of the computational 2D materials database (C2DB), 2D Mater. 8, 044002 (2021).

[2] D. Geng and H. Y. Yang, Recent advances in growth of novel 2D materials: Beyond graphene and transition metal dichalcogenides, Adv. Mater. 30, 1800865 (2018).

[3] S. K. Pandey, R. Das, and P. Mahadevan, Layer-dependent electronic structure changes in transition metal dichalcogenides: The microscopic origin, ACS Omega 5, 15169 (2020).

[4] R. K. Barik and L. M. Woods, High throughput calculations for a dataset of bilayer materials, Sci. Data 10, 232 (2023).

[5] Z. G. Yu, B. I. Yakobson, and Y.-W. Zhang, Realizing indirect-to-direct band gap transition in few-layer two-dimensional MX$_2$ (M = Mo, W; X = S, Se), ACS Appl. Energy Mater. 1, 4115 (2018).

[6] A. Ambrosetti, Jr., N. Ferri, R. A. Distasio, and A. Tkatchenko, Wavelike charge density fluctuations and van der Waals interactions at the nanoscale, Science 351, 1171 (2016).

[7] R. J. Maurer, C. Freysoldt, A. M. Reilly, J. G. Brandenburg, O. T. Hofmann, T. Björkman, S. Lebègue, and A. Tkatchenko, Advances in density-functional calculations for materials modeling, Annu. Rev. Mater. Res. 49, 1 (2019).

[8] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6, 15 (1996).

[9] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).

[10] P. Giannozzi Jr, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo et al., Advanced capabilities for mate-

rials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter 29, 465901 (2017).

[11] Y. Li, Y. Rao, K. F. Mak, Y. You, S. Wang, C. R. Dean, and T. F. Heinz, Probing symmetry properties of few-layer MoS$_2$ and h-BN by optical second-harmonic generation, Nano Lett. 13, 3329 (2013).

[12] C.-Y. Wang and G.-Y. Guo, Nonlinear optical properties of transition-metal dichalcogenide MX$_2$ (M = Mo, W; X = S, Se) monolayers and trilayers from first-principles calculations, J. Phys. Chem. C 119, 13268 (2015).

[13] A. Autere, H. Jussila, A. Marini, J. R. M. Saavedra, Y. Dai, A. Säynätjoki, L. Karvonen, H. Yang, B. Amirsolaimani, R. A. Norwood, N. Peyghambarian, H. Lipsanen, K. Kieu, F. J. G. de Abajo, and Z. Sun, Optical harmonic generation in monolayer group-VI transition metal dichalcogenides, Phys. Rev. B 98, 115426 (2018).

[14] X. Wen, Z. Gong, and D. Li, Nonlinear optics of two-dimensional transition metal dichalcogenides, InfoMat 1, 317 (2019).

[15] S. Patankar, L. Wu, B. Lu, M. Rai, J. D. Tran, T. Morimoto, D. E. Parker, A. G. Grushin, N. L. Nair, J. G. Analytis, J. E. Moore, J. Orenstein, and D. H. Torchinsky, Resonance-enhanced optical nonlinearity in the Weyl semimetal TaAs, Phys. Rev. B 98, 165113 (2018).

[16] K. Takasan, T. Morimoto, J. Orenstein, and J. E. Moore, Current-induced second harmonic generation in inversion-symmetric Dirac and Weyl semimetals, Phys. Rev. B 104, L161202 (2021).

[17] E. Drueke, J. Yang, and L. Zhao, Observation of strong and anisotropic nonlinear optical effects through polarization-resolved optical spectroscopy in the type-II Weyl semimetal $T_d-WTe_2$, Phys. Rev. B 104, 064304 (2021).

[18] Q. Xu, Y. Zhang, K. Koepernik, W. Shi, J. van den Brink, C. Felser, and Y. Sun, Comprehensive scan for nonmagnetic Weyl semimetals with nonlinear optical response, npj Comput. Mater. 6, 32 (2020).
[19] J. You, S. Bongu, Q. Bao, and N. Panoiu, Nonlinear optical properties and applications of 2D materials: Theoretical and experimental aspects, *Nanophotonics* **8**, 63 (2019).

[20] L. Zhou, H. Fu, T. Lv, C. Wang, H. Gao, D. Li, L. Deng, and W. Xiong, Nonlinear optical characterization of 2D materials, *Nanomaterials* **10**, 2263 (2020).

[21] H. Soo and M. Krüger, Fluctuational electrodynamics for non-linear media, *Europhys. Lett.* **115**, 41002 (2016).

[22] H. Soo and M. Krüger, Fluctuational electrodynamics for non-linear materials in and out of thermal equilibrium, *Phys. Rev. B* **97**, 045412 (2018).

[23] H. Soo, D. S. Dean, and M. Krüger, Particles with nonlinear electric response: Suppressing van der Waals forces by an external field, *Phys. Rev. E* **95**, 012151 (2017).

[24] D. Massa, A. Ambrosetti, and P. L. Silvestrelli, Many-body van der Waals interactions beyond the dipole approximation, *J. Chem. Phys.* **154**, 224115 (2021).

[25] T. Stedman, D. Drosdoff, and L. M. Woods, van der Waals interactions between nanostructures: Some analytic results from series expansions, *Phys. Rev. A* **89**, 012509 (2014).

[26] J. C. Hopkins, D. M. Dryden, W.-Y. Ching, R. H. French, V. A. Parsegian, and R. Podgornik, Dielectric response variation and the strength of van der Waals interactions, *J. Colloid Interface Sci.* **417**, 278 (2014).

[27] M. W. Cole, D. Velegol, H.-Y. Kim, and A. A. Lucas, Nanoscale van der Waals interactions, *Mol. Simul.* **35**, 849 (2009).

[28] V. V. Gobre and A. Tkatchenko, Scaling laws for van der Waals interactions in nanostructured materials, *Nat. Commun.* **4**, 2341 (2013).

[29] J. Sarabadani, A. Naji, R. Asgari, and R. Podgornik, Many-body effects in the van der Waals–Casimir interaction between graphene layers, *Phys. Rev. B* **84**, 155407 (2011).

[30] D. Drosdoff and L. M. Woods, Quantum and thermal dispersion forces: Application to graphene nanoribbons, *Phys. Rev. Lett.* **112**, 025501 (2014).

[31] V. Z. Lozovskii and B. I. Khudik, The new mechanism of physical adsorption on solid surface. I. Adsorption of nonpolar molecules, *Phys. Status Solidi* (b) **158**, 511 (1990).

[32] V. Z. Lozovskii and B. I. Khudik, The new mechanism of physical adsorption on solid surface. II. Adsorption of polar molecules, *Physica Status Solidi* (b) **160**, 137 (1990).

[33] R. W. Boyd, *Nonlinear Optics* (Academic Press, Rochester, NY, 2007).

[34] N. K. Balla, P. T. C. So, and C. J. R. Sheppard, Second harmonic scattering from small particles using Discrete Dipole Approximation, *Opt. Express* **18**, 21603 (2010).

[35] N. K. Balla, E. Y. S. Yew, C. J. R. Sheppard, and P. T. C. So, Coupled and uncoupled dipole models of nonlinear scattering, *Opt. Express* **20**, 25834 (2012).

[36] $E_0 = \frac{1}{2} A_{mn} P_{0,m} P_{0,n} + \frac{1}{3} B_{mnq} P_{0,m} P_{0,n} P_{0,q} + \frac{1}{4} G_{mnqp} P_{0,m} P_{0,n} P_{0,q} P_{0,p}$, $\hat{H}_h = \frac{1}{2}[A_{mn} + 2B_{mnq}P_{0,q} + 3G_{mnqp}P_{0,q}P_{0,p}]\hat{Q}_m\hat{Q}_n$ and $\hat{H}'_{anh} = \frac{1}{3}(B_{mnq} + 3G_{mnqp}P_{0,p})\hat{Q}_m\hat{Q}_n\hat{Q}_q + \frac{1}{4}G_{mnqp}\hat{Q}_m\hat{Q}_n\hat{Q}_q\hat{Q}_p$ is the anharmonic Hamiltonian.

[37] Y. V. Shtogun and L. M. Woods, Many-body van der Waals interactions between graphitic nanostructures, *J. Phys. Chem. Lett.* **1**, 1356 (2010).

[38] H.-Y. Kim, J. O. Sofo, D. Velegol, M. W. Cole, and A. A. Lucas, van der Waals forces between nanoclusters: Importance of many-body effects, *J. Chem. Phys.* **124**, 074504 (2006).

[39] H.-Y. Kim, J. O. Sofo, D. Velegol, M. W. Cole, and A. A. Lucas, Van der Waals dispersion forces between dielectric nanoclusters, *Langmuir* **23**, 1735 (2007).

[40] P. L. Silvestrelli, Van der Waals interactions in DFT made easy by Wannier functions, *Phys. Rev. Lett.* **100**, 053002 (2008).

[41] J. F. Dobson and T. Gould, Calculation of dispersion energies, *J. Phys. Condens. Matter.* **24**, 073201 (2012).

[42] M. Stöhr, T. Van Voorhis, and A. Tkatchenko, Theory and practice of modeling van der Waals interactions in electronic-structure calculations, *Chem. Soc. Rev.* **48**, 4118 (2019).

[43] A. Ambrosetti, P. Umari, P. L. Silvestrelli, J. Elliott, and A. Tkatchenko, Optical van-der-Waals forces in molecules: From electronic Bethe-Salpeter calculations to the many-body dispersion model, *Nat. Commun.* **13**, 813 (2022).

[44] F. London, The general theory of molecular forces, *Trans. Faraday Soc.* **33**, 8b (1937).

[45] A. Schweig, Calculation of static electric polarizabilities of closed shell organic π-electron systems using a variation method, *Chem. Phys. Lett.* **1**, 163 (1967).

[46] A. Schweig, Calculation of static electric higher polarizabilities of closed shell organic π-electron systems using a variation method, *Chem. Phys. Lett.* **1**, 195 (1967).

[47] D. Rayane, I. Compagnon, R. Antoine, M. Broyer, P. Dugourd, P. Labastie, J. M. L'Hermite, A. Le Padellec, G. Durand, F. Calvo, F. Spiegelman, and A. R. Allouche, Electric dipole moments and polarizabilities of single excess electron sodium fluoride clusters: Experiment and theory, *J. Chem. Phys.* **116**, 10730 (2002).

[48] T. Cysne, W. J. M. Kort-Kamp, D. Oliver, F. A. Pinheiro, F. S. S. Rosa, and C. Farina, Tuning the Casimir-Polder interaction via magneto-optical effects in graphene, *Phys. Rev. A* **90**, 052511 (2014).

[49] P. Rodriguez-Lopez, W. J. M. Kort-Kamp, D. A. R. Dalvit, and L. M. Woods, Casimir force phase transitions in the graphene family, *Nat. Commun.* **8**, 14699 (2017).

[50] N. L. Manakov, S. I. Marmo, and V. D. Ovsyannikov, Magnetoelectric susceptibilities of atoms, Zh. Eksp. Teor. Fiz. **91**, 404 (1986) [Sov. Phys. JETP **64**, 236 (1986)].

[51] T. G. Pedersen, Coulomb-Zeeman-Stark problem in two dimensions, *Phys. Rev. A* **107**, 022804 (2023).

[52] B. W. Kwaadgras, M. Verdult, M. Dijkstra, and R. van Roij, Polarizability and alignment of dielectric nanoparticles in an external electric field: Bowls, dumbbells, and cuboids, *J. Chem. Phys.* **135**, 134105 (2011).

[53] K. A. Emelyanenko, A. M. Emelyanenko, and L. B. Boinovich, Van der Waals forces in free and wetting liquid films, *Adv. Colloid Interface Sci.* **269**, 357 (2019).

[54] J. D. Jackson, *Classical Electrodynamics* (John Wiley & Sons, New York, 1999).

[55] I. Feranchuk, A. Ivanov, V.-H. Le, and A. Ulyanenkov, *Non-perturbative Description of Quantum Systems* (Springer, Switzerland, 2015).