\section{Introduction}
Rational design of energy conversion materials are indispensable to develop cutting edge technologies to exploit the renewable energy resources for future thermal energy management systems. Thermoelectric (TE) technology offers an environmentally benign route for efficient conversion of waste heat into useful electricity. However, low performance TE materials turns out to be a bottleneck for designing high performance TE devices. Performance of a thermoelectric material is determined by a dimensionless quantity, figure of merit, ZT = S$\sigma^2$T/$\kappa$, where S, $\sigma$, and T are the Seebeck coefficient, electrical conductivity, and absolute temperature, respectively. $\kappa$ = $\kappa_e$+$\kappa_l$; is total thermal conductivity which includes contribution from electronic ($\kappa_e$) and lattice ($\kappa_l$) thermal conductivities. Therefore, discovery of advanced functional materials through various design strategies or fine-tuning the existing TE materials properties by an external parameter such as doping, pressure, temperature, strain, etc., to improve the performance of a material. In addition, lone pair can be used as a tool for designing many classes of functional materials\cite{Seshadri2022} in general and low $\kappa_l$ materials in particular because of the lone pair induces strong anharmonicity in the system. On the other hand, (layered) structure causes strong bonding hierarchy $i.e.,$ having rigid and fluctuating sub-lattices in a material has important implications to enhance complex interdependent thermoelectric properties, for instance, high power factor (PF=S$^2\sigma$) and low lattice thermal conductivity ($\kappa_l$).\cite{Pal2018,pandey2020ultralow} 

Bismuth oxychalcogenides are emerging class of materials due to their outstanding thermeoelectric and opto-electronic properties.\cite{Tippireddy2021} In particular, dibismuth dioxychalcogenides, Bi$_2$O$_2$Ch (Ch = S, Se, Te), their 2D counterparts,\cite{Yang2019,Li2021,Liang2021} and heterostructures\cite{Bo2022} gained tremendous research interest due to their potential applications in thermoelectrics,\cite{Luu2015,Wang2021,Tippireddy2021} ferroelectrics,\cite{Wu2017} and optoelectronics.\cite{Yang2022,Tippireddy2021,Li2021,Zhang2015,Guo2019} The ternary Bi$_2$O$_2$Ch systems are synthesized by a solid state reaction of well-known binary thermoelectric bismuth chalcogenides Bi$_2$Ch$_3$ (Ch = S, Se, Te) with Bi$_2$O$_3$. Bi$_2$O$_2$S crystallizes in the primitive orthorhombic structure (SG: $Pnmn$, Z=2 formula units (f.u.) per primitive cell), while Bi$_2$O$_2$(Se/Te) crystallizes in the body centred tetragonal (SG : $I4/mmm$, Z = 2 f.u./cell) anti-ThCr$_2$Si$_2$-type structure by having the Bi$_2$O$_2$ (PbO-type structural units in a square planar unit of oxygen atoms mono-caped with Bi atom) and  alternating layers of Se, which are held together by weak electrostatic forces in contrast to van der Waals (vdW) interactions that are generally observed in layered materials.\cite{Wei2019,Guo2021} The stereo-chemical activity of 6$s^2$ lone-pair of Bi$^{3+}$ cation is the main reason for different crystal symmetry of Bi$_2$O$_2$S over heavy chalcogen based Bi$_2$O$_2$Se and Bi$_2$O$_2$Te systems. Infra-red and/or Raman spectra\cite{Xu2019} and strain effects on bulk and monolayer Bi$_2$O$_2$Ch crystals are systematically investigated to explore their possible applications in nano-electronics.\cite{Cheng2018} Structural stability and topological properties of ternary Bi$_2$X$_2$Y (X, Y = O, S, Se, and Te) compounds are studied using dispersion corrected DFT-D2 method.\cite{Ma2018} A high pressure X-ray diffraction study\cite{Bu2020} revealed that Bi$_2$O$_2$S undergoes a structural transition from $Pnnm$ $\rightarrow$ $I4/mmm$ at 6.4 GPa due to suppression of the Bi$^{3+}$ cation 6$s^2$ lone-pair and then 2D $\rightarrow$ 3D structure up on further compression above 13.2 GPa owing to disappearance the lone-pair. Bi$_2$O$_2$Se is one of the most stable Sillen-type compounds under pressure with no structural transition being observed until 30 GPa\cite{Pereira2018}. Bi$_2$O$_2$Te is found to be stable till 50 GPa \cite{Wang2020-1} and it has promising applications in optical fields as a tuning material due to large birefringence.\cite{Chen2019} The intrinsic ultralow $\kappa_l$ behavior\cite{Wang2018,Song2019,Song2020,Guo2021} and thermoelectric properties\cite{Wang2018,Song2020} of Bi$_2$O$_2$Ch (Ch = S, Se, Te) compounds are studied using first principles calculations and Boltzmann transport theory. Effect of hydrostatic pressure on stereochemically active lone-pair and its implications on lattice dynamics and phonon transport of Bi$_2$O$_2$S is scarce in the literature, which is essential, as it provides a fundamental understanding to enhance the performance of Bi$_2$O$_2$S in thermoelectrics and/or opto-electronic applications.

