Abstract

The mixed spin-(3/2, 9/2) Blume-Capel ferrimagnetic model is investigated on the Bethe lattice using exact recursion relations (ERR). Ground-state phase diagrams in the $\left(D_A/q\left|J\right|,D_B/q\left|J\right|\right)$ plane reveal ten competing ferrimagnetic configurations and seven multicritical points—two more ground-state regions than in the previously studied (3/2, 7/2) case on the same lattice. Temperature-dependent phase diagrams are constructed in the $\left(D_A/\left|J\right|,kT/\left|J\right|\right)$ , $\left(D_B/\left|J\right|,kT/\left|J\right|\right)$ planes and in the isotropic case DA = DB = D, for coordination numbers q = 3, 4, 5, 6. The system exhibits both first- and second-order phase transitions and compensation temperatures for appropriate crystal-field values. Under an external magnetic field, the model displays single, double, and triple hysteresis loops; triple loops occur for $?0.6\le J/\left|J\right|<?0.5$ . The ten competing spin states of sublattice B (versus eight for spin-7/2) generate nine independent recursion ratios, a richer fixed-point structure, systematically higher critical temperatures, and a broader multi-hysteresis parameter range—features that are not recoverable by simple extrapolation from the lower-spin case. Results are compared with mean-field studies of the same spin pair on the simple cubic lattice and with the (3/2, 7/2) Bethe lattice.

Keywords

Blume-Capel Model, Mixed Spin, Bethe Lattice, Exact Recursion Relations, Compensation Temperature, Multi-Hysteresis, Phase Transition

Share and Cite:

1. Introduction

The study of mixed-spin Ising models has attracted increasing interest over the past decades due to their much richer magnetic properties compared to single-spin systems [1] [2]. Mixed-spin ferrimagnetic models, in particular, provide a powerful theoretical framework for understanding ferrimagnetism in insulating materials, a field of both fundamental and technological importance, especially for applications in spintronics, data storage, and molecular magnetism [3]-[7].

Numerous theoretical works have investigated binary systems combining half-integer spins (1/2, 3/2, 5/2, 7/2, 9/2) and sometimes integer spins (1, 2, 3), using various approaches such as Monte Carlo simulations [1] [8]-[12], mean-field approximations [13]-[17], effective-field theory [18], renormalization group techniques [17], and exact recursion relations on the Bethe lattice [19]-[26]. Single-ion crystal fields, exchange interactions, and lattice geometry strongly influence phase transitions, compensation temperatures, hysteresis behaviors, tri criticality, and re-entrant phenomena.

Among these systems, mixed half-integer spin ferrimagnetic models have received particular attention due to the occurrence of compensation temperatures, where the total magnetization vanishes while sublattice magnetizations remain nonzero, and the presence of multiple hysteresis loops under an external magnetic field [10] [11] [13]. The simplest combinations, such as spin-1/2 and spin-3/2, have been extensively studied [2] [23], but they generally exhibit only second-order phase transitions. More complex pairs, such as spin-3/2 and spin-5/2 [21] or spin-5/2 and spin-3 [14], have revealed re-entrant compensation lines and first-order phase transitions. More recently, the mixed-spin (3/2, 7/2) Blume-Capel ferrimagnetic model was analyzed in detail using exact recursion relations on the Bethe lattice [27], showing multiple compensation temperatures, first- and second-order phase transitions, as well as single, double, and triple hysteresis cycles.

In parallel, mean-field studies on the simple cubic lattice have explored the mixed-spin (3/2, 9/2) Blume-Capel ferrimagnetic system, revealing compensation phenomena and multiple hysteresis loops, including triple-loop cycles [28] [29]. However, to the best of our knowledge, this spin pair has not yet been investigated using exact recursion relations on the Bethe lattice, which is exact on this tree-like topology and provides a rigorous description of local spin correlations not captured by mean-field approximations; we note that quantitative comparisons with results on the simple cubic lattice remain qualitative due to the different lattice geometry and coordination number.

Non-triviality of the spin-9/2 extension. Table 1 summarizes the key structural differences between the (3/2, 7/2) case [27] and the present (3/2, 9/2) model. The spin-9/2 value introduces ten distinct spin states on sublattice B (compared to eight for spin-7/2), yielding ten competing ferrimagnetic ground-state configurations, nine independent recursion ratios Y_l, and a qualitatively richer fixed-point structure. As we demonstrate below, this added complexity translates into physically observable differences: a higher number of multicritical points in the ground-state diagram, systematically higher critical temperatures for equivalent reduced parameters, and a broader range of ferrimagnetic coupling values that support triple hysteresis loops. These features cannot be inferred from the lower-spin case by simple extrapolation.

Table 1. Structural comparison between the (3/2, 7/2) [27] and present (3/2, 9/2) Blume-Capel models on the Bethe lattice.

