Derivation of Miller’s rule for the nonlinear optical susceptibility of a quantum anharmonic oscillator

Miller’s rule is an empirical relation between the nonlinear and linear optical coefficients that applies to a large class of materials but has only been rigorously derived for the classical Lorentz model with a weak anharmonic perturbation. In this work, we extend the proof and present a detailed derivation of Miller’s rule for an equivalent quantum-mechanical anharmonic oscillator. For this purpose, the classical concept of velocity-dependent damping inherent to the Lorentz model is replaced by an adiabatic switch-on of the external electric field, which allows a unified treatment of the classical and quantum-mechanical systems using identical potentials and fields. Although the dynamics of the resulting charge oscillations, and hence the induced polarizations, deviate due to the finite zero-point motion in the quantum-mechanical framework, we find that Miller’s rule is nevertheless identical in both cases up to terms of first order in the anharmonicity. With a view to practical applications, especially in the context of ab initio calculations for the optical response where adiabatically switched-on fields are widely assumed, we demonstrate that a correct treatment of finite broadening parameters is essential to avoid spurious errors that may falsely suggest a violation of Miller’s rule, and we illustrate this point by means of a numerical example.


Introduction
Nonlinear optical effects comprise a large number of physical phenomena that arise because the response of atoms, molecules, or solids is, in general, not simply proportional to the incident electric field [1].Instead, the induced polarization features additional Fourier components with other frequencies than the external perturbation.The resulting charge oscillations in turn radiate electromagnetic waves, effectively converting part of the incident light to different frequencies.The most prominent manifestation is harmonic generation, where two or more low-energy photons from a monochromatic field are converted into one high-energy photon with an integer multiple of the original frequency, in a parametric process that preserves both energy and momentum [2,3].Further nonlinearities include the optical Kerr effect and two-photon absorption; the latter is an example of a nonparametric process that changes the quantum state of the sample due to the energy absorption and the resulting electronic excitation.The various forms of optical nonlinearities are exploited in many technological applications, such as tunable coherent light sources [4], electro-optic modulators [5], ultrafast all-optical switching [6] and signal processing [7], as well as advanced imaging techniques for chemical [8] and biological [9] systems.
In order to select the optimal combination of materials and wavelengths for a given photonic application and to further maximize the desired properties by tuning the sample geometry, external doping, or other factors, computer simulations are an effective tool [10,11].While linear optical spectra can be calculated with high accuracy from the Bethe-Salpeter equation [12], which includes excitonic and local-field effects without any empirical parameters, a comparable treatment of the nonlinear response is usually prohibitive due to the much higher numerical cost [13].Despite some theoretical progress [14][15][16][17], most first-principles simulations are hence performed at a much simpler level, the independent-particle approximation, which ignores excitonic contributions as well as localfield effects altogether [18][19][20][21][22][23][24][25].As an alternative, Miller's rule [26] provides an empirical relation between the linear and nonlinear optical coefficients of materials, which can be used to circumvent the explicit numerical computation of higher-order susceptibilities in the first place.
The historical origin of Miller's rule, first published in 1964 [26], was the observation that the ratio between the nonlinear coefficients for second-harmonic generation, measured at a fixed laser wavelength of 1.06 µm, and a certain product of first-order susceptibilities turned out very similar for a set of eleven inorganic crystalline solids, although the coefficients themselves varied by several orders of magnitude.Soon afterwards, measurements over a larger spectral range confirmed that this ratio is nearly independent of frequency [27], and similar relations were found for higher-order nonlinear coefficients, such as third-harmonic generation [28].By now, it is established that Miller's rule extends to materials as diverse as organic molecular crystals [29], ferroelectrics [30], and chalcogenide glasses [31], while optical metamaterials were found to exhibit more significant deviations [32,33].
Due to the enormous simplification that results from deriving nonlinear optical spectra purely on the basis of easily obtainable linear susceptibilities, Miller's rule is now widely used in practical materials research, not only to avoid the expensive computation of higher-order susceptibilities in theoretical simulations [34][35][36], but also to sidestep technically challenging direct measurements of the nonlinear coefficients in experimental studies and to predict them instead from the linear optical absorption or refractive index [37][38][39].Finally, Miller's rule underlies common software packages for modelling the nonlinear optical materials properties [40,41] and plays an important role in modern computer-aided materials discovery, because it facilitates a simple quantitative estimate of various nonlinearities from a small set of tabulated results for the linear optical functions [42].
In light of the broad scope of applications and obvious practical importance of Miller's empirical rule, a proper theoretical justification would be highly desirable, but there has been little progress so far despite substantial efforts.It has long been known that Miller's rule can be rigorously derived for the classical Lorentz model [43], which describes the electrons as independent point particles in a confining potential well with a periodic driving force and velocity-dependent damping.As there are no nonlinearities in a purely harmonic confinement, the derivation involves a weak anharmonic perturbation and considers only terms of first order in the anharmonicity; all higher orders are ignored.In this form, it is featured in pertinent textbooks [1] and can be generalized to any order of the nonlinearity [44].Although it is widely believed that Miller's rule applies equally to quantum-mechanical systems, an analogous mathematical derivation is not straightforward, because the Lorentz model contains some distinctly classical concepts, such as the velocity-dependent damping term, that have no quantum-mechanical counterpart.As a consequence, no full proof has been reported so far, as previous treatments provided few details and consistently ignored the damping [45].However, this simplification gives rise to a purely real dielectric function, which fails to describe absorption and is thus unsuitable for practical use in optics research.Attempts to derive Miller's rule for realistic interacting electron systems in solids require even more drastic uncontrolled approximations [46].On the other hand, despite the lack of a proper mathematical justification, the results for real materials are often analyzed based on the premises of the Lorentz model.For example, a recent theoretical study of the polaronic enhancement of second-harmonic generation in lithium niobate interpreted the large violation of Miller's rule for some elements of the calculated nonlinear susceptibility as evidence that the polaron charge in these defect configurations is trapped in a strongly anharmonic double-well potential instead of a weakly anharmonic potential well as assumed in the classical Lorentz model, although the shape could not be verified directly [47].
In order to provide a better theoretical foundation for Miller's rule, we go beyond earlier studies and present a detailed derivation for the quantum-mechanical anharmonic oscillator in this work, taking the complex nature of the linear and nonlinear susceptibilities fully into account.Following the established treatment of the Lorentz model, we only focus on terms of first order in the anharmonicity.A key feature of our approach is that we do not treat the imaginary part as originating from a classical damping force, but from the adiabatic switch-on of the external electric field, in line with its usual role in ab initio electronic-structure calculations [48].This modification obviates the need for classical particle velocities and allows us to treat classical and quantum-mechanical systems on the same footing.Despite inherent deviations in the induced dynamic polarization, we find that the mathematical expression for Miller's rule is nevertheless identical in both cases up to terms of first order in the anharmonicity.With a view to practical applications, another important result is that the linear susceptibilities that enter Miller's rule must be evaluated not only at different frequencies, as had already been known, but also with different broadening parameters in this scenario.
This paper is organized as follows.In section 2, we explain the theoretical background; besides the definition of linear and nonlinear susceptibilities, we briefly present the known results for the classical Lorentz model together with the previously established expressions for Miller's rule.In section 3, we introduce the modified model used in this work and derive perturbative solutions up to first order in the anharmonicity both in a classical and a quantum-mechanical framework.On this basis, we then compare the resulting expressions for the dynamic polarizations and discuss the implications for Miller's rule.The importance of treating the broadening parameters correctly is further illustrated by means of a numerical example.Finally, we summarize our conclusions in section 4.