Here, we systematically investigate the crystal and electronic structure, elastic properties, lattice dynamics, and phonon transport properties for the ambient ($Pnnm$) and high pressure ($I4/mmm$) phases of Bi$_2$O$_2$S. By analyzing pressure-dependent static enthalpy, lattice, bond parameters, and elastic constants, we show that Bi$_2$O$_2$S undergoes a structural phase transition from low symmetry ($Pnmn$) to a high symmetry ($I4/mmm$) structure around 4 GPa.  By analysis of electron localization function (ELF) and electronic structure, we further show that under pressure the lone-pair electrons of Bi$^{3+}$ cation are suppressed.  This suppression of lone-pair electrons reduces anharmonicity at high pressure and therefore $\kappa_l$ is enhanced.


\section{Computational details and Methodology}
The ab-initio calculations have been performed using Vienna Ab-initio Simulation Package (VASP)\cite{Kresse1996}. The electron-electron interactions $i.e.$ exchange-correlation are treated by Generalized Gradient Approximation (GGA) within the parametrization of Perdew-Burke-Ernzerhof (PBE). The electron-ion interactions are treated with pseudopotential approach. We have used DFT-D3 method\cite{Grimme2011} to capture the weak vdWs interactions. The following plane wave basis orbitals were considered as valence electrons; Bi: 5d$^{10}$,6s$^2$,6p$^3$; O: 2s$^2$,2p$^4$; and S: 3s$^2$,3p$^4$. A kinetic energy cutoff of 600 eV was used for plane wave basis set expansion and a spacing of 2$\pi$ $\times$ 0.025 $\AA^{-1}$ for k-mesh in the irreducible Brillouin zone. The plane wave cutoff energy as well as k-spacing were varied first to ensure the total energy convergence. Self consistency criteria is chosen in such way that the total energy is converged to 1e$^{-8}$ eV/atom and the maximal force on each atom is less than 1e$^{-3}$ eV/$\AA$. Electronic structure is calculated using Tran Blaha modified Becke Johnson (TB-mBJ) potential\cite{Tran2009} as implemented in WIEN2k package\cite{Blaha2002} including spin-orbit coupling. Crystal structures and electron localization functions are visualized using the VESTA software\cite{VESTA2008}. 


