Dome-shaped phase diagram in the spin-1 XY ferromagnet with biquadratic exchange and longitudinal easy-axis crystal field

We investigate the phase diagram of a spin-1 ferromagnetic XY model in the presence of a longitudinal easy-axis crystal field assuming bilinear (J) and biquadratic (I) exchange interactions between nearest neighbors spins and using the two-time Green Functions framework at the level of the Devlin strategy. Employing both analytical estimates and numerical calculations, we find that the structure of the crystal-field-induced phase boundary changes sensibly as the ratio α=I/J\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha = I/J$$\end{document} increases. In particular, when α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} overcomes a characteristic value α∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha ^*$$\end{document}, two quantum critical points appear which are connected by a dome-shaped critical line. Due to the paradigmatic nature of the anisotropic spin model here considered, we believe that our findings may provide useful insights into the physical origin of recent experimental results found for some innovative materials which exhibit two quantum critical points and dome-shaped phase diagrams induced by non-thermal control parameters driving a non-conventional quantum criticality.


Introduction
Quantum anisotropic spin models, in several situations where strong quantum fluctuations and anisotropic parameters compete, are candidates to exhibit new quantum phases of different nature and hence possible quantum phase transitions (QPTs) driven by variation of appropriate non-thermal parameters [1][2][3][4][5][6][7]. Here, quantum fluctuations are generally introduced by the non-commutativity of quantum spin operators and their competitions with non-thermal parameters arise from different interactions in the Hamiltonian, as ferromagnetic bilinear and biquadratic exchange (BQE) couplings, external magnetic fields and crystal-field anisotropies [8][9][10][11].
The recent discovery of spin-1 magnetic compounds [12,13] has led to a lot of experimental and theoretical studies involving also the effects of single-ion anisotropy (SIA). This feature has produced a strong interest in the study of quantum criticality in spin-1 models with easy-plane and easy-axis SIA (see, for instance, refs. [14] and [15]). Up to now, extensive studies have been devoted to the conventional situations where a single QCP takes place [16][17][18][19][20][21][22][23][24].
This relevant feature suggests that it may be very interesting to study analytically and numerically the thermodynamic of appropriate spin-1 models to provide insights into the basic physics for systems showing two QCPs. Bearing this in mind, in the present paper, we focus on a spin-1 ferromagnetic XY model described by the Hamiltonian where <i,j> denotes a sum over nearest neighbor lattice sites, J and I are the ferromagnetic bilinear and biquadratic exchange interactions between nearest neighbor spins, h is the longitudinal magnetic field and D measures the intensity of a longitudinal crystal field that gives rise to an easy-axis (D > 0) or easyplane (D < 0) SIA in the z-direction. In Eq. ined in their full aspects by variation of all the competing parameters (J, I, h, D). Several investigations have been performed only for particular values of the interaction parameters providing partial information about its thermodynamics. The simplest case with I = 0 in the absence of SIA, i.e., the conventional transverse XY model [16,18,43] has been extensively analyzed in the context of magnetic-field-induced quantum criticality, where QPTs or QCPs are driven by the longitudinal magnetic field h, which activates quantum fluctuations in the low-temperature regime. In the last decade, a lot of interest has been attracted by the regimes with h ≥ 0 and different signs of the SIA parameter D. In this situation, it is found [15,44,45] that the nature of SIA has a significant effect on the h-induced phase diagram and quantum criticality in the (h, T )-plane. Specifically, in the easy-plane case the quantum critical lines have the conventional shape merging into a QCP (h c (D), T = 0) for appropriate values of D with a standard quantum critical thermodynamics. The situation appears quite different when easy-axis SIA is present. The easy-plane case (D < 0) is not so surprising, except for the existence of a hinduced QCP in a limited range of the SIA parameter, producing a quantum critical scenario qualitatively similar to that found in the absence of the crystal field anisotropy (of course, the non-universal critical quantities depend on the control parameter D). The easy-axis regime is certainly more interesting due to the emergent reentrant critical lines merging into a QCP for selected values of the SIA parameter [16][17][18]. The peculiar reentrant behavior in the phase diagram has been observed for superconducting vortices [46,47], liquid crystals [48], polymeric materials [49], ferromagnetism in semiconductors [50], solid hydrogen [51], denaturation of DNA [52] and has stimulated several theoretical investigations [15][16][17]20,44,[53][54][55][56][57][58][59]. Therefore, a complete study of the paradigmatic model (1.1) with all competing interactions may be very interesting not only from a theoretical perspective.
The aim of the present paper is to investigate the structure of the phase diagram of the model (1.1) taking into account the competing effects of the crystal field and the BQE interaction by using the two-time Green functions framework according to the Devlin scheme [60] which is a reliable approach to treat the effects of SIA.
The paper is organized as follows. In Sect. 2, we present the basic equations obtained following the Devlin strategy for the spin-1 model (1.1) in the absence of the longitudinal magnetic field. Section 3 presents the stability conditions starting from the transverse susceptibility for paramagnetic phase. Analytical estimates for dimensionalities d ≥ 3 allow us to select two regimes of the BQE ratio α = I/J, where one or two D-induced QCPs occur in the (D, T )-plane. In the two-QCPs regime, low-temperature predictions suggest the possible occurrence of a dome-shaped phase diagram. Numerical results are reported in Section 4 which support the expected existence of two D-induced QCPs and related dome-shaped phase diagrams increasing the ratio α above a well defined crossover value α * ≈ 1.8 for cubic lattices. Section 5 closes the paper with concluding remarks.