Theoretical background
In this section, we introduce key quantities and establish the notation used in this work.We start with the definition of the linear and nonlinear susceptibilities and then briefly review the results for the classical Lorentz model.

Linear and nonlinear susceptibilities
An external electric field causes a displacement of charge carriers inside a sample and thus gives rise to polarization.The induced polarization density, which corresponds to the density of the induced electric dipole moments, can be expanded by orders of the external field E α (r, t) according to where the index α labels the spatial dimensions and each term has the form [1] ..αj (r, r 1 , . . ., r j ; t − t 1 , . . ., t − t j ) × E α1 (r 1 , t 1 ) . . .E αj (r j , t j ) . ( In general, the polarization density at the time t depends on the electric field at all earlier times t i ⩽ t, because the response is not instantaneous.Furthermore, the convolution in real space reflects the fact that charge rearrangements elsewhere in the system influence the effective local field at the point r due to the long-range Coulomb interaction.The key quantity that determines the response of the system is the microscopic dynamic susceptibility χ ( j) , which is written as a tensor of rank j + 1.The first-order susceptibility χ (1) is called linear, because it yields a contribution directly proportional to the external electric field, whereas all higher-order susceptibilities are called nonlinear.To respect causality, the susceptibility in (2) can only be nonzero if all of the time differences t − t i are nonnegative.The symbol ε 0 denotes the vacuum permittivity.
If the wavelength of the incident electric field is much larger than the extent of the sample, as is typically the case in optical experiments with infrared or visible light, then the spatial variation can be neglected, and the electric field may be approximated as E α (r, t) = E α (t).In this case, the total integrated dipole moment is given by can be interpreted as multidimensional Fourier transforms with the wavevector k = 0. Furthermore, if the electric field can be represented as a sum of discrete frequency contributions according to then (3) simplifies to where the Fourier transform of the macroscopic susceptibility from real times to frequencies is defined as with τ i = t − t i .Therefore, the induced polarization is also a sum of discrete frequency components The oscillations of the linear part are always identical to the external electric field, but the nonlinear part exhibits a larger spectrum with additional frequency contributions and corresponding amplitudes where the first summation is over all ordered tuples with the fixed sum ω 1 + . . .+ ω j = ω.Depending on the sign of the individual terms, this phenomenon is called either sumfrequency or difference-frequency generation.
In the special case of a monochromatic field with frequency ω, the complex Fourier series E α (t) = Êα (ω) e −iωt + Êα (−ω) e iωt (10) contains only two terms corresponding to +ω and −ω.For a real-valued physical quantity, the Fourier coefficients must satisfy Êα (−ω) = Ê * α (ω).The most prominent nonlinear phenomena in this scenario, which we also examine in this work, are second-harmonic and third-harmonic generation.The former originates from χ (2)  αα1α2 (2ω; ω, ω) Êα1 (ω) Êα2 (ω) (11) but is only observable in noncentrosymmetric samples, while the latter originates from and occurs in all systems.Another notable second-order process is optical rectification, where an oscillating, timedependent electric field gives rise to a constant, timeindependent polarization described by In line with established conventions, we use the short-hand notation χ there is no ambiguity.

Classical Lorentz model
The Lorentz model is a phenomenological model that predates quantum mechanics and treats electrons in a material as classical particles bound to stationary atomic nuclei by springs, together with a damper to ensure a finite response at the resonance frequency.The restoring force must include an anharmonic contribution in order to describe optical nonlinearities.Under these assumptions, the response of electrons with mass m and charge −e 0 to a time-dependent electric field E(t) in one spatial dimension is determined by the classical equation of motion where ω 0 is the resonance frequency of the harmonic oscillation and τ is the relaxation time.The external electric field is assumed to be monochromatic and can be written as where c.c. denotes the complex conjugate.For a weak anharmonicity, the equation of motion ( 14) can be solved by perturbation theory.For this purpose, the displacement is expanded by orders of the parameters β 2 , β 3 according to The zeroth-order term corresponds to a damped driven harmonic oscillator and has the form For a pure harmonic potential, the period of the electronic oscillation is hence identical to the external perturbation, and the amplitude of the displacement is proportional to the electric field.In contrast, the first-order term contains a number of additional discrete frequency components, whose Fourier coefficients are given by As x0 (ω) ∝ Ê(ω) according to (18), the components with frequencies 0 and 2ω are of second order in the external electric field, while those with ω and 3ω are of third order, although all are linear in the anharmonicity parameter β 2 or β 3 .Terms of second or higher order in the anharmonicity parameters are not considered in this work.

Miller's rule for the Lorentz model
For an ensemble of N independent electrons, each carrying the charge −e 0 , the total polarization equals The linear part stems exclusively from the zeroth-order term x 0 (t), the only component of the displacement directly proportional to the electric field.By inserting the expression (18) into the Fourier representation we obtain the explicit formula for the linear susceptibility.Likewise, by exploiting (25) in order to substitute x0 (ω) in ( 20)-( 23), the higher-order nonlinear susceptibilities can be written in terms of the linear susceptibility as 27) 28) 29) These relations are the essence of Miller's rule, which, in its original empirical formulation, states that the ratio between the nonlinear susceptibility χ (2) (2ω) for secondharmonic generation and the product χ (1) (2ω)χ (1) (ω) 2 of three first-order susceptibilities evaluated at different frequencies is essentially independent of ω and, moreover, varies little between different materials.For the Lorentz model, the above results confirm that is a constant function, and Miller's rule is therefore satisfied exactly in this case.Similar rules can be derived from ( 27)- (30) for the other Fourier components of the nonlinear susceptibility, such as third-harmonic generation.

