Unlocking surface octahedral tilt in two-dimensional Ruddlesden-Popper perovskites

Molecularly soft organic-inorganic hybrid perovskites are susceptible to dynamic instabilities of the lattice called octahedral tilt, which directly impacts their carrier transport and exciton-phonon coupling. Although the structural phase transitions associated with octahedral tilt has been extensively studied in 3D hybrid halide perovskites, its impact in hybrid 2D perovskites is not well understood. Here, we used scanning tunneling microscopy (STM) to directly visualize surface octahedral tilt in freshly exfoliated 2D Ruddlesden-Popper perovskites (RPPs) across the homologous series, whereby the steric hindrance imposed by long organic cations is unlocked by exfoliation. The experimentally determined octahedral tilts from n = 1 to n = 4 RPPs from STM images are found to agree very well with out-of-plane surface octahedral tilts predicted by density functional theory calculations. The surface-enhanced octahedral tilt is correlated to excitonic redshift observed in photoluminescence (PL), and it enhances inversion asymmetry normal to the direction of quantum well and promotes Rashba spin splitting for n > 1.

O rganic-inorganic hybrid halide perovskites are prone to dynamic instabilities called octahedral tilt, which involves a rigid rotation of the inorganic octahedral cages, and can appear around any of the three Cartesian directions in the crystal with either in-phase or out-of-phase ordering 1,2 . Thermally induced octahedral tilt accompanies phase change in perovskites, and has a pronounced influence on electronic band structure including the band gap and band edge effective masses [3][4][5][6][7][8] . In conventional inorganic perovskites (with the general formula of ABX 3 , A is a monovalent cation, B is a divalent cation, and X is a halide/oxide anion), it is commonly accepted that the octahedral distortion is dominated by steric effect of the A site ion occupying the interstitial central cation site of BX 6 octahedrons, in which the unmatched A and B ionic sizes induce octahedral tilt due to steric repulsion 9,10 . In hybrid organic-inorganic 3D perovskites, however, hydrogen bonding interactions between short organic A cation and the inorganic frame plays a central role in octahedral tilt [11][12][13] . For example, out-of-phase octahedra tilt occurs in the prototype MAPbI 3 (MA = CH 3 NH 3 + ) hybrid perovskite, while FAPbI 3 (FA = (NH 2 ) 2 CH + ) hybrid perovskite shows almost no octahedral tilt due to the different hydrogen bonding interactions 14 . Octahedral tilt in 3D hybrid perovskites impacts the Jahn-Teller distortion of inorganic cages thus affecting the band gap 15,16 . In addition, the structural distortions are often coupled to vibrations, affecting the Raman mode frequencies and atomic motions related to the inorganic framework 17,18 . Octahedral tilt also induces exciton-phonon coupling and generates gap states which causes a broadening of the photoluminescence spectrum [19][20][21][22][23][24][25] .
Although Goldschmidt tolerance factors can semiquantitatively predict structural distortion on the basis of steric hindrance in 3D perovskites, it is not straightforward to apply to 2D perovskites, because the layered segregated structures in 2D perovskites allow the incorporation of larger organic cations into the inorganic cages [26][27][28][29][30] . The apparent relaxation of Goldschmidt tolerance factor originates from the mitigation of strain accumulation and self-adjustable strain-balancing in structures of 2D perovskites. This has enabled a significantly expanded library of 2D RP lead halide perovskites with larger and more exotic A cations, leading to unique properties and novel applications [31][32][33][34][35][36][37][38][39] . Kepenakian et al. performed DFT calculations to understand how RPPs relieve mechanical strain at lattice-mismatched interfaces, and found that the internal elastic energy due to lattice strain is relieved by octahedral tilt of the PbX 6 octahedra, which is more pronounced for RPPs above a certain critical thickness 40 .
X-ray diffraction (XRD) had been widely used to investigate octahedral tilt or lattice distortion in both 2D and 3D perovskites (see the summary of relevant works listed in Supplementary  Table 1 and Table 2). However, XRD provides bulk averaged information with no surface sensitivity, thus it cannot be used to probe the presence of surface reconstruction that occurs on exfoliated RPPs. Here, we directly imaged octahedral tilts on the surface of 2D RPPs (BA) 2 (MA) n−1 Pb n I 3n+1 (n = 1-4; BA: CH 3 (CH 2 ) 3 NH 3 , n refers to the number of inorganic slabs per unit cell of the homologous series) using scanning tunneling microscopy. We found that the octahedral tilt on the outermost layer of RPP was enhanced significantly compared to that of the bulk after "unlocking" the two interlocked layers of BA chains by mechanical exfoliation. In particular, the influence of MA ions on the octahedral tilt was investigated from n = 1 (without MA ions) to n = 4 (n > 1 has MA ions) RPPs, from which we identified a clear trend between n and the strength of the octahedral tilt. Interestingly, the out-of-plane octahedral tilt of the exfoliated surface is correlated to the redshifted emission peak alongside the primary exciton in the photoluminescence spectra, and also contributes to Rashba spin splitting.

