{'fraction_non_alphanumeric': 0.044072900655604, 'fraction_numerical': 0.020075912042410365, 'mean_word_length': 4.074009231258953, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Complexity and quenches in models with three and four spin interactions

We study information theoretic quantities in models with three and four spin interactions. These models show distinctive characteristics compared to their nearest neighbour counterparts. Here, we quantify these in terms of the Nielsen complexity in static and quench scenarios, the Fubini-Study complexity, and the entanglement entropy. The models that we study have a rich phase structure, and we show how the difference in the nature of phase transitions in these, compared to ones with nearest neighbour interactions, result in different behaviour of information theoretic quantities, from ones known in the literature. For example, the derivative of the Nielsen complexity does not diverge but shows a discontinuity near continuous phase transitions, and the Fubini-Study complexity may be regular and continuous across such transitions. The entanglement entropy shows a novel discontinuity both at first and second order quantum phase transitions. We also study multiple quench scenarios in these models and contrast these with quenches in the transverse XY model.


I. INTRODUCTION
Over the last few years, study of information theory in quantum many body systems have fast gained popularity. A primary reason for this is that apart from being the foundations of what is perhaps the technology of the future, such studies point to the deep connection between diverse areas of physics. Indeed, various information theoretic measures such as the quantum information metric (QIM) [1][2][3] and its associated Fubini-Study complexity (FSC) [4], Nielsen complexity (NC) [5,6], Loschmidt echo (LE) [7,8] and entanglement entropy (EE) [9][10][11][12], find extensive applications in different areas of physics, from statistical systems to black holes, the connections often being realised in quantum field theory via the gaugegravity duality [13][14][15][16][17][18][19][20][21][22].
Spin systems provide an ideal laboratory for studying information theory in a quantum mechanical setting, primarily because these are often amenable to analytical results and provide deep physical insights regarding the behaviour of information geometry near zero temperature quantum phase transitions (QPTs). As of now, these are well studied in, for example, the transverse XY model, a quasi-free fermionic model with nearest neighbour (NN) interactions [23,24], and the Lipkin-Meshkov-Glick model, which is an exactly solvable model with infinite range interactions [25,26], in both time-independent and time-dependent (quench) scenarios. The purpose of the present paper is to go beyond the physics of nearest neighbour interactions and study these information theoretic quantities in spin models in the presence of three and four spin interactions. Fortunately, an exactly solv-able model with such interactions is known, and was discussed in [27] (see also [28,29]). Here, we study this model with some specific simplifying choice of parameters (that nonetheless preserves a rich phase structure), and show that there are some remarkable differences (as well as similarities) in the behaviour of information theoretic quantities when short range interactions are added to the NN one.
The motivation for this study is twofold. Firstly, to the best of our knowledge, short range interactions apart from the ubiquitous NN one and the infinite range Lipkin-Meshkov-Glick model have not been studied in the current literature on spin models. As we show in sequel, our model shows several distinctive features as far as information theoretic quantities are concerned. Secondly, the nature of the phase transitions in the model of [27] are very different from the ones in say, the transverse XY model, and it is of interest to understand the behaviour of information theory in this context, and we show in sequel that some of the well known features of models with NN interactions are modified here.
With these motivations, in this paper, we first compute the QIM and the structure of geodesics for all the phases of the theory in the presence of four spin interactions. Our analysis indicates that contrary to some known examples, the Ricci scalar of the QIM and the FSC do not capture the entire phase structure of the theory. Interestingly, these show their typical behaviours only at a particular type of second order phase transition, namely the gapped ferromagnetic to the gapless ferrimagnetic phases, and not between two ferrimagnetic phases (with different magnetisations). Next, we compute the NC for this model and show that its derivative does not diverge at second order phase transitions, but rather shows a discontinuity, again in contrast to known examples in the literature. A quench scenario is considered thereafter, where we couple our model to a spin half system, and compute the LE. Thereafter, using the density matrix renormalisation group methods [2,[30][31][32], we compute arXiv:2207.14090v1 [quant-ph] 28 Jul 2022 the entanglement entropy, and find sharp jumps across both first order and second order phase transitions, a result that is in contradiction to the ones known in the literature, namely that the EE shows discontinuity at a first-order QPT, while a cusp or a kink in the EE indicates the second-order QPT [33,34]. Finally, we also consider multiple global quench scenarios in the model and contrast this with such quenches in the transverse XY model.
In the following sections II to VI, we elaborate upon the above points, and the paper ends with our conclusions in section VII. This paper also includes two appendices A and B where some mathematical details are relegated to. Units with = 1 are used throughout.

II. THE MODEL
A quantum spin models with alternating NN couplings, and three spin and four spin exchange interactions was studied by [27]. This will be used extensively in the paper, and we first recall their main results. More details can be found in [28,29]. The model Hamiltonian is given by Here, the indices 1 and 2 indicate two sub-lattices with nearest neighbour spins occupying two different sublattices (the total number of cells will be denoted by N in what follows, where we will assume N to be odd) and S x,y,z n,1,2 are the spin-1 2 operators for sub-lattices 1 and 2 at the n-th site. Also, µ 1,2 are the Bohr magnetons in sublattices 1 and 2, H is the external magnetic field along the z axis, J 1,2 are the alternating coupling constants between nearest-neighbour spins, J 13,23 are the alternating coupling constants for three-spin interaction, and similarly, J 14,24 are the alternating coupling constants for four-spin interactions.
The large number of interaction parameters make a general analysis of the model cumbersome even if one chooses an overall scale, and we will make some simplifying choices. For ease of computation, we will in this paper, set J 24 = 0, and it is further convenient to reparametrise the Hamiltonian in terms of new couplings, h, J, and J 3 to get rid of some unwanted fractions, as : Then the Hamiltonian Eq. (1) in terms of these new couplings is The Hamiltonian Eq. (2) can be diagonalised, and dispersion relations can be obtained by using the Jordan-Wigner, Fourier and Bogoliubov transformations detailed in Appendix A, and after these transformations, we finally obtain with the dispersion relation,  Consider the nonnegative values of the magnetic field parameter h ≥ 0 (from here and onwards, the magnetic field is specified by h instead of H). There are four critical values of magnetic field, denoted by h 1,2,3 and h 13 , and the explicit expressions of h 1,2,3 can be evaluated by putting E k,1,2 = 0 for k = 0 and π, while for h 13 , we need to maximise the solution of E k,1,2 = 0 with respect to k. We find and where the expression of k m is lengthy, and is given in Appendix B. First, consider the case with a small magnetic field, 0 < h < h 3 , and values of three spin-exchange interactions, 0 < J 3 < 2/ √ 5. Here, the critical line h 13 touches the line h 1 at J 3 = 0.76 for J = 1, and h 2,3 intersect the J 3 -axis at J 3 = (2 √ J 2 )/ √ 5. In the first branch, E k,1 is negative, and in the second branch, E k,2 is positive for any k. This domain is denoted by I in the phase diagram and is referred to as the ferrimagnetic phase, as can be gleaned from the ground state magnetisation per cell in the thermodynamic limit given by [28,29] In the region h 2,3 < h < h 1 , in the first branch, E k,1 < 0 for |k| < k c1 and E k,1 > 0 for |k| > k c1 , while for the second branch, E k,2 > 0 for any k ∈ (−π, π). This region is denoted as II in the phase diagram, and is a ferrimagnetic phase with gapless excitations. Similarly, the region between h 1 < h < h 13 for 0 < J 3 < 0.76 denoted by III in the phase diagram is a ferrimagnetic phase in which the first branch E k,1 < 0, for k ∈ (−k c1 , −k c2 )∪(k c2 , k c1 ), while in the second branch, E k,2 > 0 for any k ∈ (−π, π) with Consider now higher values of three spin-exchange interaction J 3 > 2/ √ 5, and magnetic field in the range 0 < h < h 2 . Here, for the first branch, E k,1 < 0 for |k| < k c1 , and for the second branch E k,2 < 0 for |k| < k c2 . The system is in the ferrimagnetic state with gapless excitations and is denoted by IV in the phase diagram. Finally, for small values of J 3 < 0.76, both the branches are positive for h > h 13 , while for higher values of J 3 > 0.76, both branches are positive for h > h 1 . This region is denoted by V in the phase diagram, and the system will show ferromagnetic behaviour with gapped excitations.
The ground state of the model of Eq. (2) varies with parameters h and J 3 , and they are classified with distinct ranges of k values in different regions of the phase diagram. The ground state of phase I can be described as a vacuum state with the first branch fully occupied and the second branch empty, i.e., b k,1 |Ψ with k ∈ (−π, π). In region II of the phase diagram, b k,1 |Ψ II = 0 for |k| < k c1 and b † k,2 |Ψ II = 0 for k ∈ (−π, π). This implies that the ground state in phase II is |Ψ II = k |Ψ k,1,2 I for |k| < k c1 and is trivial for other values of k in this region. Similarly, in the region III and IV of the phase diagram, the ground state is given by Note that the limits of k in the ground states of the regions III and IV are the same, but the values of h and J 3 are different. The ground state of region V i.e., the ferromagnetic phase of the phase diagram is trivial since both the energy eigenvalue branches are positive for k ∈ (−π, π), and is given by By the same method, for the sake of completeness, we also draw the phase diagram for the simpler case up to three spin interactions with parameters µ 1 = µ 2 = 1, J 1 = 2, J 2 = 1, J 13 = J 23 = J 3 ≥ 0, and J 14 = J 24 = 0. The ground-state phase diagram is shown in Fig. (2), where F/AF denotes the ferromagnetic/antiferromagnetic phases, and 1 and 2 denote ferrimagnetic phases. The lines H c1,2 and H s are the lines of second-order QPT with critical values Also, the line H = 0, J 3 > 1 is the line of first order QPT, i.e., the spontaneous magnetisation changes from zero to a finite value.

III. THE QIM AND COMPLEXITY
The ground state wavefunction Eq. (9) is defined on the parameter space h and J 3 for a fixed value of J = 1. The QIM measures the distance between two quantum states |Ψ( λ) and |Ψ( λ + dλ) , where λ = {h, J 3 } and dλ is the infinitesimal separation in parameter space. The distance can be expressed as: with i = 1, 2, ...., m, and m being the dimensions of parameter space, which in our case is 2. We now find the QIM in the regions of the ground state phase diagram.
For small values of J 3 in region I, we have calculated the QIM numerically in the thermodynamic limit, and found that the numerical results are in good agreement with the analytical expressions. In all other regions, we only rely on the numerical results due to the complicated expressions of k c1,2 . We now turn to the QIM of region I, which takes the following form: where k ∈ (−π, π). In the thermodynamic limit, the procedure to evaluate these is simple : the summations in the above metric components are converted into integrals, i.e., we replace k = N 2π π −π dk, and then do a standard computation of residues in the complex plane. The analytic expressions of the metric components are too lengthy to be presented here. However, we have expanded these metric components by assuming a small value of J 3 and then analysed the Ricci scalar, R. Even then the analytic expression for R is cumbersome, and will not be presented here. Numerical values of R are shown in Fig. (3). Importantly, there is no divergence on the critical line h 3 . In region II, the components of the QIM take the same form as given in Eq. (12), but here k ∈ (−k c1 , k c1 ). We have computed the Ricci Scalar numerically, which shows oscillations in this region and is divergent on approaching the critical line h 1 . These oscillations in R are not particularly interesting, as these are due to the discontinuity in the metric tensor components, since k is not summed over the entire range (−π, π) in the metric tensor : the summation here runs for k in the range (− The integer values of λ c1 are different for distinct values of h and J 3 . For example, if (J 3 , h) = (2, 1.2), then λ c1 = 14.05, and the summation runs over 0 to 14, while for (J 3 , h) = (2, 1.4), λ c1 = 13.45, and summation runs over 0 to 13. This creates discontinuities in the nature of the metric tensor, which is reflected as oscillatory nature in the Ricci scalar plots.
Similarly, in region III, the summation in QIM tensor in Eq. (12) runs over (−k c1 , −k c2 ) ∪ (k c2 , k c1 ). In this region also, the Ricci scalar shows an apparently oscillatory nature, and has a clear divergence on approaching the critical line h 13 . The behaviour of R as a function of h is shown in Fig. (4). In the same way, we have computed R in region IV . Here, it turns out that there is no divergence on the critical line h 2 , which indicates that the geometry of the ground state is regular everywhere except on the critical lines h 13 , and h 1 for J 3 > 0.76, i.e., transition lines which separate the ferrimagnetic from ferromagnetic phases.

A. The geodesics and FSC
Now we numerically compute the geodesics and FSC in the ground state phase diagram. The method is standard and has been described in details for spin systems in [23]. Namely, to obtain the geodesics, we parametrise h and J 3 as a function of an affine parameter τ , and then we solve the two coupled second order differential equations given from where x i = {h(τ ), J 3 (τ )}, and Γ i jl are Christoffel symbols. We solve these two equations numerically in each region of the phase diagram by specifying the values of four initial conditions h(τ = 0), J 3 (τ = 0),ḣ(τ = 0), J 3 (τ = 0), and one normalisation condition g ijẋ iẋj = 1, with dot implies a derivative with respect to τ . This then determines h(τ ) and J 3 (τ ). In regions I and IV , the geodesics do not show any special behaviour since the metric and the Ricci scalar are regular throughout these regions. We have inverted the solution of geodesics to find the FSC, i.e., τ (h), using a standard root finding procedure in Mathematica. These do not have any special behaviour up to the phase boundary h 2,3 , reflecting that the ground state manifold is regular on the critical line h 2,3 .
In region II, the behavior of a typical geodesic is shown in Fig. (5). In our numerical analysis, we have chosen the initial conditions h(τ = 0) = 1, J 3 (τ = 0) = 1, J 3 (τ = 0) = 0.04, andḣ(τ = 0) = 6.78 is determined from the normalisation condition. The geodesics that start from any coordinate (J 3 , h) where J 3 < 0.76 and h 3 < h < h 1 approach the phase boundary h 1 but do not show any interesting behaviour. On the other hand, the geodesics that start from the J 3 > 0.76 and h > h 2 keeps approaching the critical line h 1 , attains a minimum distance with the critical line, but will never cross the phase boundary, as shown in Fig. (5). The FSC, and its derivative are numerically evaluated, and the derivative, ∂C F S /∂h is plotted in Fig. (6). The derivative of FSC approaches zero when the distance between the geodesic and phase boundary line is minimal. This can be explained by a simple relation that we have obtained previously in [23], namely, Here, g hh → 0 at h → h 1 for J 3 > 0.76 and J = 1 is fixed.

B. The Nielsen Complexity
For a fixed value of J = 1, the ground state of the model can be expressed in the form: where e −iφ k = (1 − i sin k)/ 1 + sin 2 k. By following [35], the NC of our model can be computed straightforwardly, i.e., C N = k |∆θ k | 2 , where ∆θ k = (θ T k − θ R k )/2, and Hereafter, the superscripts R and T on a variable are used to denote the reference and the target values respectively, of that variable. Note that unlike the Kitaev model [35] or the XY spin chain [23], the derivative of the NC does not have any divergences on any of the critical lines h T 1,2,3,13 . This derivative of NC with respect to h T is plotted in Fig. (7) as a function of h T . It can be seen that the derivative approaches zero on the critical line h T 13 , in region III of the phase diagram. This behaviour can be explained using the relation [23] Since g h T h T → 0 on h T → h T 13 , so the derivative of NC goes to zero on this critical line. A similar behaviour is observed as h T → h T 1 for J T 3 > 0.76 in region II.

IV. THE TIME-DEPENDENT NIELSEN COMPLEXITY AND LOSCHMIDT ECHO
We now consider the model Eq.
(2) with a quantum quench, assuming the sudden quench approximation. Following the work of [36], we consider our model coupled to a two-level central spin-1/2 system. Here, we will take the interaction Hamiltonian to be H I , so that the total Hamiltonian is H F = H + H I , with H given in Eq. (2), and where δ is the interaction strength, and |e is the excited state of the two-level system. By a standard method, we write the ground state |Ψ h,J3 of H in terms of that of H F , labeled |Ψ k,1,2 h+δ,J3 for the k-th Fourier mode. This gives where we have defined with the operators d k,1 and d k,2 being the Fourier operators, and θ k , u k , and v k are defined in Eq. (A5) with the arguments appropriately shifted.
The computation of the NC is now standard (see, e.g., [35]). The final result after the quench is where with the single particle excitations Λ k (h + δ, J 3 ) = h+δ 2 − J 3 cos k 2 + 1 + sin 2 k. A particularly interesting quantity in a system after quench is the Loschmidt Echo. The system is initially prepared in the ground state |Ψ h,J3 of the Hamiltonian H, and it is then quenched to evolve according to the Hamiltonian H F . The LE is defined as As discussed in [24], the relation between LE and NC for small times, L e −C N is also valid for our spin chain model. In fact, Fig. (8) shows that this relation continues to hold for large times as well. We have plotted the dynamical behavior of the time-dependent NC per system size, C N (t)/N along with negative of logarithm of LE per system size, − log(L)/N for N = 101, δ = 0.1, J 3 = 0.2, and J = 1. The solid red (scattered black squares), solid green (scattered purple triangles), and solid blue (scattered brown stars) shows the numerical results for C N (t)/N (− log(L)/N ) corresponding to h = 0.5, 1, and 1.4, for regions I, II, and III of the phase diagrams respectively.
As anticipated, the time-dependent NC first increases linearly and then starts oscillating [24,35]. However, in our model, the oscillations do not die out to approach a time-independent value, contrary to what is observed in those works. This oscillatory behaviour is observed here in all the regions of the phase diagram for any finite or large times and is well explained by the results in [24]. That is, the time-dependent NC and LE are oscillatory functions of k due to the nature of Ω k and Λ k . The oscillations die out rapidly in the cases considered in [24,35] because a large number of Fourier modes start contributing with the increase in time, and they interfere destructively. But in our four-spin interaction model, the Fourier modes are asymmetrical in k-space. As a result, even if there are a large number of maxima and minima interfering with one another with increase in time, some parts of the contributing modes still cannot interfere in a completely destructive manner. Hence, temporal oscillations are always present in the system. Secondly, the relation L = e −C N is valid for all times and in all the regions of the phase diagram. This is because the term sin 2 (2Ω k ) << 1 everywhere in the phase diagram and even on the critical lines h 1,2,3,13 , unlike in [24], where the relation is found to be invalid for the critical lines or when the quench protocol is on the critical lines.
We have also performed the analysis regarding the role of time-dependent NC and LE to identify the critical lines. For our model of Eq. (2), we find that there is no non-analyticity or discontinuity in these quantities, whenever the phase change occurs. This reflects the fact that, like the NC, these two dynamical quantities cannot signal any particular characteristics of critical behaviour of the model.

V. MULTIPLE QUENCH SCENARIOS
In this section, we consider multiple quench scenarios in our model, extending the analysis of section IV. We will consider here a series of four sudden quenches in the model. The multiple quench protocol is as follows. At time t = 0, the external magnetic field h is quenched from h 0 to h 0 + δ for a time t 1 = 15s, where δ is the interaction strength. At time t 1 , the external magnetic field after quench i.e., h 0 + δ is changed back to h 0 for a time t 2 = 15s, which is defined as no quench in the system. This sequence has been repeated four times. The multiple quench protocol involves the interaction Hamiltonian to be H multi I , so that similar to Eq. (18), the total Hamiltonian is H multi where we define Here, the quench time is taken as T = 15s, Θ is the Heaviside function, and |e is the excited state of the twolevel system as discussed in the single quench scenario. The ground state of the system and the time-dependent NC after time t = T have been discussed in details in the single quench case. The ground state after time t = 2T i.e., when h 0 + δ is changed back to h 0 is where we have defined with the operators d k,1 and d k,2 being the Fourier operators, and θ k , u k , and v k are defined earlier with the arguments appropriately shifted.
To compute the complexity, the reference and the target states are chosen to be |Ψ h,J3 and |Ψ e (t 1 , t 2 ) respectively. The final result of NC after t = 2T is to handle analytically, and we have performed numerical analysis for the same. The behaviour of the NC with time after four quenches is shown in Fig. (9), for system size N = 501, J 3 = 1.5, and T = 15s. The time evolution of the NC exhibits an oscillating behaviour whenever the system is quenched, and the magnitude of oscillations decreases with time. Also, the magnitude of NC is found to be constant whenever h 0 + δ is changed back to h 0 . In Fig. (9), the blue curve corresponds to h 0 = 1 and δ = −0.4; the orange curve corresponds to h 0 = 1 and δ = −0.2, and the green curve corresponds to h 0 = 3 and δ = 0.2. Note that the green curve is plotted by choosing parameter values such that the system comes on the critical line h 1 after the quench. As concluded in the single quench case, the NC in multiple quench scenarios also does not indicate any critical behaviour of the system. In Fig. (10), we plot the NC at large The blue and orange lines correspond to system sizes N = 501 and 1001, respectively. times for h 0 = 1, and δ = −0.2, for system size N = 501 (blue) and N = 1001 (orange). We observe here that for late times and for large system sizes, in our four spin interaction model, oscillations are regular, and their maximum amplitude decreases initially with time and later acquires nearly a constant value. These results are similar to what was observed in [39,40] for integrable models.
We have also performed the same analysis of multiple quenches in the transverse spin-1/2 XY model [1], where the quench parameters are the anisotropic coupling and transverse magnetic field. When we quench the value of the magnetic field, the behaviour of the NC is similar to the one we encountered before -like in the four-spin exchange interaction model, the NC here first increases linearly and then starts oscillating over time. The oscillations are more rapid when the system is quenched on one of the critical lines. However, no other signal of QPT is observed from the complexity versus time curve, even for the transverse spin-1/2 XY model.

VI. THE DMRG AND ENTANGLEMENT ENTROPY
We now compute the entanglement entropy of the system using DMRG methods, which we first recall. The DMRG is a variational method that uses the matrix product states (MPS) and matrix product operator (MPO) representations to find the lowest energy state of a quantum many-body system. An MPS represents a quantum state with a product of matrices. The Hilbert space of the one-dimensional spin-1/2 chain has 2 N degrees of freedom for N spins. The MPS can be manipulated to form a subspace of the larger 2 N dimensional Hilbert space which captures the essential physics. The search for the ground state using DMRG is an iterative solution of the Schrodinger equation H |ψ = E |ψ . In the first step, the Hamiltonian of the system has to be constructed as an MPO. To solve the Schrodinger equation for |ψ , we start with a guess of the MPS. The DMRG algorithm minimises the ground state energy in the space of the MPS, where matrix elements are treated as variational parameters. The minimisation is performed by a sweeping procedure where the matrix state at a site is optimised at one time, keeping all others fixed, then optimising the next matrix, and so on. When all matrices are optimised, we again sweep back and forth until convergence is achieved. The convergence is guaranteed as the energy goes down at each iteration step. The DMRG uses the entanglement properties of a bipartite system to produce a state with minimal energy. In order to quantify the entanglement or quantum information, the most common measure is the Von Neumann entanglement entropy denoted by S. It is defined in a standard fashion, by bipartition of a composite system into two subsystems, A and B, where ρ A and ρ B are the reduced density matrices of subsystems A and B respectively.

A. EE with three-spin interactions
We will numerically calculate the ground state EE in the simplest case µ 1 = µ 2 = 1, J 13 = J 23 = J 3 , and J 14 = J 24 = 0 using DMRG method. We split the spin chain of N cells at the centre into two subsystems. Using singular value decomposition, the MPS is divided into left and right bipartition. Then, the DMRG calculations are performed using the library [37] on N = 51 cells with maximum bond dimension χ = 300. The numerical results of EE as a function of J 3 is plotted in Fig. (11), for H = 0, J 1 = 1.2, and J 2 = 0.8. When H = 0, there are two QPTs. One occurs at J 3 > (J 1 − J 2 ) = 0.4, which is the onset of the first-order QPT. Another occurs at J 3 = J 1 + J 2 = 2, when the line of second-order QPT, i.e., H c2 intersects the first-order QPT line (J 3 > 0.4, H = 0). From Fig. (11), we see that the EE shows a sharp jump whenever the critical point is approached, irrespective of the order of the QPT.

B. EE with four-spin interactions
Similar to the case of three-spin interactions, we have also plotted the numerical results of EE of our model Eq.
(2) as a function of J 3 in Fig. (12). The blue and red curves in Fig. (12) correspond to h = 0.25 and h = 2.5, respectively, for a fixed value of J = 1. When h = 0.25, there are two critical points at J 3 = 0.6958 and J 3 = 1.0958. When h = 2.5, the critical point is observed at J 3 = 0.9753. The EE clearly shows sharp jumps to indicate these second-order QPTs.
We conclude that for both the first and second-order QPTs, the EE shows a discontinuous nature when the parameters takes their critical values. For the first-order QPT, the result is in agreement with [38], in which the first-order QPT is associated with the discontinuity in the density matrix elements. However, the discontinuity in EE for the second-order QPT in our case is contrary to what is observed in [38], which suggests that the discontinuity or divergence should be present in the first derivative of EE for second-order QPTs, instead of EE itself. The reason for this result is that our model behaves differently from the other known exactly solvable spin chains. As we have mentioned before, there are two branches of the energy spectrum, E k,1,2 . Depending on the magnetic field h and the three spin interaction parameter J 3 , these branches completely or partially lie above or below the zero level, resulting in different magnetisations m z . The different phases in the ground state are thus determined by the values of m z = 0/1 (antiferromagnetic/ferromagnetic) and 0 < m z < 1 (ferrimagnetic). Hence all the phase diagram characteristics of the model are determined by the signs of E k,1,2 , which can be positive, negative, or zero. It is important to note that this behaviour is different from spin models with NN interactions, in which QPTs usually occur whenever the energy gap between the branches of energy spectra vanishes. For those models, the discontinuity is present in the derivative of EE for second-order QPTs, as concluded in [38]. But in our case, for different phases in the phase diagram, the ground state wavefunction |Ψ takes distinct values, as the summation over the Fourier modes k runs over specific ranges. The reduced density operator ρ A,B in Eq. (29) is the partial trace of the outer product of the ground state wavefunction |Ψ . So it becomes discontinuous whenever the phase change occurs, irrespective of the order of QPTs. The DMRG, a numerical technique for finding the ground state energy and wavefunction, confirms our analytical results of section II by showing the discontinuous nature of EE for both the first and second-order QPTs in Figs. (11) and (12).

VII. CONCLUSIONS AND DISCUSSIONS
In this paper, we have carried out a comprehensive analysis of information theoretic geometry in a spin model with three and four spin interactions along with an alternating NN coupling. In particular, we computed the quantum information metric and the related Fubini-Study complexity, the Nielsen complexity in static as well as quench scenarios and the entanglement entropy in these models. The motivation and importance of this study lies in the fact that the nature of phase transitions here are very different from the usual ones in models with NN interactions, such as the transverse field XY model, and it is important and interesting to quantify information theoretic measures and contrast these with the ones arising out of NN interactions.
Here, we find a number of distinctive features, which are very different from ones in models with NN interactions. The derivative of the NC does not diverge across a continuous phase transition, and the Fubini-Study complexity also does not show any special behaviour across two ferrimagnetic phases. The entanglement entropy shows a discontinuity across both first and second order phase transitions. These are in sharp contrast to the ones reported in the literature, in models with only NN interactions. Our analysis points to the fact that several features of information geometry that are valid in NN interaction models cease to be valid when three and four spin interactions are included.

Appendix A
In this appendix, we list the transformations to diagonalise the Hamiltonian of Eq. (2).