Harmonic and anharmonic interatomic force constants (IFCs) are needed to calculate $\kappa_l$ from first principle calculations. For both the systems, harmonic IFCs were obtained using the finite displacement method within the Phonopy package~\cite{togo2008first,togo2015first} with 5$\times$5$\times$2 (500 atom) supercell and $\Gamma$-point only k-grid. To obtain well converged phonon frequencies a high plane wave cutoff of 600 eV and strict energy convergence criteria of 1e$^{-8}$ were used. Long range corrections to phonon frequencies were included by calculation of Born effective charges and dielectric constants within the method proposed by Wang et. al.~\cite{wang2010mixed}. Anharmonic IFCs were calculated using 3$\times$3$\times$2 supercell. For the displaced supercells configurations phonon-phonon interactions were considered out to the 11$^{th}$ nearest neighbors. This requires 972 and 500 DFT calculations for $Pnnm$ and $I4/mmm$ phases, respectively. To reduce the computational cost here, only $\Gamma$-point k-mesh along with plane wave energy cutoff of 520 eV and energy convergence criteria of 1e$^{-8}$ eV was used for Brillouin zone integration. The calculated harmonic and anharmonic force constants were used to solve the phonon Boltzmann transport equation iteratively as implemented in the ShengBTE code~\cite{li2012thermal1,li2012thermal,li2014shengbte}. A Gaussian smearing width of 0.1 was used to approximate the Dirac  distribution, and a q-point mesh of 21$\times$21$\times$9 was used for Brillozin zone integration. In these calculations three-phonon scattering and the isotopic disorder scattering~\cite{tamura1983isotope} from natural isotope mass variation are considered. 

\section{Results and Discussion}

\subsection{Crystal structure and structural phase transition}
As discussed in section I, at ambient conditions, Bi$_2$O$_2$S crystallizes in lower symmetry ($Pnmn$) structure (see Figures \ref{fig:str}(a-b)) than heavy chalcogen based systems namely Bi$_2$O$_2$Se and Bi$_2$O$_2$Te. This is due to the presence of stereochemically active $6s^2$ lone-pair of Bi$^{3+}$ cation in Bi$_2$O$_2$S and its activity might be suppressed by chemical pre-compression of heavy chalcogens (Se/Te) in Bi$_2$O$_2$Se and Bi$_2$O$_2$Te systems, which drives them to crystallize in relatively high symmetry ($I4/mmm$) structure (see Figures \ref{fig:str}(c-d)) and the $Pnmn$ structure is distorted analogue of $I4/mmm$ crystal symmetry. To investigate the bonding differences between these two phases further, we examine the ELF~\cite{becke1990simple, seshadri2001visualizing} as depicted in Figures~\ref{fig:str}(e)-(g). The ELF is projected onto the [110] plane, which contains atoms of Bi and S. The red color represents complete electron localization, while the blue color represents almost no localization. The ELF of the two phases appear similar at first glance, but significant differences can be seen around Bi-S bonds. In the $Pnmn$ phase, the Bi lone-pair can be seen by yellow regions of electron localization in Figure~\ref{fig:str}(e). When pressure is applied (Figures ~\ref{fig:str}(f)-(g)), the intensity of the yellow lobes surrounding the Bi atoms is significantly reduced, indicating that the lone-pair of Bi$^{3+}$ cation is suppressed. Due to presence of lone-pair in $Pnmn$ phase, the coordination environment around Bi$^{3+}$ cation differs greatly between these two phases. In $Pnmn$ phase, Bi and S have two in-equivalent bond lengths and these are listed as Bi-S1 (3.00 $\AA$) and Bi-S2 (3.51 $\AA$) in Table \ref{table2} and Figure~\ref{fig:str}(a). While in the $I4/mmm$ phase, Bi and S bonds have four-fold coordination with only one in-equivalent bond with a length of 3.14 $\AA$. Similarly, the bond length between Bi and O atoms in the $Pnmn$ phase ranges from 2.24 to 2.43 $\AA$, whereas in the $I4/mmm$ phase, all Bi-O bonds are identical, with a bond length of 2.19 $\AA$. Overall, the bonding is highly heterogeneous in the $Pnmn$ phase over $I4/mmm$ phase. 