Results
Unlocking surface octahedral tilt. The mechanical exfoliation of van der Waals stacked 2D material (e.g., Graphene or transition metal dichalcogenide) usually does not induce large surface structural reorganization with respect to the bulk structure because of the lack of covalent bonds in the out-of-plane directions. However, the case is different for organic-inorganic hybrid 2D perovskites owing to its molecularly soft nature. As shown in Fig. 1a and b, delaminating one of the two interlocked layers of BA chains removes the steric constraints for the remaining, thus causing it to tilt at a larger angle relative to its bulk position. To gain an insight into the surface structure relaxation, we conducted DFT simulations on n = 4 RPP to model both exfoliated monolayer surface and bulk interface (Fig. 1c). The surface of relaxed monolayer structure has only one layer of BA on top, thus it is free from the steric constraints of interlocked BA bilayer and readily undergoes relaxation. The tilt angle (Ω) of the surface BA molecule with respect to the c-axis increases from 14.4°of the bulk, unrelaxed structure to 27.1°for the relaxed structure as illustrated in Table 1. This surface-enhanced octahedral tilt consists of a rigid rotation of the anion cage which can be specifically divided into the out-of-plane tilt (versus c-axis) and in-plane tilt (in ab plane). As shown in Fig. 1d, the Pb 1 -I E -Pb 2 bond angle (θ) and I E1 -I E2 -I E3 bond angle (γ) are used to demonstrate out-of-plane and in-plane tilt, respectively. The degree of octahedral tilt is judged from the tilt angles of Δθ (180°−θ) and Δγ (|90°−γ|). In the out-of-plane tilt, the tilt of two octahedral cages toward each other narrows the θ angle by 5.2°from bulk interface to relaxed structure (Table 1). In the in-plane tilt, a stronger in-plane tilt degree is present in the monolayer surface (Δγ = 3.2°) compared to the bulk interface (Δγ = 1.6°). The corresponding bond length changes are summarized in Supplementary Figs. 1-4 and Tables 3-6.
Atomic-scale STM measurement verifying DFT simulations. The visual surface structure of RPPs was imaged using STM measurement (see the detailed sample preparation process in the Supplementary Discussion and Supplementary Fig. 5). Large-sized, pure phase single crystals of (CH 3 (CH 2 ) 3 NH 3 ) 2 (CH 3 NH 3 ) n−1 Pb n I 3n+1 of the homologous series n = 1-4 were synthesized, and thin flakes exfoliated from the bulk crystals were used for STM studies. According to the proposed model of exfoliated surface in monolayer n = 4 RPP, the stronger out-of-plane tilt (decreasing θ angle) of the two adjacent Pb-I cages toward each other will naturally bring their apical I atoms (large yellow spheres) closer, forming a "dimer-like" structure ( Fig. 2a). This is indeed verified by STM realspace atomic-scale imaging of the surface structure of the exfoliated RPP flakes and cleaved CH 3 NH 3 PbI 3 single crystals 41,42 .  Supplementary  Fig. 6, clearly verifying the orthogonal symmetry. Each unit cell contains two atoms (bright protrusions) that form a "dimer-like" structure, the distance of two nearest atoms is 4.58 Å, showing a very good agreement with the DFT relaxed model of monolayer structure of n = 4 (4.55 Å). The charge density plot is simulated as shown in Fig. 2c, most of the charge density is contributed by the terminal I atoms, thus the charge density plot reproduced the "dimer" structure under +2.3 V bias voltage.
The different orbitals of the 2D RPPs could be imaged by STM measurement by applying different bias voltages, as illustrated in Fig. 2d with the direction of electron tunneling marked by arrows.
The energies of different atomic orbitals are different with respect to the Fermi level, thus they are accessed at different bias voltages (see the calculated local density of states in Supplementary Fig. 7). When applied positive bias is larger than +2.0 V (electrons tunneling from tip to sample), apical I atoms are imaged (Fig. 2e). Reducing the voltage to +1.9 V allows the Pb atoms of n = 4 RPP to be imaged (Fig. 2f), where a nearly square lattice with dimensions of 5.89 Å was observed, consistent with the Pb-Pb bond length in the calculation model (6.09 Å). A negative bias of −2.3 V (electrons tunneling from sample to tip) is required to image the organic BA cations, as shown in Fig. 2g. Applying bias during the STM measurement did not distort the lattice, this has been verified by bias-variable STM measurements and STM simulations shown in the Supplementary Figs. 8 and 9. The experimental STM images are consistent with simulated STM images at different bias voltages. The different topologies of STM images at +1.9 V and +2.3 V are attributed to contributions from Pb and I orbitals, respectively.
Layer-dependent octahedral tilt DFT studies. To investigate the surface octahedral tilt as a function of n in 2D RPPs, DFT calculations were carried out on both monolayer and bulk RPPs from n = 1 to n = 4. The simulated structures include fully relaxed ground state (GS) of orthorhombic structures, and a "high symmetry" (HS) counterpart where all octahedrons are fixed and un-tilted ( Supplementary  Fig. 10). The degree of absolute octahedral tilt is determined from the tilt angles (Δθ n and Δγ n , n is the thickness of octahedron) with respect to the reference HS structure, the result is summarized in Fig. 3a Fig. 3c, with the shortest I-I distance in n = 4 and the longest in n = 1 RPP. Interestingly, the apical I-I distance in n = 2 is slightly shorter than that in n = 3, thus the I-I distance does not change monotonically, but depends on the interplay of various Coulombic interactions with other ions (e.g., methylammonium ions) in the structure. We should note that besides the absolute tilt angle measured with respect to a non-tilted imaginary reference HS structure, it is more important to look at the relative tilt angle (marked as blue rectangles in Fig. 3a and b) because in the relaxed bulk system, it already has an intrinsic octahedral tilt relative to un-tilted HS reference. Here, the relative tilt angle denotes the difference value between exfoliated monolayer and non-exfoliated bulk system. In the following section, we correlate the experimentally observed trends to the relative tilt.
STM. The experimental atomic arrangements of the apical I atoms from n = 1 to n = 4 ( Fig. 3d-g) exfoliated RPPs were imaged by STM. All the bright protrusions are due to apical I atoms from n = 1 to n = 4 RPPs, which has been verified by the DFT calculation of LDOS and corresponding STM image simulations ( Supplementary Figs. 11 and 12). To our delight, the STM results verified the DFT-predicted trend of apical I-I distance of monolayer RPP structures from n = 1 to n = 4. STM shows that n = 4 RPP (Fig. 3g) has the shortest apical I-I distance of 4.55 Å, while n = 2 RPP has shorter I-I distance of 4.85 Å (Fig. 3e) than n = 3 RPP of 5.02 Å (Fig. 3f). Meanwhile, n = 1 RPP (Fig. 3d) has the largest I-I distance of 5.69 Å, which fits well with the simulated model that almost no out-of-plane tilt occurred in the surface octahedron of n = 1 RPP. These results clearly validate that the out-of-plane surface octahedral tilt is directly correlated to the I-I distance, while in-plane octahedral tilt has less correlations with it.
Excitonic redshift in the exfoliated RPP crystals. A single narrow photoluminescence (PL) peak due to the recombination of free exciton characterizes the emission of single-crystalline RPP crystals of n = 1 to 4 homologs series 43 . However, the PL spectrum of exfoliated RPP flakes of~100-150 nm thickness (Supplementary Fig. 13) reveals a sharp redshifted emission peak in addition to the main exciton peak, as shown in Fig. 4a, which is distinct from defect-related broad emissions (BE) 44 . Using n = 4 RPP as an example, the redshifted narrow emission peak (NE2, full width at half maximum, FWHM~19 meV, Fig. 4b inset) at 665 nm appears just alongside the main exciton peak (NE1, FWHM~14 meV, Fig. 4b inset). Temperature-dependent PL of n = 4 ( Fig. 4b) shows that peak NE2 becomes slightly weaker in relation to the primary exciton peak NE1 upon temperature reduction to 77 K. Power-dependent PL measurement (Fig. 4c) verifies the excitonic nature of the redshifted emission peaks, where the integrated PL intensity can be fitted by the power-law relation with respect to the excitation intensity L as I~L k with k NE1 = 1.05 and k NE2 = 1.11 (Fig. 4d). Similar data for n = 2 and n = 3 are shown in Supplementary Fig. 14 and Fig. 15. Furthermore, the energy shifts between the two exciton peaks were analyzed across the RPP homologs series as shown in Fig. 4e.
Here, n = 2 RPP flake shows the largest shift with ΔE 2 = 44.7 meV, followed by n = 4 (ΔE 4 = 37.4 meV) and n = 3 (ΔE 3 = 30.4 meV) RPPs. This trend matches the relative out-of-plane octahedral tilt angle (blue rectangles in Fig. 3a) across the homologous series. In contrast to n > 1, the exciton peak blueshifts in n = 1 RPP due to phase transition of the organic cations at low temperatures 18 , and it seems that phase shift is allowed in n = 1 due to the very weak out-of-plane octahedral tilt. Importantly, the redshifted emission is uniform across the surface and is not related to edge emission 45,46 . It also does not originate from thickness-dependent photon recycling because covering the surface with hexagonal boron nitride (h-BN) removes the redshifted PL, independent of the flake thickness.
Since the out-of-plane octahedral tilt occurs after delamination, it is instructive to ask if the tilt can be suppressed if a monolayer of another material was stacked on top by van der Waals interaction. In view of the fact that h-BN has a large band gap (~6 eV) and is optically transparent and inert, we used it as a top confinement layer on freshly exfoliated RPP to examine whether it can suppress surface reconstruction and prevent the redshifted emission. Strikingly, the redshifted emission ( Fig. 5a and b) was only observed for freshly exfoliated RPPs after allowing the surface to relax in an inert atmosphere. If the exfoliated sample was immediately covered with h-BN, the redshifted emission was suppressed (Fig. 5c). To generate similar redshift in the PL for the h-BN confined sample, additional energy input in the form of laser annealing is needed as reported previously 39 . To understand the interfacial interaction, the h-BN/RPP heterostructure was modeled using DFT by stacking a monolayer of h-BN on n = 2 RPP monolayer and allowing the heterostructure to relax. As illustrated in Fig. 5d-f, we compare the changes in bond angles (out-of-plane tilt: Δθ and in-plane tilt: Δγ) and bond lengths (Supplementary Table 7) between the n = 2 bulk structure, monolayer, and h-BN covered monolayer. According to our simulations, the h-BN covered monolayer n = 2 RPP shows a very similar out-of-plane tilt angle (12.4°) as the bulk n = 2 RPP (11.5°), while it has a larger in-plane tilt (27.5°) than the bulk crystal (Table 2). Therefore, DFT simulation verifies that the outof-plane relaxation is suppressed in the presence of h-BN confinement layer compared with the uncovered monolayer perovskites.
One possibility for the redshifted PL is the coupling of the excitons to optical phonon arising from torsional motion of BA chains 22,23,47,48 ; the corresponding exciton-phonon Fröhlich coupling strength was reported to be around 30 to 60 meV, which agrees well with the redshifted PL of several tens of meV (e.g., ΔE 2 = 44.7 meV). Besides, our temperature-dependent Raman measurement is performed in exfoliated n = 2 RPP. As shown in Supplementary Fig. 16, the vibrational modes around 25 and 250 cm −1 are associated with torsional mode of the NH 3 + head of the organic spaces that is coupled to excitonic transitions 48,49 . On the other hand, our DFT simulation based on both nonadiabatic molecular dynamic and time-dependent DFT (TDDFT) methods corroborates the presence of a redshifted exitonic peak in exfoliated n = 2 RPP in addition to the main exciton peak seen in the bulk crystal ( Supplementary  Fig. 17).