Two-time Green functions framework
Although no ordering occurs in the z-axis, the Devlin scheme can be conveniently applied for a non conventional description of the magnetic properties of the model (1.1) in terms of the ensemble average of the operator (S z i ) 2 even in the absence of the longitudinal magnetic field. To this aim, we put h = 0 in what follows and introduce the retarded Green functions (GFs)

1)
where S ± i = S x i ± iS y i are the spin raising and lowering operators, θ(x) is the step function, X(t) stands for the equilibrium average at temperature T for a generic operator X, i.e., with the Boltzman constant assumed equal to unity, and Following standard procedure, we obtain the equations of motion (EMs) for the time Fourier transform and (2.6) Here It is evident that the exact two coupled EMs (2.5) and (2.6) do not constitute a closed system of equations in the required GFs G ij (ω) and Γ ij (ω) due to the appearance of higher order exchange and crystalfield anisotropy-like GFs. Thus, one is forced to introduce appropriate decoupling schemes. Here, we follow the Devlin idea [60] using Random Phase Approximation (RPA)-like decouplings for the exchange GFs and treating exactly the crystal-field-like ones. Specifically, we adopt the following RPA scheme [10,11] which neglects the correlations between the transverse and longitudinal spin components at different lattice sites. In contrast the SIA-like terms will be treated exactly. Then, with the further Fourier transforms in , the wavevector k ranging in the first Brillouin zone (1BZ), the problem reduces to a simple closed set of two algebraic equations for G( k, ω) and Γ( k, ω) with solutions defines the dispersion relation of the spin excitations.
In these equations, one must have (1 − qJ( k) D ) ≥ 0 for assuring the reality of ω( k) and to simplify the notation we have defined Finally, the spectral theorem [10,11] for two generic operators A and B . With standard algebra one has, for instance, In particular in the static case at the same site, one has From previous equations, it is evident that the unknown thermodynamic quantity q = q z , as a function of the intensive parameters (T, D), plays a key role. Indeed, all the thermodynamic properties can be determined if q = q(T, D) is known, in analogy with the role played by the longitudinal magnetization m = S z in models where it is non-zero. A self-consistent equation for q can be derived combining the kinematic identity valid for spin S = 1 and the spectral theorem. The result is the selfconsistent equation in which q must satisfy the boundary condition −2 ≤ q ≤ 1. If in this equation we express ω( k) and D in terms of q (see (2.10) and (2.11)) we obtain the required self-consistent "Devlin equation" is the "structure factor" and, for short notation, At zero temperature Eq. (2.16) reduces to since the dispersion relation ω( k) is always a nonnegative quantity.

