Designing exotic many-body states of atomic spin and motion in photonic crystals

Cold atoms coupled to photonic crystals constitute an exciting platform for exploring quantum many-body physics. For example, such systems offer the potential to realize strong photon-mediated forces between atoms, which depend on the atomic internal (spin) states, and where both the motional and spin degrees of freedom can exhibit long coherence times. An intriguing question then is whether exotic phases could arise, wherein crystalline or other spatial patterns and spin correlations are fundamentally tied together, an effect that is atypical in condensed matter systems. Here, we analyse one realistic model Hamiltonian in detail. We show that this previously unexplored system exhibits a rich phase diagram of emergent orders, including spatially dimerized spin-entangled pairs, a fluid of composite particles comprised of joint spin-phonon excitations, phonon-induced Néel ordering, and a fractional magnetization plateau associated with trimer formation.

R ich phenomena in condensed matter arise when quantum spin systems couple to phonons or orbital degrees of freedom of the underlying crystal lattice. Perhaps the most famous example is the spin-Peierls model [1][2][3][4] , wherein the spin interaction leads to a lattice instability resulting in a ground state of singlet pairs and a bond-ordered density wave. Motivated by this emergence of new physics, it is tempting then to consider the most extreme limit of coupling between spin and motionwhere the spin-carrying particles are completely free to move, and the spin-dependent forces thus dictate the properties of the emergent spatial order. Such a material would constitute a novel 'quantum crystal' that has not existed before, in which the emergent spatial patterns and spin properties are intricately locked together, and where driving one would automatically affect the properties of the other.
One possible route toward realizing such a material involves the new experimental platform consisting of cold atoms coupled to photonic crystal structures. Photonic crystals 5 are periodic dielectric structures in which the propagation of light can differ significantly from uniform media. An important feature of photonic crystals is the appearance of photonic band gaps, where strong interference in scattering from the periodic dielectric yields a complete absence of propagating modes within some bandwidth. Nominally, an excited atom whose transition frequency resides in the gap would not be able to spontaneously emit; instead, it has been shown that an atom-photon bound state can form, in which the atom becomes dressed by a localized photonic cloud [6][7][8][9][10] . The tight spatial confinement associated with this photon yields large dispersive forces on proximal atoms that depend on the atomic internal 'spin' states, which thus realizes the required spin-dependent forces to possibly observe the phenomena described above.
In this paper, we first describe the key features of atoms coupled to photonic crystals. While this interface in principle enables the realization of many different Hamiltonians 9-11 , we focus on one model where atoms are trapped in a weak one-dimensional external potential, and where a short-range spin-dependent force can be made sufficiently strong to exceed the external potential. To understand the emergent orders of this system, we begin by treating the motion of the atoms classically and their spins quantum mechanically. We find an effect reminiscent of the spin-Peierls transition, in which the atoms spatially dimerize and realize a high degree of entanglement within each dimer. We then proceed to a fully quantum model. Using density matrix renormalization group (DMRG), we find a rich variety of quantum phases beyond the spin-Peierls state, such as a state where spin and phonon excitations form composite particles, phonon-induced Néel ordering, and spatial trimers associated with magnetization plateaus. While here we study a specific model to create correlated spin-orbital quantum matter, more generally this work suggests that spin-orbital coupling can be a dominant phenomenon in all hybrid systems of atoms and photonic crystals. Similar considerations could also apply to a number of other atomic systems where spatially-dependent spin interactions can be realized, including polar molecules [12][13][14] , Rydberg atoms 15 , ion chains [16][17][18] , and atoms in high-finesse cavities 19 .