Discussion
The results collected above allow us to examine the microscopic mechanism of surface octahedral tilt and its dependence on the nth dimensionality of 2D hybrid perovskites. According to our DFT calculations, octahedral tilt is intrinsic to the PbI 6 octahedral cage and occurs in bulk 2D RPPs. Depending on the interplay of electrostatic, steric, and hydrogen bonding interactions of BA + and MA + cations, the degree of octahedral tilt varies from n = 1 to n = 4 RPPs. Our finding is that in the exfoliated structure, the surface octahedral tilt is even stronger than the bulk and detailed structure comparisons of these two are demonstrated in Supplementary Fig. 18 and Tables 8-10. Similar to 3D hybrid perovskites, MA cations play a vital role in tuning the bulk octahedral tilt: (1) in n = 1 RPP where MA cations are absent, only in-plane tilt occurs; (2) in n ≥ 2 RPPs, absolute out-of-plane octahedral tilt increases due to steric effects from MA molecules, shown as orange triangles in the Fig. 3a. The out-of-plane octahedral tilt is stronger in n > 1 RPPs homologues that have a greater number of inorganic layers, this is due to the accumulating internal charge and structural mismatch with added MA ions. Interestingly, the surface octahedral tilt generates varying configurations of the top BA molecules owing to the minimization of electrostatic energies and steric repulsions. As shown in Fig. 6a and b, the BA organic chain is a polar molecule with a positively charged "head" of NH 3 + ion (shaded as blue) and a negatively charged "tail" of C 4 H 9 − ion (shaded as white). We compare the DFT-optimized surface structures of monolayer n = 1 (Fig. 6a) and n = 4 (Fig. 6b) RPPs to reveal how the configuration of BA molecules is commensurate with the octahedral tilt on account of the N-H···I hydrogen bond connecting them. In n = 1 RPP, the BA cations interlaced with each other due to the interaction with apical I atom (shaded as red) of surface PbI 6 octahedron, where a larger in-plane octahedral tilt and a small out-of-plane tilt exist. In the case of n = 4 RPP, the BA molecules adopt a "head-to-head" configuration so that their NH 3 groups are positioned to have strong hydrogen bonding interactions with the two apical I atoms of surface PbI 6 octahedron with large outof-plane tilt, resulting in alternating long and short I-I distances at surface that could be seen as the apparent "dimer" structures in the STM image of Fig. 3g.
We compare the surface (interface) octahedral tilt between monolayer and bulk structures according to the DFT simulations. The out-of-plane octahedral tilt was enhanced in exfoliated n ≥ 2 RPPs. n = 2 has the largest relative tilt angle of 13.9°(Δθ Exfoliated − Δθ Bulk ), following by n = 4 (5.2°) and n = 3 RPPs (3.1°), as shown by the blue rectangle plot in Fig. 3a. This trend is consistent with the energy shift between the two exciton peaks in the PL measurements (Fig. 4e), where n = 2 has the largest shift. On the other hand, the in-plane tilt is less correlated with the redshifted exciton peaks observed in PL. One possible explanation is that the out-of-plane octahedral tilt results in a stronger vibronic coupling of the organic cations with the lead iodide octahedral cage compared to the in-plane octahedral tilt.
The structural relaxation at the top layers is transmitted to the subsurface layers due to the hydrogen bonding interactions that causes synergistic movement of PbI 6 cages. Although the structural distortion is expected to decrease with distances from the surface, it must also necessarily imply that structural asymmetry builds up in the direction perpendicular to the plane of the quantum well (Supplementary Figs. 19-21 and Tables 11-13). The inversion asymmetric potential may lift the spin degeneracy of two-dimensional (2D) electronic bands-the so-called Rashba effect, causing it to split into two spin-polarized bands (Fig. 6c), where the momentum offset (k 0 ) refers to the distance between the extrema of the splitting band and the high symmetry point in k-space, and the splitting energy (ΔE) is the energy difference between them. Here, the Rashba splitting parameter is defined as: Redshifted PL peak appears for n ≥ 2 RPPs, whereas n = 1 RPP shows blue-shifted PL peak caused by phase transition instead.
DFT calculations were performed to see if the strength of the Rashba band splitting is correlated to the octahedral tilt observed for the RPP crystal 50 . The Brillouin zone and corresponding Rashba spin-polarized bands of the monolayer n = 2 structure are shown in Fig. 6d and e, in which the band minimum of the split band shifts from the high symmetry Γ point along X-Γ-Y directions. The band splitting energy and Rashba splitting coefficients are ΔE = 9 meV, α R = 1.15 eV·Å for electron and ΔE = 34 meV, α R = 1.55 eV·Å for hole along the Γ-Y band direction. The indirect band gap is 42 meV lower than the direct band gap at Γ point. The corresponding Rashba-type spin-textures of the spin-split valance bands are projected onto the 2D X-Y plane near Γ point in Fig. 6f. In the inner and outer branches, the sense of rotation (anti-clockwise and clockwise) of the spin-textures are opposite due to the spin-orbit interaction. By contrast, the band structure of the monolayer n = 1 RPP (Fig. 6g) shows no spin splitting. In view of the fact, there is minimal out-of-plane octahedral tilt in n = 1 as compared to strong out-of-plane octahedral tilt in n = 2, we can infer that the strength of the Rashba band splitting may be correlated to the out-of-plane octahedral tilt, among other factors.
Our studies show that the surface pairing of iodine atoms observed in the STM images of exfoliated 2D RPP perovskites is correlated to the out-of-plane octahedral tilt of the Pb-I cage at the surface, which is accentuated compared with the bulk. The octahedral tilt is a result of balancing steric forces and hydrogen  bonding interactions with methylammonium ions. The surfaceenhanced out-of-plane octahedral tilt is correlated to the emergence of a redshifted peak in the photoluminescence spectra, where the extent of the redshifts can be correlated to the relative octahedral tilt (surface versus bulk) from n = 2 to n = 4. Nonadiabatic molecular dynamic and TDDFT methods suggest that as a result of the surface octahedral tilt, sub-gap excitonic states are created and these can give rise to the redshifted PL. In addition, out-of-plane octahedral tilt at the surface enhances inversion asymmetry, and DFT calculations show that Rashba spin splitting is pronounced for n > 1 RPP perovskites that exhibit enhanced out-of-plate octahedral tilt.