The present work therefore examines the mixed half- integer spin (3/2, 9/2) Blume-Capel ferrimagnetic model on the Bethe lattice using exact recursion relations. We focus on ground-state and finite-temperature phase diagrams, compensation temperatures, and multiple hysteresis phenomena, with systematic comparisons to mean-field results [28] [29] and to the (3/2, 7/2) case [27].

The paper is organized as follows. Section 2 presents the model and the ERR method. Section 3 discusses phase diagrams, compensation temperatures, and hysteresis cycles. Section 4 provides conclusions and perspectives.

2. The Model and Hamiltonian

We consider a ferrimagnetic Blume-Capel Ising model with mixed spins $s_i=\pm 3/2,\pm 1/2$ on sublattice A and ${\sigma }_j=\pm 9/2,\pm 7/2,\pm 5/2,\pm 3/2,\pm 1/2$ on sublattice B, arranged on a Bethe lattice with coordination number q.

The Hamiltonian reads

$\mathscr{H}=-J{\displaystyle \sum }_{i,j}{s}_i{\sigma }_j-D_A{\displaystyle \sum }_{i\in A}{s}_i^2-D_B{\displaystyle \sum }_{j\in B}{\sigma }_j^2-h\left({\displaystyle \sum }_{i\in A}{s}_i+{\displaystyle \sum }_{j\in B}{\sigma }_j\right),$ (2.1)

where J < 0 is the antiferromagnetic exchange interaction, D_A and D_B are single-ion crystal-field anisotropies, and h is the external magnetic field. The Bethe lattice is a tree-like structure with no loops, where each site has q nearest neighbors; this topology allows exact calculations via recursion relations without approximations for q = 3, 4, 5, 6.

2.1. Partition Function and Probability Distribution

The partition function is $Z={\text{Σ}}_{conf}{\text{e}}^{-\beta H}={\text{Σ}}_{Spc}P\left(Spc\right)$ , where $\beta =1/\left(k_{B}T\right)$ . Following the standard ERR procedure [27], the unnormalize probability for a central spin $s_o$ on sublattice A is

$P\left(\left\{s_0\right\}\right)={\text{e}}^{\beta \left(D_A{s}_0^2+h{s}_0\right)}{\displaystyle {\prod }_{k=1}^q{Q}_n}\left(s_0|{\sigma }_1^{\left(k\right)}\right),$ (2.2)

where $Q_n\left(s_0\text{|}{\sigma }_1^{\left(k\right)}\right)$ is the partition function of the k-th branch of generation n starting from the nearest-neighbor spin ${\sigma }_1^{\left(k\right)}$ on sublattice B, computed recursively as

$Q_n\left(s_0{\sigma }_1\right)={\text{e}}^{\beta \left(J{s}_0{\sigma }_1+D_B{\sigma }_1^2+h{\sigma }_1\right)}{\displaystyle \prod }_{l=1}^p{Q}_{n-1}\left({\sigma }_1{s}_2^{\left(l\right)}\right),$ (2.3)

with $p=q-1$ branches per neighbor spin. Defining the branch partition function on sublattice A as $f_n\left(\left\{s_0\right\}\right)={\displaystyle {\sum }_{\left\{{\sigma }_1\right\}}{Q}_n\left(s_0|{\sigma }_1\right)}$ , Equation (2.3) yields

$f_n\left(s_0\right)={\displaystyle \sum }_{{\sigma }_1}{\text{e}}^{\beta \left(J{s}_0{\sigma }_1+D_B{\sigma }_1^2+h{\sigma }_1\right)}{\left[f_{n-1}\left({\sigma }_1\right)\right]}^p.$

Considering all ten spin values of sublattice B (${\sigma }_1\in \left\{\pm 9/2,\pm 7/2,\pm 5/2,\pm 3/2,\pm 1/2\right\}$ ), Equation (2.4) expands to