Results
Atoms coupled to photonic crystal waveguides. Photonic crystals 5 are periodic dielectric structures in which the propagation of light can differ significantly from uniform media (Fig. 1a). The dispersion relation in such structures consists in general of different bands, between which can appear bandgapsfrequency regions in which the light cannot propagate inside the crystal (Fig. 1b). Particularly rich phenomena are predicted to arise when an atomic transition is driven by a laser at a frequency within the bandgap. A specific example is illustrated in Fig. 1a, where two identical atoms are coupled to an 'alligator' photonic crystal waveguide (PCW), which consists of two separate waveguides whose modes hybridize with one another. Atoms have recently been coupled to such a structure in experiments described in refs 20-23. Here, the atoms are assumed to have three relevant electronic levels, with two ground (or metastable) states # j i, " j i connected by a common excited state e j i. The transition between ground state " j i and excited state e j i is globally driven by an external laser with frequency o L and Rabi frequency O L , while the transition between ground state # j i and e j i is coupled to the guided modes of the waveguide. In principle, an atom in state " j i could Raman scatter a pump photon into the waveguide and flip to state # j i. However, when the frequency o sc ¼ o L þ o m À o k of that scattered photon lies within the bandgap (see inset of Fig. 1b), it is unable to propagate and instead forms a bound state of length L around the atomic position. A second atom nearby in state # j i can subsequently absorb that photon, resulting in an effective spin flip interaction between the two atoms 9 . The effective spin Hamiltonian, generalized to many atoms, takes the form h j denotes the spin lowering operator from " j i to # j i, and conversely for s þ . J and L are the strength and characteristic length of the interaction respectively, which are tunable through the laser where a is the lattice constant of the PCW), is the Bloch function associated with the electric field at position x i . A more detailed microscopic derivation of the effective Hamiltonian (1) based on ref. 9 is provided in Supplementary Note 1 (also see Supplementary  Fig. 1). Here, we will assume that the atoms are tightly trapped in the transverse directions such that the position along x is the only dynamical variable. Absent any motional effects (i.e., if f is constant), equation (1) corresponds to the 'XX' spin model in 1D 24 .
The possibility to tune J, L and even the type of spin interaction makes the atom-photonic crystal interface a promising candidate for the simulation of many-body spin models with long-range interactions, when atoms are trapped at fixed positions [9][10][11] . Beyond photonic crystals, there are also a number of other proposed approaches to realize spin models with atoms 14,16-18,25-27 . Here, our goal is to investigate phenomena that can occur when motion is included as well. In particular, if one treats the position variables x i as dynamical degrees of freedom, the Hamiltonian in equation (1) should be regarded as a spin-dependent potential, wherein the forces experienced by the atoms can depend on the spin correlations. As these forces originate from the dispersive forces associated with photons confined to the nanoscale, their magnitude can be comparable to or much larger than those associated with conventional optical trapping forces 9 , implying that the physics of spin-motion coupling can become prominent in such systems. For example, using an electronic transition in a typical alkali atom, J can approach the GHz scale, as compared to BMHz scales for the excited state spontaneous emission rate and external trap frequencies (see Supplementary Note 1 for details). In emphasizing the role of spin correlations on motion, we also extend previous ideas involving self-organization of atoms in cavities or waveguides due to optical forces, where the atoms are treated essentially as classical dielectric particles with no internal degrees of freedom [28][29][30][31][32][33] .
Classical motion. We propose a realistic experimental setup, which highlights the interplay of spin and motion. Atoms interact via the Hamiltonian of equation (1), and are separately trapped by an external, spin-independent optical lattice H trap ¼V L P i sin 2 k tr x i (this could originate from optical fields in another guided band far from the atomic resonance). Peculiarly, this lattice traps atoms at the nodes of the Bloch function, and thus nominally hides the atoms from the PCW interaction. Despite not being a fundamental requirement to see spin-motion coupling, we assume that the trapping wavelength is such that atoms are localized around even nodes of the Bloch wave functions, i.e., k tr ¼ k/2 ¼ p/(2a), where a is the length of the unit cell of the PCW, as pictured in Fig. 2a. It can be readily shown that within our model, trapping atoms at every site would yield a phase transition with discontinuous change in the atomic positions.
We consider the Hamiltonian in the case of one atom per trapping site and an external magnetic field that can polarize the atoms with energy h along z: where d i denotes the displacement of atom i from the bottom of its external well. In the present section we treat the atomic position classically, while investigating the case of quantum motion in the next section. We assume that the coupling strength J is positive. For simplicity, in equation (2) we also ignore the self-interaction term (i ¼ j), which can be compensated by an external potential.
To study the many-body ground state of Hamiltonian (2) without any assumption about the spatial configuration is very difficult. Furthermore, for L=a ) 1 the long-range character of the interaction makes the spin model relatively difficult, even for fixed positions. As a consequence, we restrict our attention to the case LBa, for which we can make a nearest-neighbor approximation. We can get an intuition of the possible ground state configuration of a system of many atoms by considering how just two atoms in neighboring sites interact. If the atoms remain at the bottom of their trapping wells, the function f(x 1 , x 2 ) ¼ 0 as these positions coincide with nodes of the Bloch functions. However, the PCW interaction energy would become negative, if the two atoms were to form a triplet state, (or a singlet for Jo0), and simultaneously displace toward each other to form a spatial dimer. Such a process would become energetically favourable overall for a certain ratio of J/V L . Motivated by this simple case we make an ansatz that the spatial configuration of the many-body ground state consists of dimerized pairs. In particular, we assume that where d represents the displacement from the trap center, as pictured in Fig. 2a. This is reminiscent of the lattice instability that creates entangled dimers in the spin-Peierls model 1 , but with the substantial difference that our system becomes noninteracting in the absence of dimerization (as the atoms are at the nodes). In the following, we treat d as a variational parameter and proceed to solve the spin ground state exactly.
The nearest-neighbor spin Hamiltonian can be mapped to a chain of spinless fermions through standard Jordan-Wigner transformation 34 , with the presence/absence of a fermion on a site corresponding to spin up/down, respectively. Because of the staggered spatial configuration, it is natural to define a unit cell j consisting of a pair of dimerized atoms (labelled L,R). Two different spin couplings J S;W d ð Þ¼J sin 2 kd e À 2a Ç 2d ð Þ =L then characterize the interaction between atoms within the same dimer, and between consecutive atoms R,L in neighboring dimers, respectively (see Fig. 2a). The Hamiltonian then reads where c (L,R),j are fermion annihilation operators for site j. Just as in the standard Jordan-Wigner transformation, this two-spin per-site Hamiltonian can be exactly or numerically diagonalized by moving to Fourier space, as we describe in Supplementary Note 2.
By minimizing the ground-state energy of H(d) with respect to d we find the optimal spatial configuration (within the ansatz). In Fig. 2b we plot the resulting value of d as a function of the interaction strength J and of the magnetic field h (in units of V L ). In the J À h plane one can clearly distinguish a critical value of the spin interaction strength, J crit (h), above which a phase transition occurs from a non-interacting to a dimerized state. The increase in spin entanglement with dimerization can be quantified by taking the two-particle reduced density matrix r 2S of atoms within a dimer, and calculating its overlap with the triplet state, T S ðdÞ¼ T h jr 2S T j i. We plot T S (d) in Fig. 2c for zero magnetic field. For d ¼ 0 this quantity tends to the value in the conventional XX spin model, T S (0) ¼ (1/2 þ 1/p) 2 E0.67, while for large values of d and small L it tends to 1. Similarly, defining an analogous quantity T W (d) between consecutive atoms in neighboring dimers, we find a decrease in correlation with increasing dimerization.
Quantum motion. We now consider a quantum description of motion and spins, which is relevant, e.g., if the motion is initially cooled to its ground state. We keep the assumption of tight trapping of the atoms around the minima of the external potential, such that tunneling of atoms between sites can be neglected. We then proceed by projecting the Hamiltonian of equation (2) onto the lowest two motional bands, and denote by a j i i and b j i i the associated Wannier functions localized around site i, as shown in Fig. 3a (see Supplementary Note 3). This represents the minimal model in which spin and motion can couple, since superpositions of states a j i and b j i yield spatial wave-functions that are displaced from the site centers. We have also performed calculations involving a third band to verify that the conclusions made from the two-band approximation do not qualitatively change (see Supplementary Note 4 and Supplementary Fig. 2). Within the two-band approximation, it will be convenient to introduce a set of pseudo-spin operators on each site, The overall Hamiltonian can thus be expressed in terms of these operators as The terms proportional to D and h describe the energy arising from the band and magnetic field, respectively, while the remainder describes the spin-motion coupling. Here (w a,b being the Wannier functions of states a j i and b j i). In the following we take L ¼ 2a and the ratio between the trapping lattice depth V L and the recoil energy E R to be 20, for which numerical evaluation of the Wannier functions yields Z 0 E0.54, Z a E0.06 and Z b E0. 16. The terms between brackets contain the dependence of the spin interaction on the motional state of the atoms and have a simple physical explanation. The dominant s x is x i þ 1 term has largest amplitude when both atoms sit in an equal superposition of states a j i and b j i (i.e., the wavefunction is maximally displaced from the center), which reflects that the atoms are trapped at nodes of the PCW. The other terms, which are smaller, originate from the exponentially decaying interaction and are responsible for spatial dimerization. It is interesting to note that this Hamiltonian constitutes an extreme case of spin-orbit coupled systems, as neither an orbital kinetic energy nor a motion-independent spin interaction appear.
We study the phase diagram of Hamiltonian (4) in the g À h plane by means of a finite-size density matrix renormalization group (DMRG) algorithm 35 . The resulting phase diagram for 0rg,hr2D is shown in Fig. 3b for N ¼ 62 atoms, where we can clearly distinguish at least six phases. The procedure by which these phases are numerically demarcated is discussed in detail in Supplementary Note 3, while in the main text we describe the salient physical properties of each phase. First, for sufficiently large magnetic fields h4h crit (g), with h crit (0) ¼ 0, the spins are fully polarized and thus the spin-motion coupling has no effect. The many-body state is thus separable, with each atom residing in the lowest motional band, c j i¼ a; # j i N ('P' phase in Fig. 3b).
Along the g-axis up to g crit we have a Néel ordered phase 'N', where the magnetization per atom M z ¼1=ð2NÞ P i hs z i i is zero and the Néel order parameter F¼ð1=NÞ P i ð À 1Þ i hs z i i has a finite value, as shown in Fig. 3c. This phase also extends to finite values of h with a lobe-like shape. The existence of this phase can be predicted analytically in the weak coupling regime, i.e., for g/ 2D small, such that the high-energy excitations associated with populating the upper band can be effectively integrated out. In particular, through a Schrieffer-Wolff transformation 36 on equation (4) one obtains the following effective Hamiltonian acting only on the spin degrees of freedom: Here J 1 ¼ g 2 (1 þ 4w 2 )/2D, J 2 ¼ g 2 w 2 /D and w ¼ Z a /(Z 0 L). Hamiltonian (5) describes a nearest neighbor anti-ferromagnetic (AF) Ising model with an additional XX term coupling next-nearest neighbors, with all such terms mediated by virtual phonons. For example, the spin-motion term in equation (4) proportional tos x is x i þ 1 enables a fluctuation where two consecutive atoms, anti-aligned in their spins, jump to the higher band and exchange their spins, before returning to the original state (see Fig. 3d). This process results in a lower energy for the anti-aligned configuration and produces the longitudinal s z i s z i þ 1 À 1 À Á term in (5). For zero magnetic field, given that J 1 ) J 2 the ground state exhibits AF ordering along z (FE1), as illustrated in Fig. 3c. On the other hand, for h4h crit (g) all spins are in state # j i. Intuitively, one can expect that the transition from Néel ordering to polarized occurs with all # j i spins in the Néel , which has two phase transitions to polarized phases (for subchain B) at h ¼ 2(J 1 ± J 2 ). It follows that for ho2(J 1 À J 2 ) the total system (A and B) is in the Néel phase, while for h42(J 1 þ J 2 ) it is in the P phase. In between the two phases the subchain melts under the effective XX model. Since J 2 ( J 1 , in the g À h plane this transition region is too narrow to quantitatively match the DMRG results to the XX model predictions, although the effective theory gives correctly the boundary between N and P at h crit (g)Eg 2 /D for g ( D (solid line in Fig. 3b).
The Néel order extends to values of g=D\1 where the low-energy description of (5) is no longer accurate, and decreases discontinuously to zero with the onset of a new phase of dimerized triplets (labelled 'D' in Fig. 3b). This phase is characterized by zero magnetization and a non-zero spin triplet dimer order parameter, defined as being the spin triplet state and r i,i þ 1 the two-site spin reduced density matrix (Fig. 3c). It also has a non-zero spatial dimer order parameter, defined as D x ¼ð1=NÞ P i ð À 1Þ is x i . The entangled dimerized structure is evident in Fig. 4a, where we plot the triplet fraction in the two-particle density matrix, T h jr i;i þ 1 T j i and the displacement s x i in a part of the chain for (g, h) ¼ (1.7, 0.2)D. Also, we can observe thats x i $ 1. Thus, the two-band approximation for the atomic motion is technically violated since the displacement from the trap center is saturated. However, in Supplementary Note 4, we present calculations involving a third motional band, which allows for a greater maximum displacement of atoms. These calculations exhibit a slower onset of saturation with increasing g and no appearance of new phases (at least within the range of parameters considered). Together, this suggests that an exact calculation involving all bands, although directly unfeasible, would produce a result similar to the previously discussed case of classical motion, with a steadily increasing degree of dimerization and triplet fraction with increasing g.
For simultaneously large values of g and h, there is a spinmotion fluid phase ('SMF') where the system is gapless and the magnetic field strongly polarizes the spins, such that M z is close to À 1/2. This phase corresponds with good approximation to the ground state of the XX Hamiltonian Here t z i is the Pauli matrix with eigenstates + j i¼ a; # j i and * j i¼ b; " j i, while t AE are associated raising and lowering operators. Thus, this phase corresponds to a dilute fluid of composite flips of spin and motion. The existence of this phase can be understood by noting that for large magnetic field, the system is only dilutely populated by spins pointing up. Thus the terms in equation (4) proportional to Z a,b that are responsible for dimerization can be neglected. The structure of the remaining Hamiltonian connects naturally the states + j i directly to * j i, in the form of H þ (see Supplementary Note 3 for details). The locking between spin and motional correlations can be observed in Fig. 4b, where the expectations values of s z i ands z i obtained with DMRG are plotted for a representative point in the phase. The oscillations of s z i and s z i are due to the open boundary conditions in a finite system and are observable also in a pure XX model. In Fig. 4c the magnetization curve predicted by H þ is compared with the numerical result from the DMRG study of the full Hamiltonian for g ¼ 1.6D, showing good agreement, while in Fig. 3b the predicted boundary with the 'P' phase h crit (g)E À D þ 2g is represented by a dashed line.
For À 1=4tM z o0, H þ no longer serves as a good description for the ground state. Most of this region consists of a set of phases 'U' whose origin is not completely understood yet. However, for strong interactions g=D\1, the system qualitatively appears to behave as an interacting Luttinger liquid for the t particles. Numerical evidence is shown in Fig. 5a, where the twopoint correlation functions C X ij hX þ i X À j i À hX þ i ihX À j i are plotted for various X¼t; s;s, for a representative set of values (g, h) ¼ (1.74, 1.38)D. In particular, if t behaves as a Luttinger liquid, then the long-range decay of interactions is predicted to have a power law form of C t ij $ ð À 1Þ j À i j j j À i j j À 1=2K (ref. 37). The inset of Fig. 5a plots the absolute value jC t ij j on a log-log scale, which confirms an approximate power law decay. On the other hand, correlations of the other degrees of freedom exhibit more erratic behavior. Similar observations hold for the density correlation functions D X ij hX z i X z j i À hX z i ihX z j i (Fig. 5b). We fit the Luttinger parameter K (ref. 37) from the numerical data, taking the ten (g, h) values indicated by red stars in Fig. 3b across the SMF to U boundary. These fits are performed on the points |i À j|44 of C t ij , in order to reduce the influence of short-range corrections, which exist even for an ideal Luttinger liquid 37 . The inset of Fig. 5a shows the best fit (red dashed line) for (g, h) ¼ (1.74, 1.38)D, while the fitted values of K for all ten chosen (g, h) points are plotted in Fig. 5c. We have also simultaneously plotted x, the sum of the squares of the residuals between the best linear fit on a log-log scale and the numerical data. We note that while the choice of region of exclusion of |i À j|r4 in taking the fit is somewhat arbitrary, modifying this region (or excluding no points at all) does not change the qualitative conclusions. The decrease below K ¼ 1 is indicative of the formation of a charge density wave phase with quasi-longrange order, i.e., algebraic decay of the correlation functions. We thus denote this part of the phase diagram as SMF(CDW). The precise boundary of this phase and the nature of the transition to neighboring phases is still not completely understood.
Approaching M z -À 1/6 we notice that not only the fitted value of K tends to zero but also the quality of the fit decreases rapidly, as indicated by the increase in the residual error x. This indicates a change of the decay of the correlation function from polynomial to exponential. This is in agreement with the fact that the decrease in K is also known to facilitate the possibility of phases with spontaneously broken symmetry, which is observed in our system as well. At M z ¼ À 1/6 (one third of the maximum magnetization), we observe indeed the presence of a plateau in the magnetization curve (Figs 3c and 4c), for values of g sufficiently large. In this region the ground state assumes a trimerized configuration, as shown in Fig. 5d, where s z i ;s z i and the displacements x i are plotted. While we are not able to predict the appearance of such a plateau in our model from first principles, we note that all of its features are consistent with the conditions of ref. 38. In particular, our Hamiltonian allows for a gapped phase with spontaneously broken symmetry in the ground state with spatial periodicity n ¼ 3, provided that the quantization condition n(S À |M z |) ¼ integer is satisfied (here S ¼ 1/2 is the total spin). Such a gapped phase should be accompanied by a magnetization plateau.

