text:
# Specific heat and magnetic measurements in Nd0.5Sr0.5MnO3, Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3 samples

## I Introduction

Compounds like La0.5Ca0.5MnO3 and Nd0.5Sr0.5MnO3 present a real-space ordering of Mn3+ and Mn4+ ions, named as charge ordering (CO). Close to the charge ordering temperature (TCO), these materials show anomalies in resistivity, magnetization and lattice parameters as a function of temperature, magnetic field and isotope mass,. At low temperatures both, ferromagnetic and antiferromagnetic phases, could coexist . However, a relatively small external magnetic field can destroy the CO phase and enforces a ferromagnetic orientation of the spins.

Moreover, electron microscope analysis has revealed convincing evidence that CO is accompanied by orientational ordering of the 3d3 orbitals on the Mn3+ ions, called orbital ordering (OO). The physical properties in CO manganese perovskites are thought to arise from the strong competition among a ferromagnetic double exchange interaction, an antiferromagnetic superexchange interaction, and the spin-phonon coupling. These interactions are determined by intrinsic parameters such as doping level, average cationic size, cationic disorder and oxygen stoichiometry. Microscopically, CO compounds are particularly interesting because spin, charge and orbital degrees of freedom are at play simultaneously and classical simplifications that neglect some of these interactions do not work. More detailed information on the physics of manganites can be found in a review paper by Myron B. Salamon and Marcelo Jaime.

We have shown that polycrystalline samples of La0.5Ca0.5MnO3 and Nd0.5Sr0.5MnO3 presented an unusual magnetic relaxation behavior close to each critical temperature, . However, a clear understanding of all these features has not been reached yet. An alternative to a bulk characterization like magnetization would be to perform specific heat measurements. In contrast to magnetization, which has a vector character, the specific heat is an scalar. Therefore, a comparison between both types of data could give valuable information.

J. E. Gordon et al. reported specific heat measurements for a Nd0.67Sr0.33MnO3 sample and found a Schottky-like anomaly at low temperatures. They associated this result to the magnetic ordering of Nd3+ ions and the crystal-field splitting at low temperature. F. Bartolomé et al. also found Schottky-like anomaly in a closely related compound of NdCrO3. They proposed a crystal-field energy level scheme in agreement with neutron-scattering studies in the same sample. In two papers V. N. Smolyaninova et al.,  studied the low temperature specific heat in Pr1-xCaxMnO3 (0.3$`<`$x$`<`$0.5) and La1-xCaxMnO3 (x=0.47, 0.5 and 0.53). They found an excess specific heat, C´(T), of non-magnetic origin associated with charge ordering. They also reported that a magnetic field sufficiently high to induce a transition from the charge ordered state to the ferromagnetic metallic state did not completely remove C´(T). However, no Schottky anomaly was found in any of these compounds.

Here, we report magnetic and specific heat measurements with applied magnetic fields between 0 and 9 T and temperatures between 2 and 300 K for Nd0.5Sr0.5MnO3, Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3 samples. All these compounds presented a Schottky-like anomaly at low temperatures, as well as characteristic peaks associated to the charge ordered, antiferromagnetic and/or ferromagnetic transitions. We have already reported a short version of preliminary results about these topics . However, as far as we know, detailed specific heat measurements in these compounds have not been published.

## II Experimental methods

Polycrystalline samples of Nd0.5Sr0.5MnO3, Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3 were prepared by the sol-gel method. Stoichiometric parts of Nd2O3 (Ho2O3) and MnCO3 were dissolved in HNO3 and mixed to an aqueous citric acid solution, to which SrCO3 or CaCO3 was added. The mixed metallic citrate solution presented the ratio citric acid/metal of 1/3 (in molar basis). Ethylene glycol was added to this solution, to obtain a citric acid/ethylene glycol ratio 60/40 (mass ratio). The resulting solution was neutralized to pH$``$7 with ethylenediamine. This solution was turned into a gel, and subsequently decomposed to a solid by heating at 400 oC. The resulting powder was heat-treated in vacuum at 900 oC for 24 hours, with several intermediary grindings, in order to prevent formation of impurity phases. This powder was pressed into pellets and sintered in air at 1050 oC for 12 hours. X-ray diffraction measurements indicated high quality samples in all cases.

The magnetization measurements were done with a Quantum Design MPMS-5S SQUID magnetometer. Specific heat measurements were done with a Quantum Design PPMS calorimeter. The PPMS used the two relaxation time technique, and data was always collected during sample cooling. The intensity of the heat pulse was calculated to produce a variation in the temperature bath between 0.5 % (at low temperatures) and 2% (at high temperatures). Experimental errors during the specific heat and magnetization measurements were lower than 1 % for all temperatures and samples.

## III Results and Discussion

### A Magnetization measurements

Figure 1 shows the temperature dependence of magnetization, measured with a 5 T applied magnetic field and zero field cooling conditions, in polycrystalline samples of Nd0.5Sr0.5MnO3, Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3. All three compounds present a charge ordering transition below 160, 250 and 271 K, respectively , , . It is interesting to note that the relation between the charge ordering temperature and the antiferromagnetic ordering temperature changes from one sample to the other. In the first case they are approximately coincident (T$`{}_{CO}{}^{}`$TN), in the second case the charge ordering temperature is much higher (T$`{}_{CO}{}^{}`$TN), and in the third case there is not a long range antiferromagnetic transition.

The Nd0.5Sr0.5MnO3 sample presents a ferromagnetic transition at TC=250 K and an antiferromagnetic transition at TN=160 K. The Nd0.5Ca0.5MnO3 compound presents a strong maximum near TCO, but shows an unexpected minimum close to the antiferromagnetic transition temperature TN=160 K. Usually an antiferromagnetic transition is accompanied by a maximum in the temperature dependence of the magnetization. For temperatures lower than 20 K both Nd0.5Sr0.5MnO3 and Nd0.5Ca0.5MnO3 samples show a sharp increase in the magnetization. This sharp increase in magnetization have been associated to a short range magnetic ordering of the intrinsic magnetic moment of Nd3+ ions. However, no long range ferromagnetic order of the Nd3+ ions was found in neutron diffraction measurements at these low temperatures, .

Differently from the two previous samples, the Ho0.5Ca0.5MnO3 compound do not present a strong maximum at the charge ordering temperature in the magnetization versus temperature curve. The existence of charge ordering in Ho0.5Ca0.5MnO3 was suggested by T. Terai et al.  studying simultaneously the magnetization and resistivity curves. Figure 2 shows the inverse of the DC susceptibility measured with 0.1 mT (inset) and 5 T (main figure) in the Ho0.5Ca0.5MnO3 sample. The measurement at 5 T clearly indicates a transition between two linear behaviors near 270 K. This latter result, together with the specific heat measurements, that we will discuss below, reinforces the idea of the existence of a charge ordering transition in the Ho0.5Ca0.5MnO3 sample, although the details of the interactions in this compound differ from those in Nd0.5Sr0.5MnO3 and Nd0.5Ca0.5MnO3. On the other hand, the DC susceptibility measurement in Ho0.5Ca0.5MnO3 with a probe field of 0.1 mT indicates two possible Curie-Weiss intervals: below and above 100 K (see straight lines).

We fitted the magnetization data, measured with a probe field of 0.1 mT for all samples, to a Curie-Weiss law: $`M/H\mu _{eff}^2/(TT_\mathrm{\Theta })`$, where $`M/H`$ was the DC susceptibility, $`\mu _{eff}`$ was the effective paramagnetic moment and $`T_\mathrm{\Theta }`$ was the Curie-Weiss temperature. A simplified theoretical estimate of the effective paramagnetic moment gives values of 5.1 $`\mu _B`$, 5.1 $`\mu _B`$ and 8.7 $`\mu _B`$ for Nd0.5Sr0.5MnO3, Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3 samples, respectively. The positive or negative value of $`T_\mathrm{\Theta }`$ suggests the predominance of local ferromagnetic or antiferromagnetic interactions between the spins in the corresponding interval of temperature used. However, it is important to remember that the Curie-Weiss law is only a mean field approximation, that do not consider fluctuations.

In our case the fitting interval was taken between 300 and 400 K in the Nd0.5Sr0.5MnO3 and Nd0.5Ca0.5MnO3 samples. The $`\mu _{eff}`$ and $`T_\mathrm{\Theta }`$ values were 2.2 $`\mu _B`$ and 2.1 $`\mu _B`$, and 260 K and 215 K for Nd0.5Sr0.5MnO3 and Nd0.5Ca0.5MnO3 samples, respectively. The discrepancies between our experimentally measured values of $`\mu _{eff}`$ and the theoretical ones in the samples of Nd0.5Sr0.5MnO3 and Nd0.5Ca0.5MnO3, are essentially due to the vicinity of the TC or TCO temperatures to the temperature interval that we could use to fit the curves to the Curie-Weiss law. F. Millange et al.  measured the effective paramagnetic moment, in the interval between 450 and 800 K, for a Nd0.5Ca0.5MnO3 sample and found a value roughly in agreement with theory. M. T. Causa et al. also showed that the transition to a Curie-Weiss regime in several manganites occurs approximately for T $`>`$2 TC . In the temperature interval approximately between 2 TC and TC magnetic ions start to form local clusters. This behavior is different from the one of a simple ferromagnetic sample like iron, where the transition from the ferromagnetic state to the Curie-Weiss regime occurs in a much shorter temperature interval.

For the Ho0.5Ca0.5MnO3 sample we fitted the Curie-Weiss law to two intervals: from 100 to 350 K and from 2 to 100 K. In the Ho0.5Ca0.5MnO3 sample the $`\mu _{eff}`$ and $`T_\mathrm{\Theta }`$ values were 9.7 $`\mu _B`$ and 3.4 $`\mu _B`$, and -38 K and -16 K for the high and low temperature intervals, respectively. The effective paramagnetic moment found in the high temperature interval for the Ho0.5Ca0.5MnO3 sample is roughly in agreement with the theoretical value, and the negative value of $`T_\mathrm{\Theta }`$ indicates that the main local interactions are antiferromagnetic. However, the transition approximately below 100 K to a different Curie-Weiss regime suggests that the predominance of short range antiferromagnetic interactions are giving space to short range ferromagnetic ones, with a higher value of $`T_\mathrm{\Theta }`$. This conclusion is also supported by the magnetic field dependence of the magnetization, as we will discuss below.

Figure 3 shows the magnetization hysteresis cycle at 2 K for the three studied samples. The applied magnetic field was increased from 0 to 5 T, decreased to -5 T and then increased back to 5 T again. The Nd0.5Sr0.5MnO3 curve is characteristic of a mixture of two phases: one ferromagnetic and another antiferromagnetic. The ferromagnetic part is easily oriented at low magnetic field values and shows a hysteretical behavior. The almost linear and reversible dependence for magnetic fields higher than approximately 1 T, indicates a gradual destruction of the antiferromagnetic phase. R. Mahendiram et al. reported, in a sample with the same composition at 50 K, that for magnetic fields higher than approximately 5 T, the magnetization started to increase rapidly, and for magnetic fields above 10 T, it slowly approached the ferromagnetic saturation value. Because at 2 K these transition fields are much higher than 5 T we were unable to see them. The Nd0.5Ca0.5MnO3 curve did not show any trace of a ferromagnetic component at low fields. In fact, the curve is linear and reversible for the whole magnetic field interval. As before, this linearity characterizes the gradual destruction of the antiferromagnetic phase. F. Millange et al. reported, for a sample with the same composition, a set of M vs. H curves with applied magnetic fields up to 22 T. They found sharp transitions and a large hysteresis at 130 K, between applied magnetic fields of 12 and 18 T. Their results were interpreted as evidence of the existence of a spin-flop transition.

