Michel Rérat
email: michel.rerat@univ-pau.fr

Philippe D'arco

Valentina Lacivita

Fabien Pascale

Roberto Dovesi

From anisotropy of dielectric tensors to birefringence: a quantum mechanics approach

Keywords: Refractive index, Birefringence, (Non)linear electric susceptibility tensor, Anisotropy, Quantum mechanical simulation, CRYSTAL code, Gaussian-type basis set

The way quantum mechanical ab initio computer codes allow to compute, through perturbation theory (the so-called SC-CP, self-consistent coupled-perturbed scheme), many properties resulting from the interaction of the electric field with a crystalline system is illustrated. The polarizability, which leads to the dielectric tensors as well as to the refractive indices and to he birefringence of materials, is the simplest on this list. Higher order tensors, like the first and second hyperpolarizabilities, can be obtained as well with the CRYSTAL code here used. These properties, resulting from the Taylor expansion of the total energy of the solid as a function of the electric field, belong to a large family of phenomena generated by combining in different ways the frequencies of the fields. Second-harmonic generation (SHG), Pockels effect, intensity-dependent refractive index (IDRI), and other quantities now accessible to experiment can be computed at a relatively low cost and with high accuracy.

Introduction

In this contribution, we illustrate the way modern quantum mechanical methods allow to compute the (hyper)polarizability tensors and, consequently, optical properties such as the refractive index and birefringence, through which the anisotropy of the physical properties of crystalline compounds manifests itself.

The properties mentioned above can be obtained by investigating the interaction of the electromagnetic field with a periodic infinite system (the model implies, without serious consequences, that the crystalline compound is infinite).

The equations describing this interaction can be formulated at various levels (for example: relativistic or non-relativistic quantum mechanics). As none of these equations can be solved exactly, many approximations must be introduced, whose importance should be discussed carefully, and one should possibly verify numerically how severe these approximations are.

In a very broad sense, all these equations are solved by performing at various steps series expansions, so that the differential equations transform in matrix equations.

This essentially requires: (a) to compute matrix elements (these, in turn, are the sum of integrals, many of which are This paper is the peer-reviewed version of a contribution presented at the Conference on Anisotropic Properties of Matter, organized by Giovanni Ferraris and held at Accademia Nazionale dei Lincei in Rome, October 16-17, 2019. 1 3 bielectronic four center six dimensional integrals); (b) sum, multiply, and diagonalize matrices that can easily reach very large dimensions ( 10 3 -10 6 , as typical cases).

The larger the matrices, the more accurate the calculation.

Obviously, this kind of linear algebra requires the use of computers (clusters of PC) or supercomputers containing 10 2 -10 4 processors.

The starting point for describing the crystalline system is the (stationary, or time-independent) Schrödinger's equation:

where Ĥ is an operator, called Hamiltonian in memory of classical mechanics developed by Lagrange (1736Lagrange ( -1813) ) and then Hamilton (1805Hamilton ( -1865)), and E the energy associ- ated to the wavefunction .

The real limit of this equation is that Ĥ and depend on many variables as atomic positions, that is ≡ ( 1 , 2 , 3 , … , N ) . In the case of an infinite crystalline system, N goes to infinity. But also for a molecule like, say, benzene, containing 12 atoms and 42 electrons, for a total of 162 Cartesian coordinates (plus the spin), the Schrödinger's equation cannot be solved exactly, and its approximate solution requires a huge amount of skills and computational effort. We will not dwell on the techniques, hypotheses, and approximations that bring this intractable problem to something that can be tackled. The interested reader can refer to several excellent textbooks covering quantum mechanics and computational chemistry methods. We simply mention that these approximations bring to the so-called Hartree-Fock or Kohn-Sham methods, in which a single particle (say electron) is moving in the field created by all the other electrons (whose wavefunction is unknown; mean field theory). This implies, in turn, that these equations must be solved through a self-consistent field (SCF) scheme.

We can now suppose that we are able to describe with reasonable accuracy the ground state of a crystalline system.

We are then faced with the problem of the description of the electromagnetic field, and of its interaction with the solid. This interaction is described through a well-known tool of quantum mechanics, namely the perturbation theory, that takes the form of a Taylor expansion of the system's energy in powers of the electric field. As usual, the series is truncated after a few terms due to (a) computational costs and (b) hopefully, the rapid convergence. Also in this case, the solution of the resulting equations requires an iterative scheme. At the very end, the crucial points (crucial due to numerical accuracy and computational cost) are: evaluation of multicenter integrals and multiplication of very large matrices.

We are then considering a crystalline solid, and an electric field operating on it. What is the information that can be obtained as a response of the system to this perturbation? Let (1) Ĥ𝛹 = E𝛹 , us consider the Taylor expansion of the total bulk energy of the system with respect to the field amplitude 0 , truncated to the fourth order: where 0 , 0 , 0 , and 0 are the permanent dipole moment, polarizability, and first and second hyperpolarizabilities of the free system, respectively (the conventional negative sign is such that the dipole moment is defined as the sum over the charges multiplied by their position, and such that the polarizability of the ground state is positive). The symbol ⊗ indicates the outer product of vectors. Given the electric field conversion coefficient from the atomic units (a.u.) system to the international system of units (SI): 1 a.u. = 5.14 × 10 11 Vm -1 , a large static field of say 50 kVcm -1 (maximum field amplitude that can be applied with electrodes at the surface of a slab before electric breakdown), is smaller than 10 -5 a.u. . If we apply a field of this amplitude to a molecule, say water, with a polarizability equal to ∼ 10 bohr 3 (the polarizability unit is equivalent to a volume in a.u.), the energy variation due to the polarization would be equal to 10 -9 E h . As the ratio between the terms appearing in Eq. 2 ( 0 ∕ 0 , 0 ∕ 0 ) is generally smaller than 10 3 a.u. , to be multiplied by an additional field intensity of 10 -5 , it is clear that contributions approach rapidly the numerical accuracy limit of quantum mechanical calculations. This is why, 0 and 0 coefficients have been considered in the past of low interest, and terms as 0 , corresponding to power five of 0 , have been neglected in Eq. 2. However, if the expanded quantity is not the total energy, but some higher term evaluated analytically, then the power of 0 for obtaining 0 and 0 is lower. This is the case, for example, when the polarizability 0 is computed analytically, with a reduction by 2 of the power of 0 . Moreover, lasers with much higher intensity than static electric fields are now available, which can allow to access experimentally many non-negligible second-and third-order non-linear optical (NLO) effects.

It should be stressed that Eq. 2 can provide a lot of information:

(a) as 0 is a vector with three components, E x , E y , and E z , it turns out that 0 , 0 , 0 , and 0 are tensors of rank 1, 2, 3, and 4, respectively, whose components can vary from case to case, allowing access to important specific features. (b) The electric fields appearing in Eq. 2 can be different from each other (say 1 , 2 … ), and combined in dif- ferent ways.

(2) Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

E = E 0 -0 ⋅ 0 - 1 2! 0 0 ⊗ 0 - 1 3! 0 0 ⊗ 0 ⊗ 0 - 1 4! 0 0 ⊗ 0 ⊗ 0 ⊗ 0 ,

R

(c) Obviously, the various electric fields can depend on frequency as in light radiation, ≡ ( ) = 0 cos( t) with the angular -frequency of the corresponding photon ( ℏ with ℏ = h∕2 , h being the Planck constant equal to 6.63 × 10 -34 Js ), further increasing the number of possible ways of perturbing the system. For example, in the second-harmonic generation (SHG) experiments, interaction of light with matter can provide one scattered photon of energy ℏ2 from two photons of energy ℏ , the intensity of the scattering light depending on the frequency-dependent first hyperpolarizability, (-2 ; , ) , of the material, as we will see later on. (d) As a corollary to point (c), the frequency can be used for perturbing both nuclei and electrons, or just the latter.

The following section will deepen some of the main topics introduced here.

Methods

Dipole moment and (hyper)polarizability

Definitions

The dipole moment of a finite system is a vector defined, in the atomic unit system, as:

where the absolute value of the electron charge is equal to 1 ( e = -1 a.u. ). Z N and are the nuclear charge and position of the Nth atom, and ( ) is the electron charge density in . For a given geometry, Z N and N are fixed, and the value of the dipole moment can be obtained if the electron density is known in the whole space. In quantum chemistry, this latter, which is the square of a wave function, , describing the electronic state of the system (generally its ground state) can be obtained by solving the time-independent Schrödinger's equation (see Eq. 1). Using Dirac's notation, ∫ ( )d in Eq. 3 can now be replaced by ⟨0� �0⟩ where ⟨0� and �0⟩ rep- resent the bra and ket of the ground state.