Results and discussion
Although the derivation of Miller's rule for the classical Lorentz model is straightforward, as shown in the previous section, it remains unclear whether and how these results can be generalized to quantum-mechanical systems.The principal obstacle is that the imaginary part of the susceptibility (26) vanishes unless the relaxation time τ is finite, but the origin of this parameter within the Lorentz model is a velocity-dependent friction force, a distinctly classical concept with no quantum-mechanical analogue.On the other hand, the susceptibility and the associated dielectric function become unphysical if the scattering time is neglected, because a purely real response function diverges at the resonance frequency ω 0 and fails to describe essential phenomena like the optical absorption spectrum, which corresponds to the imaginary part of the dielectric function, even at a qualitative level.In order to overcome these problems, we propose a modified model in this work, which avoids the classical notion of velocitydependent forces and instead treats the imaginary part of the dielectric function in the same fashion as in modern quantummechanical ab initio simulations.Although closely related to the Lorentz model, it can hence be treated on the same footing in a classical and a quantum-mechanical framework and thereby allows us to investigate the validity of Miller's rule in both cases.The model consists of independent electrons in an anharmonic potential that includes a cubic and a quartic term, subject to an adiabatically switched on monochromatic electric field with γ > 0, which gives rise to a time-dependent potential energy +e 0 E(t)x.In contrast to classical physics, where metastable states can well reside in local minima of unbounded potentials, the quartic term in V(x) is necessary here to ensure the existence of a proper quantum-mechanical ground state, which requires the potential to be bounded from below.We note that, in principle, it is possible to consider more general switching functions, but nonadiabatic switching necessarily involves different frequency components to achieve the desired time dependence and thus implies additional nonlinearities due to complicated sum-frequency mixing.Such processes are not considered in this work, which focuses on harmonic generation.Furthermore, adiabatic switching is ubiquitous in ab initio calculations, and the proper application of Miller's rule in this scenario is hence of special concern.