A recent HP-XRD study\cite{Bu2020} revealed that Bi$_2$O$_2$S undergoes a structural phase transition from orthorhombic ($Pnmn$) to tetragonal ($I4/mmm$) structure at 6.4 GPa due to suppression of the $6s^2$ lone-pair and then it transforms from 2D to 3D structure upon further compression above 13.2 GPa. To get further insights on structural phase transition and anharmonicity, we have systematically investigated the high pressure behavior of Bi$_2$O$_2$S up to 10 GPa. As shown in Figure~\ref{fig:lattice}(a), the enthalpy difference between ambient and high pressure phases is gradually decreasing with pressure up to 3 GPa and then after both phases exhibit iso-enthalpic nature above 4 GPa. This result suggests a continuous phase transition from $Pnmn$ $\rightarrow$ $I4/mmm$ in Bi$_2$O$_2$S. To get further insights on the structural phase transition, we have analyzed the pressure dependent lattice constants, bond parameters as a function of pressure and equation of state (EOS). The calculated lattice constants are decreasing monotonically with pressure.  Interestingly, the lattice constants 'a' and 'b' are merging into a single lattice constant above 4 GPa and lattice constant 'c' continues to be 'c' for high pressure tetragonal ($I4/mmm$) phase (see Figure \ref{fig:lattice}(b)) and this result is consistent with the variation of lattice constants as a function of pressure as observed in the HP-XRD measurements.\cite{Bu2020} The calculated pressure coefficients by fitting the pressure-lattice parameter data to a quadratic expression show that the $Pnmn$ lattice is highly compressible along a-axis.  As shown in Figure \ref{fig:lattice}(c), the calculated volumes for both $Pnmn$ and $I4/mmm$ phases are converging and become almost equal above 3 GPa until the maximum studied pressure range (0-10 GPa) in this work and it is in excellent agreement with the HP-XRD measurements.\cite{Bu2020}  As illustrated in Figures \ref{fig:BL}(a-b), the calculated in-equivalent bond lengths (Bi-S $\&$ Bi-O) and angles (S-Bi-S $\&$ O-Bi-O) of lower symmetry ambient $Pnmn$ phase are converging and merging into a single in-equivalent bond length/angle for higher symmetry $I4/mmm$ phase above 4 GPa. The static enthalpy calculations and pressure evolution of structural properties of ambient and high pressure phases clearly demonstrate that Bi$_2$O$_2$S undergoes a continuous phase transition from $Pnmn$ $\rightarrow$ $I4/mmm$ under high pressure. 


\subsection{Cross-band-gap hybridization and enhanced Born effective charges}
It is well-known fact that the standard DFT functionals underestimate the band gaps\cite{Heyd2005,Verma2017} by 30-40 $\%$ for semiconductors and insulators due to self-interaction, to overcome this problem, we have used TB-mBJ\cite{Tran2009} potential to calculate electronic structure of Bi$_2$O$_2$S. The calculated band gaps are 0.82 and 0.67 eV without and with inclusion of spin-orbit coupling (SOC) within PBE-GGA, respectively for $Pnmn$ phase and the corresponding band gaps using TB-mBJ potential are 1.52 and 1.33 eV. The obtained PBE-GGA band gap 0.82 eV without SOC is consistent with previously obtained band gap of 1.0 eV,\cite{Zhang2015} and the TB-mBJ band gap with SOC 1.33 eV is consistent with HSE06 band gap of 1.25 eV.\cite{Wu2017} As shown by Shi et al.~\cite{Shi2016}, materials with fully occupied lone-pair ($ns^2$, n = 4, 5, 6) show distinct features in their electronic structure by having an additional band (dominated by fully occupied cation $s$-states) below the valence band compared to the materials without lone-pair cation. The conduction band is mainly dominated by extended $p$-states of lone-pair cation in contrast to the cation $s$-states in materials without lone-pair cation. The calculated projected electronic density of states (PDOS) for both ambient orthorhombic and high pressure tetragonal (at 5 GPa) phases of Bi$_2$O$_2$S are presented in Figures \ref{fig:DOS}(a) and (b), respectively. As discussed above, the bottom of the valence band is mainly derived from $s$-states of Bi$^{3+}$ cation for both of these phases ($Pnmn$ and $I4/mmm$). As shown in Figure \ref{fig:DOS}, the top of the valence band is mainly dominated by the $p$-states of S$^{2-}$ and O$^{2-}$ anions and also minor contribution from $s$ $\&$ $p$-states of the Bi$^{3+}$ cation, while the conduction band is mainly dominated by extended $p$-states of the Bi$^{3+}$ cation, which hybridize with anion (S$^{2-}$ and O$^{2-}$) $p$-states of the valence band, thus, results in a significant cross-band-gap hybridization in Bi$_2$O$_2$S. The highly dispersive nature of $p$-states of S$^{2-}$ and O$^{2-}$ anions and Bi$^{3+}$ cation at the top of the valence band and minimum of the conduction band, respectively provides high mobility for light-induced e-h separation and their transport in photocatalysis.\cite{Zhang2015} Moreover, the electronic structure also gives a fundamental insight on stereochemical activity of a lone-pair.\cite{Waghmare2003} The cation $s$-states are broader (see Figure \ref{fig:DOS}(a)) and are inter-mediated by unoccupied cation $p$-states with anion $p$-states. When the activity of lone-pair is diminishing, the cation $s$-states (below the valence band) are narrow and have sharp peaks in the PDOS (see Figure \ref{fig:DOS}(b)). The observed cross-band-gap hybridization leads to a mixed ionic-covalent character, which causes enhanced lattice polarization in lone-pair cation containing compounds.\cite{Du2010,Du2010InX} This could be clearly seen from Born effective charges (BECs), as presented in Table \ref{Born}, the calculated BECs are significantly enhanced for both phases of Bi$_2$O$_2$S compared to the formal charges of cation (+3) and anions (-2). Enhanced BECs bring the lattice in the proximity of ferroelectric instability\cite{Du2010} thereby lowering the lattice thermal conductivity ($\kappa_l$)~\cite{pandey2020lattice}. Therefore, it is very intriguing to investigate elastic and phonon transport properties of lone-pair containing materials in general and in particular for Bi$_2$O$_2$S, especially to understand what is the effect of lone-pair on lattice dynamics and phonon transport at high pressure? 