In the presence of a time-dependent electromagnetic (E(r,t), B(r,t)) field, the expression of Ĥ becomes:

where and are the momentum and position of the electron, m its mass ( m = 1 a.u. ), c the speed of light ( c = 137. a.u. ), and V C the Coulomb potential. The electromagnetic field is

(3) = ∑ N Z N N -∫ ( )d , (4) Ĥ = 1 2m - e c ( , t) 2 + eU( , t) + V C ( ),
defined from the vector A(r,t) and scalar U(r,t) potentials via Maxwell's equations:

and there is an infinite number of (A(r,t),U(r,t)) couples describing one electromagnetic field. Obviously, the solutions of the time-dependent Schrödinger's equation must be independent from the choice of the so-called gauge. This is indeed the case of the energy E-eigenvalue as well as of | ( , t)| 2 . Only the phase of the wave function depends on the choice of the gauge. For a static electric field: ( ) = 0 e ⋅ with a small wave number: q = 2 ∕ compared to (i.e., for -wavelength much larger than the size of the studied system: 𝜆 ≫ 1 Å ), ( ) is generally defined from the scalar potential only as follows:

which leads for a constant field 0 to: in the electric dipole moment approximation for finite systems (0D or molecules), Ĥ0 being the Hamiltonian operator of the unperturbed system.

The energy E of the molecule in presence of the 0 -field can be developed as a Taylor series as shown in Eq. 2; similarly, for the dipole moment, we have: where 0 is the permanent dipole moment of the molecule, the linear response 0 is called the polarizability, and 0 and 0 are the first and second hyperpolarizabilities. The previous electric responses to the electric field are intrinsic properties of the molecule: they are equal to the first, second, and third derivatives of the induced dipole moment ( -0 ) with respect to the field 0 , at zero field ( 0 → ). Moreo- ver, being 0 and 0 vectors, 0 , 0 and 0 are tensors of rank 2, 3, and 4, respectively.

In the case of a frequency-dependent electric field ( ( ) = 0 cos t ), the linear term of the dipole moment induced by the oscillating field is oscillating at the same frequency of the field (forced oscillation if the field frequency is small compared to proper resonance frequencies of the system) leading to a frequency-dependent polarizability, ( ) . For a molecule in its electronic ground state ( �0⟩ ), ( 5) Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

( , t) = - ( , t) t -U( , t) (6) ( , t) = × ( , t), (7) ( ) = -U( ), (8) Ĥ = 2 2m + V C ( ) -e ⋅ 0 = Ĥ0 -e ⋅ 0 (9) = - dE d 0 = 0 + 0 0 + 1 2! 0 0 ⊗ 0 + 1 3! 0 0 ⊗ 0 ⊗ 0 + ⋯ , R E V I S
the "dynamic" polarizability is a sum of two contributions that can be separately calculated in the Born-Oppenheimer approximation:

1. The electronic ( e ) contribution, for a fixed ( N ) geometry:

where n = E n -E 0 are the allowed transition energies from the �0⟩-ground to �n⟩-excited electronic states with E 0 and E n as corresponding eigenvalues of the unper- turbed Hamiltonian operator Ĥ0 (if we are interested in the polarizability of the ground state) and n = ⟨0� �n⟩ is the corresponding transition dipole moment.

The associated mean value oscillator strengths:

are such that ∑ n f n is equal to the number of electrons involved in these transitions, and the ( n , f n ) couples reproduce the UV-visible spectrum.

See the work of [START_REF] Orr | Perturbation theory of the non-linear optical polarization of an isolated system[END_REF] for the expression of hyperpolarizabilities. 2. The ionic or nuclear relaxation ( nr ) contribution:

where i is the frequency of the normal mode ( Q i ) and i the Born charge ( d 0 ∕dQ i with 0 as defined in Eq. 3 at zero field) obtained at the equilibrium geometry. The ( i , i ) couples reproduce the infrared (IR) spectrum.

We refer to the work of [START_REF] Kirtman | On the contribution of mixed terms in response function treatment of vibrational nonlinear optical properties[END_REF] for the treatment of the vibrational hyperpolarizabilities that can be evaluated from the IR and (hyper) Raman spectra.

nr is zero in purely covalent materials as diamond or silicon, and also negligible for other ionic materials if the field frequency ( ) corresponds to the UV-visible range of energy 𝜔 ≫ 𝜔 i . Indeed, the Born charge d 0 ∕dQ i , which is actually a charge divided by the square root of ion mass, is small with respect to 1 a.u., and

2 i ∕( 2 i -2 ) ≃ -2 i ∕ 2
, which is then negative, becomes very small in absolute values compared to the static vibrational polarizability contribution: 2 i ∕ 2 i . Also the vibrational contribution to and can be neglected in the UV-visible frequency range for the same reason as for , if all the external fields are frequencydependent. The exception is the particular case of the intensity-dependent refractive index, IDRI, a process depending on (-; , -, ) , that includes terms with opposite sign (10)

e (𝜔) = ∞ ∑ n≠0 2𝜔 n n ⊗ n 𝜔 2 n -𝜔 2 , ( 11
) f n = 2 3 n ⟨0� �n⟩ 2 (12) nr (𝜔) = 3N-6 ∑ i=1 i ⊗ i 𝜔 2 i -𝜔 2 ,
phases, ± , generating then a static field and a vibrational contribution.

Calculation of the microscopic response properties

The electronic contribution to the polarizability (Eq. 10) is a second-order perturbation energy, the expression of which is a sum extended to an infinite number of excited state terms, a series which generally converges very slowly. Moreover, it also depends on the continuum. For example, if we use all the true (and well known) discrete spectral states of the H atom, its polarizability value is less than 4 bohr 3 , while its exact value is 9∕2 bohr 3 (without taking into account the relativistic effect) [START_REF] Coulson | II.-the Van der Waals force between a proton and a hydrogen atom[END_REF][START_REF] Mcdowell | Exact static dipole polarizabilities for the excited s states of the hydrogen atom[END_REF][START_REF] Traini | Electric polarizability of the hydrogen atom: a sum rule approach[END_REF][START_REF] Bishop | Polarizability and hyperpolarizability of atoms and ions[END_REF]). Fortunately, it is not necessary to know all the excited discrete and continuum states of the electronic system to obtain a "good" static polarizability value, or dynamic polarizabilities for a field frequency smaller than the first resonance one. The coupled-perturbed Hartree-Fock (CPHF) method proposed by [START_REF] Hurst | Ab initio analytic polarizability, first and second hyperpolarizabilities of large conjugated organic molecules: Applications to polyenes C 4 H 6 to C 22 H 24[END_REF], consisting in the independent-particle model using relaxed occupied and virtual orbitals via a self-consistent process in the presence of the external field, generally leads to results in good agreement with the experiments, in particular when a density functional theory (DFT) Hamiltonian with a percentage of the exact or Hartree-Fock (HF) non-local exchange potential is used, as is the case of the hybrid B3LYP [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF]) functional (20% of HF exchange). Indeed, response properties which are n-order perturbation energies depend strongly on the band gap (which determines the lowest transition energy values, n ≥ gap , on the denominator of Eq. 10 for electronic polarizability, and of hyperpolarizabilities), generally too large with HF [START_REF] Evarestov | Ab initio calculation of the ultraviolet-visible (UV-vis) absorption spectrum, electron-loss function, and reflectivity of solids[END_REF]) and too small with pure DFT [START_REF] Yakovkin | The problem of the band gap in LDA calculations[END_REF], but reasonably described by hybrids [START_REF] Garza | Predicting band gaps with hybrid density functionals[END_REF].

Infinite periodic systems

Expression 8 of the Hamiltonian, that includes the position operator, , as the perturbation operator associated to the external field, comes from the dipole moment approximation that can be used for finite systems (molecules) in the presence of a constant field, or a field with a large wavelength with respect to the size of the electronic system ( 𝜆 ≫ 1Å).

For infinite periodic systems, the electric potential ⋅ 0 becomes infinite at (±∞) , where the electronic density is finite, so that the wavefunction is not square summable in Eq. 3. Moreover, infinite periodic systems described by symmetry translated cells and for which the crystalline orbitals (CO) are combinations of Bloch functions following Born von Karman (BVK) conditions (last cell in each direction is bound with the first one) have an ill-defined dipole moment. Actually , which is a non-periodic potential, is not Hermitian in the BVK CO basis set depending of the points of the reciprocal space. Then, we must restart from the expression of an external electric field described by a scalar potential gauge: = 0 e ⋅ , where 0 is the field amplitude and q its wave number 2 ∕ , and determine dipole moment transitions between i( ) and j( � ) COs involved in the field perturbation at the → limit for a constant field. Using the momentum conservation: = � -, it follows that the dipole moment operator becomes [START_REF] Blount | Solid state physics[END_REF][START_REF] Otto | Calculation of the polarizability and hyperpolarizabilities of periodic quasi-one-dimensional systems[END_REF][START_REF] Rérat | Comparison of the polarizability of periodic systems computed by using the length and velocity operators[END_REF]:

and the perturbative Hamiltonian (multiplied by the constant 0 field) is Hermitian and remains block diagonal in the reciprocal space for 𝜆 ≫ 1 Å.

The basic equations of the self-consistent coupled-perturbed (SC-CP) computational scheme for periodic systems can be found in Ferrero et al. (2008a, b) for the CPHF calculation of (hyper)polarizabilities, adapted to Kohn-Sham Hamiltonian (CPKS) in [START_REF] Orlando | Calculation of the first static hyperpolarizability tensor of three-dimensional periodic compounds with a local basis set: A comparison of LDA, PBE, PBE0, B3LYP, and HF results[END_REF] and to frequencydependent electric fields in [START_REF] Evarestov | Ab initio calculation of the ultraviolet-visible (UV-vis) absorption spectrum, electron-loss function, and reflectivity of solids[END_REF], [START_REF] Rérat | Computation of second harmonic generation for crystalline urea and KDP. An ab Initio approach through the coupled perturbed Hartree-Fock/Kohn-Sham scheme[END_REF], [START_REF] Maschio | Calculation of the dynamic first electronic hyperpolarizability (-; 1 , 2 ) of periodic systems. Theory, validation, and application to multi-layer MoS 2[END_REF], whereas many numerical examples are reported in Sect. 3.

Refractive index

In this section, we are going to look at the refractive index of materials, which is a macroscopic property depending on the (hyper)polarizability described in the previous section, as we will see further. Let us look first at the definition of this optical property and then at the resulting birefringence of materials.

Definition

Electric induction ( ) and magnetic excitation ( ) fields are related to the frequency( )-dependent electromagnetic ( , ) field as follows [START_REF] Condon | Theories of optical rotatory power[END_REF]):

where r and r are the relative electric permittivity and magnetic permeability (matrices) of the medium, respectively; 0 and 0 are the permittivity and permeability (constants) of the vacuum, such that 0 0 c 2 = 1 ; and is the

(13) Ωk = + 𝚤 k = -𝚤e -𝚤 ⋅ k e 𝚤 ⋅ , (14) = r 0 - c t (15) = r 0 + c t ,
chirality (matrix) responsible of the rotation angle of a polarized electromagnetic field.

Then, for materials with r ≃ (low magnetic perme- ability) and = (no chirality), the refractive index is (Condon 1937):

Optical indicatrix and birefringence

As seen above, the relative dielectric tensor r is a secondrank symmetric tensor usually represented by a Hermitian matrix. This matches an ellipsoid [START_REF] Nye | Physical properties of crystals[END_REF], the optical indicatrix, whose equation is:

and with semi-axis lengths given by the square roots of the dielectric tensor eigenvalues, n i = √ i (i = X, Y, Z) , corre- sponding to the principal refractive indices of the medium. We will assume that the indices along the semi-axis X, Y, Z, are ordered by increasing value.

Hence, the phenomenon of birefringence is estimated as the difference = n Zn X . In principle, all crystals are birefringent and the specific indicatrix properties depend on the crystal symmetry. However, some special directions-the optical axes-exist which select as many circular sections within the indicatrix. For cubic minerals, which are optically isotropic, the optical indicatrix (see Fig. 1) is a sphere (null birefringence) defined by a unique refractive index n. The uniaxial optical indicatrix of tetragonal, hexagonal, or trigonal minerals is an ellipsoid of revolution characterized by two independent semi-axis of length n and n . The axis of revolution of the indicatrix parallels the c-direction of the mineral and corresponds to n . Such ellipsoid possesses a single circular section perpendicular to the c-direction which corresponds to the optical axis. The birefringence is

= |n -n | . If n 𝜖 > n 𝜔 ( n = n Z , n = n X = n Y ) the indicatrix and the cor- responding mineral are said uniaxial positive. If n 𝜖 < n 𝜔 ( n = n X , n = n Y = n Z )
, the indicatrix and the corre- sponding mineral are said uniaxial negative. Any other mineral is biaxial, having an optical ellipsoid with two circular sections (and as many optical axes) of indices n Y (see Fig. 1). n X and n Z lie on the plane of the optical axes and bisect the angles between them. Depending on whether the acute angle between the optical axes (2V) is bisected by n Z ( 2V = 2V Z ) or by n X ( 2V = 2V X ), the crystal is said to be positive or negative. 2V Z + 2V X = . The following equa- tion defines the relationship between V and n X , n Y and n Z : Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

(16) = � ( r r ) -1 2 ± � -1 ≃ √ r . X 2 n 2 X + Y 2 n 2 Y + Z 2 n 2 Z = 1 327

Non-linear optics and birefringence

Intense electric fields can also provide birefringence as demonstrated in the following. Equation 14 for = can also be written as follows:

where the polarization vector, , is the dipole moment per unit volume induced by the macroscopic field, and can be developed in a Taylor series as follows:

where (n) are the (non)linear electric susceptibility tensors of rank (n + 1) , the units of which are the inverse of field to the (n -1) th power.

Let us consider the field produced by the light wave of frequency together with an external static electric field 0 :

where is the vector amplitude of the wave. Then, a nonzero component value of (2) leads to an additive term proportional to E 0 -modulus in the expression of the polarization

cos V = n X n Y √ √ √ √ n 2 Z -n 2 Y n 2 Z -n 2 X . (17) = r 0 = 0 + , (18) = 𝜖 0 ( (1) + (2) ⊗ + (3) ⊗ ⊗ + ⋯), (19) = 0 + cos t,
vector which oscillates with the same frequency as the electromagnetic field:

and induces a linear variation of the refractive index with respect to the modulus of 0 :