Classical solution
The classical equation of motion for our model is which differs from the Lorentz model ( 14) by the absence of the damping term and the different time dependence of the external electric field, now given by (34).Analogous to the treatment in section 2.2, the equation of motion can be solved by perturbation theory.For this purpose, the displacement x(t) is expanded by powers of the parameters β 2 , β 3 as in (16).
Inserting this expression into (35) and comparing terms of equal order in the perturbation yields the set of equations All terms of second and higher order in β 2 , β 3 are again disregarded.The zeroth-order differential equation ( 36) corresponds to a driven harmonic oscillator with the solution The resulting Fourier coefficient is very similar but not identical to (18).In particular, it can be written exactly as the sum of two complex Lorentzian functions, which is not true for the classical Lorentz model and justified only as an approximation in the limit of large relaxation times.A comparison of the two expressions demonstrates that the adiabatic parameter γ fulfills a similar mathematical role as τ /2, generating an imaginary part even without classical velocity-dependent damping.After inserting x 0 (t) on the right-hand side of (37), it is straightforward to confirm that the ansatz solves the first-order differential equation with the complex Fourier coefficients

Quantum-mechanical solution
The quantum-mechanical analogue to the classical equation of motion is the Schrödinger equation.It determines the time evolution of the wavefunction ψ(x, t), from which all observable properties can be derived.The Schrödinger equation for the system considered in this work has the form with the time-dependent Hamiltonian of a driven harmonic oscillator and the time-independent anharmonic part of the potential Following the same approach as in the classical case, we employ perturbation theory and expand the wavefunction by orders of the parameters β 2 , β 3 according to Inserting this expression into (45) and comparing the terms in each order again yields a set of differential equations that can be solved in succession to obtain the individual components of the wavefunction.Although this strategy appears natural in light of the analogous classical solution, we note that it is in fact highly unusual from a technical point of view: Standard applications of time-dependent perturbation theory in quantum mechanics assume a well-defined stationary initial state and then treat time-dependent external fields in a perturbative manner.In this case, the roles are reversed, however, as we take the time-dependent Hamiltonian of the driven harmonic oscillator H (0) (x, t) as the unperturbed system and the time-independent anharmonic part of the potential H (1) (x) as the perturbation.This causes some mathematical challenges but has the advantage that it allows a nonperturbative treatment of the external electric field without truncation at a finite order.