\subsection{Elastic properties, lattice stability and phonon transport}
To explore the mechanical and dynamical stability, we calculated second order elastic constants (SOEC) and phonon dispersion curves, respectively for $Pnmn$ and $I4/mmm$ phases of Bi$_2$O$_2$S at ambient as well as at high pressure. The calculated SOEC at 0 and 5 GPa for both the phases are given in Table \ref{elastic}. At 0 GPa, the $Pnmn$ phase obeys necessary and sufficient conditions of Born stability criteria\cite{Mouhat2014} and no imaginary frequencies are observed, which indicate that the $Pnmn$ phase is mechanically and dynamically stable at ambient pressure (see Figure \ref{fig:PHDISP}(a)). The $I4/mmm$ phase doesn't satisfy the Born stability criteria (C$_{44}$ $\textgreater$ 0) and also it possesses imaginary frequencies along $\Gamma$-M direction (see Figure \ref{fig:PHDISP}(b)) at 0 GPa, which clearly demonstrate that $I4/mmm$ phase is mechanically and dynamically unstable at ambient pressure for Bi$_2$O$_2$S.

When a non-zero hydrostatic pressure is applied to a particular crystal system, the Born stability criteria must be revised with modified elastic constants of a crystal system at a given hydrostatic pressure (P) conditions. According to Sinko and Smirnov\cite{Sinko2002,Sinko2004} the modified elastic constants at a given pressure are $\widetilde{C}_{ii}$ = $C_{ii}$ - P, (for $i$ = 1-6) and $\widetilde{C}_{12}$ = $C_{12}$ + P,  $\widetilde{C}_{13}$ = $C_{13}$ + P and $\widetilde{C}_{23}$ = $C_{23}$ + P. Therefore, the modified Born stability criteria under hydrostatic pressure for orthorhombic crystal symmetry can be found in our previous work,\cite{Yedu2015} whereas for tetragonal system the same are given as follows:
\begin{equation}
\begin{split}
C_{11} - C_{12} - 2P  >  0, \\ P^2 + P (C_{11} + C_{12} + 4C_{13}) + 2C^2_{13} - C_{33}(C_{11} + C_{12}) > 0, \\  C_{44} - P  >  0, \hspace{0.2in} C_{66} - P > 0
\end{split}
\end{equation}
Under high pressure ($\textgreater$ 4 GPa), for instance at 5 GPa (see Table \ref{elastic}), the calculated SOEC for both $Pnmn$ and $I4/mmm$ phases have the following relationships: C$^{Pnmn}_{11}$ $\approx$ C$^{Pnmn}_{22}$ $\approx$ C$^{I4/mmm}_{11}$, C$^{Pnmn}_{44}$ = C$^{Pnmn}_{55}$ = C$^{I4/mmm}_{44}$ and C$^{Pnmn}_{12}$ $\approx$ C$^{Pnmn}_{13}$ = C$^{I4/mmm}_{12}$. These relationships strongly suggest a continuous phase transition from $Pnmn$ to $I4/mmm$ is occurring in Bi$_2$O$_2$S under high pressure and is in good accord with the pressure dependent structural properties such as pressure dependent lattice constants and bond parameters (see Figures \ref{fig:lattice} and \ref{fig:BL}). In addition, no imaginary frequencies are found for the computed phonon dispersion curves at 5 GPa, which strongly suggest that $I4/mmm$ phase becomes dynamically stable at 5 GPa.   