$\begin{array}{c}{f}_n\left(s_0\right)={\text{e}}^{\beta \left(J{s}_0\frac{9}{2}+D_B\frac{81}{4}+h\frac{9}{2}\right)}{\left[f_{n-1}\left(9/2\right)\right]}^p+{\text{e}}^{\beta \left(J{s}_0\frac{7}{2}+D_B\frac{49}{4}+h\frac{7}{2}\right)}{\left[f_{n-1}\left(7/2\right)\right]}^p\\ +{\text{e}}^{\beta \left(J{s}_0\frac{5}{2}+D_B\frac{25}{4}+h\frac{5}{2}\right)}{\left[f_{n-1}\left(5/2\right)\right]}^p+{\text{e}}^{\beta \left(J{s}_0\frac{3}{2}+D_B\frac{9}{4}+h\frac{3}{2}\right)}{\left[f_{n-1}\left(3/2\right)\right]}^p\\ +{\text{e}}^{\beta \left(J{s}_0\frac{1}{2}+D_B\frac{1}{4}+h\frac{1}{2}\right)}{\left[f_{n-1}\left(1/2\right)\right]}^p+{\text{e}}^{\beta \left(-J{s}_0\frac{1}{2}+D_B\frac{1}{4}-h\frac{1}{2}\right)}{\left[f_{n-1}\left(-1/2\right)\right]}^p\\ +{\text{e}}^{\beta \left(-J{s}_0\frac{3}{2}+D_B\frac{9}{4}-h\frac{3}{2}\right)}{\left[f_{n-1}\left(-3/2\right)\right]}^p+{\text{e}}^{\beta \left(-J{s}_0\frac{5}{2}+D_B\frac{25}{4}-h\frac{5}{2}\right)}{\left[f_{n-1}\left(-5/2\right)\right]}^p\\ +{\text{e}}^{\beta \left(-J{s}_0\frac{7}{2}+D_B\frac{49}{4}-h\frac{7}{2}\right)}{\left[f_{n-1}\left(-7/2\right)\right]}^p+{\text{e}}^{\beta \left(-J{s}_0\frac{9}{2}+D_B\frac{81}{4}-h\frac{9}{2}\right)}{\left[f_{n-1}\left(-9/2\right)\right]}^p.\end{array}$ (2.4)

Substituting $s_0=+3/2,-3/2,+1/2,-1/2$ into Equation (2.5) yields directly the four branch partition functions $f_n\left(\pm 3/2\right)$ and $f_n\left(\pm 1/2\right)$ used in the numerical fixed-point iteration. Their explicit forms follow immediately and are not reproduced separately, since each is a special case of Equation (2.5).

Similarly, the branch partition function for sublattice B is

$\begin{array}{c}{f}_{n-1}\left({\sigma }_1\right)={\text{e}}^{\beta \left(J{\sigma }_1\frac{3}{2}+D_A\frac{9}{4}+h\frac{3}{2}\right)}{\left[f_{n-2}\left(3/2\right)\right]}^p\\ +{\text{e}}^{\beta \left(J{\sigma }_1\frac{1}{2}+D_A\frac{1}{4}+h\frac{1}{2}\right)}{\left[f_{n-2}\left(1/2\right)\right]}^p\\ +{\text{e}}^{\beta \left(-J{\sigma }_1\frac{1}{2}+D_A\frac{1}{4}-h\frac{1}{2}\right)}{\left[f_{n-2}\left(-1/2\right)\right]}^p\\ +{\text{e}}^{\beta \left(-J{\sigma }_1\frac{3}{2}+D_A\frac{9}{4}-h\frac{3}{2}\right)}{\left[f_{n-2}\left(-3/2\right)\right]}^p.\end{array}$ (2.5)

Evaluating Eq. (eq:fn_sigmaB 2.6) for each of the ten spin states ${\sigma }_1\in \left\{\pm 9/2,\pm 7/2,\pm 5/2,\pm 3/2,\pm 1/2\right\}$ provides the ten branch partition functions $f_{n-1}\left({\sigma }_1\right)$ entering Equation (2.5). Again, each is a direct substitution into Equation (2.6) and the five resulting expressions are used directly in the iteration without being reproduced in full.

For sublattice A (spin-3/2), three independent ratios are defined with respect to the reference state $s_0=-1/2$ :

$X_1=\frac{f_n\left(+3/2\right)}{f_n\left(-1/2\right)},X_2=\frac{f_n\left(+1/2\right)}{f_n\left(-1/2\right)},X_3=\frac{f_n\left(-3/2\right)}{f_n\left(-1/2\right)}$ (2.6)

For sublattice B (spin-9/2), ten branch partition functions arise. Nine independent ratios are defined with respect to the reference state ${\sigma }_1=-1/2$ (note $Y_6\equiv 1$ by construction):

$\begin{array}{l}{Y}_1=\frac{f_{n-1}\left(+\frac{9}{2}\right)}{f_{n-1}\left(-\frac{1}{2}\right)},Y_2=\frac{f_{n-1}\left(+\frac{7}{2}\right)}{f_{n-1}\left(-\frac{1}{2}\right)},Y_3=\frac{f_{n-1}\left(+\frac{5}{2}\right)}{f_{n-1}\left(-\frac{1}{2}\right)},\\ Y_4=\frac{f_{n-1}\left(+\frac{3}{2}\right)}{f_{n-1}\left(-\frac{1}{2}\right)},Y_5=\frac{f_{n-1}\left(+\frac{1}{2}\right)}{f_{n-1}\left(-\frac{1}{2}\right)},Y_6\equiv 1,\end{array}$ (2.7)

$\begin{array}{l}{Y}_7=\frac{f_{n-1}\left(-3/2\right)}{f_{n-1}\left(-1/2\right)},Y_8=\frac{f_{n-1}\left(-5/2\right)}{f_{n-1}\left(-1/2\right)},\\ Y_9=\frac{f_{n-1}\left(-7/2\right)}{f_{n-1}\left(-1/2\right)},Y_{10}=\frac{f_{n-1}\left(-9/2\right)}{f_{n-1}\left(-1/2\right)}.\end{array}$ (2.8)