Methods
2D RPP single crystals with the chemical formula (BA) 2 (MA) n−1 Pb n I 3n+1 (n = 1-4) were synthesized via a temperature-programmed solution precipitation method which have been reported before 39  STM measurements. STM measurements were carried out in ultrahigh vacuum with a base pressure of 2 × 10 −11 mbar using a commercial Omicron LT STM system. All STM images were collected under constant-current mode with electrochemically etched tungsten tips and all given voltages refer to the sample. Before STM measurements, the tip was conditioned and calibrated on a clean Au (111) surface. Images were collected at 4.5 K.
Steady-state photoluminescence measurement. The 400-nm fs laser pulses, which were obtained by frequency doubling of the 800-nm fs laser pulses (50 fs, 1 kHz) generated from a Ti: sapphire laser (Coherent Libra-HE) via the use of a BBO crystal, were used as the excitation light for the PL measurement. The 400-nm fs laser beam was focused onto the sample by an objective lens (Olympus LMPLFLN 50X/0.5). The PL emission from the surface of the sample was collected by the same objective lens and then coupled to a spectrum analyzer (Princeton SpectraPro 2750 integrated with a Princeton EMCCD) for spectrum analysis. For low-temperature (i.e., from 77 to 350 K) PL measurement, the sample was kept inside a vacuum chamber (Linkam optical DSC600 temperature-controlled stage) with an optical window. The optical window, which was purged with N 2 gas, was used to perform PL measurement at low temperature.
First-principles simulation. Density functional theory (DFT) calculations were performed with the Vienna ab initio simulations package (VASP) 51 . The Perdew-Burke-Ernzerhof (PBE) type of the generalized gradient approximation (GGA) was used as the exchange-correlation functional 52 . The van der Walls interactions were described by the DFT-D2 method 53 . A plane wave cutoff energy of 400 eV with a 3 × 3 × 1 Monkhorst-Pack grid was adopted for structure optimization, whereas a denser 5 × 5 × 1 grid for electronic calculation. All the structures were completely relaxed until the force on the atoms is <0.01 eV/Å. The (BA) 2 (MA) n−1 Pb n I 3n+1 (010) surface was represented using a periodic slab based on the space group Cc2m containing n inorganic atomic layers with a vacuum thickness of 15 Å. The spinorbital coupling was considered for the Rashba band and spin texture calculation.

Data availability
All data generated and analyzed in this study are included in the Article and its Supplementary Information, and are also available from authors upon request.