In the M vs. H curve for the Ho0.5Ca0.5MnO3 sample, two characteristics are well noticed. One, that there is not hysteresis at all, and second, that the magnetization values at 5 T are well below the theoretical saturation value of 8.7 $`\mu _B`$. In order to try to discriminate between a paramagnetic or ferromagnetic ordering at 2 K, the experimental points were fitted by a Brillouin function:

$$M=(N/V)\gamma JB_J(x)$$
(1)
where $`\gamma =g\mu _B,\beta =1/k_BT`$ and

$$B_J(x)=[(2J+1)/2J]\mathrm{coth}[(2J+1)x/2J](1/2J)\mathrm{coth}(x/2J)$$
(2)
Here, $`x=\beta \gamma JH`$, $`N/V`$is the number of ions per unit volume, $`g`$ is an effective gyromagnetic ratio, $`\mathrm{J}`$ is the total angular momentum,$`\mathrm{T}`$ is the temperature, $`\mu _B`$ is the Bohr magneton and $`k_B`$ is the Boltzmann constant. Equation 1 is valid for a set of identical and not interacting ions of angular momentum J. Because our sample has three different magnetic ions we fixed J=8 (Ho3+), J=2 (Mn3+) or J=3/2 (Mn4+) and let the effective gyromagnetic ratio to change freely in each case. The best fitting, represented by the solid line in figure 3, was found for J=8 and g=0.35. Although it is difficult to enunciate a more conclusive statement from these results, they suggest that there is not a long range ferromagnetic order for the Ho0.5Ca0.5MnO3 sample, at 2 K.

### B Specific heat at high temperatures

Figure 4 shows specific heat measurements with a zero applied magnetic field from 2 to 300 K in the Nd0.5Sr0.5MnO3, Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3 samples. In order to facilitate the visualization, the curves for Nd0.5Sr0.5MnO3 and Ho0.5Ca0.5MnO3 were displaced 20 J/mol K upside and downside, respectively. Specific heat measurements give information about both lattice and magnetic excitations. At high temperatures the excitations from the lattice vibrations are dominant and decrease as the temperature decreases. The magnetic contribution can be obtained by subtracting the lattice part from the experimental values.

The Dulong and Petit model predicts that the thermal contribution at high temperature increases assintotically to the value $`C_{DP}=3nR`$ , where $`n`$ is the number of atoms in the compound´s formulae and $`R`$ is the ideal gas constant. Therefore, $`C_{DP}=`$125 J/(mol K) for all studied samples in this work, since their unit cell contain 5 atoms each. At 300 K the experimental specific heat values were 110 J/(mol K) (88 % of $`C_{DP}`$) for the Nd0.5Sr0.5MnO3 and Nd0.5Ca0.5MnO3 samples and 100 J/(mol K) (80 % of $`C_{DP}`$) for the Ho0.5Ca0.5MnO3 sample. These values suggest that the Debye temperature should be, whithin the experimental error, approximately equal for all three compounds.

Continuous lines in figure 4 represent the fitting of the thermal background, in the interval from 30 to 300 K, by the Einstein model given by:

$$C_{Einstein}=3nR\underset{i}{}a_i\left[\frac{x_i^2e^{x_i}}{\left(e^{x_i}1\right)^2}\right]$$
(3)
where $`x_i=T_i/T`$ . We used three optical phonons ($`i`$ =1, 2, 3) with energies $`\mathrm{T}_i`$ (in Kelvin) and relative occupations $`\mathrm{a}_i`$. The Einstein model for the specific heat considers the oscillation frequency (or energy) independently from the wave vector, which is a valid approximation for the optical part of the spectrum. The values of temperatures (energies) found were 148, 438 and 997 K for Nd0.5Sr0.5MnO3, 152, 432 and 1035 K for Nd0.5Ca0.5MnO3 and 147, 438 and 1023 K for Ho0.5Ca0.5MnO3. These values are similar to those reported, using the same model, by A. P. Ramirez et al. in a La0.37Ca0.63MnO3 sample and Raychaudhuri et al. in a Pr0.63Ca0.37MnO3 sample. However, we should point out that the values found are not unique because there are six varying parameters during the fitting (one energy and one occupation coefficient for each oscillation mode).

An alternative method to determine the oscillation mode energies of the crystalline lattice is to study the Raman spectrum of a given compound. For example, V. Dediu et al. found that, at a measured temperature of 50 K for a Pr0.65Ca0.35MnO3 sample, the strongest peaks were located at energies of 683 and 877 K. Besides, E. Granado et al. found that, at a measured temperature of 200 K for a Nd0.66Ca0.33MnO3 sample, the strongest peaks were located at energies of 388, 439, 690 and 877 K. However, these peaks were not the only ones, they could change with temperature and they are not narrow. Furthermore, the Raman spectrum could depend on variables like the frequency, the intensity of the excitation laser and the luminosity of the sample. All of these make difficult a straightforward comparison between both methods to determine the oscillation energies of a crystalline lattice.

Figure 5 represents the differences between the experimental data and the fitted curves in figure 4. There is a maximum at 231 K for the Nd0.5Sr0.5MnO3 sample, which is correlated with the ferromagnetic transition at 250 K in the corresponding magnetization curve (figure 1). A second maximum, which could be partially associated to the antiferromagnetic and charge ordering transitions at 160 K, appears at 180 K. However, lattice parameters in this compound change rapidly between approximately 110 K and 250 K. This variation in lattice parameters changes the intensity of the interactions between the atoms, and consequently the oscillation frequency of the phonons mode, contributing to the specific heat in the second peak. Figure 5 also shows that there is a maximum at 243 K for the Nd0.5Ca0.5MnO3 sample. This maximum correlates with the corresponding charge ordering temperature at 250 K. Differently from the Nd0.5Sr0.5MnO3 sample, there is not magnetic ordering in this high temperature interval for the Nd0.5Ca0.5MnO3 sample. In this latter compound lattice parameters change very rapidly between approximately 200 and 250 K . An inflection point in the C vs. T curve appears at 141 K for the Nd0.5Ca0.5MnO3 sample. Besides, there is a maximum at 141 K in figure 5, but its height is relatively small compared to the maximum at 243 K. The Neel temperature corresponding to this compound is 160 K. Therefore, these results lead us to conclude that the specific heat variations due to the antiferromagnetic order are small compared to those induced in the charge ordering and ferromagnetic transitions. For the Ho0.5Ca0.5MnO3 sample, we found a maximum at 276 K. This maximum is correlated with the transition seen at high magnetic fields in the inverse DC susceptibility curve of figure 2. These experimental results and the resistivity measurements reported by T. Terai et al., indicate the existence of a phase transition between 270 and 280 K. However, electron diffraction studies would be needed to unambiguously classify this transition as CO.

V. Hardy et al. reported similar measurements in a single crystal of Pr0.63Ca0.37MnO3. They found a specific heat peak at 229 K and resistivity and magnetization curves indicated a charge ordering transition at 227 and 228 K, respectively. They also point out that part of the peak amplitude was due to the crystalline lattice distortions. Moreover, the specific heat curve showed an inflection point at 150 K, close to the corresponding TN = 160 K. Similar characteristics were also found in a polycrystalline sample of La0.37Ca0.63MnO3.

Results in figure 5 allow us to calculate the variation in entropy ($`\mathrm{\Delta }`$S), associated to the charge ordering, ferromagnetic and antiferromagnetic transitions:

$$\mathrm{\Delta }S=\underset{Ti}{\overset{Tf}{}}\frac{\left(CC_{ph}\right)}{T}𝑑T$$
(4)
where $`\mathrm{T}_i`$ and $`\mathrm{T}_f`$ are two temperatures conveniently chosen to delimitate the interval of interest and $`\mathrm{C}_{ph}`$ is the specific heat due to the lattice oscillations.

For the Nd0.5Ca0.5MnO3 sample the entropy variation between 201 and 301 K was $`\mathrm{\Delta }S(`$T$`{}_{CO}{}^{})=`$ 2.0 J/(mol K). A. K. Raychaudhuri et al. reported an entropy variation, close to the charge ordering transition in the compound Pr0.63Ca0.37MnO3, of 1.8 J/(mol K) with zero applied magnetic field and 1.5 J/(mol K) with an 8 T magnetic field. On the other hand, Ramirez et al.  found $`\mathrm{\Delta }S(`$T$`{}_{CO}{}^{})=`$ 5 J/(mol K) in a La0.37Ca0.63MnO3 sample. All these results correspond to those expected for a charge ordering transition . On the other hand, the entropy variation (not related with phonons), calculated between 118 and 201 K for the Nd0.5Ca0.5MnO3 sample, was $`\mathrm{\Delta }S(`$T$`{}_{N}{}^{})=`$ 0.80 J/(mol K). We associate this smaller $`\mathrm{\Delta }\mathrm{S}`$ value to the antiferromagnetic order at T$`{}_{N}{}^{}=`$160 K.

The entropy variation between 133 and 274 K for the Nd0.5Sr0.5MnO3 sample was $`\mathrm{\Delta }S(`$T$`{}_{CO+FM}{}^{})=`$ 3.6 J/(mol K). Considering that the entropy variation associated to the charge ordering transition is the same as for the Nd0.5Ca0.5MnO3 sample, we found that the entropy variation associated to the ferromagnetic transition is approximately 1.6 J/(mol K). Using a model proposed by J. E. Gordon et al.  we estimated that the change in entropy needed to a full ferromagnetic transition in the Nd0.5Sr0.5MnO3 sample was 12.45 J/(mol K). Therefore, the fact that the entropy variation found is 13 % of the theoretical value, confirms that only a small part of the spins order ferromagnetically. The same report of J. E. Gordon et al.  found that the entropy variation associated to the ferromagnetic order in a Nd0.67Sr0.33MnO3 sample was approximately 10 % of the theoretical value. The entropy variation, not associated to lattice oscillations, between 216 and 294 K for the Ho0.5Ca0.5MnO3 sample was $`\mathrm{\Delta }S(`$T$`{}_{CO}{}^{})=`$ 0.78 J/(mol K). This $`\mathrm{\Delta }S(`$T$`{}_{CO}{}^{})`$ for the Ho0.5Ca0.5MnO3 sample is smaller than the corresponding value for the Nd0.5Ca0.5MnO3 sample, suggesting a different nature of the charge ordering transition.

### C Specific heat at low temperatures

For temperatures lower than 30 K, the phonon contribution to the specific heat starts to take a smaller part, and the C vs. T curve shows a slow variation. Figure 6 reproduces the specific heat measurements for temperatures between 2 and 30 K in the samples of (a) Nd0.5Sr0.5MnO3, (b) Nd0.5Ca0.5MnO3 and (c) Ho0.5Ca0.5MnO3. Measurements were made in the presence of applied magnetic fields of 0 T (open squares), 5 T (closed circles), 7 T (open up triangles) and 9 T (closed down triangles). Note that close to 5 K all curves show a Schottky-like anomaly.

It is important to stress that high values of specific heat were found at low temperatures for the three samples presented here. These values are similar to those reported by J. E. Gordon et. al. in a sample of Nd0.67Sr0.33MnO3. However, the absolute values of specific heat at 2 K, reported by J. J. Hamilton et al. in samples of La0.67Ba0.33MnO3 and La0.80Ca0.20MnO3, were more than 100 times smaller. Similarly, V. Hardy et al., in a single crystal of Pr0.63Ca0.37MnO3, found values approximately equal to those reported by J. J. Hamilton et al.. The high values of specific heat found in our work could be interpreted as an increase in the effective mass of the electrons due to localization, which is also consistent with the insulating behavior revealed by electrical resistivity measurements, , .

Continuous lines in figure 6 indicate the fitting of the experimental data between 15 and 30 K by the following expression:

$$C=\beta _{2n+1}T^{2n+1}$$
(5)
Here, $`\mathrm{C}`$ is the specific heat, $`\mathrm{T}`$ is the temperature and the parameters $`\beta _{2n+1}`$ represent the contribution of the phonon modes. Notice that we did not include the lowest temperature interval to avoid the Schottky anomaly. To be able to fit the whole temperature interval we have chosen values of $`n`$ from 1 to 4. Nonetheless, this large number of free parameters difficult an unique determination of each one. Since from resistivity measurements, ,  all the studied samples show an insulating behavior at low temperatures, and the applied magnetic fields are not strong enough to destroy this characteristic, the expected linear contribution, from the free electrons to the specific heat, is zero. However, other kind of excitations could lead to a linear contribution. This can imply an implicit error of the fitting model. Moreover, we could not resolve in our data a term of type T3/2. This term is usually interpreted as an evidence of the existence of ferromagnetic interactions. However, previous studies in samples clearly identified by other techniques as ferromagnetic, have not found this term in the specific heat either, .

The values of $`\beta _3`$ change between 0.28 mJ/(mol K4) at H=0 T in Nd0.5Sr0.5MnO3 to 1.57 mJ/(mol K4) at H=5 T in Ho0.5Ca0.5MnO3. The corresponding Debye temperatures (TD), obtained from $`\beta _3`$, will be plotted in figure 7a. The graphs in the insets of figure 6 show the differences between the specific heat experimental data and the phonon contribution to the specific heat, extrapolated to low temperatures from the fitting in the temperature interval between 15 and 30 K.

Figure 7 shows (a) the magnetic field dependence of the Debye temperature, (b) the variation of magnetic entropy between 2 and 20 K ($`\mathrm{\Delta }S`$) and (c) the Schottky temperature (TS) in the three studied samples. The entropy variation was calculated from equation 4 and the Schottky temperature was determined from the maxima

in the insets of figure 6. The Debye temperature was calculated using the values of $`\beta _3`$ and the following equation:

$$T_D=\left(\frac{12\pi ^4nR}{5\beta _3}\right)^{1/3}$$
(6)
where $`\mathrm{n}`$ is the number of atoms in the unit cell and $`\mathrm{R}`$ is the ideal gas constant. We should point out that the Debye temperature was estimated from the low temperature data. This procedure lead to a Debye temperature that is slightly different than the actual value of temperature for which the specific heat saturates. The Debye temperature from the low temperature data decreases with the increase of the applied magnetic field. We have also made specific heat measurements with a 9 T magnetic field, at high temperatures, for several charge ordered compounds, and they show an increase of the specific heat when compared to the 0 T case, in agreement with the magnetic field dependence of the Debye temperatures. Other authors have made an initial assumption that the Debye temperature is magnetic field independent, which is not supported by our experimental results.

Figure 7 also shows that TS grows with the increase of the external magnetic field in all cases. However, the growth of TS seems to be saturated for a magnetic field of 5 T in the sample of Ho0.5Ca0.5MnO3, while there is no sign of TS saturation in the other two compounds. It is also interesting to note here the relative low TS values. For the reagent compounds of Nd2O3 and Ho2O3 the peak in the specific heat, measured in zero magnetic field, were found at approximately 10 and 9 K, respectively. The fact that TS is lower in the manganese compounds suggests that the collective charge ordered phase could be determining a smaller splitting in the energy levels.

The entropy variation, associated to the Schottky anomaly, grows as a function of magnetic field in the sample of Nd0.5Sr0.5MnO3. The same result is clearly visualized from the height of the Schottky anomaly in the inset of figure 6a. J. E. Gordon et al. also reported a similar increase in a sample of Nd0.67Sr0.33MnO3. However, the Schottky entropy variation decreases with the increase of magnetic field in the Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3 samples. In these two cases, the C vs. T curves in figure 6b and 6c, indicate that the local minimum, at a temperature above the Schottky anomaly, disappears with the increase of the external magnetic field. This is also reflected in the decrease of the height of the peak with the increase of the applied magnetic field (insets of figures 6b and 6c). These results seem to indicate that the magnetic field dependence of the Schottky anomaly is more correlated with the presence (Nd0.5Sr0.5MnO3) or not (Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3) of ferromagnetic interactions than with the particular type of magnetic ion (Nd3+ or Ho3+).

The expected entropy variation from the magnetic ordering of Nd3+ or Ho3+ ions could be estimated as $`\mathrm{\Delta }S`$=0.5 R ln(2), where $`\mathrm{R}`$ is the ideal gas constant. The actually found variation correspond to values from 63 to 77 % of the expected ones in the Nd0.5Sr0.5MnO3 sample, 80 to 62 % in the Nd0.5Ca0.5MnO3 sample and 68 to 25 % in the Ho0.5Ca0.5MnO3 sample, for magnetic fields between 0 and 9 T, respectively. J. E. Gordon et al. found that the entropy variation associated to the ordering of Nd3+ ions in Nd0.67Sr0.33MnO3 was approximately 85 % of the expected value.

Let us consider that Nd (Ho) ions be oriented by a molecular field interaction (Hmf), and not by the exchange interaction between pairs of Nd-Nd ions (Ho-Ho). Assuming that Hmf does not change with the external magnetic field, it is possible to estimate it, using a mean field model and the peak temperature in the specific heat. Considering a two level energy splitting $`\mathrm{\Delta }(\mathrm{H})`$, due to a magnetic moment $`\mathrm{m}`$ in an external magnetic field $`\mathrm{H}`$, one finds in zero applied field $`\mathrm{\Delta }(0)=2mH_{mf}`$, and in $`\mathrm{H}=9`$ T the value changes to $`\mathrm{\Delta }(9`$ T$`)=2m[H_{mf}+9`$ T$`]`$. One can also use that the energy splitting could be related to the peak temperature in the specific heat by $`\mathrm{\Delta }=\mathrm{k}_B\mathrm{T}_S/0.418`$, a relation valid for a two level Schottky function. Solving this system of two linear equations we found that $`\mathrm{H}_{mf}=11.4`$ T and $`\mathrm{m}=0.43\mu _B`$ in Nd0.5Sr0.5MnO3, $`\mathrm{H}_{mf}=20.6`$ T and $`\mathrm{m}=0.44\mu _B`$ in Nd0.5Ca0.5MnO3 and $`\mathrm{H}_{mf}=13.6`$ T and $`\mathrm{m}=0.58\mu _B`$ in Ho0.5Ca0.5MnO3. These values of $`\mathrm{m}`$ are smaller than those obtained from susceptibility measurements at high temperatures. However, they are similar to the one found by J. E. Gordon et al. using the same method ($`H_{mf}=`$10 T and $`m=`$ 0.8 $`\mu _B`$ in Nd0.67Sr0.33MnO3).

The ground state of the Nd3+ and Ho3+ ions are usually denoted as 4I9/2 and 5I8, where I stands for an orbital angular momentum L=6, the superprefix specify the total spin as $`2\mathrm{S}+1`$ and the subscript the total angular momentum $`\mathrm{J}`$. The number of the lowest energy levels is given by $`2\mathrm{J}+1`$, which leads to 5 Kramers doublets in the ground state of the first ion and a singlet and 8 Kramers doublets in the second ion. F. Bartolomé et. al. showed that the second doublet in the Nd3+ ion was approximately 120 K (in energy) above the lowest doublet. As this temperature is about 10 times higher than the temperature where the Schottky anomaly appears, the contribution of the second doublet is expected to be small. In a previous report we showed that a two level Schottky function (only one doublet) did not fit properly our experimental data at low temperatures. The same result was verified for the new experimental data presented here. One alternative, justified by the existence of several different grains in polycrystalline samples, is to consider a distribution of energy splitting around the value that would correspond to a single crystal in the same two level Schottky model. Although the fitting results using this second approach improved a little bit, we found that they still remained unsatisfactory.

At first sight someone might be tempted to correlate the existence of the Schottky anomaly with the presence of an intrinsic magnetic moment in Nd3+ and Ho3+ ions (in contrast with La3+ ions without magnetic moment and no Schottky anomaly in the manganite). However, specific heat measurements reported by V. Hardy et al. in a compound of Pr0.63Ca0.37MnO3 (Pr3+ ions have approximately the same magnetic moment as Nd3+ ions) did not show any Schottky anomaly. Moreover, Ho3+ ions have an intrinsic magnetic moment almost 3 times bigger than Nd3+ ions, but the Schottky temperature at zero magnetic field were 2.73 K in Nd0.5Sr0.5MnO3, 5.08 K in Nd0.5Ca0.5MnO3, and 4.39 K in Ho0.5Ca0.5MnO3. Probably the existence of the Schottky anomaly is related with the Kramers theorem. It states that an ion possessing an odd number of electrons, no matter how unsymmetrical the crystal field, must have a ground state that is at least doubly degenerate. This could lead to the thermal depopulation that produces the Schottky anomaly in the specific heat. Ions of Ce, Nd, Sm, Gd, Dy, Er and Yb all have an odd number of electrons and their respective oxides present a Schottky anomaly in the specific heat. However, the Kramers theorem does not exclude that ions with an even number of electrons might also have a doubly degenerate ground state. This is the case of the Ho ions, although it is noteworthy that the magnetic entropy associated with the Schottky transition in the compound with Ho had the lowest value (see above). Besides, our results indicate that the higher TS are found in the compounds with lower tolerance factors: Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3. The tolerance factor compares the Mn-O separation with the separation of the oxygen atom and A-site occupant in AMnO3. In this way the tolerance factor characterizes the angle of the Mn-O-Mn bond. Compounds with higher tolerance factors present wider conduction bands, which lead to a stronger ferromagnetic interactions and lower electrical resistivity.

## IV Conclusions

We have made a magnetic characterization of Nd0.5Sr0.5MnO3, Nd0.5Ca0.5MnO3 and Ho0.5Ca0.5MnO3 polycrystalline samples. It allowed us to identify the ferromagnetic, antiferromagnetic and charge ordering phases in each case. There are two important differences among these samples. The first one is that the tolerance factor, which influences the width of the conduction band, decreases from Nd0.5Sr0.5MnO3 to Ho0.5Ca0.5MnO3. The second one is that the rare earth ions have intrinsic magnetic moments: J=9/2 for Nd3+ and J=8 for Ho3+. These intrinsic magnetic moments experience a short range order at low temperatures. Moreover, the relation between the charge ordering temperature and the antiferromagnetic ordering temperature changes from one sample to the other. In the Nd0.5Sr0.5MnO3 sample they are approximately coincident (T$`{}_{CO}{}^{}`$TN), in the Nd0.5Ca0.5MnO3 sample the charge ordering temperature is much higher (T$`{}_{CO}{}^{}`$TN), and in the Ho0.5Ca0.5MnO3 sample there is no a long range antiferromagnetic transition.

We also reported, to our knowledge for the first time, specific heat measurements with applied magnetic fields between 0 and 9 T and temperatures between 2 and 300 K in all these three samples. Each curve was successfully fitted at high temperatures by an Einstein model with three optical phonon modes. Close to the charge ordering and ferromagnetic transition temperatures the specific heat curves showed peaks superposed to the characteristic response of the lattice oscillations. The entropy variation corresponding to the charge ordering transition was higher than the one corresponding to the ferromagnetic transition. Near 160 K in the two compounds with Nd3+ ions the specific heat curve showed an abrupt change in slope, which were correlated to the corresponding antiferromagnetic transition.