Transverse susceptibility, stability condition and crystal-field-induced quantum critical points
In our scheme, the physical solution q of the Devlin equation and the magnetic properties of the model can be obtained through the transverse susceptibility χ ⊥ = −G( k = 0, ω = 0) which is given by where provides the energy gap of the spin excitations. By inspection of Eq. (3.1), it is evident that the thermodynamic stability condition χ −1 ⊥ ≥ 0 is assured by the simple requirement 0 < Q ≤ 1. From this condition the following basic features emerge: i) being Q a definite positive quantity, q and D must have always the same sign; ii) the reality condition for the energy spectrum is automatically satisfied through the stability condition, since γ( k) ≤ 1; iii) the equation Q = 1 defines the phase boundary where χ ⊥ → +∞ and the energy gap ω(0) = |D|(1 − Q) 1/2 vanishes. Thus, this criticality condition locates possible critical points for appropriate space dimensionalities.
It is important to emphasize that, in contrast to q and D, the physical SIA parameters D and q do not have always the same sign and one can have a negative q even in the easy axis regime (0 < D < I(0) 6 (4 − q)). Exploiting the above feature (iii), the physical solution of the Devlin equation (2.16) at the criticality is determined setting Q = 1, thus yielding the selfconsistent equation that, at T = 0, reduces to Here F (−1/2) = 1 N k 1/(1 − γ( k)) 1/2 is the particular case n = −1/2 of the so called "structure sums" In the limit of infinite lattice (for N → ∞) this sum reduce to an integral which converges for any space-dimensionality d > 1. For three-dimensional lattice: in the simple cubic (sc), F (−1/2)) = 1.11536, in the body-centered cubic (bcc) it is F (−1/2)) = 1.08731 and in the face-centered cubic (fcc), F (−1/2)) = 1.07341.
Equation (3.4) can be simply solved yielding the global representation from which we see that q depends on the SIA parameter D and on the BQE interaction I(0) only through signD.
The solution corresponding to signD = +1 is that is a positive solution giving rise to the transverse autocorrelation function which measures the out-plane quantum spin fluctuations. Correspondingly, the T = 0 critical value of D is having defined α ≡ I J as the relative intensity of the BQE parameter. The other solution for q, corresponding to signD = −1, is the out-plane quantum spin fluctuation now is and the related quantum D-coordinate is Notice that D c2 can be positive or negative depending on whether α is greater or less than α * ≡ 2F −1 (−1/2), respectively. Previous T = 0 analysis suggests that, for α < 2F −1 (−1/2), our spin model exhibits a D-induced QCP (D c1 ) located in the easy-axis region D c1 > 0 and another QCP (D c2 ) located in the easy plane region D c2 < 0. On the contrary, if α > 2F −1 (−1/2) both critical points are present in the easy axis region. Observe that D c2 < D c1 for all the allowed values of α and at the crossover point α * one has D c2 = 0 and D c1 = The structure of the phase diagram of the model (1.1) at our level of approximation can be deduced inserting the low temperature solution q(T ) of the Devlin equation (2.16) in the critical condition Q = 1 obtaining in this way the locus of critical points ending in D c1 and D c2 . We find in this way a combination of two temperature depending branches, one driven by q(T ) > 0 and the other by q(T ) < 0, which merge at q = 0, corresponding to D = 2 3 I(0) for T c /J(0) = (2/3)F (−1). A conventional critical scenario is expected to occur for BQE ratio α < 2F −1 (−1/2), when a single QCP exists in the D, T -plane. In this low-BQE regime, the solution q = q(T ) close to the QCP (D c1 , T = 0), with q c1 and D c1 given by Eqs. (3.6) and (3.8), can be determined setting signD = +1 in Eq. (3.3). At low temperature, by iteration with initial point q c1 > 0, we get where the integral can be estimated in a standard way for d > 2 assuming in the integrand a small-k expansion depending on the coordination number z of the lattice γ( k) 1 − 1 z k 2 , so the longitudinal quadrupolar parameter can be expressed as follows (3.13) where K d = 2 1−d π −d/2 /Γ(d/2) and Γ(x) and ζ(x) denote the Euler gamma and the Riemann zeta functions, respectively. With this ingredients, it is immediate to obtain the required asymptotic equation of the critical line using the criticality condition. We get with D c1 − 2 3 I(0) > 0. Let us consider now the high-BQE regime α > 2F −1 (−1/2) where two quantum critical points exist. Their (T = 0)-coordinates D c1 and D c2 are given by Eqs. (3.8) and (3.11) with correlated q values (3.6) and (3.9), respectively. For the branch of the critical line ending in the first QCP (D c1 , T = 0), we obtain the expression (3.14), as in the previous regime. Similar to this case, for the second QCP, the low-T longitudinal quadrupolar parameter reads Then, the critical line equation close to the QCP (D c2 , T = 0) is We observe that both the critical lines are characterized by the same shift exponent ψ = d−1 typical of the continuous O(2)-vector TIM-like models [4,5,61] with dynamical critical exponent z = 1 as for the easy-plane and easy-axis SIA cases when α < 2F −1 (−1/2). More interestingly, they exhibit a different orientation of the curvature in the (D, T )-plane (on the right for D c2 (T ), on the left for D c1 (T )). Previous estimates provide preliminary indications towards a peculiar phase diagram of the model when α > α * .
To have additional information about the structure of the phase diagram we provide below analytical estimates of the critical temperature T c as function of D for different values of the BQE coupling in some asymptotical regimes. As a first step, we solve by iteration the Devlin equation at criticality for D = 0 and therefore we must assume in (3.3) signD = −1, hence q(T ) < 0. So the equation to solve reads Focusing on the low-T regime a low-wave-vector expansion in the integrand of (3.17) yields Then, the related equation of the critical line follows from the critical condition Incidentally, we observe that, since it must be q(T ) < 0 for internal consistency, this equation requires that α/6 < 1 at our level of approximation. Finally, the low-T estimate of the critical temperature as function of the biquadratic parameter is (

3.21)
Notice that this is a decreasing function that goes to zero as α → α * , with α * the crossover value between the two BQE regimes. It is possible to estimate T c (α) also at finite temperature by applying to Eq. (3.3) the Laurent expansion in terms of the Bernoulli numbers B 2n . Bearing in mind the previous iterative procedure with starting point where a n = B2n (2n)! F (n − 1)F (−1) 2n−1 . In particular to second order in α, we have (3.24) i.e., T c first increases with α (from α = 0) and then decreases. The behavior is confirmed by exact numerical calculations (see next section). Since T c (α) is the initial critical temperature of any critical line in the (D, T )plane, its behavior with α will play an important role to determine the structure of the phase diagram. Finally, in this section, we obtain an estimate of the finite temperature T c (D) when a single QCP exists. One first solves approximatively the Devlin equation (3.3) to obtain q(T ) then, using the critical condition on q(T ), a conventional iterative procedure leads to with 0 ≤ D/J(0) < 4. Notice that at D = 0 this equation reproduces T c (α), i.e., the starting critical points in the absence of SIA. Previous analytical estimates do predict a structure of a phase diagram appearing as a sequence of "dome shaped" critical lines increasing α and this is confirmed by numerical calculation.

