Materials with more than two switchable ferroic properties are of great interest due to potential use in four state memory devices. Materials that exhibit spontaneous ordering of magnetization, electric polarization and elastic strain are called multiferroics [1] [2] [3] [4] [5] [6] [7] . These materials find several technological applications such as refractory electrodes, thermistors, and thermoelectric materials. The coupling of the ferroelectric and ferromagnetic order parameters is particularly attractive if found in the same phase, both from fundamental science and technology point of view [8] . These compounds find use in spintronics and as data storage devices. It is interesting to note that most of the multiferroics, currently of interest are perovskite (ABO3) structured. These structures have a capacity to accommodate multitude of structural distortions and are able to incorporate almost every element. Ideal perovskite structure is a framework of corner sharing octrahedra hosting B cations, while A cations are placed in the resulting dodecahedral sites. Structural instabilities of these biferroic oxides (ABO3) depend on the A cation size (with same B cation) and change in B cation. The effects of external parameters like temperature, pressure, chemical composition drive these distortions, which lead to several rich and unique physical properties. These physical properties are derived from tilt of BO6 octahedra, polar cation displacement leading to ferroelectricity, etc. Multiferroic materials with coupled Ferromagnetic (FM) and Ferroelectric (FE) order parameters are promising for developing new generation of electrically and magnetically controlled multifunctional devices. A large number of studies on chromate based perovskites have been performed in recent years due to its rich physics and simple structures [7, [9] [10] [11] [12] .
There are two types of multiferroics: i) The materials that exhibit magnetic ordering and ferroelectric ordering up to a very high temperature, but the ordering temperature as well as origin of ferroic properties are very different, hence coupling is very weak; ii) On the other hand, there are compounds with the same ordering temperature and origin of ferroelectric and magnetic ordering with strong coupling, but unfortunately the transition temperature as well as the magnitude of electric dipole moments are very small [13, 14] . Hence, across the globe intense research is going on to achieve both the requirement of strong coupling as well as high transition temperature.
The compounds RCrO3 (R=Y, La) are magnetic but do not exhibit ferroelectricity in pure single phase.
However, mixing them together results in enhanced ferroelectric properties [15] [16] [17] [18] [19] [20] [21] [22] [23] . The mixing leads to strain, which could cause the strain mediated ferroelectricity in the mixed phase. Thermal expansion could also act as a source of strain in the mixed compound. Hence, it is important to understand the thermal expansion behavior and its anisotropy that would be useful to design the multiferroic material in future.
RCrO3 (R=Y, La) are found in orthorhombic structure (space group Pnma) at ambient condition ( Fig 1) . With respect to the ideal cubic perovskite structure (space group Pm-3m), this orthorhombic phase is obtained by an anti-phase tilt of adjacent CrO6 octahedra; the R cations is usually in 8 coordination (Fig 1) . Although Pbnm is a centrosymmetric space-group and non-compatible to ferroelectricity, but local non-centrosymmetric nano regions give rise to ferroelectricity in YCrO3 [18, [23] [24] [25] [26] . Previous work on YCrO3 indicated [27] that it is ferrimagnet with weak ferromagnetism (TN~140 K). Subsequently it has been shown that YCrO3 exhibits [28] canted antiferromagnetic behaviour below 140 K and a ferroelectric transition around 473 K. More recent X-ray and neutron powder diffraction, and Mössbauer spectroscopy measurements on YCrO3 also indicate [29] that it is antiferromagnetic and the direction of the moment on Cr 3+ ion is along the c-axis. LaCrO3 is also reported [29, 30] to be G-type antiferromagnetic, and the direction of the magnetic moment is found along [30] the crystallographic aaxis.
Ferroelectricity is absent in LaCrO3 due to the large size of La 3+ ion which might prevent some structural instability necessary for ferroelectric transition (in comparison to Y 3+ ). In order to understand the structural instabilities in these oxides, several theoretical and experimental studies have been carried out. Understanding the driving force behind ferroelectric instability in oxides like LaCrO3 has warranted several electronic structure studies, and experimental Raman studies [8, [19] [20] [21] [22] 31, 32] . Understanding the correlation between individual phonon modes and structural distortions, make interesting ongoing studies.
Temperature dependence of Raman measurements [33] on YCrO3 show that the A1g and B2g Raman active modes about ~ 560 cm -1 undergo significant changes across the magnetic transition temperature. Therefore, these modes are likely to show strong coupling with magnetic ordering. There have been studies to understand the changes in the physical properties of YCrO3 with doping; it is found that its electrical conductivity increased with Ca-doping. High pressure synchrotron powder diffraction experiments have been carried on YCrO3 up to 60 GPa to study its evolution with increasing pressure [34] .
High resolution neutron diffraction studies have been reported to understand the average and local structure in YCrO3 in order to explain the electric polarization seen above 430 K despite the centrosymmetric phase [18] .
In this paper, we report neutron inelastic scattering measurements over a wide range of momentum transfers, which enables to separately identify the contributions from magnetic excitations (at low momentum transfers) and phonons (at high momentum transfers). The neutron measurements have been performed as a function of temperature, which reveal significant changes across the magnetic phase transition; however, the magnetic scattering persists in the paramagnetic phase due to paramagnetic relaxations. We have also performed extensive ab-initio calculations of the phonon spectra and identified the effect of magnetism. Further, we use the ab-initio calculations to derive the anisotropic elasticity and the thermal expansion, including the negative thermal expansion along the a-axis of the orthorhombic structure. The results are in very good agreement with the available data. The theory helps to identify the anharmonic phonon modes behind the anomalous thermal expansion behaviour.
Polycrystalline samples of LaCrO3 and YCrO3 were prepared by solid state reaction method. The extending up to 7 Å -1 . About 10 grams of polycrystalline sample has been used for the measurements over 100-550 K. The polycrystalline sample was put inside a cylindrical niobium sample holder and mounted in a cryoloop. The data analysis was performed by averaging the data collected over the angular range of scattering using ILL software tools [35] to get neutron cross-section weighted phonon densities of states. The phonon density of states g (n)( E) in the incoherent [36, 37] one-phonon approximation is extracted from the measured dynamical structure factor S(Q,E) as follows for exchange correlation following the parameterization by Perdew, Burke and Ernzerhof. [41, 42] . The plane wave pseudo-potential with a plane wave kinetic energy cutoff of 900 eV was adopted. The integration over the Brillouin zone is sampled using a k-point grid of 2×2×2, generated automatically using the Monkhorst-Pack method [43] . The above parameters were found to be enough to obtain a total energy convergence of less than 0.1 meV for the fully relaxed (lattice constants & atomic positions) geometries. The total energy is minimized with respect to structural parameters. The Hellman-Feynman forces are calculated by the finite displacement method. The total energies and force calculations are performed for the 17 distinct atomic configurations resulting from symmetrical displacements of the inequivalent atoms along the three Cartesian directions (±x, ±y and ±z). The convergence criteria for the total energy and ionic forces were set to 10 -8 eV and 10 -5 eV A -1 , respectively. The phonon energies (and the dispersion curves and density of states) were extracted from subsequent calculations using the phonopy-1.14 software [44] . The phonon calculation has been done considering the crystal acoustic sum rule. The calculated structures (Table I) of LaCrO3 and YCrO3 agree well with the available experimental data.
The thermodynamic and transport properties of material with magnetic ions also contributed significantly by magnetic interaction between these ions. These additional magnetic interactions between controls translational and rotational flexibility of atoms in the crystal, hence, control the dynamics of the material. In addition to that the associated spins of the magnetic ions form collective spin excitation known as magnons in magnetic ordered phase. While in the paramagnetic phase they do contribute in the free energy in the form of spin configurational entropy. Hence to understand the effect of magnetic interaction on thermodynamics properties, it is important to include the magnetic interaction in the system. The standard DFT formalism allow the inclusion of magnetic interaction in two ways, the collinear spin polarized magnetic calculation where the spin is considered as a one-dimensional variable and the Hamiltonian in rotationally invariant with respect to lattice. The other way to include the magnetic interaction is known as noncollinear spin polarized calculation, this is very computationally expensive method and the spin is treated as a three-dimensional object.
Both the compounds YCrO3 and LaCrO3 are reported [29, 30] to be G-type antiferromagnetic. The structure relaxation for both the compounds was performed in the non-magnetic as well as magnetic configurations. For LaCrO3, on-site Hubbard correction is applied within the Dudarev approach [45] using Ueff = U − J = 7.12 eV [46] . The same value of Ueff has been used in the calculations of YCrO3. For YCrO3, we have performed the phonon calculation considering the G-type magnetic configuration and moment on Cr atoms along a-axis (non-collinear calculation) as well along c-axis (collinear calculation).
As discussed below, the calculations show that the change in magnetic moment direction does not result in significant changes in the overall phonon density of states. However, we have seen that for some zone centre modes their energy changes significantly (~5%) due to change in the magnetic moment direction.
For LaCrO3, we have performed only the collinear calculation. The non-collinear calculation does not converge for required accuracy for phonon calculation, hence we have not performed the non-collinear phonon calculation for LaCrO3.
Thermal expansion behavior of any material is an important thermodynamic property essential for material design for specific application. Here, we have computed the linear thermal expansion behavior in both the compounds using quasi-harmonic approximation that has been found useful in many similar compounds. The thermal expansion calculation is done using pressure dependence of phonon frequencies in the entire Brillouin zone [47] . Each phonon mode of energy Eqj (j th phonon mode at wavevector q in the Brillouin zone) contributes to the thermal expansion coefficient, which is given by the following relation for an orthorhombic system [48] :
Where sij are elements of the elastic compliance matrix, s=C -1 at 0 K (C is the elastic constant matrix), V0 is volume at 0 K and CV(q, j, T) is the specific heat at constant volume due to j th phonon mode at point q in the Brillouin zone. The mode Grüneisen parameter of phonon energy Eq,j is given as [49] ,
The volume thermal expansion coefficient for an orthorhombic system is given by:
The elastic constants (Cij) and elastic compliance matrix elements (sij) are given in Table III . The elastic compliances show negative values of shear components, which implies that elongation along one axis leads to contraction to the conjugate axis.
As described above, the distribution of phonon energies averaged over Brillouin zone does not show significant change with magnetic configuration; hence, we do expect that the thermodynamic quantities, which are functional of phonon density of states, might not change significantly with magnetic orientation.
A. Temperature Dependence of Phonon Spectra Therefore, two Q-domains were considered; i.e., high-Q (4 to 7 Å -1 ) and low-Q (1 to 4 Å -1 ) in order to extract the magnetic contribution in the INS data at low-Q and the phonon contribution at high-Q.
The temperature dependence of the Bose-factor corrected S(Q, E) plots of YCrO3 and LaCrO3 are shown in Fig 2. At low temperatures (up to 315 K), the low-Q data show a larger elastic line as compared to the high-Q data. We could see that the energy spectrum close to the elastic line is also very sensitive to Q, i.e. at lower Q (1 to 4 Å -1 ), there is larger intensity in the low energy spectrum and at higher momentum transfer (Q~4-7 Å -1 ), it reduces significantly. As noted above, the magnetic signal is more pronounced at low Q, and vanishes at higher Q, following the magnetic form factor. We speculate that this quasi-elastic scattering originates from paramagnetic spin fluctuations.
The low-Q inelastic spectrum in YCrO3 shows significant intensity at energy around 20 meV at that there is significant contribution by magnon excitations at about 20 meV. Further, the intensity around zero energy is contributed by paramagnetic scattering of neutrons by Cr ions. The inelastic spectrum averaged over high Q does not show any prominent change with temperature. However, in case of LaCrO3, we do not observe any significant difference between the low-and high-Q data at 175 K (TN~291 K) (Fig 2) , and so, the contribution from the magnetic excitations seems to be insignificant.
As mentioned above, the neutron inelastic scattering spectrum from powder sample contains rich In the inelastic experiment, we measure the dynamical structure factor S(Q, E), which contains the information of both the magnon and phonon density of states. Hence, the dynamical structure factor contributed from the magnon dynamics would also exhibit distinct von-Hove singularities as described above. To probe these features in our inelastic neutron spectrum, we have further processed our data to extract the magnon structure factor Smag(Q,E) by subtracting the high-Q inelastic neutron spectrum from the low-Q inelastic spectrum. We obtained the Smag(Q,E) for YCrO3 ( Fig. 3(a) ). The magnon dynamical structure factor consists of one von Hove singularity at around ~20 meV energy, which is a strong signature of the G-type AFM structure. This result supplements the previous diffraction observation and justifies our calculation with the G-type AFM structure.
Here, our aim was to describe the richness of the information contained in the INS spectrum, which can be further utilized to characterize the material properties. In order to confirm that the peak observed at around 20 meV in Smag(Q,E) has magnetic origin and hence should follow the magnetic form factor of associated magnetic ion, we have plotted ( Fig. 3(b) ). This Figure gives the integrated intensity in the Smag(Q,E) over 17-22 meV as a function of momentum transfer from 1-7 Å -1 . We have also calculated the Q dependence of magnetic form factor of Cr +3 using analytical method [52] and compared it with the variation of the integrated peak intensity (Smag(Q,E)) with Q. We found a very good agreement between the two, which unambiguously confirms that the peak ~20meV in Smag(Q, E) is contributed by magnetic excitation (magnons).
Further, the spectra derived from S(Q,E) data using Eq. (2), within the incoherent approximation in both the Q regimes, are shown in Fig. 4 . We find that for both the compound the low-Q data show large variation in the intensity as a function of temperature. However, for the high-Q data we do not observe any significant change in the spectra except at the lowest temperature. Further, in the low-Q data of YCrO3 at 100 K, which is below the magnetic transition temperature (~140 K), there is a large intensity of the low energy inelastic spectra (~20 meV) as compared to the data collected at higher temperatures.
This is expected to be due to a strong magnetic signal. At 175 K, it is found that there is considerable decrease of the intensity of the low energy peaks around 20 meV indicating loss of the magnetic signal.
We have calculated the relaxed structure for both the compounds in the non-magnetic as well as magnetic configurations, and the results are given in Table I . Without the magnetism the a-lattice constant is significantly overestimated by about 5 %. We find that the magnetism substantially reduces the overestimate to about 2 %. Further, in order to see the effect of the direction of magnetic moment on the structure we performed both collinear (moment along c-axis) and non-collinear calculations (moment along a-axis) for YCrO3. We find that the relaxed lattice parameters are quite similar in the two cases and differ from the experiments by about 2% (Table I) , which is acceptable in standard DFT calculation. The calculated magnetic moment of Cr ion in both the compounds is found to be 3.0 µB, which is due to the fact that Cr atom occurs in +3 ionic states in the compounds and all the remaining d 3 electron occupy the t3g orbital. It may be noted that above the magnetic transition temperature the compounds become paramagnetic, where the magnetic ordering is lost which results in marginally reduced effective interaction between magnetic ions. Hence, in order to reproduce the phonon spectrum in paramagnetic phase we cannot use nonmagnetic calculation (magnetic moment of Cr= 0 µB, which is equivalent to magnetic quenching). The calculation of phonon spectrum in paramagnetic phase can be obtained by calculating the phonon spectrum over a large set of randomly oriented spin configurations which is not a feasible option but may be in near future one can use machine learning methods to achieve the same.. It has been reported that the Raman modes [33] in YCrO3 around ~ 560 cm -1 (~70 meV) show significant changes across the magnetic transition temperature, and therefore are expected to be influenced by magnetic interactions, which is consistent with the above ab-initio results.
The phonon density of states of LaCrO3 and YCrO3 differ significantly (Fig 5) below 50 meV, which are dominated by Y and La dynamics. The differences are mainly attributed to difference in various bond lengths and strengths in both the compounds. However, the high energy phonon modes above 50 meV, which mainly arise from Cr-O stretching, do not change much since in the both the compound the Cr-O bond length is nearly same.
We compare in Fig. 6 the inelastic neutron scattering measurements at T=310 K with phonon calculations that includes the magnetic interaction (spin polarised phonon calculation). The phonon density of states obtained from measurements is neutron weighted. Hence, in order to compare with the data we have done the neutron weighting of the calculated phonon data using equation (2) . The calculated phonon density of states for YCrO3 shows good agreement with the measurements. We observe slightly underestimated phonon energies with respect to measurements. That could be understood better in terms of overestimated lattice parameters or under-binding effect of GGA exchange correlation function in DFT approach. In case of LaCrO3, the calculated phonon density of states is found to be in good agreement with the neutron data below 55 meV, while there is significant underestimation of high energy phonon density of states. The significant underestimation of high energy modes might be attributed to strong correlation effect of La-f electrons, which is not very well accounted in DFT methods. However, qualitatively overall spectrum matches very well with the measurements.
Further, we have shown (Fig. 6 ) the contribution of individual atoms to the neutron inelastic scattering spectrum. We found that for both the compounds the contribution from oxygen atoms is dominated across the entire spectral range. The Cr and R (=Y, La) atomic dynamics are significantly contributed in the spectrum below 70 meV and 30 meV respectively.
The phonon density of states is an integrated spectrum of phonons in the entire Brillouin zone. To give a sense of quantitative changes in specific phonon modes due to collinear and non-collinear magnetic configuration in YCrO3, we have shown the calculated zone-center phonon frequencies in Fig 7 and Table   II in the G-type magnetic configuration with magnetic moments along the a-and c-axis. One may note that a few zone-center phonon modes show significant change with magnetic configuration, indicating that these modes exhibit strong spin-phonon coupling. We find that most of the zone-center modes in YCrO3 above 45 meV show significantly change in phonon energies (~ 5%) with magnetic ordering along the a-and c-axis. This also reflects strongly in the density of states calculations shown in Fig. 5 . Further, we have shown the calculated zone-center phonon modes of LaCrO3 in Table II . We may observe that most of the zone-center phonon modes in LaCrO3 are lower in energy than that in YCrO3, which might be due to larger La mass and weaker bonding between La and O. Experimental data of the zone-center modes for YCrO3 and LaCrO3 available from Refs [26, 53] and [54, 55] respectively, are in fair agreement with our ab-initio calculations.
In Fig. 8 , we have shown the eigenvectors of some of the representative modes in YCrO3, which show significant change in phonon energy with magnetic configurations. The most striking feature of these modes is that these modes are mainly dominated by oxygen dynamics. The dynamics of oxygen will cause variations in the O-Cr bond and O-Cr-O bond angle, which are critical for antiferromagnetic interaction. Hence, different types of magnetic configurations result in difference in strength of the force constants for these bonds and bond angles, which in turn reflect in changes in the phonon frequencies. As the stretching and bending modes in these compounds occur at high energies, we expect that significant changes in the high energy Raman and IR modes would be observed in presence of strong magnetic field.
Neutron powder diffraction has been used to study the high temperature behaviour of RCrO3 (R = Y, La). An orthorhombic to a rhombohedral structural phase transition has been reported [56] in LaCrO3 at about 533 K, while in YCrO3 measurements show [57] no transition upto 1200 K. The calculated anisotropic Grüneisen parameters averaged over Brillouin zone are shown in Fig. 9 , 
One may note that if the Grüneisen parameters were isotropic ( = = ) in these compounds, values; hence, in Eq. (12) , all the three terms lead to the NTE along the a-axis. We find that in Eqs. (13) and ( The calculated linear thermal expansion coefficients are used to calculate the lattice parameters and volume as a function of temperature. In Fig 11, we have compared the calculated values with the available experimental data [56, 57] on RCrO3 (R=Y, La). The experimental measurements are not available in the NTE region at low temperature in YCrO3. The experimental observation of high temperature thermal expansion behavior of YCrO3 is well reproduced (Fig 11) with our calculation. As noted above, the overall phonon spectrum in YCrO3 is not much sensitive to the orientation of the magnetic moments. Here also we find that calculated thermal expansion behavior in YCrO3 using both the magnetic configurations (Fig 10 and 11 To investigate the specific anharmonic phonon energies responsible for the NTE along the a-axis in these compounds, we have calculated the contribution to the linear thermal expansion coefficients ( Fig   12) from phonons of various energies averaged over the Brillouin zone at 300 K in YCrO3 and LaCrO3 respectively. We find that phonons of energy below 20 meV (Fig 9) are largely behind the NTE along the a-axis. Below 20 meV, the calculated phonon density of states of LaCrO3 is significantly larger in comparison to that of YCrO3, resulting in larger contribution to the NTE in LaCrO3 along the a-axis.
It is important to mention here that a large number of phonons in the entire Brillouin zone contribute [47, 58, 59] to the thermal expansion behaviour. However, to understand the anisotropic expansion behavior and large NTE along the a-axis in LaCrO3, we have identified a few representative low-energy phonon modes which contribute negative expansion along the a-axis while positive expansion in other directions. We have selected the modes at the zone-center and zone-boundary since they are easier to visualize compared to the modes at general wave-vectors in the middle of the Brillouin zone.
These selected modes have been plotted in Fig 13 and involve displacements of the oxygen and La atoms.
The La displacements are larger in the modes of 13.9, 15.6 and 15.7 meV than that in the 19.7 meV mode, which is consistent with the partial density of states (Fig. 5) . The oxygen movements mainly distort the CrO6 octahedra, along with small overall librations. Among these modes, the largest contraction along the a-axis is produced by the mode at 13.9 meV at (0 1/2 0), which involves large displacements of La and O atoms along the b-axis. We need to caution that this may be an oversimplified picture, but it is illustrative of the NTE behaviour.
In this article, we have presented the ab-initio lattice dynamics calculations and inelastic neutron scattering measurements of YCrO3 and LaCrO3. The scattering with low momentum transfer (Q=1-4 Ǻ - [29, 30] in the orthorhombic phase (space group Pnma) of RCrO3 (R=Y, La). Non-collinear and Collinear correspond to the ab-initio calculations performed with magnetic moment on Cr atoms along a-and c-axis respectively. The experimental values of the magnetic moment for YCrO3 and LaCrO3 are from Ref [29] and [50] respectively. The low-Q and high-Q, unity-normalized, excitation spectra, g (n) (E), inferred from the neutron energy gain S(Q,E) data, within the incoherent approximation. [56, 57] on RCrO3 (R=Y, La) (solid circles) are from literature. Calculations for LaCrO3 are shown in the orthorhombic phase [56] , which is stable up to about 533 K.