Entropies and IPR as Markers for a Phase Transition in a Two-Level Model for Atom–Diatomic Molecule Coexistence

A quantum phase transition (QPT) in a simple model that describes the coexistence of atoms and diatomic molecules is studied. The model, which is briefly discussed, presents a second-order ground state phase transition in the thermodynamic (or large particle number) limit, changing from a molecular condensate in one phase to an equilibrium of diatomic molecules–atoms in coexistence in the other one. The usual markers for this phase transition are the ground state energy and the expected value of the number of atoms (alternatively, the number of molecules) in the ground state. In this work, other markers for the QPT, such as the inverse participation ratio (IPR), and particularly, the Rényi entropy, are analyzed and proposed as QPT markers. Both magnitudes present abrupt changes at the critical point of the QPT.


Introduction
The study of phase transitions in quantum systems is a topic of present interest, usually referred to as quantum phase transitions (QPT) [1][2][3]. Since the seminal Gilmore and collaborators works [4][5][6], there have been numerous papers characterizing QPTs in two-level quantum systems of different dimensionalities used to model nuclear and molecular systems, such as the interacting boson model (IBM) and the vibron model (see [7][8][9] and references therein).
In connection with this, particular Hamiltonians based on algebraic structures that are exactly solvable have been proposed, such as the Lipkin [10], the Bose-Hubbard [11], the Jaynes-Cummings, the Tavis-Cummings and the Dicke models [12][13][14], just to cite a few of them. These models present specific dynamical symmetries that correspond to different equilibrium configurations of the system in the ground state. The algebraic structure of these models allows for simple solutions in some cases, which provide important references for more complex systems In this work, a solvable two-level model that represents the coexistence of atoms and homo-nuclear diatomic molecules is used to study QPTs [15][16][17][18]. The model is briefly presented in Sect. II, where the matrix elements relevant in the model Hamiltonian are given explicitly on a basis with two labels: the number of molecules and the number of atoms. The eigenvalues and eigenvectors of the Hamiltonian are easily obtained by diagonalizing the corresponding matrix. The thermodynamic or large particle number limit of the model is also presented so as to classify the QPT and analyze the critical transition point. The model has one control parameter that drives the system from a molecular condensate, in one phase, to a new phase in which atoms and molecules coexist. Usual markers for the critical point are the ground state energy and the behavior of an order parameter that is zero in one phase and different from zero in the other one. Usually this order parameter is the expected value in the ground state of the number of atoms (or the number of molecules). In this work, we propose the use of the inverse participation ratio (IPR) and the Rényi entropy as other good markers for the critical point. They are presented in Section 3. Section 4 is for conclusions.

The Model for the Atom-Diatomic Molecule Coexistence
A simple two-level model designed to describe a system of two coexisting components, individual atoms and diatomic homo-nuclear molecules, is worked out. In Figure 1, the model is represented schematically. ω 0 − ω a b Figure 1. Schematic representation of the model used in this work for the atom-diatomic coexistence. This is a two-level model. Diatomic molecules (b) are in the lower level, while single atoms (a) are in the upper level. The quantity ω 0 − ω represents the energy needed for separating the molecule into its two single atoms. This figure was taken from [16].
Each component in the model is represented in terms of bosons. There are two boson types: a and b. Boson type a represents individual atoms of energyhω 0 /2, whereas b−bosons represent diatomic molecules of energyhω. Atoms and molecules interact, and the proposed Hamiltonian is (h = 1 is used in this work) [16] where is the total number of atoms and is a conserved quantity. This magnitude gives the size of the system. Moreover,n a = a † a is the particle number operator of bosons of type a (atoms) andn b = b † b is the particle number operator of b−bosons (diatomic molecules).
The expected values of these operators are the numbers of particles of the two boson types (n a and n b ). To make everything simpler, only even M-values are considered in this work. Furthermore, λ is a control parameter that drives the system from one phase to the other. Given that ω 0 > ω, for λ = 0 the ground state is just a molecular condensate without any single atom. However, as λ increases, the interaction produces a more balanced atommolecule distribution. Thus, depending on the control parameter λ, the system has two phases: one with just molecules and another with molecules and atoms.