Zeroth order: driven harmonic oscillator.
The timedependent Schrödinger equation of a driven harmonic oscillator in one dimension is Its exact analytic solution, first derived by Schrödinger [49], involves a sequence of transformations that reduce it to a simpler solvable form.In the first step, the coordinate system is transformed by substituting x = y + ζ(t), where ζ(t) initially denotes an arbitrary time-dependent function.This defines a new shifted wavefunction With this change of variables, the Schrödinger equation becomes where where we have introduced the classical Lagrangian for the driven harmonic oscillator as a short-hand notation.The latter only depends on t but not on y.Therefore, it does not affect the spatial variation of the wavefunction but only its phase.In the third step, another unitary transformation ϕ (y, t) =: e i ´t 0 L(x0(t ′ ),ẋ0(t ′ ),t ′ )dt ′ /h ξ (y, t) (56) This is the time-dependent Schrödinger equation of a simple harmonic oscillator, whose stationary solutions are well known and can be written as in terms of the normalized eigenfunctions and the corresponding energy eigenvalues E n = hω 0 (n + 1/2) of the time-independent Hamiltonian, where n is a nonnegative integer.The symbol H n denotes the Hermite polynomial of order n.Finally, by reverting the transformations (56), (52), and (50), we obtain the desired solutions of the original Schrödinger equation ( 49) of the driven harmonic oscillator The classical trajectory (38) and the external electric field (34) are both known exactly.Therefore, it is straightforward to evaluate the integral over the Lagrangian L(x 0 (t ′ ), ẋ0 (t ′ ), t ′ ) analytically, but this is in fact not necessary in the present context, because the phase shift does not contribute to the physical observables and cancels out in the subsequent steps.Besides the canonical coherent states considered here, other solutions of (57) with distinct physical properties and no classical analogues, such as pulsating states [50], may also be derived.They correspond to different initial conditions and will be ignored in the following.As an example, the pulsating states are time dependent because the width of the wave packet varies periodically, but as it remains symmetric around the origin, there is no observable polarization.The quantities of central interest in this work are the induced polarization and the associated susceptibilities.The dipole moment P(t) = −Ne 0 ⟨x(t)⟩ is proportional to the spatial displacement of the electron charge, which is defined as the expectation value of the position operator x.The zeroth-order contribution is obtained from the wavefunctions of the driven harmonic oscillator, taking the full strength of the external electric field but no anharmonic perturbation into account.By inserting the previously derived solution (60) and substituting x = y + x 0 (t), this simplifies to The first integral vanishes, because the integrand is an odd function of y, while the second integral equals the norm of the wavefunction and yields 1.As a consequence, the quantummechanical expectation value coincides exactly with the solution of the classical equation of motion and hence with the classical trajectory derived in section 3.1, where the Fourier coefficient is given by (39).Notably, it does not depend on the quantum number n and is therefore identical for all eigenstates.