As discussed above, under the application of pressure the lone-pair of Bi$^{3+}$ cation is suppressed, which can have important consequences on phonon transport. To investigate this further we analyze the lattice dynamics and phonon transport in both $Pnmn$ and $I4/mmm$ phases of Bi$_2$O$_2$S. The calculated phonon dispersion and phonon density of states (PHDOS) are shown in Figures~\ref{fig:PHDISP}(a) and (b) for $Pnmn$ and $I4/mmm$ phases, respectively. In general, for both $Pnmn$ and $I4/mmm$ phases, the three acoustic modes are comprised of very low frequencies ($\leq$ 1.5 THz) and they overlap with optical branches. In both phases, the optical phonons show very low frequency modes due to presence of heavy metal Bi and they hybridize with the acoustic phonon modes. The optical phonon branches in high-frequency regions ($>$ 4THz) are very dispersive and exhibit large phonon group velocity, indicating the possibility of contribution to $\kappa_l$. The PHDOS in these compounds can be divided into three main regions: Bi dominated region below 1 THz, mixed Bi and S region above 1 THz, and O dominated region above 6 THz. Despite the huge mass contrast between Bi and S (m$_{Bi}$/m$_S$ = 6.5), there is a significant overlap between Bi and S in the PHDOS. This overlap is larger for the $Pnmn$ phase (1 THz - 5.5 THz) than in the $I4/mmm$ phase (3 THz - 5.5 THz). The sulfur PHDOS in the $Pnmn$ phase shows several small peaks in 1 THz - 5.5 THz and a large peak around 5 THz. The $I4/mmm$ phase also shows sulfur peaks in the PHDOS, but these are less pronounced, reflecting more dispersive optical branches. The large overlap between acoustic and optical phonons can lead to enhanced scattering of the heat-carrying acoustic phonons --- suggesting the possibility of relatively lower $\kappa_l$ for the $Pnmn$ phase.

The temperature-dependent $\kappa_l$ of the two phases are compared in Figure~\ref{fig:Kappa}(a) along in-plane and cross-plane directions. As seen, $\kappa_l$ of the ambient $Pnmn$ phase is much lower than the high pressure $I4/mmm$ phase. At room temperature, the average $\kappa_l$ values are 1.47 and 3.63 W/m-K for $Pnmn$ and $I4/mmm$ phase, respectively. Since both $Pnmn$ and $I4/mmm$ phases have an anisotropic crystal structure, the same is reflected in the calculated $\kappa_l$. Due to the orthorhombic structure, the $Pnmn$ phase exhibits different $\kappa_l$ along in-plane x and y directions and are become comparable at higher temperatures. As can be seen from Figure~\ref{fig:Kappa}(a) both phases show anisotropy in $\kappa_l$ with $\kappa^{\parallel}_l$/$\kappa^{\perp}_l$ ratio as 2.14 and 2.58 in $Pnmn$ and $I4/mmm$ phase, respectively. For instance, the calculated thermal anisotropy 2.14 of ambient phase of Bi$_2$O$_2$S is comparable with thermal anisotropy 2.0 for Bi$_2$O$_2$Se in $I4/mmm$ symmetry.\cite{Guo2021} As expected for layered materials with weak electrostatic interactions,\cite{Guo2021} the $\kappa_l$ along the cross-plane direction is smaller than the in-plane direction. Interestingly the calculated thermal transport anisotropy is much smaller than what is typically observed in vdW layered materials such as, MoS$_2$ where $\kappa^{\parallel}_l$/$\kappa^{\perp}_l$  of nearly 30 was found~\cite{meng2019thermal}. This clearly distinguishes between the role of interlayer weak electrostatic and vdW bonding in determining the thermal anisotrpy in layered materials. 