Numerical implementation and reproducibility. All fixed-point equations are iterated in the log-ratio representation $\ell X_k=\ln X_k$ ($k=1,2,3$ ) and $\ell Y_l=\ln Y_l$ ($l=1,\cdots ,5,7,\cdots ,10$ , with $Y_6\equiv 1$ excluded), which prevents arithmetic overflow at low temperatures. All log-sum-exp operations are stabilised as $\ln {\displaystyle {\sum }_i{\text{e}}^{a_i}}={\max }_i{a}_i+\ln {\displaystyle {\sum }_i{\text{e}}^{a_i-{\max }_i{a}_i}}$ .

Ground-state initialisation. For each parameter set $\left(D_A,D_B,q,J\right)$ , initial log-ratios are provided by a dedicated ground-state routine that runs the ERR at $T_{\text{gs}}=0.05\left|J\right|/k_B$ from all $4\times 10=40$ pure-spin initialisations (every combination of the four $s_A$ states and ten ${\sigma }_B$ states), using a strict convergence tolerance ${\epsilon }_{\text{gs}}={10}^{-14}$ and a maximum of $N_{\text{gs}}=20000$ iterations. The internal energy is computed at each resulting fixed point, and the log-ratios of the minimum-energy fixed point are taken as the starting values for the subsequent temperature (or field) sweep. Among fixed points of comparable energy, preference is given to those with $M_A>0$ to ensure consistency of sign conventions across the phase diagram.

Main temperature and field sweeps. During the sweep, the converged solution at step $n$ is used directly as the initial guess for step $n+1$ (continuation method), so that the iteration tracks a given solution branch continuously rather than always relaxing to the global free-energy minimum. This is essential for following metastable branches (first-order transitions, hysteresis). The iteration is stopped when ${\text{max}}_k\left|\text{Δ}\ell X_k\right|<\epsilon ={10}^{-10}$ , with a hard limit of $N_{\text{max}}=5000$ iterations per step.

The sublattice magnetization $M_A$ is the thermal average of $s_0$ :

$M_A=\frac{3/2A_1{X}_1^q+1/2A_2{X}_2^q-1/2A_3-3/2A_4{X}_3^q}{A_1{X}_1^q+A_2{X}_2^q+A_3+A_4{X}_3^q},$

with $A_1={\text{e}}^{\beta \left(D_A\cdot 9/4+h\cdot 3/2\right)}$ , $A_2={\text{e}}^{\beta \left(D_A\cdot 1/4+h\cdot 1/2\right)}$ , $A_3={\text{e}}^{\beta \left(D_A\cdot 1/4-h\cdot 1/2\right)}$ , $A_4={\text{e}}^{\beta \left(D_A\cdot 9/4-h\cdot 3/2\right)}$ .

Similarly, $M_B$ is

$\begin{array}{c}{M}_B=[9/2{\text{e}}^{\beta \left(D_B\cdot 81/4+h\cdot 9/2\right)}{Y}_1^q+7/2{\text{e}}^{\beta \left(D_B\cdot 49/4+h\cdot 7/2\right)}{Y}_2^q\\ +5/2{\text{e}}^{\beta \left(D_B\cdot 25/4+h\cdot 5/2\right)}{Y}_3^q+3/2{\text{e}}^{\beta \left(D_B\cdot 9/4+h\cdot 3/2\right)}{Y}_4^q\\ +1/2{\text{e}}^{\beta \left(D_B\cdot 1/4+h\cdot 1/2\right)}{Y}_5^q-1/2{\text{e}}^{\beta \left(D_B\cdot 1/4-h\cdot 1/2\right)}\\ -3/2{\text{e}}^{\beta \left(D_B\cdot 9/4-h\cdot 3/2\right)}{Y}_7^q-5/2{\text{e}}^{\beta \left(D_B\cdot 25/4-h\cdot 5/2\right)}{Y}_8^q\\ -7/2{\text{e}}^{\beta \left(D_B\cdot 49/4-h\cdot 7/2\right)}{Y}_9^q-9/2{\text{e}}^{\beta \left(D_B\cdot 81/4-h\cdot 9/2\right)}{Y}_{10}^q]/{\mathcal{Z}}_B,\end{array}$

where ${\mathcal{Z}}_B$ is the corresponding normalization sum (denominator with all positive signs and the same exponential pre-factors).

2.2. Conditions for Second-Order Phase Transitions

Second-order transitions occur when both sublattice magnetizations vanish continuously. At the critical temperature $T_c$ , the conditions are