Absolute values of specific heat close to 2 K were about 100 times higher in our samples than in other charge ordering samples like Pr0.63Ca0.37MnO3. At low temperatures the specific heat curve, in all three studied samples and measured magnetic fields, showed a Schottky-like anomaly. In all cases an increase in the applied magnetic field moves the Schottky peak to higher temperatures. However, the height of the peak and the entropy, behave differently in Nd0.5Sr0.5MnO3 than in the other two samples. We could not successfully fit the curves by either assuming a singlet or a distribution of two-level-Schottky anomaly. Besides, we did not find a straightforward correlation between the maximum temperature of the Schottky anomaly and the magnetic moment of the rare earth ions. More experiments are clearly necessary to unambiguously identify the origin of the Schottky anomaly and its possible correlation with the charge ordered phase.

We thank the Brazilian science agencies FAPESP and CNPq for the financial support.

FIGURE CAPTIONS

Figure 1. Temperature dependence of the magnetization, with a 5 T applied magnetic field, in field cooling–warming condition for the three polycrystalline samples studied. Magnetization is given in Bohr magnetons per manganese ion. It is indicated TC, TN and TCO for in each sample. Note that charge ordering can occur in the presence of a ferromagnetic, antiferromagnetic or a paramagnetic phase.

Figure 2. Inverse of DC susceptibility measurements as a function of temperature in the Ho0.5Ca0.5MnO3 sample. The main graph shows the measurement made with a 5 T applied magnetic field, while the inset with a probe magnetic field of 0.1 mT. Lines are only guides to the eye.

Figure 3. Magnetization versus applied magnetic field at 2 K for the three samples studied. After a zero field cooling the magnetic field was increased from 0 to 5 T, decreased from 5 T to -5 T and increased again from -5 T to 5 T. The solid line represents a fitting by a Brillouin function (See the text for details). Magnetization is given in Bohr magnetons per manganese ion.

Figure 4. Specific heat measurements with zero magnetic field from 2 to 300 K for Nd0.5Sr0.5MnO3 (open squares), Nd0.5Ca0.5MnO3 (closed circles) and Ho0.5Ca0.5MnO3 (open triangles). Continuous lines represent the fitting of the phonon background to the Einstein model (see text). The Nd0.5Sr0.5MnO3 and Ho0.5Ca0.5MnO3 curves were displaced 20 J/mol K upside and downside, respectively.

Figure 5. Differences between the experimental specific heat data and the corresponding phonon background curves in the Nd0.5Sr0.5MnO3 sample (open squares), Nd0.5Ca0.5MnO3 sample (closed circles) and Ho0.5Ca0.5MnO3 sample (open triangles).

Figure 6. Specific heat measurements between 2 and 30 K in the samples of (a) Nd0.5Sr0.5MnO3, (b) Nd0.5Ca0.5MnO3 and (c) Ho0.5Ca0.5MnO3. Measurements were made in the presence of applied magnetic fields of 0 T (open squares), 5 T (closed circles), 7 T (open up triangles) and 9 T (closed down triangles). Continuous lines represent the fitting of the 15 to 30 K temperature interval data to the phonon contribution, as it is explained in the text. The graphs in the insets show the difference between the experimental values and the extrapolation of the phonon contribution to temperatures lower than 15 K.

Figure 7. Debye temperature (TD), entropy variation between 2 and 20 K ($`\mathrm{\Delta }S`$) and Schottky temperature (TS) as a function of the applied magnetic field in Nd0.5Sr0.5MnO3 (open squares), Nd0.5Ca0.5MnO3 (closed circles) and Ho0.5Ca0.5MnO3 (open triangles).


format: text

text:
# Renormalization for cubic frequency invariant tori in Hamiltonian systems with two degrees of freedom

## 1. Introduction

We consider the dynamics of a particle in a potential $`V`$ described by the following Hamiltonian :

$$H(p,x,t)=\frac{1}{2}p^2+V(x,t),$$
(1)
where $`V`$ is $`2\pi `$-periodic in $`x`$ and $`t`$. For $`V=0`$, the phase space is foliated with invariant tori. For $`V`$ sufficiently small and regular, the KAM theorem states that most of these invariant tori persist.
The purpose of this paper is to construct renormalization transformations in order to study the persistence of invariant tori with a given cubic frequency $`\omega `$ smaller than one, satisfying

$$\omega ^3=a\omega ^2+b\omega +1,$$

where $`(a,b)^2`$.

For Hamiltonian systems with two and three degrees of freedom, renormalization-group transformations in the framework of Ref.  have been successfully applied for the determination of the threshold of the break-up of invariant tori . They have also been used for the analysis of the properties of critical invariant tori (at the threshold of the break-up) .
For two degrees of freedom and for reduced quadratic irrational frequencies, renormalization transformations have been defined with a fixed frequency vector, i.e. at each iteration of the transformation, the frequency vector of the torus was unchanged. This feature was based on the fact that the continued fraction expansion of the frequency was periodic. For cubic (and more generally for non-quadratic irrational) frequencies, this is not the case, and thus ergodic renormalization was constructed from the entries of the continued fraction of the frequency. This transformation allows the accurate determination of critical couplings by comparison with Greene’s criterion  and Laskar’s frequency map analysis  (see Ref. ).
In order to define a renormalization with a fixed frequency vector, one has to add one degree of freedom to the system as explained below. For degenerate Hamiltonians with three degrees of freedom, a renormalization transformation has been defined for the spiral mean frequency vector. This transformation allows one to determine accurately the thresholds of break-up of tori by comparison with Laskar’s Frequency Map Analysis .
In this article, we propose to define and study numerically two renormalization transformations for invariant tori with cubic frequency : one based on the continued fraction of the frequency acting within a space of Hamiltonians with two degrees of freedom, and the other one with a fixed frequency vector acting within a space of degenerate Hamiltonians with three degrees of freedom.
We apply both renormalization methods to the numerical determination of the threshold of invariant tori of a given one-parameter family of Hamiltonians, Escande’s paradigm model  which is a forced pendulum model. We consider three cases : the spiral mean frequency, the Tribonacci mean frequency, and the $`\tau `$-mean frequency.

## 2. Renormalization from the continued fraction expansion

In order to investigate the torus with frequency $`\omega `$, we map the time-dependent Hamiltonian (1) into a time-independent Hamiltonian with two degrees of freedom, written in terms of actions $`𝑨=(A_1,A_2)^2`$ and angles $`𝝋=(\phi _1,\phi _2)𝕋^2`$ :

$$H(A_1,A_2,\phi _1,\phi _2)=\frac{1}{2}A_1^2A_2+V(\phi _1,\phi _2).$$

The two-dimensional invariant torus is now located around $`A_1=\omega `$. We shift the action $`A_1`$ in order to locate this torus around $`A_1=0`$. The Hamiltonian becomes :

$$H(𝑨,𝝋)=𝝎𝑨+\frac{1}{2}A_1^2+V(\phi _1,\phi _2),$$

where the frequency vector of the considered two dimensional torus is now $`𝝎=(\omega ,1)`$.

The renormalization we define relies upon the continued fraction expansion of $`\omega `$ :

$$\omega =\frac{1}{a_0+\frac{1}{a_1+\mathrm{}}}[a_0,a_1,\mathrm{}].$$

This transformation will act within a space of Hamiltonians $`H`$ of the form

$$H(𝑨,𝝋)=𝝎𝑨+V(𝛀𝑨,𝝋),$$
(2)
where $`𝛀=(1,\alpha )`$ is a vector not parallel to the frequency vector $`𝝎`$. We assume that Hamiltonian (2) satisfies a non-degeneracy condition : If we expand $`V`$ into

$$V(𝛀𝑨,𝝋)=\underset{\genfrac{}{}{0pt}{}{𝝂^2}{k0}}{}V_𝝂^{(k)}(𝛀𝑨)^ke^{i𝝂𝝋},$$

we assume that $`V_\mathrm{𝟎}^{(2)}`$ is non-zero.
The transformation contains essentially two parts : a rescaling and an elimination procedure. We follow the renormalization scheme defined in Ref. . Similar renormalization transformations have been defined in Refs.  for vector fields.

(1) Rescaling : The first part of the transformation is composed by a shift of the resonances, a rescaling of time and a rescaling of the actions. It acts on a Hamiltonian $`H`$ as $`H^{}=H𝒯`$ (see Ref.  for details) :

$$H^{}(𝑨,𝝋)=\lambda \omega ^1H(\frac{1}{\lambda }N_a𝑨,N_a^1𝝋),$$
(3)
with

$$\lambda =2\omega ^1(a+\alpha )^2V_\mathrm{𝟎}^{(2)},$$
(4)
and

$$N_a=\left(\begin{array}{cc}a& 1\\ 1& 0\end{array}\right),$$
(5)
and $`a`$ is the integer part of $`\omega ^1`$ (the first entry in the continued fraction expansion). This change of coordinates is a generalized (far from identity) canonical transformation and the rescaling $`\lambda `$ is chosen in order to ensure a normalization condition (the quadratic term in the actions has a mean value equal to 1/2). For $`H`$ given by Eq. (2), this expression becomes

$$H^{}(𝑨,𝝋)=𝝎^{}𝑨+\underset{𝝂,k}{}V_𝝂^{{}_{}{}^{}(k)}(𝛀^{}𝑨)^ke^{i𝝂𝝋},$$
(6)
where

$`𝝎^{}=(\omega ^{},1)\text{ with }\omega ^{}=\omega ^1a,`$ (7)
$`𝛀^{}=(1,\alpha )\text{ with }\alpha ^{}=1/(a+\alpha ),`$ (8)
$`V_𝝂^{{}_{}{}^{}(k)}=r_kV_{N𝝂}^{(k)}\text{ with }r_k=(1)^k2^{1k}\omega ^{k2}(a+\alpha )^{2k}\left(V_\mathrm{𝟎}^{(2)}\right)^{1k}.`$ (9)
We notice that the frequency of the torus is changed according to the Gauss map :

$$\omega \omega ^{}=\omega ^1\left[\omega ^1\right],$$
(10)
where $`\left[\omega ^1\right]`$ denotes the integer part of $`\omega ^1`$. Expressed in terms of the continued fraction expansion, it corresponds to a shift to the left of the entries

$$\omega =[a_0,a_1,a_2,\mathrm{}]\omega ^{}=[a_1,a_2,a_3,\mathrm{}].$$

Remark: After renormalization (one or more iterations) the frequency vector $`𝝎=(\omega ,1)`$ is changed into $`𝝎^{}=\mu M𝝎`$ with $`MGL_2()`$. In order to define a renormalization with a fixed frequency vector $`𝝎^{}=𝝎`$, the frequency $`\omega `$ has to be a quadratic irrational. For a cubic frequency (e.g., satisfying $`\omega ^3=a\omega ^2+b\omega +1`$), such renormalization cannot be defined with a fixed frequency vector. The idea is to add one more dimension and consider frequency vectors of the form $`(\omega ,1,\omega _1)`$. If $`\omega _1`$ is carefully chosen, there exists $`\mu `$ and $`MGL_3()`$ such that $`𝝎=\mu M𝝎`$, and then one can define a similar renormalization with the exception that now there is an additional degree of freedom corresponding to the additional frequency $`\omega _1`$.

(2) Elimination : The second step is a (connected to identity) canonical transformation $`𝒰`$ that eliminates the non-resonant modes of the perturbation in $`H^{}`$.
Following Ref. , we consider the set $`I^{}^2`$ of non-resonant modes as the set of integer vectors $`𝝂=(\nu _1,\nu _2)`$ such that $`|\nu _2|>|\nu _1|`$. Other choices of resonant modes can be done without affecting the convergence and speed of the algorithm (see Ref.  for a detailed discussion). A mode which is not an element of $`I^{}`$, will be called resonant. The canonical transformation $`𝒰`$ is such that $`H^{\prime \prime }=H^{}𝒰`$ does not have any non-resonant mode, i.e. it is defined by the following equation :