To understand the thermal transport behavior further we calculate the cumulative $\kappa_l$ ($\kappa_{ac}$) as a function of phonon frequency, which is shown in Figure~\ref{fig:Kappa}(b) and (c) for $Pnmn$ and $I4/mmm$ phases.  As can be seen there, in both the phases, the optical phonons contribute significantly to $\kappa_l$ particularly along the in-plane direction. The phonon frequency range contributing to cross-plane and in-plane $\kappa_l$ is also quite different. The cross-plane $\kappa_l$ is mostly contributed phonons with frequency $<$ 5 THz, whereas phonons with frequencies up to 12 THz contribute to in-plane $\kappa_l$ in both the phases. This can be understood by analyzing phonon group velocities which are shown in the Figure~\ref{fig:Kappa}(d) and (e) for $Pnmn$ and $I4/mmm$ phases, respectively. In both the phases for the majority of the frequency range $v_g$ is higher along in-plane direction than the cross-plane direction. This is more prominent for optical phonons in the frequency window of 6-12 THz, where in-plane $v_g$ is more than two times larger than the cross-plane $v_g$. These higher in-plane $v_g$ are the main reason for higher $\kappa_l$ along the in-plane direction. For the $Pnmn$ phase, the $v_g$ along in-plane x and y directions are comparable which is consistent with the comparable $\kappa_l$ along these directions. 

We then focus on exploring the origin of the relatively low $\kappa_l$ in the $Pnmn$ phase. To explain the calculated differences in $\kappa_l$, we look into three factors, phonon group velocity, specific heat, and three phonon scattering rates, that determine $\kappa_l$. As explained above the phonon group velocities are comparable in both phases --- therefore phonon velocities are not the reason for the lower $\kappa_l$ in the $Pnmn$ phase. As for the second-factor specific heat, the $Pnmn$ phase has only slightly higher ($\sim$ 1.3 \%) specific heat than the $I4/mmm$ phase. This indicates that the difference in phonon anharmonicity should be the main reason for different $\kappa_l$. The phonon anharmonicity governs  the strength of phonon-phonon interactions which determines the overall intrinsic thermal resistance. Greater anharmonicity results in shorter phonon lifetimes and lower $\kappa_l$. Indeed, as shown in Figure~\ref{fig:Kappa}(f), the three phonon scattering rates of the $Pnmn$ phase are significantly larger than the $I4/mmm$ phase particularly in the frequency range 6-15 THz. The optical phonons in $Pnmn$ phase are more dispersive than in $I4/mmm$ phase, which enables it to satisfy the conservation laws required for three phonon scattering resulting in higher scattering phase space and therefore higher scattering rates. 