Exact Solution of the Eigenvalue Problem
An obvious basis for studying the Hamiltonian (1) is obtained by giving the number of molecules n b and the number of individual atoms n a : |n a , n b . Since M = 2n b + n a is conserved, one can use alternatively the notation |M, n b with M being fixed and defining the system.
The matrix elements of (1) in the mentioned basis are trivial and produce a tridiagonal matrix that can be easily diagonalized for each selected M-value. The relevant matrix elements are: For a given M, the matrix to be diagonalized is of dimensions (M/2 + 1), since one can have from zero molecules (only M atoms) to M/2 molecules (no atoms). A simple diagonalization of the corresponding tridiagonal matrix will provide one with all Hamiltonian eigenvalues and eigenfunctions. In particular, given an M-number, this diagonalization allows one to obtain the ground state energy and the corresponding wavefunction as a function of the control parameter λ. This can be used to study the ground state phase transition of the system as a function of λ.
In particular, once the ground state wavefunction is obtained, |gs(λ) , one can use it to calculate the expected value of the number of atoms gs(λ)|n a |gs(λ) . We show in the next subsection that this magnitude behaves as an order parameter. It is zero in one phase and different from zero in the other one. In a later section we will show other observables that can be used as markers for the critical point of the phase transition. In order to study a reference for the phase transition, a mean field study of the model is presented next.