with This Pockels effect responsible of the birefringence: = 1 0 , depends on the odd rank tensor, (2) , which is null for materials with an inversion symmetry. However, the birefringence of centrosymmetric materials can be seen if they own a large third-order non-linear (3) susceptibility (tensor of rank four) and for intense electric fields. Indeed, it comes from the quadratic term in in the polarization vector expression (Eq. 18). In the intensity-dependent refractive index (IDRI) Kerr effect (see [START_REF] Boyd | Nonlinear optics. Chapter 4: the intensity-dependent refractive index[END_REF], an intense beam of light in the crystal can itself provide the modulating electric field: without the need for an external field to be applied. The product of frequency-dependent fields leads again to a polarization vector oscillating with the same frequency as the laser field:

and then, we have: with the light intensity. In that case, the birefringence is given by 2 , i.e., by the IDRI-Kerr non-linear electric (3) (-; ,-, ) susceptibility. A recent application referring to yttria-stabilized zirconia can be found in [START_REF] Marcaud | Third order nonlinear optical susceptibility of crystalline oxide yttria-stabilized zirconia[END_REF]. Several other NLO processes also appear when combining the Taylor development of the polarization vector shown in Eq. 18 with the field expression given in Eq. 19 (see [START_REF] Saleh | Fundamentals of photonics: electrooptics[END_REF][START_REF] Boulon | Nonlinear optics: molecular engineering[END_REF], as SHG linked to (2) (-2 ; , ) and third harmonic generation (THG) linked to (3) (-3 ; , , ) , but they do not affect the refractive index at -frequency.

(20) ( ) = 0 ( (1) (-; ) + 2 (2) (-; ,0) 0 ) cos t, ; ) and 1 ( ) = (2) (-; ,0) ∕ 0 ( ). Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

(21) ( ) = √ + (1) (-; ) + 2 (2) (-; ,0) 0 ≃ 0 ( ) + 1 ( ) 0 0 ( ) = √ + (1) (-
(22) = cos t = 1 2 e t + e -t (23) ( ) = 0 (1) (-; ) + 3 4 (3) (-; ,-, ) 2 cos t, (24) ( ) ≃ 0 ( ) + 3 8n 0 ( ) (3) (-; ,-, ) 2 ≃ 0 ( ) + 2 ( )

Relation between (non)linear susceptibilities and (hyper)polarizabilities

The macroscopic (non)linear optical properties, (n) , in the expression (Eq. 18) of the polarization vector are linked to the microscopic , , … (hyper)polarizability properties of a unit cell of the material, recalling that is the dipole moment per unit cell volume V induced by the mean (macroscopic) field ( ) felt by the cell, and that can be also developed as follows:

The relative electric permittivity (or dielectric matrix) r at zero field is then linked to the polarizability of the unit cell:

where V is the unit cell volume. Similarly, we have:

(2) = 1 2! 1 0 ∕V and (3) = 1 3! 1 0 ∕V , with 1∕ 0 = 4 in atomic units.
It follows that the , , and tensors obtained from the SC-CP calculation discussed in Sect. 2.1.3 for infinite periodic systems allow to determine the (non)linear susceptibilities, (n) , and the refractive index and birefringence.

Applications and comparison with experimental data

The examples should provide an idea of the quality of the results that can be obtained by simulation, and in particular with the CRYSTAL code [START_REF] Dovesi | Crystal17 User's Manual[END_REF][START_REF] Dovesi | Quantum-mechanical condensed matter simulations with CRYSTAL[END_REF][START_REF] Dovesi | The CRYSTAL code, 1976-2020 and beyond, a long story[END_REF]. Before illustrating these examples, it should be mentioned that (obviously) the results depend deeply on the computational parameters; the most important are listed below:

-The variational basis set: in the following examples, a localized Gaussian basis set in split valence or multiple zeta forms, usually including polarization functions, is used. This is in general sufficient to provide accurate results. -The level of the theory: all calculations are performed at the DFT level, and the most reliable, in our opinion, are obtained with hybrid functionals, containing a fraction of the exact Hartree-Fock exchange. However, as the superiority of one functional with respect to the others is frequently a matter of discussion, in some cases we will compute the same property with various functionals.

(

) = 1 V + 1 2! ⊗ + 1 3! ⊗ ⊗ + ⋯ . (26) r = + (1) = + 1 0 ∕V, 25
-The intrinsic numerical accuracy of the implemented algorithms, that in the CRYSTAL code is very high, so that its influence on the final results can be considered negligible.

One additional point must be underlined, when comparing simulation and experiments: the experimental determinations to be compared with are, in many cases, scarce, or affected by large uncertainty, so that error bars on the two sides (simulation and experiment) should be considered. Just to mention an example: the complete determination of all components of a tensor (and many of the ones mentioned above are third or fourth-order tensors, with many components) requires that the experiment is repeated with different orientations of the crystal, possibly with various polarization of the light. The obtained results are often not directly the specific constant to be inserted in the tensor table, but a linear combination of them generating a system of (linear) equations, whose solution (in particular when small and large numbers are involved) can produce strongly correlated final values. The uncertainty is increased in some cases by the fact that, intrinsically, the experiment is unable to determine the signs of the constants.

One advantage of simulation is that all components of the tensors are determined in a single shot, so that the accuracy of the various terms is the same.

The drawback of simulation is that quantum mechanical calculations refer to T = 0 K, so that the effect of tempera- ture can just be guessed, or evaluated a posteriori with simple and, in general, not very accurate tools.

A few more words about anisotropy: each one of the properties listed above, when referred to a gas or a liquid, reduces to the trace (for order two tensors) or to a few invariants (for higher order tensors). In the solid state, on the contrary, all specificity related to orientation is contained in the tensor. Tensors are defined with reference to a cartesian frame (there are rules for defining the orientation of the lattice vectors with respect to the cartesian frame), so that the constants take the following form, for a fourth-order tensor: T ijkl , where i, j, k, l can be x, or y or z. The first manifestation of anisotropy is that some of these components are null for symmetry reasons. Or, viceversa, if the components that should be null are not, some sort of deformation of the lattice with respect to the ideal situation must be taken into account.

The examples refer to the (hyper)polarizabilities of polyacetylene (PA) for which both electronic and ionic contributions are important, the NLO properties of benchmarks as urea ( CH 4 N 2 O ) and potassium di-hydrogen phosphate (KDP, KH 2 PO 4 ), and the birefringence of a series of uni- axal and biaxal minerals. Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

The case of polyacetylene: the importance of the ionic contribution

The electronic (hyper)polarizabilities (with zero-point averaging included) which are calculated from the SC-CP method described above do not take into account the socalled pure vibrational effects, which can be quite important.

In the CRYSTAL code, these ionic effects are taken into account analytically for the polarizability, (Eq. 12), and for the first hyperpolarizability, (Eq. 6 of [START_REF] Rérat | Computation of second harmonic generation for crystalline urea and KDP. An ab Initio approach through the coupled perturbed Hartree-Fock/Kohn-Sham scheme[END_REF], from the IR and Raman spectra, but not for , the second hyperpolarizability (hyper-Raman data are also required; see [START_REF] Champagne | Analysis of the vibrational static and dynamic second hyperpolarizabilities of polyacetylene chains[END_REF].

Also the dynamic e ( ) (or e (-; ) ) and e (-; 1 , 2 ) (but not the corresponding e ) can be obtained from CRYSTAL.

The finite field (FF) scheme was implemented in the CRYSTAL code before SC-CP, to compute numerically the static response properties of molecules and periodic systems (see [START_REF] Darrigan | Implementation of the finite field perturbation method in the crystal program for calculating the dielectric constant of periodic systems[END_REF]. Using then, the finite field nuclear relaxation (FF-NR) scheme of [START_REF] Bishop | A simple method for determining approximate static and dynamic vibrational hyperpolarizabilities[END_REF] that mixes the FF and SC-CP methods, several NLO processes due to the second hyperpolarizabilities can be studied by fitting the Taylor developments of the dipole moment, polarizability, and first hyperpolarizability with respect to the static finite field , at the equilibrium geometry optimized in the presence of the field or not, E and 0 respectively. For example, the Taylor development of the SC-CP (electronic) polarizability is:

where:

The superscript nr indicates the nuclear relaxation approximation for the (field-free) equilibrium vibrational contribution, and the (circular) frequencies of the applied fields are given (as usual) in parentheses, e.g.:

(0;0, 0) = (-; 1 , 2 ) with static applied fields i = 0 and = 1 + 2 .

( Note that the fits on the dipole moment with respect to the static field should lead, in principle, to linear terms identical to the electronic and total static polarizabilities, the latter being the sum of the electronic (Eq. 10) and nuclear relaxation (Eq. 12) contributions to the polarizability: since, in either case, only harmonic vibrational terms are included. However, the static hyperpolarizabilities also contain contributions due to anharmonic force constants and anharmonic electrical property derivatives (see, for example, [START_REF] Torrent-Sucarrat | Basis set and electron correlation effects on initial convergence for vibrational nonlinear optical properties of conjugated organic molecules[END_REF]). To isolate the nuclear relaxation term, one can either subtract the analytically determined electronic term or calculate the difference between numerical values from the Taylor developments with and without geometry optimization in the presence of the static field. When the geometry is not specified, it is 0 ; the omitted frequencies are zero. The subscript → ∞ in Eqs. 29 and 30 refers to the infinite optical frequency (high-frequency or UV-visible frequency) approximation. In addition to harmonic terms, first-order anharmonic contributions are also included for nr (-; , 0, 0)| →∞ with this FF-NR method (see [START_REF] Bishop | A simple method for determining approximate static and dynamic vibrational hyperpolarizabilities[END_REF]; for the other two NLO processes, nr (-; , 0)| →∞ and nr (-2 ; , , 0)| →∞ (the latter being obtained from fits of e ), the first-order anharmonicity terms vanish.

) e tu ( 0 , ) = e tu ( 0 , ) + ∑ v e tuv E v + 1 2 ∑ v,w e tuvw E v E w + ⋯ (28) e tu ( E , ) = e tu ( 0 , ) + ∑ v b tuv E v + 1 2 ∑ v,w g tuvw E v E w + ⋯ , 27
The measured values of non-linear optical properties ordinarily correspond to the sum of vibrational and electronic contributions. In principle, the two may be separated experimentally as well as computationally. For the former, this requires frequency-dependent measurements, as discussed above for , and implied when passing for hyperpolarizabilities as suggested by [START_REF] Shelton | Hyperpolarizability dispersion measured for Kr and Xe[END_REF].

Let us consider now the case of all-trans polyacetylene (PA) [START_REF] Lacivita | Static and dynamic coupled perturbed Hartree-Fock vibrational (hyper)polarizabilities of polyacetylene calculated by the finite field nuclear relaxation method[END_REF], lying in the xy plane, with alternating double and single C-C bonds along the x periodic direction (a double C-C bond is included within each unit cell). The center of the unit cell is an inversion point (which annihilates odd order energy perturbation terms, i.e., and in our case) and lies on a C 2 -axis perpendicular to the xy h mirror plane, which relates z and -z directions, so that all the components of the and tensors containing an odd number of z indices vanish. A finite field is applied along the non-periodic directions (y, z, and mixed yz) to obtain the various non-zero independent tensor components of the vibrational (hyper)polarizabilities. is a fourth-order tensor consisting, in principle, of 3 4 = 81 components, tuvw . Several components are, however, null or equivalent either by point symmetry or permutation of indices, the latter depending on the number of static field indices. For example, in the case of nr tuVW (-; , 0, 0)| →∞ , only the permutations P t,u (associated with the → ∞ Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

limit) and P V,W (between two static fields) leave the prop- erty invariant. According to the number ( m = 4, 2 or 1) of static fields, there are three nuclear relaxation contributions to the second hyperpolarizability:

(a) nr 4,tuvw = nr TUVW (0;0, 0, 0), (b) nr 2,tuvw = nr tuVW (-; , 0, 0)| →∞ and (c) nr 1,tuvw = nr tuvW (-2 ; , , 0)| →∞ .

Seven different finite fields were applied along each direction (y, z and y = z ) in Lacivita et al. (2012), namely | | = 0.1 , 0.5, 1, 2.5, 5, 7.5, 10 × 10 -3 a.u. . Overall, 21 field-dependent geometry optimizations, , followed by the SC-CP1 calculations at the first order of perturbation were used to generate ( ) , e ( ) (left hand side of Eq. 28) and e ( ) . To extract the nuclear relaxation (hyper)polarizabilities from the Taylor expansions of the dipole moment and (hyper) polarizabilities (right-hand side of Eqs. 29 and 30 for the polarizability), an additional set of SC-CP2 calculations at the second-order of perturbation was performed at the fieldfree optimized geometry, 0 . Total and nuclear relaxation contributions to are reported in Table 1.

Let us consider first the transverse (in-plane and non-periodic) yyyy-component of . The vibrational contribution to the static value, nr 4 (0;0, 0, 0) = 1169 a.u. for four static fields, is almost 50% of the total value, g (0;0, 0, 0) = 2419 a.u. , while this percentage decreases to 25% and -2% when only two or one fields are static, nr 2 (-; , 0, 0) and nr 4 (-2 ; , , 0) , respectively. The same comment applies to components including the out-of-plane z-direction and off-diagonal yz indices.

In the longitudinal x-direction of PA, no components of can be calculated from the fit of the dipole moment with respect to static fields, since its x -component is ill-defined. Then, only components of having 1, 2, or 3 x-indices can be obtained from fits of e xy(z) and e xyz , e xx and e xxy(z) , and e xxx with respect to static E y and E z fields. This means that at least one field must be frequency-dependent, leading then to a small nuclear relaxation contribution nr 1 to g with respect to the electronic one. However, the vibrational contribution is larger than the electronic one for the two-static field g non-linear optic property, particularly when both in-plane periodic x and non-periodic y components are involved. In the case of xxyy (-; , 0, 0) , the vibrational contribution, nr 2 , is equal to +2.332 × 10 4 a.u. , while the total value is smaller: g = +9877 a.u. , showing that the electronic contribution has an opposite sign with respect to the vibrational contribution, being equal to -1.334 × 10 4 a.u..

It is worth noting that nr xxxx for the intensity-dependent refractive index (IDRI) process depends on the Raman intensities only, as nr xx depends on the IR intensities, and can be directly obtained as follows [START_REF] Champagne | Analysis of the vibrational static and dynamic second hyperpolarizabilities of polyacetylene chains[END_REF]):

The CRYSTAL code computes the Raman intensities, and then allows to analytically obtain this parallel IDRI vibrational contribution too, which is of the same order of magnitude as the electronic one, around 6 × 10 6 a.u. at the HF level of calculation with a 6-31G basis set including ghost atoms (basis set B in [START_REF] Lacivita | Static and dynamic coupled perturbed Hartree-Fock vibrational (hyper)polarizabilities of polyacetylene calculated by the finite field nuclear relaxation method[END_REF]).

(32 Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

) nr xxxx (-; , -, ) = 2 3N-6 ∑ i=1 d e xx dQ i 2 2 i -2 ,
The static nr xxxx (0;0, 0, 0) term, instead, needs also the determination of the product of IR and hyper-Raman amplitudes, (d x ∕dQ i )(d e xxx ∕dQ i ) , to be evaluated at the equilib- rium geometry. The latter is not yet available in the CRYS-TAL code. Nevertheless, in the polyacetylene case, the IR intensity is small, leading to 𝛼 nr xx ≪ 𝛼 e xx [0.2 and 171.5 a.u., respectively (Lacivita et al. 2012)] and then most probably to a small contribution, such that the very large nr xxxx (0;0, 0, 0) value can be evaluated from Eq. 32 with = 0.

3.2 KDP and urea: dielectric tensor e , second-harmonic generation e , and the effect of the field wavelength

In the previous example of polyacetylene, the subscript → ∞ in Eqs. 29 and 30 means that only the electronic transitions, the energies of which are much larger than the IR mode ones, are involved in the corresponding high-frequency ("infinite frequency") electric field perturbation, but the frequency was set to zero in the calculation. In the present section, the effect of the field wavelength is studied on the electronic contribution to the so-called high-frequency or optical dielectric tensors.

Let us consider the examples of tetragonal KDP and urea. In Table 2, the non-null components of the optical dielectric tensor and SHG susceptibility (electronic contribution only, e and e ) at zero frequency, as obtained at the HF level and with various DFT functionals, are reported, and compared with the experimental determinations at = 1064 nm wavelength. The ionic contribution to electric properties at this wavenumber is negligible, but the corresponding photon energy ( ∼ 1 eV ) is by far smaller than the gap value (and the UV-visible absorption edge), confirming that it can be considered as null in the electronic contribution as we will see further: e+nr =1064nm ∼ e →∞ . is a symmetric second-order tensor; only two components, xx and zz , are independent and non-null for symmetry reasons. As regards the third-rank SHG = 1 2 (2) , only one non-equivalent component (xyz) survives. It should be underlined that part of the reduction of the number of independent terms is due to the intrinsic symmetry of the physical property (for , its symmetric character reduces the constants from 9 to 6; in a similar way, the intrinsic symmetry of SHG reduces, for a triclinic compound, the non-null and non-equivalent terms from 27 to 10). In the last line of Table 2, the band gap is also reported, due to its relevance in determining the numerical values of the various quantities.

Tables 2 and3 permit to discuss two points characterizing the simulation of the reported properties, and of the SC-CP approach:

(a) the effect of the adopted functional; (b) the effect of the self-consistent treatment.

As regards point (a), in the tables, the results obtained with five of the most popular approaches, namely HF, LDA (Perdew and Zunger 1981), PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF], PBE0 [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: the PBE0 model[END_REF], and B3LYP [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF], are reported. The HF Hamiltonian is known to be affected by a universal overestimation of the band gap E g (15.99 vs 7.12 eV, + 125%). This leads to a systematic underestimation of the dielectric properties ( -10% for e xx , -49% for d e xyz ): remember that CPHF and CPKS are SC-CP perturbative schemes in which the gap appears in the denominator in the analytical definition of the optical properties, see, for example, Eq. 10 where n ≥ gap.

At the other extreme, the LDA energy gap is underestimated (5.72 vs 7.12 eV, -20% ), because of the self-interac- tion error; as a consequence, also the (hyper)polarizabilities are generally overestimated (+ 4% for e xx , + 25% for d e xyz

). Gradient corrections (e.g., PBE) provide only small improvements for the gap ( -16% ), and for the dielectric properties (+ 4% for e xx , + 20% for d e xyz

). When the two hybrid functionals, B3LYP (+ 12% for the gap) or PBE0 (+ 20%, about the same error, with opposite sign, of LDA), are used, the difference with respect to experiment is usually smaller than when using LDA or PBE (0% for e xx , -4% for d e xyz for B3LYP, -1% for e xx , Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

-8% for d e xyz for PBE0). Note, however, that if the other component of the dielectric tensor is considered, e zz , the LDA or PBE results are closer to experiment than the ones of hybrids (+ 1%, 0%, -4% , -4% and -12% , from left to right in the table).

Let us consider now the effect of the SCF process, which permits to the system to respond to the electric field perturbation [point (b) above]. The self-consistent coupled-perturbed, SC-CP, results can be compared with the uncoupled SOS (sum over states; iteration 0 of the SC-CP process) data shown in parentheses in any second row of the table.

A few comments concerning e :

-In all cases, the coupled-perturbed scheme improves the SOS results. -The larger the distance from the experiment at the SOS level, the larger the SC-CP correction (SC-CP minus SOS): for HF, e xx varies by 0.32 (from 1.71 to 2.03) and from -24 to -10% ; for LDA, at the opposite side of the table, e xx decreases from 2.56 to 2.34 and from + 14 to + 4%.

It is interesting to notice that for hybrids, and in particular for B3LYP, the difference between SOS and SC-CP is quite small (0.01 and 0.04 for e xx and e zz for B3LYP, and 0.03 and 0.00 for PBE0 for the same components, with always smaller than 4%). In summary, hybrid functionals seem to require a much smaller correction from the coupling than LDA, PBE, and HF.

The above comments apply also to the SHG d e xyz data; the effects are, however, much larger in percentage. Therefore, for HF, increases from -83% to -49% ; for LDA, it decreases from + 66 to + 25%; for B3LYP from -13 to -4% , with a relatively modest change in absolute value from 0.34 to 0.37 pm/V.

The data shown in Table 2 are static, but they refer to electronic calculations or measurements in which a high (UV-visible) field frequency has been used. In Table 3, we explore, for both urea and KDP, the effect of the field wavelength, for the electronic dielectric tensor [START_REF] Rérat | Computation of second harmonic generation for crystalline urea and KDP. An ab Initio approach through the coupled perturbed Hartree-Fock/Kohn-Sham scheme[END_REF]. Computed SHG results for urea and KDP (see Table 4) show that B3LYP reproduces rather well the experimental values measured at 1064 and 600 nm. For urea, the quasi-isotropic electronic contribution is slightly smaller than the experimental value at = 1064 nm [START_REF] Halbout | Efficient phasematched second-harmonic generation and sum-frequency mixing in urea[END_REF]. At this wavelength, the vibrational contribution may not be completely negligible. For that reason, the double harmonic vibrational d nr xyz(zxy) (-2 ; , ) components were calculated at = 1064 nm (see Eq. 6 in [START_REF] Rérat | Computation of second harmonic generation for crystalline urea and KDP. An ab Initio approach through the coupled perturbed Hartree-Fock/Kohn-Sham scheme[END_REF]; their value is 0.041 (0.035) pm V -1 for B3LYP, and has the same sign as the electronic contribution d e xyz(zxy) (-2 ; , ) . The total B3LYP d e+nr xyz(zxy) value is, then, equal to 1.027 and 1.018 pm V -1 which falls essentially at the outer limit of the error bars for the experimental value: Halbout et al. (1979) ( d e+nr xyz ≈ d e+nr zxy for wavelength larger than 600 nm). Vibrational anharmonicity and/or temperature effects, which would increase the magnitude of this term, could be among the reasons of the small discrepancy. At 600 nm, the calculated vibrational contribution is four times smaller than at 1064 nm and, thus, can be neglected. The B3LYP value (1.371 and 1.361 pm V -1 ) in this case is well within the experimental window: [START_REF] Bäuerle | Phasematched second harmonic generation in urea[END_REF]. Finally, for KDP, the B3LYP electronic value of 0.41 pm V -1 is in perfect agreement with the available experimental reference, i.e., 0.41 pm V -1 [START_REF] Singh | CRC handbook of laser science and technology, supplement 2: optical materials[END_REF]). Again, Table 4 clearly shows the well-known tendency of LDA and GGA functionals to grossly overestimate highorder electric susceptibilities as the value of the wavelength approaches the resonance. It is noteworthy that d e xyz increases at each frequency when the percentage of HF exchange Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

d 14 = 1.2 ± 0.1 pm V -1 of
1.3 ± 0.3 pm V -1 of
decreases. This correlates with the predicted band gap for the different Hamiltonians. For urea, E g decreases from HF (14.0 eV) to PBE0 (7.4 eV) and B3LYP (6.9 eV) to PBE (5.2 eV) and LDA (4.8 eV); see Table 2 for KDP.

Figure 2 shows the variation of SHG d e xyz as a function of , the field wavelength, for three levels of theory: HF, LDA, and B3LYP. The dots represent the experimental determination. Its error bar is also indicated. The above discussion concerning the effect becomes here very clear: above 600 (HF), 700 (B3LYP), or 1100 (LDA) nm d e xyz remains essentially constant when the wavelength of the field is varying. As for other properties, B3LYP performs best.

Refractive index n and birefringence ı of minerals

The refractive index and birefringence are directly related by simple equations to the components of the dielectric tensor, as shown in Sect. 2.2.

Here, we will focus on three aspects:

(a) the relative importance of the electronic and ionic contributions. We remind that the former requires a CPHF calculation at fixed geometry (experimental or calculated); the ionic contribution is more expensive, as it requires to build the full Hessian matrix for generating the vibrational frequencies. (b) the dependence of these contributions on the wavelength (c) the anisotropy, that is the different response to the electric field applied in different directions.

As discussed previously, the relative importance of the electronic e (Eq. 10) and ionic nr (Eq. 12) contributions to polarizability depends on the frequency. In the UV-visible region, the electronic contribution is much larger than the ionic one. Therefore, we will first consider results obtained at the sodium source wavelength, D = 589.3 nm , taking into account the electronic contribution alone.

A large set of minerals, namely fluorite ( CaF 2 ), periclase (MgO), corundum ( Al 2 O 3 ), quartz ( SiO 2 ), rutile ( TiO 2 ), anatase ( TiO 2 ), calcite ( CaCO 3 ), aragonite ( CaCO 3 ), anda- lusite, sillimanite, and kyanite (three Al 2 SiO 5 polymorphs), forsterite ( Mg 2 SiO 4 ), topaz ( Al 2 SiO 4 F 2 ), and perovskite ( CaTiO 3 ), plus CaO have been selected to cover a large range of refractive index and birefringence values. Table 5 reports data calculated for this set of systems, using the B3LYP hybrid functional, relatively severe computational (in pm/V) of urea and KDP computed using various Hamiltonians at → ∞ limit and equal to 1064 and 600 nm Data as in [START_REF] Rérat | Computation of second harmonic generation for crystalline urea and KDP. An ab Initio approach through the coupled perturbed Hartree-Fock/Kohn-Sham scheme[END_REF]. Experimental data from [START_REF] Levine | Large local-field effects in the second-harmonic susceptibility of crystalline urea[END_REF], [START_REF] Halbout | Efficient phasematched second-harmonic generation and sum-frequency mixing in urea[END_REF], [START_REF] Bäuerle | Phasematched second harmonic generation in urea[END_REF] for urea and from Refs. [START_REF] Eckardt | Measurement of nonlinear optical coefficients by phase-matched harmonic generation[END_REF], [START_REF] Singh | CRC handbook of laser science and technology, supplement 2: optical materials[END_REF] Fig. 2 Variation of the SHG d e xyz tensor component of bulk urea with respect to the field wavelength, , at the HF, B3LYP, and LDA levels of calculation. Experimental data from [START_REF] Levine | Large local-field effects in the second-harmonic susceptibility of crystalline urea[END_REF]Allan (1993), Halbout et al. (1979), [START_REF] Bäuerle | Phasematched second harmonic generation in urea[END_REF] Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

conditions, and the basis sets developed by [START_REF] Peintinger | Consistent gaussian basis sets of triple-zeta valence with polarization quality for solidstate calculations[END_REF]. In all but two cases, the experimental geometry has been used. The refractive indices obtained at the optimized geometry are extremely close to the ones computed at the experimental geometry, the exceptions being topaz and perovskite, for which both data are reported in the table.

Figure 3 shows the birefringence ( ) as a function of the refringence expressed by the proxy ñ . In general, being the birefringence a small fraction of the refractive indices, the intermediate refractive index is a good approximation of the refringence. Therefore, ñ = n Y , ñ = n , and ñ = n for biaxial, uniaxial, and cubic crystals, respectively. Consider first the cubic systems, on the zero axis. For fluorite and CaO, simulation and experiment coincide, whereas, for periclase, the experimental value is slightly larger (1.735 vs 1.725). This is always the case for all systems: when the two circles do not overlap, the experimental ñ value is always slightly larger than the computed one, the difference being of the order of 2-3%.

Three sets of polymorphs have been considered: (rutile, anatase), (calcite, aragonite), and (andalusite, sillimanite, kyanite). As already observed for the andradite-grossular solid solution [START_REF] Lacivita | Anomalous birefringence in andradite-grossular solid solutions: a quantum-mechanical approach[END_REF], polymorphs with close density (for example, andalusite, sillimanite) present similar indices. When the density is different, the denser system has the largest ñ value, as is the case of anatase-rutile and of andalusite-kyanite, according to the Gladstone-Dale or Drude law [START_REF] Anderson | The relation between refractive index and density of minerals related to the Earth's mantle[END_REF]. In summary, the ñ experimental data are well reproduced, and the residual error is small, and always with the same sign.

The agreement between the experimental and calculated birefringence (the difference between the smallest and largest refractive indices) is better than the one for ñ by about one order of magnitude. This is due to the fact that inaccuracies due to basis set limitations, use of a specific functional, definition of the equilibrium geometry, and numerical approximations are to a large amount the same for the different components of the dielectric tensor, and then cancel when computing birefringence.

The comparison of calculated and experimental optical sign and angle is more delicate. Figure 3 confirms that for uniaxial crystals, the optical sign is correctly predicted, also for quartz or corundum which are characterized by a small birefringence. For corundum, the discrepancy is the largest in the set; note, however, that the experimental birefringence is the smallest in the set. This indicates that the evaluation of the optical sign of weakly birefringent uniaxial crystals is delicate.

For the biaxial crystals, the optical sign and angle are connected. In general, the agreement between calculated and experimental optical signs is good. A closer inspection reveals, however, the difficulty of obtaining the precise shape of the indicatrix, as the interesting case of topaz, orthorhombic, shows. The "module" of the calculated birefringence compares satisfactorily to the experimental one (0.0102) either at the experimental (0.0066) or optimized (0.0078) For each mineral (col. 1: Flu: Fluorite, Per: Periclase, CaO, Cor: corundum, Qtz: quartz, Ru: rutile, Ana: anatase, Cal: calcite, Ara: aragonite, And: andalusite, Sil: sillimanite, Ky: kyanite, Fo: forsterite, Top: topaz, Pv: perovskite) the crystalline system is reported (col. 2 : Cub: cubic, Tet: tetragonal, Trig: trigonal, Rho: rhomboedrical, Ort: orthorhombic, Tric: triclinic). Columns 3-5, 6, and 7 report the experimental refractive indices, the birefringence, and the 2V Z angle, respectively. The corresponding calculated values are given in columns 8-10, 11, and 12, respectively. The sign associated with the birefringence is the optical sign. For uniaxial crystals, n Y = n and n Z = n for the positive ones and n X = n and n Z = n for the negative ones. 2V Z is smaller (larger) than 90 degrees for positive (negative) biaxial compounds. Opening braces indicate sets of polymorphs. Experimental refractive indices from the compilation at [START_REF] Shannon | Refractive index and dispersion of fluorides and oxides[END_REF]; when several data are available, preference has been given to data agreeing with [START_REF] Fleischer | Microscopic determination of the nonopaque minerals[END_REF]. The accuracy on the refractive experimental indices is ±(0.0001 -0.0010) . The variations between different measurements are usually of the order of 0.001. Opening parentheses associate calculations performed at the experimental (first line) and optimized (second line) geometry *See text about birefringence of perovskite Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

geometry. However, the sign of the calculated indicatrix reverses according to the considered geometry.

We mention here also the CaTiO 3 perovskite case, whose crystals are finely twinned. We are not aware of any experimental determination of the refractive indices, as Table 5 shows. We found, however, a single experimental determination of the birefringence [START_REF] El-Mallah | Birefringence of CaTiO 3 and CdTiO 3 single crystals as fonction of temperature[END_REF]. We indicate this situation with an ellipse, rather than with a circle. The value computed at the experimental geometry is quite close to the experimental value (0.019 vs 0.018), whereas the value at the optimized geometry is slightly larger (0.030).

A few words now concern the optical angle 2V, shown in Table 5, for which we consider again topaz and perovskite. In both cases, the calculated refractive indices at the two geometries are quite close, but the optical angles are different, and strongly dependent on these small differences. The reason is that the size of the indicatrix depends on the refractive index values, whereas the shape is mainly controlled by the two partial birefringences: 1 = n Yn X and 2 = n Zn Y . Note that = 1 + 2 (Wright 1951). For 𝛿 < 0.05 and 1.400 < n X < 2.000 , the main contribution to cos 2V Z is given by 2 -1 . The partial birefringence being smaller than the birefringence, small changes of the partial birefringence can change the optical sign. Therefore, at this stage, the calculated sign of the indicatrix of biaxial materials has to be considered as poorly constrained.

As anticipated, at wavelengths shorter than the IR region, the polarizability, and, consequently, the dielectric matrix (Eq. 26) and refractive index (Eq. 16), can be approximated by its electronic part only, since nr ∼ -∑ i 2 i ∕ 2 ∼ 0 -(see Eq. 12 with 𝜔 ≫ 𝜔 i ). However, defining a precise limit for IR is not simple. It depends on the frequency modes ( i ) and corresponding Born charges ( i ) of the considered systems. Practically, 4000 nm (equivalent to 2500 cm -1 ) is a provi- sional limit that should be applied to most of the studied systems. Hence, experimental indices measured at wavelength shorter than about 4000 nm should be larger than refractive indices evaluated at infinity, when considering only the electronic part ( n e →∞ ), since e = ∑ n f n ∕( 2 n -2 ) (see Eq. 10) increases with respect to till the first resonance = n in the UV-visible spectrum ( 𝜆 < 600 nm ) for many minerals.

Then, we have performed a literature search to find experimental dispersion of the refractive index with respect to the photon energy between the IR and UV-visible absorption spectra ( ∈ [200 -10,000] nm ). Available data are represented in Fig. 4. Table 6 reports the largest wavelength ( max ) at which the refractive index of the Fig. 3 Birefringence ( ) versus refringence ( ñ ). See text for the definition of ñ , that depends on the crystalline system. Positive and negative uniaxial or biaxial minerals are reported in the upper (blue circles) and lower (red circles) part of the figure, respectively. Open and closed circles correspond to calculated and experimental values, respectively. The ellipse corresponds to perovskite and is centered on the reported interval of indices and the largest measured birefringence (see text); the optical sign is unknown. Abbreviations are as follows: Flu: fluorite, Per: periclase, Cor: corundum, Ru: rutile, Ana: anatase, Cal: calcite, Ara: aragonite, Fo: forsterite, And: andalusite, Sil: sillimanite, Ky: kyanite, Qtz: quartz, Top: topaz, Pv: perovskite. When circles overlap, in all cases, they refer to the same compound, that is indicated only in one of the circles Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

considered materials have been measured, and the corresponding measured indices ( n max ), as well as n D , corresponding to the yellow line of Na. For quartz, calcite and rutile, max is in the near IR and significantly shorter than 4000 nm. For the first two systems, n max is close to the n e →∞ limit value represented by horizontal lines in Fig. 4. For rutile, whose max is close to the visible spectrum represented by the vertical grey strip, n max is smaller than n e →∞ . The difference n e →∞n max is comparable to its equivalent at D . For fluorite, corundum, and periclase, max largely exceeds 4000 nm and the n max values are much smaller than n e →∞ , due to the negative ionic contribution, see Eq. 12 when 𝜔 > 𝜔 i . The larger n max , the larger n e →∞n max (Fig. 4). The difference reaches almost 10% in the cases of corundum and fluorite and the electronic polarizability alone fails to reproduce indices measured at large wavelengths.

Frequency calculations were performed to evaluate nr ( ) (Eq. 12) and the polarizability was evaluated considering the two contributions. In the case of fluorite, corundum and periclase, including the ionic contribution, yield calculated refractive indices ( n e+nr max ) in very good agreement with the experimental values (Table 6). The difference between calculated and measured indices is now of the order of the difference at D . For rutile, calcite, and quartz, addition of the ionic contribution slightly reduces the calculated indices as expected.

The relative e and e+nr contributions to the refractive index in the various spectral regions can be appreciated from Fig. 4 where experimental dispersion data and calculated n D , n e →∞ and n e+nr max are reported for the materials considered. As anticipated above, the electronic contribution yields a fairly good value of the index n D (large open symbols in the grey vertical strip). However, n e →∞ (horizontal lines) deviate significantly from indices measured at 𝜆 > 4000 nm . In fact, this "limit" depends on materials. From our data, a better provisional "limit" should be 2500 nm. Above 4000 (or 2500) nm, the ionic contribution cannot be neglected and brings the calculated indices ( n e+nr max ) in very good agreement with experimental indices.

Conclusions

In this document, it has been shown that ab initio quantum mechanical simulation can be used for the calculation of a large set of properties related to the perturbation of an electric field on a crystalline system. The simplest ones, like the dielectric tensor, the refractive indices, the birefringence, correspond to the second order terms in the Taylor expansion of the total energy of the system as a function of the field strength. The laser technology permits nowadays to access experimentally also to the third-and fourthorder terms in the expansion. The combination in various ways of the frequencies of the involved fields multiplies the number of phenomena and physical features that can be accessed experimentally. A large set of these (first and second hyperpolarizability, second-harmonic generation, intensity-dependent refractive indices, and many others) can be obtained from the CRYSTAL code used here, at relatively low cost for small-medium-size periodic systems (the most interesting, due to the tensor nature of many of these quantities).

One of the big advantages of simulation is that the full set of constants defining the tensor are obtained with a single calculation, whereas many experiments must be performed for obtaining the same result.

The combination with other perturbations, for example the strain of the unit cell, generates a new set of important →∞ allow to appreciate the role of the ionic contribution to refractive index. For uniaxial crystals, solid and dashed lines refer to n and n , respectively. max is indicated along the vertical line (see Table 6). Color refers to the compound. The grey strip indicates the visible region. The vertical scale of the upper part is half that of the lower part. Notice that no dispersion data have been found for n of corundum Rendiconti Lincei. Scienze Fisiche e Naturali 1 3

properties, like the piezoelectric (third order) or photoelastic (fourth order) tensor, available as well in a very simple way from CRYSTAL.

Quantum mechanical simulation appears then an essential tool for the accurate exploration of many tensorial properties of crystalline compounds.

  is the dipole moment of a polymer (AB) n→∞ , that of AB: -(A-B)-(A-B)-(A-B)-or the one of BA: -A)-(B-A)-(B-A)-(B-?

Fig. 1

 1 Fig. 1 Representation of the indicatrix of a positive biaxial crystals. The directions of the principal axes of the ellipsoid X, Y, and Z are represented by a thick green line, and the associated red thin arrows indicate the corresponding refractive indices n X , n Y , and n Z . The grey surfaces are the circular sections of the indicatrix, and their radius is constant and equals n Y . The normals to the circular sections or optical axis (OA) form the acute 2V Z angle bisected by the Z semi-axis. For uniaxial crystals n X = n Y , usually indicated as n , and n Z is indicated as n ; the two circular sections merge in a single one orthogonal to Z

Fig. 4

 4 Fig.4Refractive index as function of wavelength. Large and small symbols correspond to calculated and measured indices, respectively. Square, circle, downward triangle, upward triangle, diamond, and pentagon correspond to fluorite, corundum, periclase, rutile, calcite, and quartz, respectively. Empty large symbols within and outside the vertical grey strip indicate n D and n e+nr max , respectively. The horizontal lines at height n e →∞ allow to appreciate the role of the ionic contribution to refractive index. For uniaxial crystals, solid and dashed lines refer to n and n , respectively. max is indicated along the vertical line (see Table6). Color refers to the compound. The grey strip indicates the visible region. The vertical scale of the upper part is half that of the lower part. Notice that no dispersion data have been found for n of corundum

Table 1

 1

	FF-NR static and dynamic vibrational (nuclear	Fitted values
	relaxation) contributions to the second hyperpolarizability (in a.u.) of PA obtained by fitting		(a) g	nr 4	(b) g	nr 2	(c) g	nr 1
	(a) the dipole moment, (b) the polarizability and (c) the first hyperpolarizability versus the	yyyy 2419 ± 44 zzzz 2843 ± 283 371 1169	1672 ± 26 2780 ± 16	422 308	1223 ± 8 2514 ± 6	-27 42
	finite field (according to Eq. 28 in the polarizability case)	xxxy -xyyy -	--	-1675 ± 80	-5247	(1.373 ± 0.004) × 10 5 -3614 ± 5	2900 -42
		xxyy -	-	9877 ± 1009	2.332 × 10 4 (-1.358 ± 0.002) × 10 4 -140
		xxzz	-	-	(1.158 ± 0.241) × 10 4 5874	5766 ± 7	60
		yyzz 1541 ± 17	933	861.3 ± 15.7	253.2	658.1 ± 3.1	50
		zzyy	= yyzz	= yyzz 671.7 ± 7.3	63.7	602.6 ± 0.9	-5.5
		yzyz	= yyzz	= yyzz 1855 ± 31	1247	= yyzz	= yyzz
		xyzz	-	-	879.4 ± 202.8	659.5	250.2 ± 1.9	30.3
		xzyz	-	-	2740 ± 86	2520	= xyzz	= xyzz
		HF Hamiltonian and 6-31G(d) basis set including ghost atoms

Table 2

 2

	Coupled-perturbed optical dielectric constants, e xx and e zz , and SHG susceptibility, d e xyz (in pm/V), of tetragonal KDP (space group I 4d2 ) at zero frequency and different levels	e xx e zz	HF 2.025 (1.712) 1.868	PBE0 2.223 (2.195) 2.039	B3LYP 2.230 (2.244) 2.046	PBE 2.328 (2.526) 2.132	LDA 2.340 (2.562) 2.159	Exp. 2.23 a 2.13 a
	of theory		(1.639)	(2.043)	(2.086)	(2.317)	(2.369)
		d e xyz	0.197	0.355	0.373	0.467	0.488	0.38 b
			(0.066)	(0.302)	(0.341)	(0.609)	(0.647)
		E g	15.99	8.51	7.99	6.00	5.72

A split valence basis set was used with d functions on H and f-functions on K, P, and O. SOS (sum over states) values in parentheses. E g is the energy gap (in eV). Calculated data from

[START_REF] Lacivita | Calculation of the dielectric constant and first nonlinear susceptibility (2) of crystalline potassium dihydrogen phosphate by the coupled perturbed Hartree-Fock and coupled perturbed Kohn-Sham schemes as implemented in the CRYSTAL code[END_REF]

. Experimental values from a Polyanskiy (2020) and b Eckardt and Byer (1991) at = 1064 nm

Table 3

 3 High-frequency dielectric tensor components

			(nm)		HF	PBE0	B3LYP	LC-BLYP	PBE	LDA	Exp.
	e xx and e zz of urea (three top lines) and KDP (three bottom lines) computed using various Hamiltonians at → ∞ limit and equal to 1064 and 600 nm	Urea	∞ 1064	e xx e zz e xx e zz	1.901 2.187 1.907 2.195	2.059 2.433 2.070 2.448	2.070 2.451 2.081 2.467	2.057 2.430 2.067 2.446	2.149 2.567 2.163 2.588	2.186 2.599 2.202 2.621	2.194 2.529
			600	e xx	1.920	2.093	2.106	2.091	2.195	2.237	2.220
				e zz	2.214	2.483	2.504	2.481	2.635	2.671	2.577
		KDP	∞	e xx	2.025	2.223	2.230	2.219	2.327	2.341
				e zz	1.868	2.039	2.046	2.049	2.132	2.160
			1064	e xx	2.030	2.231	2.239	2.227	2.339	2.352	2.231
				e zz	1.872	2.046	2.053	2.055	2.141	2.169	2.131
			600	e xx	2.041	2.250	2.258	2.245	2.363	2.377	2.277
				e zz	1.880	2.060	2.068	2.069	2.160	2.189	2.155
		Data from Rérat et al. (2015). Experimental values from Rosker et al. (1985) for urea and Polyanskiy
		(2020) for KDP

Table 4

 4

	Calculated SHG
	high-frequency electric susceptibilities d e xyz and d e zxy

Table 5

 5 Calculated and experimental refractive indices both at ( D ) for various minerals characterized by different symmetry

Table 6

 6 Comparison between calculated and measured refractive indices n D and n max max in nm in the last column. For uniaxial minerals, superscripts and underline that the lines report n and n , respectively

	d Ghosh (1999)
	e Devore (1951)

a Malitson (1963)

b Malitson (1962)

c

Stephens and

[START_REF] Stephens | Index of refraction of magnesium oxide[END_REF]

Acknowledgements

Access to the HPC resources of CINES/IDRIS/ TGCC obtained thanks to the Grants 2018-[A0040810471] (Fabien Pascale) and 2018-[A0050810537] (Philippe D'Arco) made by GENCI are warmly acknowledged. High-Performance Computing resources were partially provided by the EXPLOR centre hosted by the University de Lorraine. This project has received funding from the ANR (Agence Nationale de la Recherche)-FOIST (project number 18-CE24-0030-03). Part of this work was granted access to the HPC resources of [CCRT/CINES/IDRIS] under the allocation 2018-2019 and 2019-2020 [A0040807031] made by GENCI (Grand Equipement National de Calcul Intensif). We also acknowledge the Direction du Numérique de l'Université de Pau et des Pays de l'Adour for the computing facilities provided.

Rendiconti Lincei. Scienze Fisiche e Naturali 1 3