As introduced before the main difference between orthorhombic $Pnmn$ and tetragonal $I4/mmm$ phases, is the absence of Bi 6$s^2$ lone-pair electrons in the $I4/mmm$ phase, which may explain its relatively less anharmonicity.  To gain further insights into the role of lone-pairs electrons in bonding and anharmonicity, we next calculate the potential energy by off-centering the Bi, O, and S atoms along Cartesian x and z directions. The resulting potential energy plots are presented in Figures~\ref{fig:PES}(a) and (b) for $Pnmn$ and $I4/mmm$ phases, respectively. As can be seen there, the displacement potential energy of Bi and S atoms are significantly different along x and z directions in both $Pnmn$ and $I4/mmm$ phases. We find that for Bi and S atomic displacements in the x or y directions, the potential energy curve is relatively flat, whereas displacement along z-direction requires more energy. In both phases the energy required for O displacement is comparable. The potential energy analysis stipulates that Bi and S vibrational modes could be easily excited by low energy. This is consistent with the calculated PHDOS, where despite the huge mass difference, PHDOS of Bi and S atoms show significant overlap.  Furthermore, for Bi vibration in $Pnmn$ phase potential energy curve shows asymmetric features and deviates from parabolic (harmonic) behavior. The polynomial fit of this potential energy well gives anharmonic coefficient as 2.99 and 1.68 along x and z directions, respectively. Such large deviation from typical parabolic behavior corroborates the strong anharmonicity of this phase. For comparison,  in $I4/mmm$ phase the anharmonic coefficients for Bi atomic vibrations along z direction is 1.56, and for Bi displacement along x direction the potential is harmonic (coefficient of anharmonic fit is zero). The above displacement potential energy analysis confirms that the $Pnmn$ phase is more anharmonic than $I4/mmm$ phase.  This is further confirmed by plotting the mean square displacements (MSD) for Bi and S atoms which are shown in Figure~\ref{fig:PES}(c) for both $Pnmn$ and $I4/mmm$ phases. The MSD for Bi and S atoms are higher in $Pnmn$ phase  (Bi = 0.012 $\AA^2$, and S= 0.024 $\AA^2$ at 300 K), than in the $I4/mmm$ phase (Bi = 0.009 $\AA^2$, and S= 0.011 $\AA^2$ at 300 K). This is expected because the longer Bi-S2 bond lengths in $Pnmn$ phase make it easier to move the atoms from their equilibrium positions during the vibration than in the $I4/mmm$ phase. Overall, the anharmonicity of the $6s^2$ lone-pair of Bi$^{3+}$ cation plays a pivotal role in determining the phonon transport in Bi$_2$O$_2$S and it is suppressed (in $I4/mmm$ phase) upon compression as can be seen from chemical pre-compression of Se/Te atoms in Bi$_2$O$_2$Se and Bi$_2$O$_2$Te compounds. 

\section{Conclusions}
In summary, we have systematically investigated the effect of hydrostatic pressure on stereo-chemically active 6$s^2$ lone-pair containing Bi$^{3+}$ cation with weak electrostatic interlayer bonding compound, Bi$_2$O$_2$S and their implications on crystal structure, lattice dynamics and phonon transport properties using first principles calculations and Boltzmann transport theory. We predict that Bi$_2$O$_2$S undergoes a continuous structural phase transition from low symmetry distorted ($Pnmn$) structure to a high symmetry ordered ($I4/mmm$) structure at around 4 GPa, which is comparable to the transition pressure 6.4 GPa determined from HP-XRD measurements\cite{Bu2020}. The mechanical and dynamical stability of both $Pnmn$ and $I4/mmm$ phases are confirmed from the calculated elastic constants and phonon dispersion curves, respectively. The obtained enhanced Born effective charges (BECs) clearly demonstrate that both the phases show significant cross-band-gap hybridization from their electronic structure suggesting near ferroelectric instability and is favorable to achieve low $\kappa_l$. The obtained low $\kappa_l$ value 1.47 W-m/K for $Pnmn$ phase over 3.63 W-m/K for $I4/mmm$ phase is due to suppression of lone-pair of Bi$^{3+}$ cation at high pressure which reduces anharmonicity in $I4/mmm$ phase. The highly anharmonic behavior of Bi and S atoms of $Pnmn$ phase demonstrated from PES and MSD at ambient pressure is mainly responsible for low $\kappa_l$. Overall, the present work provides an insight on how hydrostatic pressure effects the phonon transport in  stereo-chemically active 6$s^2$ lone-pair containing Bi$_2$O$_2$S. The present study could stimulate exploring pressure dependent phonon transport properties in materials especially with lone-pair cation, which will provide a new avenue of pressure modulated low $\kappa_l$ materials for future thermal energy applications.


\section{Acknowledgments}
NYK and TP contributed equally to this manuscript. NYK would like to thank Institute for Advanced Computational Science, Stony Brook University for providing computational resources (Seawulf cluster). TP is supported by Research Foundation Flanders (FWO-Vl). SCRR would like thank RGUKT Basar for providing computational facilities. 