Numerical data and dome-shaped structure of phase diagram
This section reports numerical data regarding the phase diagram of the spin-1 XY model with BQE interactions and easy-axis SIA obtained by means of the solution of the Devlin equation at criticality for the out-plane quadrupolar parameter q. Our primary interest is to derive a reliable numerical scenario about the structure of the phase boundary in the (D, T )-plane by variation of the reduced parameters D/J and α = I/J. In particular, we will focus on the appearance of partial and full dome-shaped critical lines increasing the BQE parameter α, as suggested by previous analytical estimates. With this plan in mind, the starting step is to solve numerically the general Devlin equation (2.16) under the stability condition 0 < Q = J(0)q/[D − I(0)(4 − q)/6] ≤ 1 which is the basic ingredient to explore the thermodynamics of our spin model in the disordered phase.
In Fig. 1, we have plotted the parameter q as a function of reduced temperature T /J by variation of D/J for α = 0.5 (panel a) and α = 2 (panel b), which are two representative values in the BQE regimes α < α * and α > α * , where a single QCP and two QCPs are expected to occur, respectively. In this context, we have estimated analytically α * = 2F −1 (−1/2) = 1.793143 for the simple cubic (sc) lattice. In both panels, the black solid line provides the parameter q evaluated on the phase boundary (where Q = 1(χ −1 = 0)), which allows to analyze numerically the phase diagram of the system. From this figure it is evident that, consistently   Figure 2 presents the critical lines in the (D, T )-plane by variation of the BQE ratio α for simple cubic lattice.
Here it is shown that, increasing α, a different structure of the phase boundaries takes place. For 0 ≤ α ≤ α * (α * depends on the coordination number) each critical line terminates in a single D-induced QCP with a incomplete dome-shaped behavior (see, for instance, panels a-d). Further increasing α, a non-conventional regime occurs with dome-shaped critical lines bounded by two QCPs (panels e-f). More specifically, for the BQE ratio in the interval 0 ≤ α ≤ 1, usually studied in literature in the absence of SIA [23,[62][63][64], the cubic lattice exhibits a single QCP with conventional critical lines. However, our analysis indicates that this feature persists also for 1 ≤ α ≤ α * . The situation gradually changes in the BQE regime with α ≥ α * displaying dome-structure critical lines bounded by two QCPs (panel f). It is also observed that the dome-amplitude of different critical lines decreases for increasing values of the BQE ratio α.
It is worth emphasizing that the reduction of T c (α) for increasing α ≥ 0 is crucially involved in determining the structure of the phase diagram. In Fig. 3, we have plotted T c as a function of α for a sc lattice, for four values of the SIA parameter D. Here it is shown that, for D = 0, T c (α) vanishes exactly at α = α * , which marks the crossover between the regimes with presence of one and two QCPs.
A scenario quite similar to the phase diagram found before for a sc lattice takes place also for bcc (z=8) and fcc (z=12) lattices.
To analyze the physical mechanism producing the dome-shaped phase diagram scenario in our anisotropic XY model beyond the Devlin-like framework is a prohibitive issue especially in the absence of reliable information about its ground state. However, at our level of approximation, the emerging dome-shaped phenomenon can be qualitatively traced back to the conflicting competition between the easy-plane BQE coupling and the transverse (out-plane) easy-axis SIA, where the latter tends to inhibit the in-plane ordering. Hence, some qualitative insights can be obtained analyzing the effects of the aforementioned competition. Increasing the crystal field, the z-component of the spins is enhanced implying a reduction of the in-plane order and hence a lowering of the transition temperature to disordered phase. The physical situation is drastically different from the case of an easy-plane crystal field which favors the in-plane order [18,43,44]. Here, the z-component of the spins decreases as the crystal field is enhanced implying a stronger ordering in the XY-plane and hence a higher transition temperature to disorder. Previous mechanism may provide a simple explanation of the appearance of disordered regions which are not present in the conventional phase diagrams in the absence of SIA [42,[65][66][67].
Bearing all this in mind, some insights can be provided on the possible physical origin of the domes in the phase diagram by observing the evolution in Fig. 2 from panel (a) to panel (f). Panel (a) depicts a typical critical line in the absence of BQE interaction (α = 0), which merges into a QCP arising from the quantum fluctuations produced by the presence of SIA. The BQE coupling starts to play a sensible role increasing the ratio α as shown through the panels (b-e). As expected from the analytical estimates, for fixed α < α * ≈ 1.8 the critical lines exhibit a single D-induced QCP and a partial dome-shaped behavior which becomes more and more marked for small D from the panel (b) to panel (e). We argue that this feature is imputed to the more stabilized in-plane ordering produced by the increased BQE ratio up to the crossover value α * = 1.8 between the one-QCP regime and the two-QCPs one (panel e). When α = α * , with T c (D = 0, α * ) = 0, the system is in a state with an incipient QCP (D c2 = 0) and an effective one (D c1 = 0). In the panel (f), we present a full dome-shaped critical line typical of the regime α > α * , where two D-induced QCPs occur. Here, the critical line shows a peculiar behavior which can be interpreted as follows. Since the easy-axis crystal field fluctuations tend to inhibit the in-plane ordering, it is reasonable to conjecture that for D c2 < D < D * (q < 0 ) the BQE coupling dominates and the in-plane ordering becomes active producing an increased critical temperature up to the a maximum T * c (D = D * , q = 0) signaling the balance of the two effects. For D * < D < D c1 (q < 0) quantum SIA fluctuations start to produce a sensible reduction of the XY order and increasing D the critical temperature decreases and vanishes at D = D c1 . Finally, no critical singularity occurs for D > D c1 .
We conclude this section showing some numerical data for thermodynamics of paramagnetic phase obtained using the solution of the general Devlin-like equation under stability condition 0 < Q < 1. In this context, the basic quantity is the transverse susceptibility χ ⊥ (T, D; α).
In Fig. 4, we have plotted the inverse susceptibility in the paramagnetic phase as a function of the easy-axis SIA parameter D for different isotherms and three distinct values of the BQE ratio α. These figures show that, when two D-induced QCPs are present (see Fig. 4c) the transverse susceptibility has a nonstandard behavior: at fixed temperature T it increases with D < D c2 (T ) by approaching the low-D branch of the critical line and diverges as D → D − c2 (T ) while along the high-D branch it diverges as D → D + c1 (T ) and is a decreasing function of D for D > D c1 (T ). Besides, along the D-axis it intercepts in two distinct critical values of D for T < T * c , similarly to a typical reentrant phase diagram [15][16][17]44]. Additional information can be extracted from Fig. 5, where the transverse susceptibility is plotted as a function of T at fixed D for selected values of α in the paramagnetic phase around a typical critical line. It is evident that its behavior, as compared to the corresponding one in Fig. 4, appears more complex (panels b-c) especially in the regime where two QCPs are present (panel c). Notice that, in all cases (panels a-c), the susceptibility tends to zero increasing T and diverges at the QCP as T −1 , as expected. Previous numerical data suggest also that, in all the regimes χ ⊥ decreases increasing the temperature at fixed D. Fig. 4 Inverse susceptibility χ −1 as a function of the single-ion-anisotropy parameter D for different isotherms at fixed values of the biquadratic ratio α. In a, we have α = 0, showing a conventional situation; in b, it is α < α * , so we are in a regime with one-QCP, but the inverse susceptibility have a non-monotonic behavior; in c, we have α > α * ; hence, there is a characteristic dome-shaped critical line with two-QCPs At this stage, using the general solution of the Devlinlike equation (2.16) and the numerical results for transverse susceptibility it is possible to derive information about other thermodynamic quantities in the paramagnetic phase (as internal energy, specific heat, and so on) by employing the basic formalism of the powerful two-time GF method and, in particular, the spectral theorem.

Concluding remarks
In this paper we have studied analytically and numerically the phase diagram of the spin-1 ferromagnetic XY model with BQE coupling (I) and longitudinal easyaxis crystal field (D) on cubic lattice. Following the Devlin strategy [60], in the EMs we have employed RPA decouplings for exchange higher-order GFs and treated exactly the crystal-field anisotropy terms. The main results are as follow: (i) for BQE ratio α = I/J below a certain crossover value α * a single D-induced QCP exists in the (D, T )-plane with a conventional phase diagram quite similar to that found in the absence of biquadratic interactions; (ii) increasing the BQE parameter the phase diagram is sensibly modified, with the appearance of critical lines showing a partial domelike structure with a single QCP. Remarkably, with a further increase of the BQE ratio, above a thresh-old α * , two D-induced QCPs appear delimiting a full dome-shaped critical line. Notice that, although explicit numerical results were obtained for simple cubic lattices (z = 6), a similar scenario can be easily derived for bcc (z = 8) and fcc (z = 12) lattices. To our knowledge, no theoretical studies have been performed in the past for investigating the non-conventional conflicting competition which may occur in simple anisotropic XY model between BQE and the quantum fluctuations induced by an easy-axis longitudinal SIA. At this stage, one can object that the above emerging non-conventional two-QCPs scenario is simply a product of the Devlin-like framework adopting the simple decoupling RPA scheme used in the EMs for the exchange higher-order GFs.
Theoretical improvements for testing this criticism can be obtained using Callen-like decouplings [68] for the exchange terms [21,22,24] providing an estimate of the correlation effects between spins in different sites. However, as suggested by previous general frameworks, this higher-order level of approximation, especially when in-plane BQE interaction and longitudinal easy-axis crystal field are simultaneously taken into account in the EMs, should complicate sensibly the mathematical treatment without a substantial change of the predictions arising from the simplest and controllable RPA scheme used through the present paper. So, we speculate that our approach is able to capture, at least qualitatively, the basic physics which governs the two-QCPs scenario with the related dome- shaped phase diagrams. In any case, the study of the competition between the contrasting BQE interaction and the longitudinal easy-axis SIA in our paradigmatic anisotropic spin-1 XY model (1.1) in h = 0 becomes, of course, a subject of experimental interest and may stimulate further theoretical, numerical and experimental researches.