First order: anharmonic perturbation.
The terms of first order in the parameters β 2 , β 3 from the Schrödinger equation ( 45) fulfill the relation The first-order wavefunction is obtained from time-dependent perturbation theory as which exploits the fact that the eigenfunctions (60) of the harmonic oscillator at a fixed time t, even in the presence of an external driving force, form a complete orthonormal set that satisfies the orthogonality condition as well as the completeness relation where δ n,n ′ denotes the Kronecker delta and δ(x − x ′ ) the Dirac delta distribution.The corresponding first-order contribution to the electron displacement, i.e. the expectation value of the position operator, is given by After inserting (65), this can be rearranged to The term n ′ = n may be omitted, because all integrals are real in this case.Together with the prefactor 1/(ih), this makes for a purely imaginary contribution to the sum, which is exactly canceled by its complex conjugate.Next, we insert the previously derived solution (60) for the wavefunctions of the driven harmonic oscillator.After the substitutions x = y + x 0 (t) and x ′ = y ′ + x 0 (t), we thus obtain In order to evaluate the integrals, we define the matrix elements which can be calculated from a recursion formula as shown in the appendix, where the relevant values for j ⩽ 4 are also listed.The integral over y in (70) equals Ω nn ′ (1) + x 0 (t)Ω nn ′ (0) and yields As n ′ = n is excluded in (70), only the terms with n ′ = n + 1 or n ′ = n − 1 contribute to the sum, which leaves The integral over y ′ can be calculated in the same way, because the Hamiltonian H (1) (y ′ + x 0 (t)) is a polynomial of degree 4 in y ′ .Using the formulas for Ω n±1,n (j) listed in the appendix, the relevant contributions are Inserting these expressions into (73) yields where we have introduced the symbol as a short-hand notation for the integrals over t.By exploiting the fact that O + j (t) and O − j (t) are, by definition, complex conjugates of each other, this can be simplified to which depends only on the real-valued sums O + j (t) + O − j (t).As the trajectory x 0 (t) of the classical driven harmonic oscillator is known, the integrals can be evaluated analytically; the details of the derivation and explicit expressions for j ⩽ 3 are again listed in the appendix.When these are inserted, we obtain the final result Compared to the classical solution (40), the quantummechanical expectation value contains two additional terms, a constant displacement that is independent of the external electric field and a firstorder correction to the linear response ⟨x (1) both of which depend on the quantum number n through the common factor 2n + 1.The remaining terms on the right-hand side of (81), which are of second and third order in the external electric field, coincide with the classical trajectory x 1 (t) and do not depend on n.

Miller's rule
For the driven harmonic oscillator, the quantum-mechanical expectation value ⟨x 0 (t)⟩ is identical to the classical trajectory x 0 (t), but this is no longer true if the potential contains an anharmonic part, as the first-order correction ⟨x 1 (t)⟩ to the quantum-mechanical solution differs from that of the classical trajectory x 1 (t).The origin of the two additional terms in (81) can be easily understood, as illustrated in figure 1 for a weakly anharmonic potential with the parameter values V(x) = 1 2 x 2 + 1 10 x 3 + 1 100 x 4 in units where lengths are measured in multiples of √ h/(mω 0 ) and energies in multiples of hω 0 .In the absence of external driving forces, the equilibrium position of a classical particle is at x = 0 for the harmonic potential 1  2 x 2 as well as the full anharmonic potential V(x), because V(0) = 0 is the global minimum in both cases.On the other hand, quantum-mechanical particles possess a finite zero-point energy of E 0 = 1/2 in these units, and the groundstate wavefunction extends into the barriers on either side.In the harmonic case, the probability density is an even function of x, like the potential itself, so that the expectation value of the position operator is x = 0 by symmetry and equals the classical equilibrium position.If anharmonic terms are included, the distribution becomes asymmetric, however, and the centroid shifts towards the flatter slope.As a consequence, the expectation value is offset from the classical equilibrium position as described by the constant ⟨x (0) 1 ⟩.If the external electric field is turned on, then the modified shape of the potential away from the origin also influences the dynamic behavior of delocalized quantum-mechanical particles more strongly than that of classical particles, whose motion is confined to the vicinity of x = 0, leading to the additional term ⟨x 1 (t)⟩ that depends on the field strength.As the wavefunction becomes more delocalized and extends further into the barriers for higher energies, the magnitude of these terms increases with the quantum number n.
Although the dynamics of the classical and quantummechanical particles deviate, they are still closely related.In order to explore the implications for Miller's rule, we now sort the results, which were grouped by orders of the anharmonicity parameters β 2 and β 3 so far, according to powers of the electric field.The zeroth-order component ⟨x 0 (t)⟩ = x 0 (t) is linear in Ê(ω), whereas the first-order 1 (t)⟩ contains quadratic and cubic terms and, only in the quantum-mechanical case, additional constant and linear contributions ⟨x (0)

⟩ and ⟨x
(1) 1 (t)⟩, respectively.From the linear displacement ⟨x 0 (t)⟩ + ⟨x (1) 1 (t)⟩ = x 0 (t) + O(β 3 ), we thus obtain the Fourier coefficients of the induced polarization where the deviation between the response of classical and the quantum-mechanical particles affects only the terms of first order in β 3 .In both cases, the linear susceptibility is hence given by χ (1) From now on, to keep track of the adiabatic switching parameter γ that determines the linewidth of the Lorentzian functions, we write it explicitly in the argument of the first-order susceptibility, which may thus be interpreted as a function of a complex frequency variable.Likewise, the relevant component for second-harmonic generation is According to the definition, this entails a nonlinear susceptibility To keep the notation simple, we omit the parameter γ in the argument of the nonlinear susceptibility, because there is no ambiguity in this case.By comparing the formulas above, one sees that the ratio χ (1) (2ω + i2γ) χ (1) is frequency independent if terms of second order in β 2 and β 3 are ignored, as done throughout this work.Despite the different dynamics, we thus find that the classical and the quantummechanical system fulfill Miller's rule and, moreover, that its mathematical formulation is indeed identical in both cases up to first order in the anharmonicity.