Mean Field for the Model Hamiltonian
In order to develop a mean field study for this model, it is useful to introduce the operators: which close under the su(1, 1) commutation relations. Using the Holstein-Primakoff expansion [19], a new c−boson can be introduced as In terms of bosons b and c, the Hamiltonian (1) can be written as To perform a semiclassical analysis of the system, the usual relation with atom coordinates and momenta (x, p) and diatomic molecule coordinates and momenta (y, q) from the harmonic oscillator are introduced: These are canonical transformations, and in the thermodynamic limit, i.e., M → ∞, the operator's position and momentum commute. In addition, in this limit the factor 1/2 can be negligible in comparison with a term multiplied by M. Then, by introducing these relations in the Hamiltonian, it is written as with To analyze the properties of the ground state of the system in the thermodynamic limit, it is useful to rewrite the Hamiltonian (16) using polar coordinates: x = r cos α ; p = r sin α ; y = s cos β ; q = s sin β.
Then it can be shown that the Hamiltonian (16) can be written as It seems clear from this Hamiltonian that the minimum energy corresponds to cos(α − β) = −1 (it is the value that makes the second term, and therefore, the H minimum, since the other terms are positive, r 2 and s 2 ). This corresponds to α − β = π. Any choices of α and β such that they differ by π give the minimum energy. A possible combination is α = 0 and β = −π, which correspond to p = 0 and q = 0, but any other selection of α and β (and, correspondingly, of p and q) that satisfies α − β = π will give the same minimum energy surface per particle: which is Equation (22) in terms of x and y, taking α = 0 and β = −π (or equivalently This equation can be obtained in a more straightforward way from Equation (2) using coherent boson states. However, it is interesting to illustrate some tools, such as those presented above, that potentially can be used to extract finite size effects in the system (expanding the potential energy surface in 1/M powers), thereby going beyond the mean field description. Anyway, Equation (22) gives the classical energy surface associated with this model.
On the other hand, the M conservation leads to the condition which for p, q → 0 gives x 2 + y 2 = 1. This allows us to reduce the original two dimensional problem to another one with only one effective degree of freedom, x ∈ [−1, 1]. Taking into account that the sign selection y = − √ 1 − x 2 produces lower energy than the positive sign, the energy surface per particle can be written as where ∆ω = ω 0 − ω.
The condition for the minimum is and provides two solutions, This last solution provides energy lower than x 1 = 0 when λ is larger than a critical value that we call λ c . It is also a solution of the problem of the expression of x 2 with a positive sign in front of the square root, but the written expression, with the minus sign, always gives lower energy. A value for λ c can be obtained to make x 2 = 0: which gives the critical point for the transition, For given values of ω and ω 0 (this fixes λ c ), the minimum energy per particle as a function of λ is obtained: Equation (29), with x 2 from Equation (26), gives an analytic expression for the minimum of the energy surface per particle as a function of the control parameter λ. In Figure 2, the large-M limit of the ground state energy per particle (panel a), its first derivative (panel b) and its second derivative (panel c) are represented, respectively, for the cases ω 0 = 2 and ω = 1. In the three plots it is clear that at λ = 0.5 there is a structural change in the system. This value is the λ c given in Equation (28). Furthermore, the discontinuity of the second derivative indicates that this is a second-order (or continuous) phase transition. From Figure 2 it is clear that the system undergoes a second-order QPT at λ c .
Since we can solve the problem exactly for a finite M, in Figure 3 the mean field result for the ground state energy per particle is represented, together with the exact numerical calculations with M = 50 and M = 700 for the cases ω 0 = 2 and ω = 1 that produce λ c = 0.5. The mean field calculation is depicted by a black line; the exact numerical results for M = 50 are shown as a green line; and for M = 700, the exact numerical results are shown as a dashed red line. For a system of size M = 50, the exact numerical result fits quite well with the analytical mean-field, except in a small region close to the critical point (finite-size effects). Nevertheless, the bigger the system is, the better the agreement with the mean field calculation. This is shown in Figure 3 for M = 700, which is basically indistinguishable from the mean-field result. In order to show better the convergence, an inset is included in Figure 3 representing a function ε defined as: as a function of λ for different M sizes.

Figure 2.
Large-M limit (mean-field) results of the system as a function of the control parameter λ for the cases of ω 0 = 2 and ω = 1. In panel (a) the ground state energy per particle is represented.
In panel (b), its first derivative with respect to λ is plotted. Finally, in panel (c), the second derivative of the ground state energy per particle is given as a function of λ. The system undergoes a QPT for λ c = 0.5. In addition to the energy, one can calculate analytically at the mean field level (large M limit) the expected value for the number of atoms of type a. From Equations (9) and (12), one gets the relation n a = 2n c . Using the definitions of c † and c as a function of x and p, and by taking the classical limit (p → 0 and [x, p] = 0), one obtains easily that the number of individual atoms per particle, n a /M, is x 2 . The expected value of this observable in the system ground state is then with x 2 being given in Equation (26). This expression can be compared with real finite-M calculations to check how fast is the convergence to the large-M limit, and consequently how large are the finite-M effects.
In Figure 4, the large M limit of the expected value ofn a /M in the ground state is plotted as a function of λ. It is clear that this observable is an order parameter for the phase transition, since it is zero in one phase and different from zero in the other one. When λ → ∞, this order parameter tends to 2/3, as given by Equations (31) and (26). The latter means means that there would be a coexistence phase of atoms and molecular particles in which it is equally likely for an atom to either be chemically bonded or to remain unbound. In Figure 5, the exact expected value forn a /M is presented for M = 700, along with three different selections for ∆ω: 1, 2 and 3, that lead to λ c : 0.5, 1.0 and 1.5, respectively. For this large M value, the plots match the mean field result given by Equations (31) and (26). It is clearly seen in Figures 4 and 5 that this order parameter marks the critical point (represented in Figure 5 with filled dots for each ω selection).  In all cases, we have checked that the numerical results tend to the mean field approximation expressions as M is increased, and that the critical point corresponds to Equation (28).
In our model, as in the Tavis-Cummings and Jaynes-Cummings models [20], quantum fluctuations are zero, and consequently, these fluctuations cannot be responsible for the corresponding vacuum instability. In this respect, some researchers consider that this is not a quantum phase transition. However, this model possesses non-analyticity in the ground state, in agreement with a continuous quantum phase transition. As such, it is a matter of taste whether the transition should be termed quantum or not.

Other Markers for the QPT
In this section, we propose other markers for the critical point in the QPT.

Inverse Participation Ratio
The inverse participation ratio (IPR) is defined as This magnitude measures the degree of delocalization of a quantum state within a specific basis. The coefficients c (k) i are the coefficients of the state k in the used basis. On one hand, in case of full localization, the k state is one of the basis states, so only one c i = 1 and the IPR will be close to 1. On the other hand, if the state k is equally distributed among all basis states, then the normalization condition is with D being the dimension of the matrix diagonalized. Then In this case, the maximum IPR is obtained IPR max = M/2 + 1. Consequently, any values of IPR between 1 and M/2 + 1 are expected in general. For the model discussed here, an IPR = 1 is expected for λ = 0, since in this case our Hamiltonian eigenstates are those of the harmonic oscillator. For other λ−values, the Hamiltonian eigenstates will be a mixture of harmonic oscillator states and the IPR will increase. However, not every state of the harmonic oscillator will "participate" in the eigenstate of the coexistence phase. Only a linear combination of states in which the expected number of atoms is 2/3 of M will contribute. Thus, IPR will reach a constant but smaller size than the maximum possible value.
The numerical results from the exact diagonalization of the system Hamiltonian have already been presented, and these were compared to the mean-field results in the preceding section. For a given M, the exact diagonalization produces the ground state, and consequently provides the coefficients c gs M,n b . With these, one can calculate the IPR (32). In Figure 6, the ground state IPR values for M = 700 and for different ∆ω choices as a function of the control parameter λ are presented. The IPR marks clearly the transition of the system at the corresponding λ c . The ground state is well localized (small IPR) in the harmonic oscillator basis for λ below to the QPT critical point, whereas it tends to be delocalized for values of λ above the critical value. Indeed, an abrupt change of IPR occurs at λ c = ∆ω/2 in the QPT.
A natural question in relation to Figure 6 is: what are the asymptotic values for λ → 0 and λ → ∞? In order to show the λ → ∞, we plot in Figure 7 the IPR for the cases ∆ω = 1 and M = 700. It is seen that the IPR for large λ is around 30.
Whilst an IPR = 1, or approximately 1, is expected for λ ≤ λ c , for which the ground state is close to a molecular condensate (the state is basically |M, n a = 0, n b = M/2 ). For larger λ−values the Hamiltonian eigenstates will be mixtures of harmonic oscillator states, in which there will be more than one relevant state, and therefore, the IPR will increase. The limit of M/2+1 is obtained when all the basis states are contributing with equal weight. However, this is not reasonable, and states in which the number of atoms is n a = 2M/3, and consequently n b = M/6 (we notice that n a + 2n b = M), are expected to have larger weights. In fact, if one assumes for the wavefunction coefficients a binomial distribution with D = M/2 components, |M, n a , n b , whose probability of n b is p = 1/3, the corresponding IPR would be around 31. Although the distribution in our ground state is not exactly binomial, something similar is expected. In that case, the IPR will not reach the maximum possible value, and an IPR around 31 is expected for M = 700. In Figure 8, the binomial distribution for D = 350 that corresponds to M = 700 (basis dimension 351) and p = 1/3 (which corresponds to n b = M/6 = D/3) is represented against n b (dashed red line). Superimposed is the plot for the calculated ground state wavefunction components squared for the case M = 700 and λ = 1000 (full blue line). The similarity of the distributions is clearly shown, and that is why the IPR value for the large λ limit is close to the corresponding binomial distribution (around 30 in the case of Figure 7). In order to show the behavior of the IPR as a function of the system size, we present in Figure 9 the IPR for different M-values. From this figure, we can observe that the bigger the size of the system, the sharper the change in the value of the IPR at the critical point. A final comment on the IPR maximum observed right after the critical point: This is seen in Figure 7. This exact same behavior is confirmed to exist for all sizes. It is not more or less accentuated depending on M. We have already established which states are relevant in both the molecular condensate phase and the coexistence phase. However, right after the critical value is reached, the state of minimum energy is given by a linear combination of a set of states wider than the one observed for large λ values. It is a sort of transition region in which more components (fluctuations) are participating in the ground state wavefunction. As a consequence, the greatest value of the IPR is observed right there.

Renyi Entropy
Information was first defined rigorously by Claude Shannon [21]. It is a magnitude that measures how much communication it takes to transmit a message. If one has a discrete list of possible messages (events) with different probabilities, that wishes to transmit, the information value of every message depends on that probability. For instance, if one were to repeat the same message over and over, the information transmitted is measured with lower units of information. Conversely, if within this list of repeated messages something different is suddenly communicated only once, it is considered to give much more information. In other words, information measures how surprising, how unlikely, an event is. Thus, information theory does not account for content or usefulness, rather it measures only the quantity of information. The later is measured by a magnitude called entropy.
Different entropies can be defined. The most popular entropy was defined by Shannon [21], where Q are the generalized coordinates (x and y for our model), and ρ(Q) is the probability density (|Φ(x, y)| 2 , in our case). Then, Here we propose to use the Rényi entropy [22,23] that depends on one parameter α, for characterising the phase transition in our system. The Rényi entropy has as a limit situation the Shannon entropy (α → 1). The Rényi entropy is defined as that for the model discussed is, For the model under study, the ground state is a combination of harmonic oscillator states in the coordinates (x, y), where N are normalization constants and H n are Hermite polynomials. The ground state coefficients c n a ,n b are obtained from the Hamiltonian diagonalization. Consequently, the entropies can be calculated with Equation (39). However, this is computationally inefficient since for large M values it implies factorials of large numbers and make the calculation very heavy and inaccurate. As of that, we prefer to go to Shannon's original idea. The entropy [21] of a state describing a physical system is a quantity expressing the diversity, uncertainty or randomness of the system. Shannon viewed this uncertainty attached to the system as the amount of information carried by its state. His idea was based on the following consideration. If a physical system has a large uncertainty and one receives information on the system, then so-obtained information is more valuable (because it is less likely) than received from a system having less uncertainty. This is why entropy is measured in units of information. Shannon also drafted in A mathematical theory of communication [21] what is one of the most popular definitions of entropy. Let n b be a discrete random variable with probability distribution {p i } of N elements. That is then Shannon entropy is given by When one takes the binary logarithm, entropy is expressed in shannons (Sh), also known as bits. Moreover, when taking the natural logarithm, as we do in this work, entropy is expressed in the natural unit of information or nat. It is merely a difference in scale (1 Sh ≈ 0.693 nat). Note that if one event is much more likely than the others, that is p j → 1 and p i =j → 0, then entropy tends to 0. In the opposite case, if all events were equally likely, then p i = 1/N, ∀i and S = log N, which is a function that increases with N. Furthermore, take notice of the fact that both the maximum and minimum possible values of Shannon entropy correspond to maximum and minimum values of IPR.
Should one need to measure the information provided by events giving it greater or lesser differences between likely and unlikely ones, a different definition of entropy would have to be used.
A generalization of Shannon entropy was made by Alfred Rényi [22]. Classical Rényi entropy for a parameter α ≥ 0 and α = 1 is defined for the same discrete random variable as The same minimum and maximum possible values of entropy Rényi are reached, independently of α. In fact, the limiting value of Rényi entropy as α → 1, that can be calculated using L'Hôpital's rule, is the Shannon entropy S = lim α→1 R (α) .
In the context of quantum theory of information, for a density matrix in a Hilbert space, ρ ∈ N (H), we can define quantum Rényi entropy [24] as If {p i } are the diagonal elements of ρ in the basis of eigenfunctions, then the quantum Rényi entropy reduces to a Rényi entropy of a random variable n b as defined in (41). This means that for the ground state of our system we can define the probabilities p i = |c i | 2 where c i are the coefficients of the ground state wavefunction. Note that we already took the dimension of the Hamiltonian matrix N as the number of elements in the discrete random distribution.
For α < 1, all random events are weighted more equally, resulting in a smaller change in entropy from one state to another. As α tends to zero, the entropy is just the logarithm of the size of the support of n b , no matter the phase.
For α > 1, all random events are weighted more differently. As α grows, more likely events make larger contributions to entropy, whereas less likely events are disregarded. This tends to give bigger differences between quantum phases.
For α → 1, we have Shannon entropy, which results in something in the middle of both cases.
Besides α dependency, there is another factor which affects entropy values. As happens with thermodynamic entropy, quantum entropy is an extensive property, meaning that it scales with the size of the system. This behavior has already been hinted at by substituting into Shannon entropy a set of values equally likely.
On account of the above, results are expressed according to the following criteria.
• Entropy dependency with M All calculations have been done for ∆ω = 1, which gives λ c = 0.5. In Figures 10-12 we can observe the dependency of different entropies with M. The transition is sharper with increasing M for all values of α. Note that entropy is independent of M in one phase but is increasingly different with larger M-values in the other phase. The reason behind this phenomenon lies in the characteristics of both phases.
In the first one, the possibility of measuring the lowest eigenvalue of the harmonic oscillator (p 0 = |c 0 | 2 ), is almost 1 and the rest are almost zero (which is why the IPR is approximately 1). Since the dimension N is irrelevant (to a certain point), because it would not really matter how many p i there are, entropy values will be very similar and will mostly depend on α. In the second phase, we need to reach a certain proportion of particles, given by a number of relevant coefficients that is proportional to M. This is the reason why the IPR also increases with M.

•
Entropy dependence with α In Figure 13, for a system of M = 700, the values of R α (λ) are represented for a set of α values both under and over the unit as well as the Shannon entropy, which is given by the limit α → 1. It is confirmed that bigger values of α make the difference between both phases more evident since it distinguishes more abruptly between likely and unlikely events. It can be also confirmed that Shannon entropy is indeed between the entropy for α < 1 and α > 1.
Perhaps plenty more examples and evaluations could be made by toying with different values of M and α. However, the most important conclusion one should make is the following. The entropy, when set to an adequate α for it to be a good marker, is yet another quantum magnitude that experiences an abrupt (but continuous, as seen for lower α values) change from one phase to another, evidencing the existence of a second-order QPT.

Conclusions
We have studied a two-level model for the coexistence of atoms and diatomic molecules. This model has been studied using mean-field techniques and shows a ground state secondorder quantum phase transition. The critical point has been obtained for the large Mnumber limit. Analytic expressions for the ground state energy per particle and for the number of atoms per particle, as a function of the control parameter λ, have been worked out. This last observable was shown to be a good order parameter. We have proposed as additional markers for the phase transition the inverse participation ratio (IPR) and different types of entropies. Both observables clearly mark the critical phase transition point. Funding: This work is part of the I+D+i projects PID2019-104002GB-C22 and PID2020-114687GB-I00 funded by MCIN/AEI/10.13039/501100011033. This work is also part of grant Group FQM-160, EU FEDER funds US-1380840 and the project PAIDI 2020 with reference P20_01247, funded by the Consejería de Economía, Conocimiento, Empresas y Universidad, Junta de Andalucía (Spain) and "ERDF-A Way of Making Europe".