${M_A|}_{T_c}=3/2{\text{e}}^{{\beta }_c{D}_A\cdot 9/4}\left(X_1^{\ast q}-X_3^{\ast q}\right)+1/2{\text{e}}^{{\beta }_c{D}_A\cdot 1/4}\left(X_2^{\ast q}-1\right)=0,$

$\begin{array}{c}{M_B|}_{T_c}=9/2{\text{e}}^{{\beta }_c{D}_B\cdot 81/4}\left(Y_1^{\ast q}-Y_{10}^{\ast q}\right)+7/2{\text{e}}^{{\beta }_c{D}_B\cdot 49/4}\left(Y_2^{\ast q}-Y_9^{\ast q}\right)\\ +5/2{\text{e}}^{{\beta }_c{D}_B\cdot 25/4}\left(Y_3^{\ast q}-Y_8^{\ast q}\right)+3/2{\text{e}}^{{\beta }_c{D}_B\cdot 9/4}\left(Y_4^{\ast q}-Y_7^{\ast q}\right)\\ +9/2{\text{e}}^{{\beta }_c{D}_B\cdot 1/4}\left(Y_5^{\ast q}-1\right)\\ =0,\end{array}$

where ${\beta }_c=1/\left(k_B{T}_c\right)$ . These conditions are satisfied in the paramagnetic phase through the spin-reversal symmetry relations

$X_1^{\text{*}}=X_3^{\text{*}},X_2^{\text{*}}=1,$

$Y_1^{\text{*}}=Y_{10}^{\text{*}},Y_2^{\text{*}}=Y_9^{\text{*}},Y_3^{\text{*}}=Y_8^{\text{*}},Y_4^{\text{*}}=Y_7^{\text{*}},Y_5^{\text{*}}=1.$

Substituting these conditions into the fixed-point equations yields the self-consistent relations at $T_c$ :

$\begin{array}{c}{X}_1^{*}=X_3^{*}=\frac{{{\displaystyle \sum }}_{\ell =1}^5{\text{e}}^{{\beta }_c{D}_B{\alpha }_{\ell }}\cosh \left({\beta }_{c}J{\gamma }_{\ell }\right)Y_{\ell }^{*p}}{{{\displaystyle \sum }}_{\ell =1}^5{\text{e}}^{{\beta }_c{D}_B{\alpha }_{\ell }}\cosh \left({\beta }_{c}J{\gamma }_{\ell }\right)Y_{11-\ell }^{*p}},\end{array}$

where $\left({\alpha }_1,\cdots ,{\alpha }_5\right)=\left(81/4,49/4,25/4,9/4,1/4\right)$ and $\left({\gamma }_1,\cdots ,{\gamma }_5\right)=\left(27/4,24/4,15/4,9/4,3/4\right)$ , with $Y_6^{\text{*}}=1$ . The symmetric Y-ratios satisfy

$Y_k^{\text{*}}=Y_{11-k}^{\text{*}}=\frac{{\text{e}}^{{\beta }_c{D}_A\cdot 9/4}{C}_1{X}_1^{\ast p}+{\text{e}}^{{\beta }_c{D}_A\cdot 1/4}{C}_2{X}_2^{\ast p}}{{\text{e}}^{{\beta }_c{D}_A\cdot 9/4}{C}_1{X}_3^{\ast p}+{\text{e}}^{{\beta }_c{D}_A\cdot 1/4}{C}_2},$