Numerical example
The mathematical relations between the linear and nonlinear susceptibilities derived in the previous section are very similar but not identical to those for the classical Lorentz model in section 2.3 due to the differing roles of the adiabatic switching parameter γ and the relaxation time τ .In particular, Miller's rule for the Lorentz model (31) relates the second-order susceptibility χ (2) (2ω) to a product of three first-order susceptibilities, which are evaluated at different real frequencies ω and 2ω but with the same value of τ .In contrast, the three first-order susceptibilities in the equivalent formula (88) for our model feature distinct switching parameters γ and 2γ.To avoid any spurious errors, it is important to recognize this discrepancy and to apply Miller's rule correctly in practical situations.In particular, ab initio electronic-structure calculations of optical response functions for solids and their surfaces typically assume adiabatically switched-on fields [48].Although the distinction between γ and 2γ might at first appear irrelevant, since the theoretical concept of adiabatic switching usually implies the notional limit γ → 0 at the end of the calculation, this is not true in actual numerical implementations.On the contrary, most material-specific simulations deliberately use a finite, not too small value of γ in order to broaden the spectral peaks and smoothen the generated functions, which reduces the computational cost because fewer data points are needed.Under these circumstances, a proper treatment of the broadening parameters in Miller's rule is clearly mandatory.
In the following, we illustrate this point numerically, taking the weakly anharmonic potential displayed in figure 1 as an example.The upper panel of figure 2 shows the real and imaginary part of the first-order susceptibility χ (1) in units of Linear susceptibility χ (1) , nonlinear susceptibility χ (2)  for second-harmonic generation, and Miller's δ (2) for the anharmonic potential displayed in figure 1. Miller's rule is only fulfilled exactly if χ (2) (2ω) is calculated from the product χ (1) (2ω + i2γ)χ (1) (ω + iγ) 2 of first-order susceptibilities evaluated at different frequencies and broadening parameters, corresponding to a combination of points on the blue and red curves in the upper panel.If the same broadening γ is used throughout (blue curves), then δ (2) has a spurious artificial frequency dependence.
Ne 2 0 /(ε 0 m) as a function of the real frequency ω, measured in multiples of ω 0 , for two different broadening parameters γ and 2γ; in this example, we choose γ = ω 0 /4.In the central panel, we display the second-order susceptibility in units of β 2 Ne 3 0 /(ε 0 m 2 ).According to (90), the exact solution χ (2) (2ω) is proportional to the product χ (1) (2ω + i2γ)χ (1) (ω + iγ) 2 .The construction of each data point thus requires a value from the blue curve and another from the red curve in the upper panel, which correspond to a broadening of γ and 2γ, respectively.Under these conditions, Miller's rule in the formulation of (88) is satisfied exactly, and δ (2) (2ω) is indeed a constant, frequency-independent function, as indicated in the bottom panel of figure 2. In contrast, if the import of the finite value of γ is ignored and χ (1) (ω + iγ) is merely treated as a numerical approximation for χ (1) (ω) in the formulas ( 28) and (31) derived for the classical Lorentz model, then one would take all input data just from the blue curves in the upper panel of figure 2. In this scenario, the product χ (1) (2ω + iγ)χ (1) (ω + iγ) 2 still resembles the exact χ (2) (2ω), but deviations in the peak heights, peak positions, and lineshapes are clearly discernible.Furthermore, the ratio between the two, displayed as blue lines in the bottom panel, exhibits a substantial frequency variation, amounting to a dip by 50% at ω = ω 0 /2, and thereby falsely suggests a violation of Miller's rule that is really just an artifact of the incorrect treatment of the finite broadening.