$$𝕀^{}(H^{}𝒰)=0,$$
(11)
where $`𝕀^{}`$ is the projection operator on the non-resonant modes; it acts on a Hamiltonian (2) as :

$$𝕀^{}H=\underset{\genfrac{}{}{0pt}{}{𝝂I^{}}{k0}}{}V_𝝂^{(k)}(𝛀𝑨)^ke^{i𝝂𝝋}.$$

We solve Eq. (11) by a Newton’s method following Ref. . We notice that in Refs.  a renormalization procedure acting within a space of quadratic Hamiltonians in the actions has been constructed. This renormalization is the one used in this paper for numerical purposes. This transformation gave the same accuracy than the renormalization defined for Hamiltonians in power series in the actions for the numerical computation of the parameters that characterize the break-up of invariant tori.

These rescaling and elimination procedures define the renormalization transformation acting on a space of Hamiltonians (2) with two degrees of freedom as $`H^{\prime \prime }=(H)=H𝒯𝒰`$. The main point is that the frequency $`\omega `$ is changed according to the Gauss map (10) each time we iterate the transformation.

## 3. Renormalization with a fixed frequency vector

In order to study a torus with frequency $`\omega `$, we consider another frequency $`\omega _1`$ such that $`(1,\omega ,\omega _1)`$ is a complete integral basis of the cubic field to which $`\omega `$ belongs to. We consider the invariant torus with frequency vector $`𝝎=(\omega ,1,\omega _1)`$. The Hamiltonian (1) is mapped into the time independent Hamiltonian with three degrees of freedom written in terms of actions $`𝑨=(A_1,A_2,A_3)^3`$ and angles $`𝝋=(\phi _1,\phi _2,\phi _3)𝕋^3`$ :

$$H(A_1,A_2,A_3,\phi _1,\phi _2,\phi _3)=\frac{1}{2}A_1^2A_2+V(\phi _1,\phi _2).$$

We shift the action $`A_1`$ by $`\omega `$ and we add $`\omega _1A_3`$ since $`A_3(t)`$ is constant ($`H`$ does not depend on $`\phi _3`$). The Hamiltonian becomes :

$$H(𝑨,𝝋)=𝝎𝑨+\frac{1}{2}A_1^2+V(\phi _1,\phi _2).$$

We assume that $`𝝎`$ is an eigenvector of a matrix $`N`$ with integer coefficients, determinant $`\pm 1`$, and a real eigenvalue $`\mu `$ of modulus smaller than one. For instance, if $`\omega _1=\omega ^2`$,

$$N=\left(\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ b& 1& a\end{array}\right),$$

and $`\mu =\omega `$.

This transformation will act within a space of Hamiltonians $`H`$ of the form

$$H(𝑨,𝝋)=𝝎𝑨+V(𝛀𝑨,𝝋),$$
(12)
where $`𝛀^3`$ is a vector of Euclidean norm one ($`𝛀=1`$), not parallel to the frequency vector $`𝝎`$. We assume that Hamiltonian (12) satisfies a non-degeneracy condition ($`V_\mathrm{𝟎}^{(2)}`$ is non-zero).
The renormalization is defined by the rescaling and elimination steps.

(1) Rescaling : The first part of the transformation is composed by a shift of the resonances, a rescaling of time and a rescaling of the actions. It acts on a Hamiltonian $`H`$ as $`H^{}=H𝒯`$ :

$$H^{}(𝑨,𝝋)=\lambda \mu ^1H(\frac{1}{\lambda }\stackrel{~}{N}𝑨,N^1𝝋),$$
(13)
where $`\stackrel{~}{N}`$ denotes the transposed matrix of $`N`$, and $`\lambda `$ is the rescaling coefficient :

$$\lambda =2\mu ^1N𝛀^2V_\mathrm{𝟎}^{(2)}.$$
(14)
For $`H`$ given by Eq. (12), this expression becomes

$$H^{}(𝑨,𝝋)=𝝎𝑨+\underset{𝝂,k}{}V_𝝂^{{}_{}{}^{}(k)}(𝛀^{}𝑨)^ke^{i𝝂𝝋},$$
(15)
where

$`𝛀^{}={\displaystyle \frac{N𝛀}{N𝛀}},`$ (16)
$`V_𝝂^{{}_{}{}^{}(k)}=r_kV_{\stackrel{~}{N}𝝂}^{(k)}\text{ with }r_k=2^{1k}\mu ^{k2}N𝛀^{2k}\left(V_\mathrm{𝟎}^{(2)}\right)^{1k}.`$ (17)
We notice that the frequency vector $`𝝎`$ of the torus is kept fixed by the transformation since from Eq. (13), we have $`\lambda \mu ^1𝝎\frac{1}{\lambda }\stackrel{~}{N}𝑨=𝝎𝑨`$.

(2) Elimination : The second step is a canonical transformation $`𝒰`$ that eliminates the non-resonant modes of the perturbation in $`H^{}`$.
Following Ref. , we consider the set $`I^{}`$ of non-resonant modes to be the set of integer vectors $`𝝂^3`$ such that $`|𝝎𝝂|\frac{1}{\sqrt{2}}𝝎𝝂`$. The canonical transformation $`𝒰`$ is such that $`H^{\prime \prime }=H^{}𝒰`$ does not have any non-resonant mode, i.e. it is defined by the following equation :

$$𝕀^{}(H^{}𝒰)=0.$$
(18)
We solve Eq. (18) by a Newton’s method following Ref. , and we apply the procedure acting on a space of Hamiltonians quadratic in the actions.

These two procedures define the renormalization transformation acting on a space of Hamiltonians (12) with three degrees of freedom as $`H^{\prime \prime }=H=H𝒯𝒰`$. The main point is that in this renormalization the frequency vector $`𝝎`$ is kept fixed at each iteration of the transformation.

## 4. Numerical results

We consider Escande’s paradigm model  which is a forced pendulum model :

$$H_\epsilon (p,x,t)=\frac{1}{2}p^2\epsilon (\mathrm{cos}x+\mathrm{cos}(xt)).$$
(19)
We map this Hamiltonian into a time-independent Hamiltonian with two degrees of freedom of the form (2), or into a time-independent Hamiltonian with three degrees of freedom of the form (12). We compute the critical coupling at which the torus is broken by the two renormalization methods $``$. We determine $`\epsilon _c`$ by iterating the renormalization transformation on $`H_\epsilon `$. If the coupling $`\epsilon `$ is smaller than a critical value $`\epsilon _c`$, the iterations converge to some integrable Hamiltonian $`H_0`$ (the coefficients $`V_𝝂`$ for $`𝝂\mathrm{𝟎}`$ vanish), and if $`\epsilon `$ is larger than $`\epsilon _c`$, the iterations diverge :

$`^nH_\epsilon H_0(𝑨)\text{ for }|\epsilon |<\epsilon _c,`$
$`^nH_\epsilon \mathrm{}\text{ for }|\epsilon |>\epsilon _c,`$
as $`n`$ tends to infinity. For numerical convenience, each renormalization $``$ is defined for quadratic Hamiltonians (all the iterations $`^nH_\epsilon `$ are quadratic in the actions), and at each step, we truncate the Fourier series by neglecting all the modes $`𝝂=(\nu _i)_{i=1,2\text{ or }i=1,2,3}`$ such that $`\mathrm{max}_i|\nu _i|>L`$. Other way of truncating the Fourier series can be chosen. The one chosen here is the simplest choice for the numerical implementation.
Since the renormalization for three degrees of freedom involves more Fourier coefficients (proportional to $`L^3`$, compared to $`L^2`$ for the renormalization for two degrees of freedom), we expect the resulting transformation to be much slower than the renormalization for two degrees of freedom. However, it turns out from numerical computations that in order to reach a given accuracy, both methods require approximately the same amount of computational time, i.e. the cut-off parameter $`L`$ for three degrees of freedom is chosen smaller than the one for two degrees of freedom. In what follows, the thresholds are computed with approximately the same amount of computational time.
We apply the above procedures for three different frequencies : the spiral mean, the Tribonacci mean, and the $`\tau `$-mean frequencies.

### 4.1. spiral mean torus

We consider an invariant torus with frequency $`\omega =\sigma ^1`$ where $`\sigma `$ is the spiral mean , i.e. it is the real solution of $`\sigma ^3=\sigma +1`$. It is approximately $`\sigma 1.32472`$. The frequency $`\sigma ^1`$ belongs to the complex cubic field with negative discriminant $`23`$. The break-up of KAM tori with spiral mean frequency of Hamiltonian systems with two degrees of freedom has been previously studied using scaling law analysis of Greene’s residues of nearby elliptic periodic orbits .

#### Renormalization for  $`𝝎=(\sigma ^1,1)`$

The first entries of the continued fraction expansion of $`\sigma ^1`$ are $`[1,3,12,1,1,3,2,3,2,4,2,141,80,2,5,\mathrm{}]`$. For a cut-off parameter $`L=25`$, the renormalization defined in Sec. 2 gives $`\epsilon _c0.0155`$.

#### Renormalization for  $`𝝎=(\sigma ^1,1,\sigma ^2)`$

The frequency vector $`𝝎=(\sigma ^1,1,\sigma ^2)`$ is an eigenvector of the matrix

$$N=\left(\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 0& 1& 1\end{array}\right).$$

For a cut-off parameter $`L=15`$, applying the renormalization defined in Sec. 3 with the above matrix $`N`$, we obtain $`\epsilon _c0.016`$. If we choose $`\omega _1=\sigma ^2+\sigma ^1+1`$ instead of $`\omega _1=\sigma ^2`$, the renormalization defined by the matrix

$$N=\left(\begin{array}{ccc}1& 1& 1\\ 1& 0& 1\\ 1& 0& 0\end{array}\right),$$

gives the same result ($`\epsilon _c0.016`$) as the renormalization defined for $`\omega _1=\sigma ^2`$.
Up to the numerical precision, the renormalizations for two and three degrees of freedom give the same result for the critical threshold $`\epsilon _c`$.

### 4.2. Tribonacci torus

We consider an invariant torus with frequency $`\omega =\varsigma ^1`$ where $`\varsigma `$ is called the Tribonacci number and is the real solution of $`\varsigma ^3=\varsigma ^2+\varsigma +1`$. Its value is approximately $`\varsigma 1.83928`$. The frequency $`\varsigma ^1`$ belongs to the complex cubic field with discriminant $`31`$.

#### Renormalization for $`𝛚=(\varsigma ^1,1)`$

The first entries of the continued fraction expansion of $`\varsigma ^1`$ are $`[1,1,5,4,2,305,1,8,2,1,4,6,17,5,1,\mathrm{}]`$. For a cut-off parameter $`L=25`$, we obtain $`\epsilon _c0.02186`$.

#### Renormalization for  $`𝝎=(\varsigma ^1,1,\varsigma ^2)`$

The frequency vector $`𝝎=(\varsigma ^1,1,\varsigma ^2)`$ is an eigenvector of the matrix

$$N=\left(\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 1& 1& 1\end{array}\right).$$

For a cut-off parameter $`L=10`$, the renormalization with the above matrix $`N`$ gives $`\epsilon _c0.0219`$. We notice that in this case, a very good agreement between both renormalizations is obtained even for small cut-off parameters $`L`$.

### 4.3. $`\tau `$-mean torus

We consider an invariant torus with frequency $`\omega =\tau ^1`$ where $`\tau `$ is the real solution of $`\tau ^3=\tau ^2+2\tau +1`$. Its value is $`\tau =2\mathrm{cos}(2\pi /7)1.24698`$. The frequency $`\tau ^1`$ belongs to the real cubic field with positive discriminant 49.

