Thermodynamics of General Heisenberg Spin Tetramers Composed of Coupled Quantum Dimers

Advances in quantum computing technology have been made in recent years due to the evolution of spin clusters. Recent studies have tended towards spin cluster subgeometries to understand more complex structures better. These molecular magnets provide a multitude of phenomena via exchange interactions that allow for advancements in spintronics and other magnetic system applications due to the possibility of increasing speed, data storage, memory, and stability of quantum computing systems. Using the Heisenberg spin–spin exchange Hamiltonian and exact diagonalization, we examine the evolution of quantum energy levels and thermodynamic properties for various spin configurations and exchange interactions. The XXYY quantum spin tetramer considered in this study consists of two coupled dimers with exchange interactions α1J and α1′J and a dimer–dimer exchange interaction α2J. By varying spin values and interaction strengths, we determine the exact energy eigenstates that are used to determine closed-form analytic solutions for the heat capacity and magnetic susceptibility of the system and further analyze the evolution of the properties of the system based on the parameter values chosen. Furthermore, this study shows that the Schottky anomaly shifts towards zero as the ground-state of the system approaches a quantum phase transition between spin states. Additionally, we investigate the development of phase transitions produced by the convergence of the Schottky anomaly with both variable exchange interactions and external magnetic field.


Introduction
In the journey towards the advancement of quantum computing, the tunability of molecular magnets has led to a better understanding of quantum tunneling applications, quantum dots, and anisotropic effects in materials in both theoretical and experimental research [1][2][3][4][5][6][7][8]. Past studies involving the utilization of spin and other variable parameters in such clusters provide evidence of improvements in quantum computation with regards to faster computers with expansive memory and data storage due to refinement of spin switching and spintronics through the use of molecular magnets as quantum gates [9][10][11][12][13][14][15][16]. A further look into how particular systems evolve in time when subjected to different scenarios may allow for optimization of various technological advancements.
A characteristic trait of molecular magnets is that they have variable spin total ground states depending on the types of interactions within the systems. As such, ferromagnetically dominated systems tend to have nonzero spin total ground states, while antiferromagnetically dominated systems typically have spin zero ground states [17]. Some of the crucial variables that allow for the exploitation of this ground state are the exchange interactions between the constituent atoms of the system. These exchange interactions entangle the atoms with various strengths and bring about correlations outside of the classical realm [12,18]. Using the Heisenberg spin-spin Hamiltonian, one can determine the properties of spin-interacting atoms based on the total spin and interaction energies. One particular point of interest that arises from these interactions are the locations of quantum phase transitions (QPTs) as system shifts between spin states and how these points manifest with variable interactions and even external magnetic fields.
At a given QPT, the system's ground state changes the total spin and changes the ground-state energy. Altering the system's energy levels impacts other properties down the line dependent on the system's partition function [8,[17][18][19]. At the low-temperature limit, thermal excitations are minimal, and these effects can be deemed to be of quantum tunneling nature [8]. The addition of an externally applied magnetic field also results in the Zeeman splitting effect and adds a term to the Hamiltonian to separate any degenerate energy levels and to produce a singular ground-state level. The ability to manipulate, create, or destroy exchange interactions allows for a multitude of desired systems for various applications. Studying how the system's bulk properties change with variable interaction energies allows for a functional design.
A major distinguishing feature of molecular magnets is that the combined smaller subgeometries of a magnetic system have individual excitations that make up the larger system as a whole [20,21]. This finding greatly simplifies the study and understanding of larger spin systems by allowing one to break down individual components of more complex systems to observe excitations as parameters of the system are varied and then characterized [10]. With this at hand, closed-form expressions for thermodynamic and magnetic properties can be determined [22]. These smaller clusters may be single molecular magnets surrounded by non-magnetic ligands and are isolated from long-range magnetic interactions or part of a larger lattice [7,17,23,24]. With this in mind, there have been a number of different studies that probe specific cluster configurations for bulk properties in thermodynamics both experimentally and theoretically to gain a better understanding of the larger picture [10,[17][18][19][20][21][22][25][26][27][28][29]. Continuing to theoretically and experimentally research spin clusters of various configurations, constituents, and parameters has the potential to lead to major improvements in device applications.
A promising realm to explore is the ability of antiferromagnetic (AFM) spin clusters to be used in quantum computing. It has been shown that spin clusters with AFM exchange interactions have little to no total magnetic moment and can be used to define qubits, which allow for initialization, gate operations, and readouts using previously known techniques for single spins [1,30]. The property of a small total magnetic moment in these AFM clusters allows for added stability in the presence of an external magnetic field as the ground state does not exhibit spin splitting. Furthermore, the systems can also be manipulated to speeds in the THz range for high-density storage and ultra-fast switching due to the typically small energy shifts between states [30]. These features offer an advantage over ferromagnetic (FM) spin clusters since FM systems tend to have larger total magnetic moments and are likely to be unstable in environments with external magnetic fields [1,31].
The focus of this paper is to provide a complete understanding of the quantum spin tetramer and its bulk thermodynamic properties (ground-state energy, heat capacity, and magnetic susceptibility). Through an examination of the evolution of these properties in various configurations (square, uncoupled, and coupled) of the system, we are able to assess the transitions in both the frustrated/non-frustrated cases. By exploring how the bulk properties evolve as the system is subjected to different exchange interactions, spin values, and an externally applied magnetic field of varying strength, we can provide a guide for experimental characterization and assessment of different spin cluster and molecular magnetic systems.
In this study, we examine variable quantum spin tetramers of general S in both the homogenous (XXXX) and heterogenous (XXYY) spin configurations. We provide examples of calculations that detail the evolution of the thermodynamic properties as spin values are varied from 1/2 to 5/2. The tetramer configurations consist of three different exchange interactions that can be changed independently (antiferromagnetic and ferromagnetic), symmetric and asymmetric systems (XXXX and XXYY spin configuration, respectively) are examined for bulk properties excluding magnetic susceptibility, one mixed-valence case (XXYY spin configuration) is covered to show how magnetic susceptibility evolves with varying exchange interaction strengths in comparison to ground-state/quantum phase transitions and heat capacity, and the systems are explored and compared with and without a varying magnetic field. A low-temperature restriction (from 0 to 1 K) is used to approach the study in the quantum regime.

The Quantum Spin Tetramer Representation
The tetramer system studied in this paper comprised four magnetically entangled atoms with various exchange interactions among them. There were two possible dimer interactions (taken between diagonal components) and four interactions around the perimeter of the system that were assumed to have the same strength. Figure 1 shows the general representation for the tetramer geometry with dimer exchanges and other possible configurations. Figure 1a shows the general representation containing one dimer exchange between S 1 and S 3 and another dimer exchange between S 2 and S 4 (referred to as S d1 and S d2 , respectively). S d1 and S d2 has exchange interactions of α 1 J and α 1 J, respectively. Dimer-dimer interactions are also presented around the edges of the tetramer (S 1 with S 2 , S 2 with S 3 , S 3 with S 4 , and S 4 with S 1 ) and have an exchange interaction denoted as α 2 J. Figure 1b,c show scenarios that are brought about by the extremes of completely turning off/on different exchange interactions within the system. Figure 1b is the result of completely turning off the α 2 J exchange to produce two uncoupled dimers. Alternatively, Figure 1c shows the result of switching off the α 1 J and α 1 J exchanges to produce the square configuration. It is important to note that while we show this structure as a square, the spatial configuration is not relevant as only the magnetic spin interactions are considered in the energy levels.
The values of S i explored throughout this study are mixtures of 1/2, 1, 3/2, 2, and 5/2, where the atoms that S d1 consisted of had an identical spin state, as does the atoms that make up S d2 , but both dimers were not necessarily in the same spin state as one another. The total spin value is the sum of all spins and is denoted as S T . Figure 1. The general representation of the tetramer. Each constituent atom has a spin value that couples to two other atoms via an exchange interaction α i J, where α i is the scaling parameter of interaction strength. (a) The coupled tetramer that contains dimer pairs S d1 and S d2 . S d1 is formed between S 1 and S 3 , whereas S d2 is formed between S 2 and S 4 . It also contains dimer-dimer interactions along the outside edges of the tetramer; (b) The uncoupled dimer configuration is formed when α 2 = 0; (c) The square configuration is formed when α 1 = α 1 = 0. It is important to note that, while we show this structure as a square, the spatial configuration is not relevant as only the magnetic spin interactions are considered in the energy levels.

General Quantum Spin Tetramer
Using a Heisenberg spin-spin exchange Hamiltonian for a magnetic system, the energy levels of multiple quantum spin tetramer configurations were derived and used to explore ground-state levels' evolution and to form the partition function for a given configuration. The partition function allows for the determination of heat capacity, magnetic susceptibility, and quantum phase transitions for various parameter conditions (temperature, magnetic field, spin, and interaction strengths). The isotropic spin-spin exchange Hamiltonian for a magnetic system is defined as: where J ij is the exchange interaction between spins S i and S j , m z is the z-component of the spin state, µ B is the Bohr magneton, g is the Landé factor, and B is an externally applied magnetic field. Applying this Hamiltonian for the general tetramer shown in Figure 1 with applied field produces where α 1 , α 1 , and α 2 are the exchange parameters and S i is the value of spin. For this study, J = 1 and the exchange parameters determine the antiferromagnetic (α i > 0) or ferromagnetic (α i < 0) interactions. With the S z basis, the Hamiltonian allows for determination of the general energy, E, as a function of S T , S d1 , S d2 , S i , m z , and B: From here, the energy eigenstates and eigenvalues are found by applying spin values to all the constituents. Setting α 2 = 0 or α 1 = α 1 = 0 produces the energy for the the uncoupled dimer or square tetramer configurations in Figure 1b,c, respectively. The individual dimer spin values and the total spin of the tetramer can be found using the spin decomposition formula for the system: which for a system of four identical spins, as in the XXXX case, will give a total spin decomposition of where the decomposition is summed over the integer of N. However, in the XXYY case, this spin decomposition is less trivial and can be shown to be With the restriction that the spins making up each respective dimer are the same, this representation is simplified to where each bracket term pertains to a dimer and contains a direct sum from 2S i to 0 in steps of 1. Equation (7) allows for total spin, S T , to be calculated from particular dimer states, and the correlated values can be plugged into Equation (3) to find the energy of the system. In the case of zero-field, each magnetic state has (2S T + 1) degenerate states and (2S 1 + 1) 2 (2S 2 + 1) 2 total states. Once the possible energy levels are determined, the systems partition function is deduced by where β = 1/(k B T). The second term of Equation (8) is summed through all N energy eigenvalues when an external magnetic field is applied. The third term is summed through all energy eigenvalues with zero-field and contains a coefficient of the number of degenerate energy states. Bulk properties such as heat capacity and magnetic susceptibility can then be determined by manipulating the systems partition function discussed in further sections.

Ground-State Energy Levels and Heat Capacities
Once the system partition function is found, it can be used to find the system heat capacity using the following: Since a change in temperature of a system leads to thermal excitations in lattice vibrations and electronic distributions, heat capacity has an overwhelming dependence on temperature [39]. Examining the heat capacity at low temperatures gives rise to the phenomenon known as the Schottky anomaly, which represents a dramatic shift in entropy as a disordered paramagnetic system becomes more ordered. In recent years, the Schottky anomaly in quantum Heisenberg spin clusters has been shown to indicate quantum phase transitions as the ground-state energy (GSE) level shifts with some tuning parameter [19]. The reason this effect is not to be expected is that the heat capacity usually increases with temperature. However, there have been sharp peaks observed in heat capacity data at low temperatures before dropping off again [40,41]. This can be explained by observing the entropy of the system during this transition as governed by where N is the total number of energy eigenstates and N 0 is the total number of ground states [17,19,20,42]. Using the spin-spin exchange Hamiltonian along with the spin decomposition, a more considerable value of S T contains more possible eigenstates (with a large portion of them being degenerated in zero-field). Nevertheless, the number of ground states does not increase as much. Given that this occurs during a small temperature change, entropy is seen to peak in correlation with the spike in heat capacity. Thus, knowing where the ground state of the system changes due to varying parameters allows one to correlate heat capacity peaks (and magnetic susceptibility steps) to particular S T values. By exploiting the Schottky anomaly and by comparing theoretical ground-state transition data with heat capacity maps, the correlation between the two can help determine spin values, exchange interactions, and other parameters in experimental systems. The following sections delve into how ground states and heat capacities theoretically evolve as multiple parameters vary. The first section involves the zero-field case with degenerate states. The second section shows how an applied magnetic field produces Zeeman splitting and separates magnetic substates to give non-degenerated ground states.

Zero-Field Ground-State Energies and Heat Capacities
Figures 2-4 cover a multitude of data for zero-field cases as spin and interaction strength are varied for the tetramer system. The first that we discuss are Figures 2 and 3, which shows how the ground-state energy/heat capacity maps transform as dimer constituent spin values and exchange parameters vary. The ground-state maps were created by determining the lowest energy state of the system over iterated parameters, and the resulting spin values were mapped for that particular energy level. Figure 2a contains the ground-state maps as functions of α 1 and α 2 , where the dimer exchange interactions, α 1 and α 1 , were assumed equal, and J = 1. The colors of the maps depict the system S T values for that particular ground state as the spin values of the individual atoms in dimer one and dimer two (S 1 from top to bottom and S 2 from left to right, respectively) varied from S i = 1/2 to S i = 5/2 in the steps of 1/2. At the same time, α 1 and α 2 varied from −2 to 2. The dark gray color is not indicated in the S T scale bar but is noted to be the |0, 0, 0, 0 > state (|S T , S d1 , S d2 , m z >), which is present in all of the ground-state maps. Due to the restriction that α 1 = α 1 , there is symmetry in the ground-state maps about the diagonal S 1 = S 2 , meaning that the energy maps for the XXYY spin configuration are the same for the YYXX spin configuration. Maintaining S 1 = 1/2 and scanning S 2 from lower to higher spin values introduces more accessible ground states that "fan" out from the symmetric spin-XXXX QPT boundaries. Restricting α 2 > 0 (AFM) and varying α 1 from −2 to 2 results in the ground-state S T value initially at |2(S 1 − S 2 )|, which then continuously decreases by one across transitions to the |0, 0, 0, 0 > state. This transition indicates that the interactions between the constituent spins of S d1 and S d2 are initially FM but have an AFM dimer-dimer interaction due to α 2 . Overall, this configuration produces the non-frustrated tetramer, with the spins of one dimer being up and the other dimer's spins being down.
As α 1 and thus α 1 are shifted to be AFM exchanges, the data show that the spins in the respective dimers do not flip immediately upon the exchanges being greater than zero. This result suggests that α 2 plays a stronger role with regards to spin-flipping in the tetramer system as, for example, it is required for α 1 = α 1 = α 2 in order for the S T transition to occur along with the S 1 = S 2 diagonal. Once this flip takes place, the tetramer is in a frustrated configuration where the dimer-dimer interaction varies from AFM to FM. Furthermore, an increased number of transitions occur as |2(S 1 − S 2 )| increases in value, and α 2 appears to lose relative strength to the individual dimer interactions. Once the transitions occur in the ground-state map along a row, the energy states' existing boundaries do not shift, but instead, new boundaries are removed/added.
Alternatively, restricting α 2 < 0 (FM) and varying α 1 from −2 to 2 results in a slightly different trend. Initially, the ground state has an S T value of 2(S 1 + S 2 ). At this point, all spins are in a single direction for an non-frustrated FM system. However, new QPTs do not appear until |S 1 − S 2 | is equal to or greater than one. Once this requirement is reached, the data show |2(S 1 − S 2 )| transitions where the initial state has S T = 2(S 1 + S 2 ) and the final transition is to the |0, 0, 0, 0 > state; this completely skips transitions to certain spin total states. The transitions between occur in decreasing steps of 1 from 2(S 1 + S 2 ), i.e., for S 1 = 1/2 and In contrast to when α 2 is an AFM interaction, when it is an FM interaction, the transition boundaries appear to shift along rows when scanning through S 2 values. Overall, α 2 appears slightly weaker in the tetramer system when it is an FM interaction regarding transitions.
The ground-state energy map in Figure 2a can be used in unison with Figure 3a,b to correlate ground-state transitions and heat capacities. The scale bar on the right side of Figure 3 denotes heat capacity (with k b factored out). The top right of the diagram in Figure 3a pertains to α 1 = α 1 = 1 and the bottom left of the diagram pertains to α 1 = α 1 = −1. As with Figure 2, the spins of the constituent atoms making up dimers S d1 and S d2 are represented by S 1 and S 2 . Due to the symmetry of the ground-state map under identical restrictions, the heat capacities are also symmetrical about the S 1 = S 2 diagonal, which is why only half of each data set is represented. Figure 2. The ground state energy levels as functions of alpha parameters zero-field and J = 1. S 1 and S 2 denote the spins of the atoms making up S d1 and S d2 , respectively. The scale on the right side of the figure denotes the total spin of the system, S T , and relates a particular value to a color in the ground-state energy diagrams. Energy increases from a low value of S T to a higher value. (a) The ground-state energy map for energy as a function of α 1 = α 1 and α 2 ; (b) The ground-state energy map for energy as a function of α 1 and α 1 , where α 2 = 1. The restriction of α 2 's value limits the attainable S T,max state through transitions. The gray areas in (a,b) denote the |0, 0, 0, 0 > state (|S T , S d1 , S d2 , m z >).
The associated heat capacities for the ground-state energy diagrams in Figure 2a as a function of α 1 and k B T/|J|. The top right of the diagram is for α 1 = α 1 and α 2 = 1, whereas the bottom left of the diagram is for α 1 = α 1 and α 2 = −1. Given the restriction that α 2 = 1 in the top right of (a), these heat capacity plots correlate with Figure 2c. Figure 3a shows heat capacity as a function of α 2 and k B T/|J|. It is apparent that the heat capacity experiences spikes for a variety of spin and α 2 values in the low-temperature limit. Due to the Schottky anomaly arising with ground-state transitions, it is expected that setting α 1 = α 1 = 1 and scanning α 2 across various spin values should correlate to ground-state transitions shown in Figure 2a. Observing the heat capacity for S 1 = 1 and S 2 = 2 shows that energy level crossings occur when α 2 = −2, −1.5, 1, 1.5, 2. Comparing this data with the respective ground-state diagram shows that these specific values for α i parameters correlate to ground-state transitions from S T = 6 → 5 → 0 → 0 → 1 → 2 with respect to the α 2 values above. The first of the two S T = 0 ground states contain dimers with individual AFM exchange interactions (|0, 0, 0, 0 >), and the second S T = 0 ground state is |0, 2, 2, 0 >, as determined by the spin decomposition and data collected. Repeating this process for the remaining top right of the diagram in Figure 3a shows the correlation of QPTs with ground-state energy transitions and heat capacity peaks. Using the scale bar for the values of heat capacity, one can see that, when the dimer exchange interactions are of AFM nature, the heat capacity spikes are relatively large throughout the diagram. The difference in heat capacity at these transitions increases with the total spin of the system, where there is a difference of approximately 2.5 (in arbitrary units) for the S 1 = S 2 = 1/2 case and a maximum difference of approximately 5.5 for the S 1 = S 2 = 5/2 when α 2 ≈ −1.25.
Using a similar process for comparing the ground-state maps in Figure 2a to the bottom left of the diagram in Figure 3a, it is apparent that a single energy level crossing occurs at α 2 = 0 for the case of ferromagnetic dimer exchanges for all XXYY spin values modeled. Furthermore, the spikes in heat capacity are smaller in value compared to the AFM cases with a difference of approximately 1.4 (in arbitrary units). Scanning α 2 values from −2 to 2 leads to a higher value of S T since all spins are in the same direction to a lower value once α 2 > 0. The system is non-frustrated in both scenarios.
The diagram in Figure 3b can also be used in unison with the ground-state maps in Figure 2a to determine where ground-state transitions and heat capacity spikes occur for the system. Figure 3b shows heat capacity (with k b factored out) as a function of α 1 and k B T/|J|. The top right of the diagram holds that α 2 = 1 and the bottom left of the diagram holds that α 2 = −1. It is also restricted that α 1 = α 1 . For the case of S 1 = S 2 = 1/2 and α 2 = 1, scanning α 1 from −2 to 2 in both Figures 2a and 3b leads to a correlated ground-state transitions and heat capacity spike, as before. Holding a constant value of S 1 and increasing S 2 produces "fanning" of ground-state transitions of the system mentioned previously. Based on the color map, the spikes in heat capacity vary as S T shifts, with the greatest difference occurring for the transition to the |0, 0, 0, 0 > state.
The bottom left of the diagram in Figure 3b shifts the value of α 2 to −1. Using this alongside Figure 2a again shows transitions forming/vanishing around the symmetric XXXX spin configuration maps. The differences in heat capacity as S T shifts is more significant for this scenario as the system transitions from an non-frustrated state (all α i < 0) with all spins aligning along the same direction to a frustrated state where the dimer exchanges are AFM. However, the dimer-dimer exchange varies from AFM to FM between nearest-neighbor spins. The sizable heat capacity difference is approximately 5.5 (in arbitrary units) for the S i = 5/2 model, which is due to the transition S T = 5 → 0.
Next, Figure 2b shows a different ground-state energy map with the restriction that α 2 = 1. This diagram shows how the ground-state energy levels transition when α 1 and α 1 do not necessarily equal each other. Due to this restriction, this diagram can only be used with the top right of the diagram in Figure 3b for heat capacity. Similar trends from the previous correlations between ground-state transitions and heat capacity can be observed by starting at the center of any ground-state map for particular spin values and then by scanning up and right to find transitions. A difference to point out with this ground-state map is the asymmetry across the XXXX spin configuration diagonal, which arises from non-identical strengths for the dimer exchange. The exchanges have the effect that the ground-state energy map for S 1 = 3/2 and S 2 = 2 can be rotated ninety degrees and mirrored to produce the map for S 1 = 2 and S 2 = 3/2. Another notable difference for this map is that it consists primarily of low spin total values with a maximum of S T = 4 in the top right and bottom left panel of the diagram. This effect is caused by the restriction of α 2 = 1. By observing Figure 2a for the case of α 2 = 1, one can see that changing the value for α 1 (and thus α 1 ) limits the accessible S T states to a maximum of S T = 4. Altering the value of α 2 to be less than zero unlocks the possibility of larger spin total states. Moving onto Figure 4: the diagram in Figure 4a is the ground-state energy map as a function of α 1 and α 2 with the restriction of α 1 = 1. The diagram in Figure 4b is the respective heat capacity as a function of α 1 and k B T/|J| with the restrictions α 1 = 1 and α 2 = 0.5. S 1 and S 2 denote the spin values of the constituent spins making up S d1 and S d2 , as before. It is assumed that J = 1. The S T scale bar is the same as the one in Figure 2. The heat capacity scale bar has a new maximum of 3.4 (in arbitrary units).
The ground-state energy map in Figure 4a again shows competing exchanges in the tetramer system. Beginning with S i = 1/2 and choosing α 2 = 1, setting α 1 = −1 produces an non-frustrated tetramer in the |0, 1, 1, 0 > state. For this to occur, α 2 must have a stronger interaction strength relative to α 1 . Since both dimer spin values are nonzero, α 2 forces the spins entangled with α 1 to have a ferromagnetic exchange to produce a nonzero S d2 value. However, the data shows that increasing α 1 allows the dimer exchanges to dominate once α 1 > 1 for this scenario. At this point, the dimer spins flip to an AFM configuration to produce a frustrated tetramer in the |0, 0, 0, 0 > state. For the cases of S 1 = 1/2 and the restrictions listed, the boundary for this transition occurs when α 1 ≥ 2(α 2 − 0.5) for α 2 ≥ 0.5. Alternatively, the same initial restrictions with the exception of setting α 2 = −1.5 produces the non-frustrated tetramer in the |2, 1, 1, 0 > state. Again, α 2 dominates α 1 and causes the spins in the system to be all in the same direction. Once the condition α 1 ≥ −(α 2 + 1) for α 2 ≤ −1 is met, the spins of each dimer will oppose each other to form the |0, 0, 0, 0 > state.
Navigating through the entirety of Figure 4a shows that increasing the spin values for the dimers results in a continuous transition boundary to the |0, 0, 0, 0 > region. By scanning across the S 1 rows, the data also show that increasing the system's spin results in new ground-state transitions for larger values of |α 2 | without shifting the previous transition boundaries. The diagram also suggests the ability to form a multitude of systems by setting restrictions to either α 2 or α 1 . For example, setting S i = 5/2 and α 1 < 0 allows for all S T states to be attained by varying α 2 from −2 to 2. The ability to tune these α i exchange parameters may help develop the tetramer system's applications via quantum tunneling, quantum dots, and anisotropic effects.
If the ground-state energy maps for this configuration were instead set to factors of α 1 and α 2 with α 1 = 1, the same ground-state maps in Figure 4a would be produced with the exception of the individual plots being in a transposed position yet the same orientation. Figure 4b confirms the Schottky anomaly for this case scenario as well. As the system's total spin is increased and the transition boundary into the |0, 0, 0, 0 > state becomes more rounded, the spike in heat capacity slowly shifts to values greater than α 1 = 0. Another feature of the heat capacity diagrams is how the left side of each one shifts at low temperature when α 2 = 0.5 lies along a transition boundary. The low left portion of the respective heat capacity plot contains a less subtle transition from C/k b ≈ 0 → 1 as in S 1 = 3/2 and S 2 = 2. Comparing this to a ground-state map where the value of α 2 is not along a transition boundary, the same transition in heat capacity occurs within a smaller temperature range, as can be seen in the S 1 = 1, 3/2 diagrams.

Heat Capacities in the Presence of a Magnetic Field
With the application of a magnetic field, B, Zeeman splitting occurs, and the separate magnetic states for degenerate energy levels split into (2S T + 1) magnetic substates and potentially create new ground states for the system [10,43,44]. The energy of the system is reduced for the m z states that are greater than zero, as shown in Equation (3). Therefore, states with higher total spin have the most dramatic decrease in energy for a given magnetic field strength. For a system that does not have all ferromagnetic exchange interactions, the applied B field aligns the spins in the system. With a sufficiently strong enough field, the maximum decrease in energy due to the splitting term occurs when m z = S T when a fully ferromagnetic system is created.
A shift in the ground-state energy levels from an applied field also correlates to shifts in QPTs and heat capacity spikes. Figure 5 shows heat capacity (with k b factored out) as a function of k B T/|J| and µ B g|B|/|J| for multiple cases where µ B is the Bohr magneton, and g is the Landé g-factor for electrons. Looking down the S 1 = S 2 diagonal in Figure 5a, there is a single transition when the magnetic field gµ B |B|/|J| = 0.5. The only change in these plots is the value of heat capacity near absolute zero, as indicated by the scale bar. The data show that the ground state is |0, 0, 0, 0 > before this transition. After the transition, the ground state is the completely ferromagnetic configuration. For S i = 1, the ground state is |4, 2, 2, 4 >. A larger difference in S T during this transition leads to greater heat capacity values residing in larger areas of the heat capacity diagrams. It was also observed that the intensity of heat capacity at low temperatures before the transition is dependent on S 1 for the parameters shown. This transition point along the diagonal for the symmetric tetramer follows where 0 = µ B g|B|/|J|. Moving away from the diagonal in Figure 5a shows that more transitions occur as the difference between S 1 and S 2 increases. It can be observed that the number of transitions for given spin values is 1 + 2|S 1 − S 2 | and the last transition occurs when = 2(S 2 α 2 + S 1 α 1 ). (12) The transitions between the first and final spikes occur in decreasing steps of one from until the total number of transitions for the spin values is reached. Transitions between different total spin states do not shift as total spin increases or decreases, but transitions are either created or removed. For the case of S 1 = S 2 = 1/2, the first transition is for S T = 0 → 2. The initial transition is always between the same S T states as S 1 is increased, and the pattern is repeatedly seen throughout the diagram.
Moving onto the top right of the diagram in Figure 5b, the restriction of α i = 1 is put in place. There are more possible transitions as a magnetic field is applied to the system since the restriction holds that all exchange interactions are AFM in nature. Each individual heat capacity diagram contains transitions that follow = α 2 , 2α 2 , 3α 2 , ..., S T α 2 .
The bottom left diagram in Figure 5b shifts the exchange restrictions to α 1 = α 1 = −1 and α 2 = 1. An easily distinguishable difference compared to the rest of the figure is removing transitions for S 1 =S 2 . Individual heat capacity plots show 4S 2 + 1 transitions, and there is always a transition when = 0. The last transition occurs when = S T α 2 , and the following transitions occur in decreasing steps of α 2 . The tetramer is in the state S T = |2(S 1 − S 2 )| immediately after = 0 and increases incrementally as each transition occurs.
Other data (not shown) were analyzed for the case of α 1 = α 1 = 1 > α 2 . It was found the number of transitions observed when an external magnetic field is applied was 0.5S T when α 2 = 0. Under this condition, the transitions in ground-state energy and heat capacity for the system occur when = nα 1 , where n is an integer value from 1 to 0.5S T . As α 2 increases but is still less than α 1 , the transitions occur when = α 1 , α 1 + α 2 , 2α 1 + α 2 , 2(α 1 + α 2 ), 3α 1 + 2α 2 , 3(α 1 + α 2 ), etc. and continues until S T transitions occur. These data suggest that creating tetramer clusters with weaker dimer interactions relative to dimer-dimer exchange results in a system that requires a less powerful magnetic field to cause spin-flipping to occur. Alternatively, a tetramer system with relatively strong dimer-dimer interactions could prove to be more stable in an environment submersed in unavoidable various magnetic field sources.

Magnetic Susceptibility
The application of a variable external magnetic field at certain temperatures causes the magnetic dipoles to align into FM or AFM configurations. Given that transition metals exhibit metallic bonding, magnetic properties arise due to partially or completely free electrons due to electrical conductivity [45]. Since magnetism is dependent on a current in a system, metals with differing numbers of unpaired electrons exhibit different responses to an applied magnetic field (hence the restriction of S i,max = 5/2 for 3d orbitals). This response is referred to as the magnetic susceptibility of the system. The Schottky anomaly in the presence of ground-state transitions has also been shown to correlate to sudden changes in the magnetization of materials composed of spin clusters in previous studies [17,19,20,22,46]. An increase in the number of magnetic substates with increasing S T can be linked to the increasing step of magnetic susceptibility in these quantum systems. Additionally, the magnetic susceptibility for the tetramer is expected to tend towards zero as the temperature approaches zero since an even-numbered cluster has the ability to have spin-0 ground states [22].
The magnetic susceptibility of a system can be determined using with M z = m z gµ B and where m z = S z T /h is the integral or half-integral magnetic quantum number [17,19]. Equation (14) makes it clear that the dependence of total spin of the system is a relatively large factor in the susceptibility.
Step increases in magnetic susceptibility data are expected to match quantum phase transitions as the total spin of the system changes about those transitions. Since the magnetic field dependence of the system is affects the energy levels of the system by shift energy states by M z B, the field dependence on the magnetic susceptibility can be easily determined.
Magnetic susceptibility data for particular mixed valence (S 1 = 3/2 and S 2 = 1) tetramer systems can be seen in Figures 6 and 7, where heat capacity (top), ground-state transitions (middle), and magnetic susceptibility (bottom) are shown tied with each other for various α i restrictions in Figures 6a,c and 7a,c. Magnetic susceptibility was multiplied by a factor of T to avoid discontinuities at T = 0. The evolution of the tetramer system is shown as particular exchange parameters varied in Figures 6b and 7b. Heat capacity (with k b factored out) and the product of magnetic susceptibility and temperature (with ((gµ B ) 2 /k B ) factored out) are shown as functions of k B T/|J|, µ B g|B|/|J| and a particular α i parameter. The ground-state maps shown between the heat capacity and magnetic susceptibility data label the ground states of the system as before and after transitions. The first four panels (left to right) contain the same ground states as labelled in the diagram, whereas the last panel contains different states for particular cases. Figure 6a is for the case of α 2 = 1 while α 1 and α 1 varied from 0 to 1 in steps of 0.250, which corresponds to the evolution diagram in Figure 6b. It is observed that the transitions in heat capacity and the ground state of the system align with the steps in magnetic susceptibility and that an increasing value of S T leads to increased susceptibility. For the case of α 1 = 0.75, the transition at = 0 shifts to a transition at = 0.25 and a new possible ground state was created; this state was found to occur once α 1 = α 1 = 0.667. The other ground states in this panel match the ground states for α 1 = 0, 0.25, 0.50 and all share the same total dimer spin states, whereas the S T = 0 state arises from a different dimer configuration. Once all α i = 1, then all possible dimer configurations of which gives S T = 0, 1, 2, 3, 4 from the spin decomposition contribute to the ground-state energy level; the energy levels are degenerated. Additionally, the maximum heat capacity (≈2.65 in arbitrary units) shifts downward in the diagrams as the tetramer evolves from the square to coupled configuration. Figure 6. The figure compares the heat capacity , ground-state energy levels, and magnetic susceptibility times temperature as functions of k B T/|J| and µ B g|B|/|J| for the mixed valence (S 1 = 3/2 and S 2 = 1) system. (a) The data for C (top), ground-state energy (GSE) (middle), and χ * T(bottom) for the case of α 1 = α 1 , which varied while α 2 = 1. (b) The evolution of the system evolving from the square configuration to the coupled tetramer configuration as α 1 = α 1 varied from 0 to 1. (c) The data for C (top), GSE (middle), and χ * T (bottom) for the case of α 1 = α 1 , which was varied while α 2 = −1. The ground states between the transitions are labelled throughout the diagrams. The first four panels (from left to right) contain the individually labelled states for those four panels. The last panel changes depending on the case. Figure 7. The figure compares the heat capacity, ground-state energy levels, and magnetic susceptibility times temperature as functions of k B T/|J| and µ B g|B|/|J| for the mixed valence (S 1 = 3/2 and S 2 = 1) system. (a) The data for C (top), ground-state energy (GSE) (middle), and χ * T (bottom) for the case of α 1 = α 1 = 1, while α 2 varied. (b) The evolution of the system evolving from the uncoupled dimer configuration to the coupled tetramer configuration as α 2 varied from 0 to 1. (c) The data for C (top), GSE (middle), and χ * T (bottom) for the case of α 1 = α 1 = −1 and varied α 2 . The ground states between the transitions are labelled throughout the diagrams. The first four panels (from left to right) contain the individually labelled states for those four panels. The last panel changes depending on the case. Figure 7a is for the case of α 1 = α 1 = 1 while α 2 varied from 0 to 1 in steps of 0.250, which corresponds to the evolution diagram in Figure 7b. Once again, the transitions in heat capacity and ground states indicate a transition in magnetic susceptibility. The ground states identified for the case of α 2 = 0 match the identical colors found throughout the cases for α 2 = 0.25, 0.50, 0.75. Once all α i = 1, the same result from Figure 6a was obtained where all possible dimer states from the spin decomposition contribute to the respective ground states. It can also be observed that the maximum heat capacity of the system (≈2.2 in arbitrary units) shifts gradually along = 0 as α 2 increased. The transitions, in this case, are in agreement with the previous statement considering that energy level crossings occur when α 2 < α 1 = α 1 (assuming both exchanges are of AFM nature), which can be found at the end of the previous section. If the individual dimer exchanges are identical, one could use these values along with to determine α 2 . Figure 6c is for the case of α 2 = −1 while α 1 and α 1 varied from 0 to 1 in steps of 0.250, which corresponds to the evolution diagram in Figure 6b. The transitions follow a similar pattern as those described for Figure 5a; there are a maximum number of transitions equal to 1 + 2|S 1 − S 2 |. The initial transition occurs at 0 (Equation (11)), and the next successive transitions occur at (Equation (12)). Furthermore, the transitions were not found to occur until α 1 = α 1 = 0.667. The data again correlates heat capacity, ground-state transitions, and magnetic susceptibility with each other. Figure 7c is for the case of α 1 = α 1 = −1 while α 2 varied from 0 to 1 in steps of 0.250, which corresponds to the evolution diagram in Figure 7b. For the case of α 2 = 0, the data obtained matches what is expected for two uncoupled dimers with FM interactions; the only transition occurs at = 0, where the state goes from |5, 3, 2, −5 >→ |5, 3, 2, 5 >. Other ground states begin to appear immediately after the dimer-dimer interaction increased, and transitions match the data obtained from Figure 5b.

Conclusions
In conclusion, we explored several in-depth aspects of the tetramer system and how its bulk properties change with respect to various parameters. We considered the XXYY quantum spin tetramer system's excitations in the coupled, uncoupled, and square configurations. Using the Heisenberg spin-spin exchange model, information about the systems energy levels, spin states, and spin decompositions were established. From there, the partition function was used to determine the thermodynamic properties such as heat capacity and magnetic susceptibility. These bulk properties were then presented in several cases as functions of temperature, applied field, exchange interactions (AFM and FM), and spin values (from S i =1/2 to 5/2). The data presented exploited the Schottky anomaly and its dependence on ground-state energy levels to correlate phase transitions between ground states, heat capacities, and magnetic susceptibility in the low-temperature range. The data were analyzed to determine transition points based on different exchange interaction values or the applied magnetic field's strength. Using the data provided, one could compare experimental results to find a best-fit match to characterize an unknown system, to gain knowledge about the possible system configurations, or to determine what constituents to use in a system to achieve the results sought after.
Further studies on the tetramer system (alongside other molecular magnet subgeometries) could help develop future applications and quantum computing advancements via spin switching. Some avenues for advancement pertain to increased speed, data storage, memory, and stability of particular systems. Furthermore, understanding how the properties of these smaller subgeometries evolve with respect to various parameters can give insight towards increasingly complex spin clusters and how to use them best to push the current technological limits.  Data Availability Statement: All data plotted was generated from the analysis of the model stated above.

Acknowledgments:
The authors acknowledge support from the Institute for Materials Science at Los Alamos National Laboratory.

Conflicts of Interest:
The authors declare no conflict of interest.