Conclusion
Miller's rule provides an empirical relation between the nonlinear and linear optical coefficients of materials, which allows a quantitative estimate of the former based on knowledge of the latter.Although it is often used to analyze data obtained from spectroscopic measurements or quantitative theoretical simulations, its formal justification is mostly based on a mathematically straightforward derivation for the classical Lorentz model.In particular, there was so far no actual proof that Miller's rule also applies to quantum-mechanical systems, as previous studies provided few details and failed to incorporate damping, which is essential to obtain a complex dielectric function and to describe optical absorption.The principal obstacle was that the Lorentz model includes damping by means of a velocity-dependent friction force, a distinctly classical concept with no quantum-mechanical analogue.
In this work, we set out to provide a unified proof of Miller's rule for the classical and the quantum-mechanical anharmonic oscillator.Compared to the Lorentz model, a key aspect of our approach was to replace the velocity-dependent damping by an adiabatic switch-on of the external electric field, a technique that is well established in electronic-structure calculations and routinely used to compute complex dielectric functions in ab initio simulations of the optical response.This modification allowed us to treat the classical and quantum-mechanical systems on the same footing, using identical potentials and external fields.Both were solved analytically in an analogous fashion, where the time-varying electric field is fully included in a nonperturbative way, while the anharmonicity is treated up to first order in time-dependent perturbation theory.Although the dynamics of the two systems, and hence the induced polarizations, deviate due to the finite zero-point energy in the quantum-mechanical framework, we found that the mathematical relation between the components of the nonlinear and the linear susceptibility is nevertheless identical in both cases up to terms of first order in the anharmonicity.Thus, the proof of Miller's rule is generalized from weakly anharmonic classical oscillators to equivalent quantum-mechanical systems, providing a better foundation and justification for its applications to real materials.
The most immediate consequence of the proper treatment of the adiabatic switch-on in our model is that the derived formulas (89)-(92) for Miller's rule, contrary to those for the friction-damped Lorentz model, feature products of firstorder susceptibilities that are not only evaluated at distinct real frequencies but also with different switch-on parameters.By means of a numerical example, we demonstrated in section 3.4 that proper implementation is indeed very important, as using the same parameter value in all factors may incur substantial errors in the peak heights, peak positions, and lineshapes that falsely suggest a violation of Miller's rule but are in fact mere artifacts.Applications of Miller's rule in the context of ab initio optical-response calculations [34,47], where adiabatically switched-on fields are ubiquitous and finite broadening parameters are routinely used to smoothen the spectral functions and reduce the numerical cost, should pay heed to this point in particular and use appropriate input data with varying broadening parameters in order to ensure correct quantitative predictions.
Finally, although the formal scope of the present proof is limited to weakly anharmonic oscillators in one dimension, the techniques employed in our derivation may be adapted for other quantum-mechanical systems and thus open a pathway for analyzing the validity of Miller's rule in more general situations.With a view to realistic nonlinear optical materials, the most relevant factors to be considered in future studies are the dimensionality, the shape of the potential, and the mutual Coulomb interaction between the electrons.leads to the recursion formula which can be iterated as necessary.The relevant matrix elements required in section 3. ) . (100) The integrals O + j and O − j defined in (79) are complex conjugates of each other.With the Fourier representation of the classical trajectory of the harmonic oscillator x 0 (t) = (x 0 (ω)e −iωt + x * 0 (ω)e iωt )e γt , they can be written as The integration is elementary and yields (106) The evaluation of (80) in fact only requires the real-valued sums O + j + O − j , which are given by

Figure 1 .
Figure 1.Potential (solid lines) and associated ground-state wavefunction (filled curves) for a harmonic (blue) and a weakly anharmonic (red) oscillator.The quantum-mechanical expectation values of the position operator are marked by vertical lines.In the harmonic case, it coincides with the classical equilibrium position, indicated by the black ball, but the anharmonic perturbation causes a finite deviation ⟨x (0) 1 ⟩ even without external electric fields.

Figure 2 .
Figure2.Linear susceptibility χ(1) , nonlinear susceptibility χ(2)  for second-harmonic generation, and Miller's δ(2) for the anharmonic potential displayed in figure 1. Miller's rule is only fulfilled exactly if χ(2) (2ω) is calculated from the product χ(1) (2ω + i2γ)χ(1) (ω + iγ) 2 of first-order susceptibilities evaluated at different frequencies and broadening parameters, corresponding to a combination of points on the blue and red curves in the upper panel.If the same broadening γ is used throughout (blue curves), then δ(2) has a spurious artificial frequency dependence.