#### Renormalization for  $`𝝎=(\tau ^1,1)`$

The first entries of the continued fraction expansion of $`\tau ^1`$ are $`[1,4,20,2,3,1,6,10,5,2,2,1,2,2,1,\mathrm{}]`$. For a cut-off parameter $`L=25`$, we obtain $`\epsilon _c0.01247`$.

#### Renormalization for  $`𝝎=(\tau ^1,1,\tau ^2)`$

The frequency vector $`𝝎=(\tau ^1,1,\tau ^2)`$ is an eigenvector of the matrix

$$N=\left(\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 1& 1& 2\end{array}\right),$$

with eigenvalue $`\mu =\tau ^1`$. For any cut-off parameter, the critical value we obtain by the renormalization procedure for this matrix $`N`$ ($`\epsilon _c0.48`$) is very far from the one obtained by the previous renormalization. The reason is that the matrix $`N`$ has not the good properties. In order to define the procedure, the renormalization transformation has to contract all the integer vectors in $`I^+`$ (the complement of $`I^{}`$ in $`^3`$) into vectors in $`I^{}`$ . The spectrum of $`N`$ is composed by $`\tau ^1,(1+\tau ^1)^1,(\tau +1)`$. The shift of the resonances is a map acting on integer vectors as $`𝝂N^1𝝂`$. Since there are two eigenvalues of modulus larger than one for $`N^1`$ ($`\tau `$ and $`1\tau ^1`$), there are two unstable directions for the shift of the resonances.

We notice that for the negative discriminant case, this problem cannot occur. If there is one eigenvalue of modulus smaller than one, then the two other complex conjugated eigenvalues are of modulus greater than one since the determinant of $`N`$ is $`\pm 1`$.

#### Renormalization for  $`𝝎=(\tau ^1,1,\tau +1)`$

The frequency vector $`𝝎=(\tau ^1,1,\tau +1)`$ is an eigenvector of the matrix

$$N=\left(\begin{array}{ccc}2& 1& 1\\ 1& 1& 1\\ 0& 1& 0\end{array}\right),$$

with eigenvalue $`\mu =(1+\tau )^1`$. Here the matrix $`N`$ has the good property of having only one eigenvalue of modulus smaller than one (the two other eigenvalues are $`\tau `$ and $`1+\tau ^1`$). The renormalization defined in Sec. 3, for a cut-off parameter $`L=7`$, gives $`\epsilon _c0.013`$. For this frequency vector, an approximate renormalization has been set up for degenerate Hamiltonians with three degrees of freedom in Ref.  in order to investigate self-similar properties of three-dimensional critical tori. It is defined from the above matrix $`N`$ and structurally stable dynamics has been found for this approximate renormalization.

Remark : For the forced pendulum model (19), the Tribonacci torus is the most robust between the three invariant tori investigated in this paper, and the spiral mean torus is more robust than the $`\tau `$-mean torus. However, this feature depends on the perturbation in a way that is not understood. The most robust invariant tori are conjectured to be noble tori in general (see Ref.  for area-preserving maps and Ref.  for the stochastic ionization of Hydrogen driven by microwaves). Therefore, tori with cubic frequency should be less important for the large-scale diffusion of trajectories. However, the results presented in this article would be useful in order to investigate the break-up of invariant tori in Hamiltonian with three degrees of freedom and its comparison with the break-up in two degrees of freedom.

## 5. Conclusion

We have constructed and studied numerically two renormalization transformations in order to compute thresholds of break-up of invariant tori with a cubic frequency in Hamiltonian systems with two degrees of freedom. The main advantage of the procedure acting in a space of degenerate Hamiltonians with three degrees of freedom, is that it is independent of the continued fraction algorithm and it is defined for a fixed frequency vector. Both methods allow an accurate computation of critical couplings.

## Acknowledgments

We acknowledge useful discussions with G. Benfatto, K. Briggs, G. Gallavotti, H.R. Jauslin, H. Koch, J. Laskar, and R.S. MacKay. Financial support from the Carnot Foundation is also acknowledged.


text:
# 1 List of major experimental results

Corrections For Two Photon Decays of $`\chi _{c0}`$ and $`\chi _{c2}`$ and Color Octet Contributions
J. P. Ma and Q. Wang

Institute of Theoretical Physics , Academia Sinica, Beijing 100080, China

## Abstract

Using the fact that the c-quark inside a charmonium moves with a small velocity $`v`$ in the charmonium rest-frame, one can employ an expansion in $`v`$ to study decays of charmonia and results at the leading order for $`\chi _{c0,2}\gamma \gamma `$ exist in the literature. We study corrections at the next-to-leading order in the framework of nonrelativistic QCD(NRQCD) factorization. The study presented here is different than previous approaches where $`\chi _{c0,2}`$ is taken as a bound-state of a $`c\overline{c}`$ pair and a nonrelativistic wave-function is used for the pair. We find that the corrections are consist not only of relativistic corrections, but also of corrections from Fock state components of $`\chi _{c0,2}`$ in which the $`c\overline{c}`$ pair is in a color-octet state. For $`\chi _{c2}`$ there is also a contribution from a Fock state component in which the pair is in a F-wave state. We determine the factorization formula for decay widths in the form of NRQCD matrix elements representing nonperturbative effects related to $`\chi _{c0,2}`$, and calculate the perturbative coefficients at tree-level. Because the NRQCD matrix elements are unknown, a detailed prediction for the decay $`\chi _{c0,2}\gamma \gamma `$ can not be made, but the effect of these corrections can be determined at certain level. Estimations show that the effect is significant and can not be neglected.

PACS numbers: 13.25.Gv, 14.40.Gx, 12.38.Bx

1. Introduction

Because of their simple final states, decays of P-wave charmonium into two photons can be good channels to determine the value of $`\alpha _s`$ and for study of properties of charmonium. In recent years various experiment groups have published their results on $`\mathrm{\Gamma }_{\gamma \gamma }(\chi _{c0})`$ and $`\mathrm{\Gamma }_{\gamma \gamma }(\chi _{c2})`$ . Some major and newest results are listed in Table 1, where errors from different sources are combined. From these data it is easy to find that although there is a tendency of consistency in latest results of CLEO and E835, values of different groups have a large discrepancy from each other. On the theoretical side, because charm quark can be taken as heavy quark and it moves with a small velocity $`v`$ inside a charmonium in its rest-frame, one may describe the $`c`$\- or $`\overline{c}`$ quark inside charmonia with a nonrelativistic wave-functions by taking charmonia as a bound-state of a $`c\overline{c}`$ pair. In most previous calculations for decays of charmonium such a nonrelativistic wave-functions is employed for $`\chi _{c0,2}\gamma \gamma `$. By expanding the small velocity $`v`$, the decay width at the leading order of $`v`$ can be expressed with the derivatives of the wave functions at origin, where the leading order is at $`v^2`$. The corrections from the next-to-leading order of $`v`$ are also studied with the nonrelativistic wave-functions, in which the corrections are only relativistic corrections. This also implies that at the next-to-leading order the charmonia $`\chi _{c0,2}`$ are still taken as a bound-state of a $`c\overline{c}`$ pair and the pair has the same quantum number of $`\chi _{c0,2}`$.

From point of view of a quantum field theory like QCD, a hadron state is a bound state of quarks and gluons and has many components composed with different number of quarks and of gluons. For charmonia it can be expanded in $`v`$ and for $`\chi _{c0,2}`$ it looks like:

$$|\chi _{c0,2}=|(c\overline{c})_1+O(v)|(c\overline{c})_8g+higherorders,$$
(1)
where the index 1 or 8 denotes color singlet or color octet respectively. From this expansion the state $`|(c\overline{c})_1`$ is a dominant component, while other components have a suppressed probability to be observed in $`\chi _{c0,2}`$, the suppression parameter is $`v`$. Hence, in a systematic study of higher-order corrections one should also include contributions from these suppressed components, i.e., a systematic expansion in $`v`$ is needed. Recently, a novel approach based on NRQCD for productions and decays of quarkonia is proposed. In this approach a decay width can be systematically expanded in the velocity $`v`$. For example, the decay width of $`\chi _{c0}`$ can be written as:

$$\mathrm{\Gamma }(\chi _{c0}\gamma \gamma )=\underset{n}{}c_n\chi _{c0}|\widehat{O}_n|\chi _{c0},$$
(2)
where $`\widehat{O}_n`$’s are operators defined in NRQCD and $`c_n`$ can be calculated with perturbative theory. There is a rule of power counting of $`v`$ for the operators $`\widehat{O}_n`$. We will follow this approach to study the next-to-leading order corrections for the decay $`\chi _{c0,2}\gamma \gamma `$. We will determine the form of relevant operators in NRQCD and calculate the corresponding coefficients. It should be noted that the role of the suppressed component $`|(c\overline{c})_8g`$ has been studied before in inclusive processes. It has been shown that only by including the contribution from $`|(c\overline{c})_8g`$ the inclusive decay of $`\chi _{cJ}`$ can be consistently predicted with QCD. It has been also shown that such components with a color octet $`c\overline{c}`$ pair of $`\psi `$ are very important for explaining the $`\psi `$ surplus at the Tevatron. In this work we will show the effect of such components in exclusive processes. The effect of $`|(c\overline{c})_8g`$ has been studied in $`\chi _{cJ}\pi \pi `$, in which a light-cone wave-function for the component was introduced and the effect was not parameterized with matrix elements of NRQCD.

At the next-to-leading order of $`v`$, the corrections, as we will see, consist of relativistic corrections and also corrections from the state $`|(c\overline{c})_8g`$. Beside these, for $`\chi _{c2}`$ there is also a correction from the state $`|(c\overline{c})_1`$ where the $`c\overline{c}`$ pair is in a state with the orbit angular momentum $`l=3`$. These corrections take a factorized form as products of pertubative coefficients and matrix elements, in which the matrix elements are defined in NRQCD and they parameterize nonperturbative effects. We will determine the form of these matrix element and calculate these coefficients at leading order of $`\alpha _s`$. Our work is organized as the following: In Sect. 2 we calculate these coefficients and give our main results. In Sect. 3 we discuss the effect of these corrections and Sect.4 is our summary.

2. The correction at the next-to leading order

We consider $`\chi _{c0,2}`$ to two photon decay in its rest frame:

$$\chi _{c0,2}(P)\gamma (k_1)+\gamma (k_2)$$
(3)
where P, $`k_1`$, $`k_2`$ denote momentum of $`\chi _{c0,2}`$ and two photons. We decompose the decay width as

$$\mathrm{\Gamma }=\mathrm{\Gamma }_1+\mathrm{\Gamma }_8,$$
(4)
where $`\mathrm{\Gamma }_{1,8}`$ denotes the contribution from a color-singlet- and color-octet $`c\overline{c}`$ pair, respectively. The color-singlet contribution can be written as:

$`\mathrm{\Gamma }_1=`$ $`16`$ $`Q_{c}^{}{}_{}{}^{4}e^4{\displaystyle 𝑑\mathrm{\Gamma }\frac{d^4p_1}{(2\pi )^4}\frac{d^4p_{1}^{}{}_{}{}^{}}{(2\pi )^4}d^4xe^{i2(k_1+k_2)x}d^4ye^{ip_1(xy)}0|\overline{c_i}(x)c_j(y)|\chi _{c0,2}}`$ (7)
$`\times `$ $`{\displaystyle d^4y^{}e^{ip_{1}^{}{}_{}{}^{}(y^{}+x)}\chi _{c0,2}|\overline{c_j^{}}(y^{})c_i^{}(x)|0\times H_{ij,i^{}j^{}}}`$
$`d\mathrm{\Gamma }={\displaystyle \frac{d^4k_1}{(2\pi )^4}}{\displaystyle \frac{d^4k_2}{(2\pi )^4}}(2\pi )^2\delta _+(k_{1}^{}{}_{}{}^{2})\delta _+(k_{2}^{}{}_{}{}^{2})`$
$`H_{ij,i^{}j^{}}=[\gamma ^\mu {\displaystyle \frac{1}{p_1/k_2/m_c}}\gamma ^\nu +\gamma ^\nu {\displaystyle \frac{1}{p_1/k_1/m_c}}\gamma ^\mu ]_{ij}\times [\gamma _\nu {\displaystyle \frac{1}{p_1/^{}k_2/m_c}}\gamma _\mu ]_{j^{}i^{}},`$
where $`c(x)`$ stands for the Dirac field of $`c`$-quark and $`i,j,i^{}`$ and $`j^{}`$ stand for color- and spin indices. This contribution can be represented by the diagrams in Fig.1. The matrix element with Dirac fields of $`c`$-quark represents the nonperturbative effect related to the initial hadron. The quantity $`H`$ is calculated with perturbative theory. An expansion in $`v`$ can be performed by expanding the Dirac fields with fields of NRQCD. This can be done by using Foldy-Wouthuysen transformation:

$`c(x)`$ $`=`$ $`e^{im_ct}[1{\displaystyle \frac{i\gamma _jD_j}{2m_c}}+{\displaystyle \frac{1}{8m_{c}^{}{}_{}{}^{2}}}(i\gamma _jD_j)^2+{\displaystyle \frac{1}{4m_{c}^{}{}_{}{}^{2}}}\gamma _j[D_0,D_j]+{\displaystyle \frac{3}{16m_{c}^{}{}_{}{}^{3}}}(i\gamma _jD_j)^3]\psi (x)`$
$`+`$ $`e^{im_ct}[1{\displaystyle \frac{i\gamma _jD_j}{2m_c}}+{\displaystyle \frac{1}{8m_{c}^{}{}_{}{}^{2}}}(i\gamma _jD_j)^2+{\displaystyle \frac{1}{4m_{c}^{}{}_{}{}^{2}}}\gamma _j[D_0,D_j]+{\displaystyle \frac{3}{16m_{c}^{}{}_{}{}^{3}}}(i\gamma _jD_j)^3]\chi (x)+O(v^4)`$
(8)
where $`\psi (x)`$ and $`\chi (x)`$ are NRQCD fields and $`\psi (\chi ^{})`$ annihilates a $`c(\overline{c})`$-quark respectively. $`m_c`$ is the pole mass of $`c`$-quark. Using this transformation the matrix elements are then of the NRQCD fields $`\chi `$, $`\psi `$ and their derivatives. Operators $`𝐃`$ and $`_i`$ are of $`v^1`$ order and $`D_0`$ and $`_0`$ are of $`v^2`$ order in NRQCD according the power counting rule. The space-time dependence of matrix element are controlled with different scales in the different directions. In the time direction it is $`m_cv^2`$, while in the spacial direction it is $`m_cv`$. The expansion in $`v`$ is then completed by doing Taylor expansion of fields $`\chi `$ and $`\psi `$ around origin. Using the power counting rule in  we collect all local matrix elements up to order of $`v^4`$, which are allowed by rotation-, charge conjugation- and party symmetries for $`\chi _{c0,2}`$, and calculate the corresponding contributions. Doing the calculation in this way is similar to the matching method by taking a free $`c\overline{c}`$ pair, as proposed in , the difference is that we only collect the contributions from operators with correct quantum numbers in the problem considered here. The calculation is rather tedious but straightforward. We obtain:

$`\mathrm{\Gamma }_1(\chi _{c0}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{6\pi Q_{c}^{}{}_{}{}^{4}\alpha ^2}{m_{c}^{}{}_{}{}^{4}}}\chi _{c0}|𝒪_{EM}(^3P_0)|\chi _{c0}{\displaystyle \frac{14\pi Q_{c}^{}{}_{}{}^{4}\alpha ^2}{m_{c}^{}{}_{}{}^{6}}}\chi _{c0}|𝒢_{EM}(^3P_0)|\chi _{c0}+𝒪(v^6)`$ (9)
$`𝒪_{EM}(^3P_0)`$ $`=`$ $`{\displaystyle \frac{1}{3}}\psi ^{}({\displaystyle \frac{i}{2}}𝐃𝝈)\chi |00|\chi ^{}({\displaystyle \frac{i}{2}}𝐃𝝈)\psi `$ (10)
$`𝒢_{EM}(^3P_0)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[{\displaystyle \frac{1}{3}}\psi ^{}({\displaystyle \frac{i}{2}}𝐃)^2({\displaystyle \frac{i}{2}}𝐃𝝈)\chi |00|\chi ^{}({\displaystyle \frac{i}{2}}𝐃𝝈)\psi +h.c.],`$ (11)
where $`\sigma _i(i=1,2,3)`$ is the Pauli matrix. With the above results we find the correction to the color-singlet contribution is only the relativistic correction. Similarly, the decay width of $`\chi _{c2}`$ is:

$`\mathrm{\Gamma }_1(\chi _{c2}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{8\pi Q_{c}^{}{}_{}{}^{4}\alpha ^2}{5m_{c}^{}{}_{}{}^{4}}}\chi _{c2}|𝒪_{EM}(^3P_2)|\chi _{c2}{\displaystyle \frac{172\pi Q_{c}^{}{}_{}{}^{4}\alpha ^2}{105m_{c}^{}{}_{}{}^{6}}}\chi _{c2}|𝒢_{EM}(^3P_2)|\chi _{c2}`$ (12)
$``$ $`{\displaystyle \frac{284\pi Q_{c}^{}{}_{}{}^{4}\alpha ^2}{105m_{c}^{}{}_{}{}^{6}}}\chi _{c2}|𝒢_{EM}^{}(^3P_2)|\chi _{c2}+𝒪(v^6)`$
$`𝒪_{EM}(^3P_2)`$ $`=`$ $`\psi ^{}({\displaystyle \frac{i}{2}}D^{(i}\sigma ^{j)})\chi |00|\chi ^{}({\displaystyle \frac{i}{2}}D{}_{}{}^{(i}\sigma _{}^{j)})\psi `$ (13)
$`𝒢_{EM}(^3P_2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[\psi ^{}({\displaystyle \frac{i}{2}}𝐃)^2({\displaystyle \frac{i}{2}}D^{(i}\sigma ^{j)})\chi |00|\chi ^{}({\displaystyle \frac{i}{2}}D^{(i}\sigma ^{j)})\psi +h.c.]`$ (14)
$`𝒢_{EM}^{}(^3P_2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[\psi ^{}({\displaystyle \frac{i}{2}}D)^{(i}({\displaystyle \frac{i}{2}}D)^{j)}({\displaystyle \frac{i}{2}}𝐃𝝈)\chi |00|\chi ^{}({\displaystyle \frac{i}{2}}D^{(i}\sigma ^{j)})\psi +h.c.],`$ (15)
where the notation $`(ij)`$ means $`T^{(ij)}=(T^{ij}+T^{ji})/2\delta _{ij}T^{kk}/3`$, i.e., the tensor is symmetric and trace-less. Unlike the case with $`\chi _{c0}`$ we find not only a relativistic correction but also a contribution from a mixing with a $`F`$-wave $`c\overline{c}`$ pair. This is allowed by symmetries, the $`c\overline{c}`$ pair is in a spin triplet state, a spin triplet state with a orbit angular momentum $`l=3`$ can have a total angular momentum $`J=2,3,4`$. The state with $`J=2`$ is just the $`\chi _{c2}`$ state.

Now we consider color octet contributions from $`|(c\overline{c})_8g`$. At the order we consider the contributions can be represented by diagrams in Fig.2. The contributions before the expansion in $`v`$ can be written as:

$`\mathrm{\Delta }\mathrm{\Gamma }_{8a}`$ $`=`$ $`Q_{c}^{}{}_{}{}^{4}e^4{\displaystyle 𝑑\mathrm{\Gamma }d^4yd^4x\frac{d^4k}{(2\pi )^4}e^{ik(yx)+i(k_1+k_2)x}0|\overline{c_i}(x)c_j(y)|\chi _{c0,2}}`$ (16)
$`\times `$ $`{\displaystyle d^4y^{}d^4x^{}\frac{d^4k^{}}{(2\pi )^4}\frac{d^4q}{(2\pi )^4}e^{ik^{}y^{}+iqx^{}}\chi _{c0,2}|\overline{c_i^{}}(0)g_sG^\mu (x^{})c_j^{}(y^{})|0\times M_{}^{a}{}_{ij,i^{}j^{}}{}^{}}`$
$`\mathrm{\Delta }\mathrm{\Gamma }_{8b}`$ $`=`$ $`Q_{c}^{}{}_{}{}^{4}e^4{\displaystyle 𝑑\mathrm{\Gamma }d^4yd^4x\frac{d^4k}{(2\pi )^4}e^{ik(yx)+i(k_1+k_2)x}0|\overline{c_i}(x)c_j(y)|\chi _{c0,2}}`$ (17)
$`\times `$ $`{\displaystyle d^4y^{}d^4x^{}\frac{d^4k^{}}{(2\pi )^4}\frac{d^4q}{(2\pi )^4}e^{ik^{}y^{}iqx^{}}\chi _{c0,2}|\overline{c_i^{}}(y^{})g_sG^\mu (x^{})c_j^{}(0)|0\times M_{}^{b}{}_{ij,i^{}j^{}}{}^{}}`$
$`\mathrm{\Delta }\mathrm{\Gamma }_{8c}`$ $`=`$ $`Q_{c}^{}{}_{}{}^{4}e^4{\displaystyle 𝑑\mathrm{\Gamma }d^4yd^4x\frac{d^4k}{(2\pi )^4}e^{ik(yx)+i(k_1+k_2)x}0|\overline{c_i}(x)c_j(y)|\chi _{c0,2}}`$ (18)
$`\times `$ $`{\displaystyle d^4y^{}d^4x^{}\frac{d^4k^{}}{(2\pi )^4}\frac{d^4q}{(2\pi )^4}e^{ik^{}y^{}iqx^{}}\chi _{c0,2}|\overline{c_i^{}}(y^{})g_sG^\mu (x^{})c_j^{}(0)|0\times M_{}^{c}{}_{ij,i^{}j^{}}{}^{}}`$
with

$`M_{}^{a}{}_{ij,i^{}j^{}}{}^{}=`$ $`[\gamma ^\alpha {\displaystyle \frac{1}{k/k_2/m_c+iϵ}}\gamma ^\beta +\gamma ^\beta {\displaystyle \frac{1}{k/k_1/m_c+iϵ}}\gamma ^\alpha ]_{ij}`$ (19)
$`\times [\gamma _\beta {\displaystyle \frac{1}{k^{}/+q/+k_1/m_c+iϵ}}\gamma _\alpha {\displaystyle \frac{1}{k^{}/+q/m_c+iϵ}}\gamma _\mu ]_{i^{}j^{}}`$
$`M_{}^{b}{}_{ij,i^{}j^{}}{}^{}=`$ $`[\gamma ^\alpha {\displaystyle \frac{1}{k/k_2/m_c+iϵ}}\gamma ^\beta +\gamma ^\beta {\displaystyle \frac{1}{k/k_1/m_c+iϵ}}\gamma ^\alpha ]_{ij}`$ (20)
$`\times [\gamma _\mu {\displaystyle \frac{1}{k^{}/+q/m_c+iϵ}}\gamma _\beta {\displaystyle \frac{1}{k^{}/+q/k_2/m_c+iϵ}}\gamma _\alpha ]_{i^{}j^{}}`$
$`M_{}^{c}{}_{ij,i^{}j^{}}{}^{}=`$ $`[\gamma ^\alpha {\displaystyle \frac{1}{k/k_2/m_c+iϵ}}\gamma ^\beta +\gamma ^\beta {\displaystyle \frac{1}{k/k_1/m_c+iϵ}}\gamma ^\alpha ]_{ij}`$ (21)
$`\times [\gamma _\beta {\displaystyle \frac{1}{k^{}/k_2/m_c+iϵ}}\gamma _\mu {\displaystyle \frac{1}{k^{}/+q/k_2/m_c+iϵ}}\gamma _\alpha ]_{i^{}j^{}}`$
Similarly we need to perform the expansion in $`v`$ for matrix element like $`\chi _{c0}|\overline{c}_i(x)g_sG^\mu (y)c_j(0)|0`$. Again with the power counting rule, one can find that the leading term for the matrix element is related to the matrix element of the operator $`\psi ^{}g_s𝐄𝝈\chi `$. This matrix element is at order of $`v^3`$ and will lead to a contribution to $`\mathrm{\Gamma }_8`$ at order of $`v^4`$. To calculate the contribution we take the gauge: $`G^0=0`$, in this gauge the electric chromofield is:

$$𝐄=\frac{}{t}𝐆.$$
(22)
Performing the expansion for matrix elements in Eq.(19), results can be obtained in a straightforward way. We obtain the contribution at $`v^4`$:

$`\mathrm{\Gamma }_8(\chi _{c0}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{4\pi \alpha ^2Q_{c}^{}{}_{}{}^{4}}{m_{c}^{}{}_{}{}^{5}}}\chi _{c0}|_{EM}(^3P_0)|\chi _{c0}+𝒪(v^6)`$ (23)
$`_{EM}(^3P_0)`$ $`=`$ $`{\displaystyle \frac{1}{6}}[\psi ^{}𝝈g_s𝐄\chi |00|\chi ^{}𝝈𝐃\psi +h.c.]`$ (24)
$`\mathrm{\Gamma }_8(\chi _{c2}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{2\pi Q_{c}^{}{}_{}{}^{4}\alpha ^2}{m_{c}^{}{}_{}{}^{5}}}\chi _{c2}|_{EM}(^3P_2)|\chi _{c2}+𝒪(v^6)`$ (25)
$`_{EM}(^3P_2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[\psi ^{}\sigma ^{(i}g_sE^{j)}\chi |00|\chi ^{}\sigma ^{(i}D^{j)}\psi +h.c.].`$ (26)
Eq.(9), Eq.(12), Eq.(23) and Eq.(25) are our main results. From these results we can see that the corrections consist not only of relativistic correction which is represented by the matrix elements $`\chi _{c0}|𝒢_{EM}(^3P_0)|\chi _{c0}`$ and $`\chi _{c2}|𝒢_{EM}(^3P_2)|\chi _{c2}`$, but also of contributions from the color octet component. For $`\chi _{c2}`$ there is also a contribution from the component in which the $`c\overline{c}`$ pair is in color single state with the orbit angular momentum $`l=3`$.

3. Numerical impact of the correction

To make theoretical predictions for decay widthes one needs to know all matrix elements in the last section. Unfortunately, there is no information available for the matrix elements appearing in the correction. There are estimations for the matrix element at the leading order of $`v`$, but for consistency we need an estimation for them at the accuracy of $`𝒪(v^4)`$. All this prevents us from detailed predictions. However, numerical impact of the correction can be studied by noting that there are several relations among these operators. Because of spin symmetry of NRQCD we can have:

$`\chi _{c2}|𝒪_{EM}(^3P_2)|\chi _{c2}`$ $`=`$ $`\chi _{c0}|𝒪_{EM}(^3P_0)|\chi _{c0}(1+𝒪(v^2)),`$
$`\chi _{c2}|𝒢_{EM}(^3P_2)|\chi _{c2}`$ $`=`$ $`\chi _{c0}|𝒢_{EM}(^3P_0)|\chi _{c0}(1+𝒪(v^2)),`$
$`\chi _{c2}|_{EM}(^3P_2)|\chi _{c2}`$ $`=`$ $`\chi _{c0}|_{EM}(^3P_0)|\chi _{c0}(1+𝒪(v^2)).`$ (27)
As discussed before, the first relation is not useful here unless we know the correction at $`𝒪(v^2)`$. Beside these relations one can also use equation of motion to obtain another type of relations. Generally there are subtleties for using equation of motion for operators. However, it is shown that one can still use equation of motion for operators which are sandwiched between physical states. For operators relevant to $`J/\mathrm{\Psi }`$ one has used this to derive this type of relations.

The equation of motion of NRQCD reads:

$$(iD_t+\frac{𝐃^2}{2m})\psi =0,(iD_t\frac{𝐃^2}{2m})\chi =0.$$
(28)
We can use these equations to replace the operator $`𝐃^2`$ with $`D_t`$. For doing this we note that the matrix element appearing in Eq.(10) can be written as:

$$0|\chi ^{}(\frac{i}{2}𝐃𝝈)(\frac{i}{2}𝐃)^2\psi |\chi _{c0}=\frac{i}{2}0|\chi ^{}𝐃\sigma 𝐃^2\psi +(𝐃^2\chi )^{}𝐃\sigma \psi |\chi _{c0}\{1+𝒪(v^2)\},$$
(29)
where the corrections terms are due to interchange of the gauge covariant derivative $`𝐃`$, and they are at order of $`𝒪(v^2)`$ and can be neglected here. Using the equation of motion we obtain:

$`0|\chi ^{}𝐃\sigma 𝐃^2\psi +(𝐃^2\chi )^{}𝐃\sigma \psi |\chi _{c0}`$ $`=`$ $`2m_ci[0|\chi ^{}𝐃\sigma D_t\psi +(D_t\chi )^{}𝐃\sigma \psi |\chi _{c0}]`$ (30)
$`=`$ $`2m_cE_{\chi _{c0}}0|\chi ^{}𝐃\sigma \psi |\chi _{c0}2m_c0|\chi ^{}g𝐄\sigma \psi |\chi _{c0},`$
where $`E_{\chi _{c0}}`$ is the binding energy. Doing the same for $`\chi _{c2}`$ we obtain the relations:

$`\chi _{c_0}|𝒢(^3P_0)|\chi _{c0}`$ $`=`$ $`m_cE_{\chi _{c0}}\chi _{c_0}|𝒪(^3P_0)|\chi _{c0}+m_c\chi _{c_0}|(^3P_0)|\chi _{c0},`$
$`\chi _{c_2}|𝒢(^3P_2)|\chi _{c2}`$ $`=`$ $`m_cE_{\chi _{c2}}\chi _{c_2}|𝒪(^3P_2)|\chi _{c2}+m_c\chi _{c_2}|(^3P_2)|\chi _{c2},`$ (31)
and the binding energy is defined as:

$$E_{\chi _{c0}}=m_{\chi _{c0}}2m_c,E_{\chi _{c2}}=m_{\chi _{c2}}2m_c.$$
(32)
It should be noted that the difference $`E_{\chi _{c2}}E_{\chi _{c0}}`$ is at order higher than $`v^2`$, hence the above relations are in consistent with those in Eq.(27). To study the numerical impact we define the following parameters:

$`a_8`$ $`=`$ $`{\displaystyle \frac{\chi _{c0}|(^3P_0)|\chi _{c0}}{m_c\chi _{c0}|𝒪(^3P_0)|\chi _{c0}}}={\displaystyle \frac{\chi _{c2}|(^3P_2)|\chi _{c2}}{m_c\chi _{c2}|𝒪(^3P_2)|\chi _{c2}}}(1+𝒪(v^2)),`$
$`a_F`$ $`=`$ $`{\displaystyle \frac{\chi _{c2}|𝒢^{}(^3P_2)|\chi _{c2}}{m_{c}^{}{}_{}{}^{2}\chi _{c2}|𝒪(^3P_2)|\chi _{c2}}}.`$ (33)
These parameters are order of $`v^2`$. With these parameters and with the relations in Eq.(31) the decay widths can be written:

$`\mathrm{\Gamma }(\chi _{c0}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{6\pi Q_c^4\alpha ^2}{m_{c}^{}{}_{}{}^{4}}}\chi _{c0}|𝒪(^3P_0)|\chi _{c0}\{1+0.1787{\displaystyle \frac{\alpha _s(m_c)}{\pi }}2.333{\displaystyle \frac{E_{\chi _{c0}}}{m_c}}3a_8\}+𝒪(v^6),`$
$`\mathrm{\Gamma }(\chi _{c2}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{8\pi Q_c^4\alpha ^2}{5m_{c}^{}{}_{}{}^{4}}}\chi _{c2}|𝒪(^3P_2)|\chi _{c2}\{15.333{\displaystyle \frac{\alpha _s(m_c)}{\pi }}1.024{\displaystyle \frac{E_{\chi _{c2}}}{m_c}}2.2738a_81.6905a_F\}`$ (34)
$`+𝒪(v^6),`$
where we have added the one-loop correction to the contribution at leading order of $`v`$. If we take $`m_c=1.5`$GeV and $`\alpha _s(m_c)0.3`$, we obtain by using experimental values of hadron masses:

$`\mathrm{\Gamma }(\chi _{c0}\gamma \gamma )`$ $``$ $`{\displaystyle \frac{6\pi Q_c^4\alpha ^2}{m_{c}^{}{}_{}{}^{4}}}\chi _{c0}|𝒪(^3P_0)|\chi _{c0}\{1+0.0170.6450.3\}+𝒪(v^6),`$
$`\mathrm{\Gamma }(\chi _{c2}\gamma \gamma )`$ $``$ $`{\displaystyle \frac{8\pi Q_c^4\alpha ^2}{5m_{c}^{}{}_{}{}^{4}}}\chi _{c2}|𝒪(^3P_2)|\chi _{c2}\{10.5100.3800.2270.169\}+𝒪(v^6),`$ (35)
where we have taken $`a_8=a_F=0.1`$. In the above equations, the numbers in brackets correspond to the terms in brackets in Eq.(34), respectively. From these numbers, we see that the effect of binding energies is very significant, also the effects from color octet contributions and from mixing with the F-wave state are not small, although their effects can not be determined precisely. However, our results show that each correction at the next-to-leading order is significant. Unless among those corrections some cancellation happens, they can not be neglected.

4. Summary and conclusion

We have presented a systematic study of high order corrections in two photon decays of $`\chi _{c0,2}`$. Using NRQCD velocity power counting rule, implicit expansion of decay width in relative velocity $`v`$ is achieved. We calculated not only first order relativistic corrections but also corrections from suppressed component $`|(c\overline{c})_8g`$ in which $`c\overline{c}`$ pair is in a color octet state. The effect of color octet has been studied extensively in inclusive processes. In this work we have studied the effect in exclusive processes. Unlike most previous theoretical calculations in which nonrelativistic wave function is employed, our calculation is based on a field theoretic framework of NRQCD. The results are expressed by products of perturbatively calculable short distance coefficients and nonperturbative matrix elements which concern with initial bound states. In the case of $`\chi _{c2}`$, a correction from orbital angular momentum $`l=3`$ is also obtained in our calculation. In previous calculations only relativistic corrections can be estimated.

Unfortunately, the matrix elements at the next-to-leading order are not known at all, although some relations among them can be obtained by symmetries and by using equation of motion. This prevents us from numerical predictions for the decays. However, numerical impact of these corrections can be estimated at certain level and with this estimation the effect of the corrections is significant and can be not neglected. It would be interesting to study these matrix elements by nonperturbative methods, especially by lattice QCD. With a nonperturbative method one may determine how important the color octet component is.