Discussion
The platform of cold atoms coupled to photonic crystals offers fascinating opportunities to create quantum materials in which spin and motion interact strongly with one another. We have analysed in detail the ground state properties of one experimentally feasible setup, but there exist many exciting avenues for future research. The field of interfacing atoms and photonic crystals is still new and rapidly developing, which makes it difficult to say precisely how the ground state or nearby states can be probed and prepared, but we briefly describe some of the possibilities here. First, it has already been demonstrated that tightly focused optical tweezers can be used to controllably position single atoms nearby nanophotonic structures and couple the atom to the optical mode 39 . Separately, there have been spectacular experiments to create arrays of up to B10 2 atoms in individual optical tweezers [40][41][42] , and demonstrated capabilities in such systems for motional ground-state cooling and spin readout 40 . An optical tweezer array applied to nanophotonic systems could then be a promising route toward both Absent single-site measurements, there are a number of global measurements that could be applied to yield signatures of the various phases. For example, it has also been theoretically and experimentally shown 28,43-45 that different atomic spatial patterns can give rise to very different global reflection and transmission spectra for a weak guided probe field. Similar to free space, a guided mode could also be used to efficiently read out global spin properties 46 . In terms of preparation of the ground state, one likely possibility would be through adiabatic evolution (given that the atomic 'spin' states are internal states that do not readily thermalize). Here, the atoms would be initially optically pumped to a separable state (such as # j i N ), which corresponds to the ground state of a single-particle Hamiltonian H s . The system could then adiabatically evolve through a Hamiltonian H(t) ¼ H s (t) þ H int (t), where the single-particle Hamiltonian is gradually turned off while the PCW interactions are turned on. Understanding the fidelity of this process requires a more thorough investigation of the excitation spectrum, which itself should exhibit non-trivial properties, including the possibility of signatures of fractional spin 47 .
The strong coupling between spin and motion more broadly invites a number of other intriguing questions. For example, it would be interesting to understand the transport properties when spin and motion strongly hybridize. Moreover, it would be highly interesting to consider models without an external lattice potential, and investigate whether the spin interaction alone can produce full spin-entangled crystallization. One might also consider models where the spin part of the interaction already exhibits non-trivial character, such as frustration or topology. Finally, in terms of applications, it would be interesting to explore whether specially engineered spin-motion Hamiltonians can give rise to useful many-body spin states (such as squeezed states for metrology), when the spin interaction alone is incapable of producing such states.

Methods
Density matrix renormalization group. The density matrix renormalization group (DMRG) algorithm is a well-established numerical method for ground state studies of one-dimensional systems 35 . It consists of approximating the exact ground state GS j i¼ P i1;::;iN c i1;::;iN i 1 ; ::; i N j iof a finite-size system of N sites with local dimension d with a matrix product state (MPS), i.e., a state of the form GS j i MPS ¼ X where the matrices A have maximum dimension D. The DMRG algorithm finds the matrices fA k ½ ik g for which hHi is minimum. This is obtained in an iterative way by optimizing a single matrix in every step, keeping constant all the others and shifting one site at every step. In 2N steps the system has been 'swept' once. The level of approximation can be determined by the relative energy difference dE after a sweep. In general for a gapped Hamiltonian, with D ¼ 100-1000 the energy converges in a few sweeps and one obtains an extremely good approximation of the ground state.
In our numerics we take D ¼ 100 and the precision on the total energy dE ¼ 10 À 7 with a maximum number of sweeps equal to 6. For the points in which convergence is not reached within six sweeps D is increased to 140 and additional sweeps are performed. We span the g À h plane with a resolution of 0.02D.
Data availability. The data that support the findings of this study are available from the corresponding author on request.