with $C_1=\text{cosh}\left(3/2{\delta }_k{\beta }_{c}J\right)$ , $C_2=\text{cosh}\left(1/2{\delta }_k{\beta }_{c}J\right)$ , for $k=1,2,3,4$ , where $\left({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4\right)=\left(9/2,7/2,5/2,3/2\right)$ are the spin values of sublattice B. These self-consistent equations are solved numerically for each parameter set $\left(D_A,D_B,q\right)$ to construct the phase diagrams of Section 3.

The total magnetization per site is defined as

$M_T=1/2\left(M_A+M_B\right).$

The compensation temperature $T_{\text{comp}}$ , when it exists, satisfies

$M_T\left(T_{\text{comp}}\right)=0\text{and}\left|M_A\left(T_{\text{comp}}\right)\right|=\left|M_B\left(T_{\text{comp}}\right)\right|\ne 0.$

3. Results and Discussion

3.1. Ground-State Phase Diagram

To map out the stability regions, we calculate the phase diagram in the $\left(D_A/q\left|J\right|,D_B/q\left|J\right|\right)$ plane at $h=0$ . The ground states minimize the energy per site

$H_0=s_i{\sigma }_j-\frac{1}{q\left|J\right|}\left(D_A{s}_i^2+D_B{\sigma }_j^2\right).$

Restricting to ferrimagnetic configurations, we compare all ten relevant spin pairs: (3/2, −9/2), (3/2, −7/2), (3/2, −5/2), (3/2, −3/2), (3/2, −1/2), (1/2, −9/2), (1/2, −7/2), (1/2, −5/2), (1/2, −3/2), (1/2, −1/2). A complete enumeration of all $4\times 10=40$ spin-pair combinations confirms that, for $J<0$ , the minimum-energy configurations are always of the ferrimagnetic type (opposite-sign spins on the two sublattices); ferromagnetic pairs (same-sign spins) and configurations with one zero-spin component have strictly higher ground-state energies throughout the entire $\left(D_A/q\left|J\right|,D_B/q\left|J\right|\right)$ plane explored here, so the restriction to the ten antiparallel pairs listed above is exact.

The resulting ground-state diagram is shown in Figure 1. Seven multicritical points (A-G) emerge where three or more phases become degenerate; their coordinates are listed in Table 2.

The diagram reveals extended stability domains for high-spin states on sublattice B, particularly (3/2; −9/2) and (3/2; −7/2), when $D_B/q\left|J\right|$ takes strongly negative values. Compared with the (3/2, 7/2) case [27] (five multicritical points), the (3/2, 9/2) model has seven such points, reflecting the additional competing ground-state configurations introduced by the spin-9/2 sublattice. This result is in excellent agreement with the ground-state diagram reported in [28] for the same spin pair on a simple cubic lattice ($z=6$ ).

3.2. Thermal and Compensation Behaviors

The magnetization curves were obtained by numerically solving the ERR on the Bethe lattice with $q=4$ and $h/\left|J\right|=0$ .

Figure 1. Ground-state phase diagram of the model in the $\left(D_A/q\left|J\right|,D_B/q\left|J\right|\right)$ plane ($q$ is the coordination number).

Table 2. Coordinate of the multicritical points in the ground-state phase diagram.

Identification of transition types and compensation points. Second-order (continuous) transitions are identified as the temperature at which the order parameter $\left|M_A\right|+\left|M_B\right|$ vanishes continuously. Numerically, $T_c$ is located at the last temperature step where $\left|M_A\right|+\left|M_B\right|>0.02$ , with linear interpolation to the threshold; the threshold value 0.02 is chosen to be well above numerical noise while small enough not to bias the estimate. First-order (discontinuous) transitions are identified by a magnetization jump exceeding $\text{Δ}M>0.08$ between two successive temperature steps along the continuous solution branch; they are further confirmed by the simultaneous existence of two distinct stable fixed points (the ground-state branch and a competing metastable branch) initialized independently for the same parameters. The precise location of first-order boundaries carries the numerical uncertainty of the temperature step used ($\delta T\approx 0.01\left|J\right|/k_B$ in the relevant panels); no additional free-energy comparison between branches is performed. Compensation temperatures are detected as sign changes of $M_T=1/2\left(M_A+M_B\right)$ occurring while $\left|M_A\right|>0.02$ and $\left|M_B\right|>0.02$ simultaneously, confirmed by linear interpolation between the bracketing temperature steps.

Results are shown in Figures 2-5.

Figure 2. Thermal variations of $M_A$ and $M_B$ for $D_A/\left|J\right|=4.0$ and selected values of $D_B/\left|J\right|$ .

Figure 2 shows results for $D_A/\left|J\right|=4.0$ . The sublattice A magnetization starts from its saturation value of 3/2 for all $D_B/\left|J\right|$ , since a large positive crystal field strongly favors $s=3/2$ . Meanwhile, $M_B$ exhibits at $T=0$ nine distinct saturation values (−1/2, −1, −3/2, −2, −5/2, −3, −7/2, −4, −9/2) for $D_B/\left|J\right|=-0.50$ through −4.00, respectively—a direct signature of the ten spin states of sublattice B. Both magnetizations vanish continuously at a common $T_c$ , indicating second-order transitions throughout.

Figure 3 displays results for $D_A/\left|J\right|=-8.0$ , organized in three panels. Panel (a): for $D_B/\left|J\right|<-0.36$ , $M_A$ starts from 1/2 (large negative $D_A$ favors $s=1/2$ ) and second-order transitions occur throughout. Panel (b): for $-0.36\le D_B/\left|J\right|<-0.33$ , the system exhibits first-order transitions, with a discontinuous jump between two competing ordered phases. Panel (c): for $D_B/\left|J\right|\ge -0.33$ , a single second-order transition is recovered.

Figure 4 shows results for $D_B/\left|J\right|=0.4$ and three values of $D_A/\left|J\right|$ . The starting values $M_A=1/2,1,3/2$ for $D_A/\left|J\right|=-12,-9,4$ are fully consistent with the ground-state diagram. The occurrence of both first- and second-order

Figure 3. Thermal behavior of $M_A$ and $M_B$ for $D_A/\left|J\right|=-8.0$ and selected values of $D_B/\left|J\right|$ (three panels).

Figure 4. Thermal variations of $M_A$ and $M_B$ for $D_B/\left|J\right|=0.4$ and different values of $D_A/\left|J\right|$ .

transitions is qualitatively consistent with the mean-field results of [28].

To clarify the compensation behavior, we examine $M_T=1/2\left(M_A+M_B\right)$ for $D_B/\left|J\right|=-2.0$ (Figure 5). All curves except $D_A/\left|J\right|=-0.8$ pass through a compensation point before vanishing at $T_c$ . The compensation behavior is only weakly affected by $D_A/\left|J\right|$ , while $T_c$ increases with $D_A/\left|J\right|$ . This is qualitatively consistent with [29].

Figure 5. Thermal behaviors of the total magnetization $M_T$ for $D_B/\left|J\right|=-2.0$ , $q=4$ , and selected $D_A/\left|J\right|$ . Compensation temperatures are visible for most parameter values.

3.3. Phase Diagrams

Figure 6 displays second-order transition lines in the $\left(D/\left|J\right|,k{T}_c/\left|J\right|\right)$ plane for the isotropic case $D_A=D_B=D$ and $q=3,4,5,6$ . The critical temperature $T_c$ increases as $\left|D\right|$ decreases (weaker crystal field favors ordering) and increases with $q$ (more nearest neighbors enhance ordering).

Figure 7 shows the $\left(D_A/\left|J\right|,k{T}_c/\left|J\right|\right)$ phase diagram. When $D_B/\left|J\right|<-0.25$ (panel a), critical curves behave analogously to Figure 6. For $D_B/\left|J\right|\ge -0.25$ (panel b), the critical lines originate at high temperatures and grow rapidly as $D_A/\left|J\right|\to \infty$ . The final phase diagram, in the $\left(D_B/\left|J\right|,k{T}_c/\left|J\right|\right)$ plane, is presented in Figure 8.

3.4. Hysteresis Properties

Hysteresis computation protocol. Hysteresis loops are computed by sweeping the external field from $h/\left|J\right|=-15$ to +15 (forward branch) and back to −15 (return branch), using $N_h=600$ equally spaced field steps in each direction

Figure 6. Finite-temperature phase diagrams in the $\left(D/\left|J\right|,k{T}_c/\left|J\right|\right)$ plane for $q=3,4,5,6$ : second-order transition lines.

Figure 7. Phase diagrams in the $\left(D_A/\left|J\right|,k{T}_c/\left|J\right|\right)$ plane for $q=4$ , with selected $D_B/\left|J\right|$ values.

Figure 8. Phase diagrams in the $\left(D_B/\left|J\right|,k{T}_c/\left|J\right|\right)$ plane for $q=4$ , with selected $D_A/\left|J\right|$ values, second-order transitions.

($\delta h/\left|J\right|\approx 0.05$ , i.e. step size 30/599). The forward sweep begins from the fixed point converged at $h/\left|J\right|=-15$ ; at each subsequent field value, the ERR iteration is initialized from the converged solution of the immediately preceding field step (sequential continuation). This procedure is essential for capturing hysteresis: in bistable regions, the forward and return branches converge to different fixed points, reproducing the physical irreversibility of the magnetization process. The same iteration parameters as the thermal sweeps are used ($N_{\text{max}}=5000$ ,$\epsilon ={10}^{-10}$ ) at every field step.

We examine the effect of temperature, uniform crystal field, and ferrimagnetic coupling $J$ on the hysteresis behavior.

For $D_A/\left|J\right|=0$ , $D_B/\left|J\right|=-1$ (Figure 9), a single hysteresis loop is present; it shrinks and disappears for $kT/\left|J\right|>6$ , reflecting the progressive loss of magnetic ordering.

Figure 10 shows the hysteresis loops for DA/J = 1 and DB/J = 0 at selected temperatures. The loop area decreases progressively with increasing temperature and vanishes near the critical temperature, indicating the transition from the ferrimagnetic phase to the paramagnetic phase.

Figure 11 shows the effect of the uniform crystal field at $kT/\left|J\right|=0.5$ . No hysteresis occurs for $D/\left|J\right|<-1.9$ ; a single loop appears for $-1.9\le D/\left|J\right|<-1.25$ and for $D/\left|J\right|\ge -1.25$ . The loop width remains nearly constant for $D/\left|J\right|\ge 0$ , so that neither the coercive field nor the remanent magnetization is significantly affected by further increases in the positive crystal field.

Figure 9. Hysteresis loops for $D_A/\left|J\right|=0$ and $D_B/\left|J\right|=-1$ at selected temperatures.

Figure 10. Hysteresis loops for $D_A/\left|J\right|=1$ and $D_B/\left|J\right|=0$ at selected temperatures.

Figure 11. Influence of the uniform crystal field on the hysteresis behavior at $kT/\left|J\right|=0.5$ .

Figure 12 shows the effect of the ferrimagnetic coupling $J$ at $kT/\left|J\right|=0.5$ , $D_A/\left|J\right|=0$ , $D_B/\left|J\right|=-1$ . For $J<-0.6$ , only a single loop is observed. In the range $-0.6\le J<-0.5$ , the system exhibits three distinct hysteresis loops—a qualitatively richer behavior than in the (3/2, 7/2) case [27]—before returning to a single loop for $J\ge -0.45$ , which eventually vanishes at $J=-0.25$ .

The broadening of the triple-loop regime in the (3/2, 9/2) system compared with (3/2, 7/2) has a clear physical origin rooted in the exchange-crystal-field competition on sublattice B. With spin-9/2, there are nine field-induced level crossings between adjacent $\sigma$ states (versus seven for spin-7/2), each corresponding to a discontinuous plateau transition in $M_T\left(h\right)$ . At low temperature and intermediate $\left|J\right|$ , the antiferromagnetic exchange energy and the crystal-field energy $D_B$ compete to stabilise adjacent spin states, so that multiple metastable plateaux remain accessible over a wider range of the ratio $J/D_B$ . This competition generates three simultaneously accessible plateaux—hence the triple-loop structure—over the window $-0.6\le J/\left|J\right|<-0.5$ , a range that is wider than the analogous window in the spin-7/2 case precisely because the higher-spin sublattice supports a larger number of competing spin-state pairs. The higher spin values $\sigma =\pm 9/2$ also carry a larger magnetic moment than $\sigma =\pm 7/2$ , which strengthens the effective exchange field experienced by sublattice A, raising $T_c$ systematically for equivalent reduced parameters. In contrast, for large $\left|J\right|$ , the exchange energy dominates and a single strongly antiparallel ground state is stabilised, suppressing the intermediate plateaux and collapsing the multi-loop structure back to a single loop.

Figure 13 quantifies the coercive field $h_c/\left|J\right|$ and remanent magnetization

Figure 12. Hysteresis loops for $D_A/\left|J\right|=0$ , $D_B/\left|J\right|=-1$ , $kT/\left|J\right|=0.5$ , and various values of J.

Figure 13. Coercive field $h_c/\left|J\right|$ and remanent magnetization $M_R$ versus $J/\left|J\right|$ at $D_A/\left|J\right|=0$ , $D_B/\left|J\right|=-1$ , $kT/\left|J\right|=0.5$ .

$M_R$ as functions of $J/\left|J\right|$ . The coercive field decreases almost linearly to zero around $J\approx -0.25$ . The remanent magnetization remains nearly constant (≈0.77) for strongly negative $J$ , then drops sharply near $J=-0.5$ , marking the transition from the multi-loop to the single-loop hysteresis regime. Strongly negative $J$ values correspond to hard magnetic behavior (high coercivity and remanence), suitable for permanent magnets and data storage, while values near zero correspond to soft magnetic behavior (low coercivity), desirable for transformers and electromagnetic cores.

4. Conclusions

We have investigated the critical, compensation, and hysteresis behaviors of the mixed spin-(3/2, 9/2) ferrimagnetic Blume-Capel Ising model on the Bethe lattice using exact recursion relations. The main findings are:

• The ground-state diagram exhibits ten competing ferrimagnetic configurations and seven multicritical points (versus five for the (3/2, 7/2) case [27]), a direct consequence of the additional competing spin-state pairs introduced by the spin-9/2 sublattice.

• The system displays both first- and second-order phase transitions, with compensation temperatures for appropriate crystal-field values. The critical temperatures are systematically higher than in the (3/2, 7/2) case for equivalent reduced parameters, owing to the stronger effective exchange field associated with the higher spin moment on sublattice B.

• Under an external magnetic field, the model shows single, double, and triple hysteresis loops. The triple-loop regime occurs over $-0.6\le J/\left|J\right|<-0.5$ , driven by nine exchange-crystal-field level crossings on sublattice B (versus seven for spin-7/2). The coercive field decreases nearly linearly to zero near $J\approx -0.25$ .

It should be emphasized that the ERR method is exact on the Bethe lattice: the absence of loops in this tree-like structure renders the recursive decomposition of the partition function rigorous, and all results presented here—phase boundaries, compensation temperatures, and hysteresis loops—are exact within this geometry. The comparisons with mean-field results [28] [29] (which neglect local correlations) and with results on the simple cubic lattice (which has loops and a different coordination number) are therefore qualitative in nature; they confirm overall consistency of the physical picture but should not be taken as quantitative benchmarks. Our findings extend those of [28] [29] by providing a rigorous treatment of local spin correlations on the Bethe lattice and demonstrate that the spin-9/2 extension produces qualitatively new physics not recoverable from the lower-spin case by extrapolation.

The tunable multi-hysteresis loops and compensation temperatures make this model a promising candidate for materials design in magnetic recording, spintronic devices, and sensors. Future work will address next-nearest-neighbor interactions, random crystal fields, and Monte Carlo simulations on finite lattices to complement the present exact results.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References