text:
# ON THE HIDDEN NUCLEAR STARBURST IN ARP 220

## 1. INTRODUCTION

Ultraluminous infrared galaxies (ULIGs) are the most luminous galaxies in the Universe and their bolometric luminosities are comparable to those of quasars ($`L_{\mathrm{bol}}>10^{12}L_{}`$). The physical origin of the large luminosity of ULIGs is unknown since the very centers of these objects, which is where the bolometric luminosity comes from (Sanders & Mirabel 1996) is very heavily obscured by dust and direct observation is not possible. Either an enshrouded starburst or an active galactic nucleus (AGN), or both, could generate the observed luminoisites.

Arp 220 is the nearest and best-studied ULIG111We use a distance of 74 Mpc, assuming $`H_0=75\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`\mathrm{\Omega }_0=1`$ with the recession velocity to the Galactic Standard of Rest, $`V_{\mathrm{GSR}}=5,531\mathrm{k}\mathrm{m}\mathrm{s}^1`$ (de Vaucouleurs et al. 1991). Recently Genzel et al. (1998) performed a mid-infrared spectroscopy program for a sample of ULIGs and demonstrated that Arp 220 has an appreciable PAH feature and so is powered by a starburst, assuming that the very core where $`L_{\mathrm{bol}}`$ originates from is not optically thick at 15 $`\mu `$m. Based on detailed radiative transfer modelling, Rowan-Robinson (2000 and references therein) argues that galaxies like Arp 220 must be star-formation powered since they emit most of their luminosity longward of 50$`\mu `$m. Additional indirect evidence that Arp 220 is powered by star formation includes (i) the lack of an X-ray detection (Iwasawa 1999), even at hard wavelengths where Compton scattering is less important; (ii) VLBI radio observations, which show many radio supernovae but no strong point radio source which could correspond to the AGN (Smith et al. 1998); (iii) a LINER222LINER = Low Ionization Nuclear Emission-line Region (Heckman 1980).-like optical emission-line spectrum which can result from a starburst-driven superwind (Veilleux et al. 1995; Taniguchi et al. 1999; Lutz, Veilleux, & Genzel 1999). This is not, however, unambiguous evidence that there is no AGN present. In any case, these measurements do not probe the very central regions, which is where most of the power is emitted. Rather they probe regions with optical depths $`\tau <1`$ at optical wavelengths.

On the other hand OH megamaser emission certainly indicates that the power source is compact (Lonsdale et al. 1994, 1998). This would be consistent with an AGN power source, but again this measurement is not a unique discriminant. The nature of the power source in cold ULIGs like Arp 220 has taken on a particular significance in recent times with the detection of a population of these objects at high redshifts which contribute a very large fraction of the extragalactic background light at far-infrared and submillimeter wavelengths (Barger, Cowie & Sanders 1999).

We now attempt to provide additional constraints on the nature of the power source in Arp 220 by investigating the so-called ionizing photon deficiency problem. Prestwich et al. (1994) performed near-infrared (NIR) spectroscopy for a sample of starburst galaxies including Arp 220 and found that the ratio of the Br$`\gamma `$ luminosity to the infrared luminosity of Arp 220 is lower by a factor of ten than what is observed in lower-luminosity (but still very luminous) starburst galaxies, suggesting a deficit of ionizing photons. This could in part be due to absorption of the ionizing photons by dust, to an AGN contributing much of $`L_{\mathrm{bol}}`$ but few ionizing photons, or to the detailed shape of the stellar initial mass function (IMF). These physical process all operate on the global properties of Arp 220 in different ways and in this paper we use data from the literature together with starburst models and an SED fitting method to give constraints on the nature of the hidden power source.

## 2. MODELS

The nuclear super-star clusters seen in optical (Shaya et al. 1994) and near-infrared (Scoville et al. 1998) $`HST`$ images are not the source of the high bolometric luminosity of Arp 220 (Shioya, Taniguchi & Trentham 2000), which is completely hidden by dust at those wavelengths. It is this hidden power source that we attempt to model. We use the population synthesis model of Leitherer et al. (1999, Starburst99) to compute bolometric luminosities and ionizing photon fluxes which we compare to the observed values derived from far-infrared and radio data respectively (where extinction is less of a problem).

The following assumptions are made:

1) we assume a model with continuous star formation. Note that there is abundant cold gas available as fuel (Scoville et al. 1991, 1998; Downes & Solomon 1998);

2) we consider the three stellar IMFs given in Table 1, where we list the IMF slope $`\alpha `$, lower-mass cutoff $`M_\mathrm{l}`$, and upper-mass cutoff $`M_\mathrm{u}`$. The values of $`M_\mathrm{l}`$ are higher than what is observed in the solar neighborhood, but seems to be more appropriate for violent starburst galaxies (e.g., Goldader et al. 1997);

3) we assume solar metallicity, appropriate for the central regions of galaxies, which is where most stars generated in this kind of dense starburst will end up (e.g. Faber 1973);

4) we assume that a hypothesized AGN coexists with the starburst and contributes a fraction $`f_{\mathrm{AGN}}`$ of the bolometric luminosity. Radio emission from this AGN is assumed to contribute negligibly to the ionizing photon flux inferred from radio recombination lines. This is reasonable since we expect an electron column density larger than $`10^{23}\mathrm{cm}^2`$ given the existence of very dense gas traced by HCN (Solomon, Downes, & Radford 1992) and CS (Solomon, Radford, & Downes 1990). Radio emission from an AGN is then suppressed because of free-free absorption.

## 3. RESULTS

### 3.1. Initial Mass Function, Age, and AGN fraction

The total bolometric luminosity of Arp 220 is $`L_{\mathrm{bol}}1.5\times 10^{12}L_{}`$ (Sanders et al. 1988). When considering the hidden power source, we need to subtract the part of the bolometric luminosity which is observed at optical and near-infrared wavelengths in the center of the galaxy. The $`K`$-magnitude of the central $`2.5^{\prime \prime }`$ diameter region of Arp 220 is 11.80 mag (Carico et al. 1990). Even with $`A_V=6.4`$ mag of extinction (Shier et al. 1996), this is only $`1\times 10^{11}L_{}`$ (assuming a bolometric correction of 0.5; Leitherer et al. 1999). This value is less than 10 % of the total bolometric luminosity of Arp 220 and we neglect it.

The ionizing photon production rate of Arp 220 is evaluated from the observation of radio recombination lines like H92$`\alpha `$: $`Q=1.3\times 10^{55}\mathrm{photons}\mathrm{s}^1`$ (Anantharamaiah et al. 2000). This is twice as large as the rate derived by Zhao et al. (1996) but Anantharamaiah et al. (2000) have wider velocity coverage; note that it is also larger than the upper limit derived by Scoville et al. (1991) from the free–free radio continuum flux density at 110 GHz, $`Q<7.5\times 10^{54}\mathrm{photons}\mathrm{s}^1`$. Again, we have to estimate the ionizing photon rate from the regions we do see in $`K`$-band and subtract it from the above values to derive the production rate for the hidden power source. The reddening-corrected ($`A_V`$ 10 mag; Larkin et al. 1995) Br$`\gamma `$ flux ($`5.9\times 10^{18}\mathrm{W}\mathrm{m}^2`$; Goldader et al. 1995) gives an ionizing photon rate, $`7\times 10^{53}\mathrm{photons}\mathrm{s}^1`$ for these visible starburst regions. Therefore the production rate for the hidden power source is: $`Q(\mathrm{hidden})1.23\times 10^{55}\mathrm{photons}\mathrm{s}^1`$.

Following Scoville et al. (1997), we adopt the stellar mass of the nucleus of Arp 220 as $`2.5\times 10^9M_{}`$. We therefore set the constraints on the stellar populations of the hidden starbursts as follows: $`8.1\times 10^{44}L_{}(\mathrm{photons}\mathrm{s}^1)^1<L_{\mathrm{bol}}/Q(\mathrm{hidden})<1.2\times 10^{43}L_{}(\mathrm{photons}\mathrm{s}^1)^1`$ and $`400L_{}/M_{}<L_{\mathrm{bol}}/M_{}(\mathrm{hidden})<600L_{}/M_{}`$. An AGN fraction $`f_{\mathrm{AGN}}>0`$ results in both of these parameters being lowered by a factor $`(1f_{\mathrm{AGN}})`$.

Now we compare the above numbers to the parameters derived for the three models described in Section 2. Figure 1 shows the evolutionary locus on a $`L_{\mathrm{bol}}/Q`$ \- $`L_{\mathrm{bol}}/M_{}`$ plane for each model. The general form of the curves is understood as follows. Low-mass stars have smaller $`L_{\mathrm{bol}}/M_{}`$ but larger $`L_{\mathrm{bol}}/Q`$ than high-mass stars. Since the fractional mass of low-mass stars increases with increasing age, as time increases, $`L_{\mathrm{bol}}/M_{}`$ becomes smaller while $`L_{\mathrm{bol}}/Q`$ becomes larger for all the models.

Let us consider the results for the three IMFs in turn.

#### 3.1.1 IMF1

For IMF1, the model overpredicts the bolometric luminosity given the observed ionizing photon flux (box A). This is the classical ionizing photon deficiency problem (e.g., Goldader et al. 1995). This could be solved in a number of ways: (i) reducing the number of massive stars which generate more bolometric luminosity per unit ionizing photon that lower mass stars. This could be achieved by either steepening the slope of the IMF (as in IMF2) or lowering the upper mass cutoff (as in IMF3). These adjustments have different effects on the properties studied here and are discussed in the next two sections. (ii) if some ionizing photons are absorbed by dust shortly after their production, the ionizing photon rate evaluated from radio recombination lines will be too low. The figure shows that IMF1 can then be consistent with observation if the hidden power source is a pure starburst with $`2\times 10^7<\mathrm{age}/\mathrm{yr}<3\times 10^7`$ and between 0.3 and 0.6 of the ionizing photons that are generated are absorbed by dust (iii) if an AGN contributes some of the power, the ionizing photon deficiency problem can be reduced. If $`f_{\mathrm{AGN}}=0.5`$, then IMF1 is marginally consistent with observation if $`4\times 10^7<\mathrm{age}/\mathrm{yr}<7\times 10^7`$ and no more than 30% of the ionizing photons are absorbed by dust. Increasing $`f_{\mathrm{AGN}}`$ much beyond 0.5 is not permitted since there are then too many ionizing photons given the bolometric luminosity and the ionizing photon deficiency problem now becomes an ionizing photon excess problem. For example, region C on Figure 1 is excluded if IMF1 is valid.

#### 3.1.2 IMF2

This steep IMF can only be consistent with observation if the age is very young ($`<5\times 10^6`$ yr) and $`f_{\mathrm{AGN}}`$ is finely tuned to have a value of about 0.5 and if the fraction of ionizing photons absorbed by dust is finely tuned to have a value of about 0.2 (but not 0). We therefore regard it as unlikely. This IMF can never by itself generate a high enough value of $`L_{\mathrm{bol}}/M_{}`$ to be consistent with observation. Increasing the AGN fraction helps a little but increasing it by too much creates an ionizing photon excess problem as described in the previous section.

#### 3.1.3 IMF3

If we adopt the IMF with a smaller upper-mass limit ($`M_u=30M_{}`$: IMF3), the evolutionary loci shift to the upper right on Figure 1. In the absence of either absorption of ionizing photons by dust or an AGN, the observational constraints are satisfied if $`6.3\times 10^6<\mathrm{age}/\mathrm{yr}<1.1\times 10^7`$. An AGN fraction greater than zero or an appreciable fraction of ionizing photons being absorbed by dust would generate an ionizing photon excess problem.

### 3.2. Current star formation rate in Arp 220

Since we adopt the constant star formation model, the star formation rate (SFR) is coupled uniquely to the bolometric luminosity for each model. Assuming the SFR of $`1M_{}\mathrm{yr}^1`$, the ionizing photon production rate is about $`Q=4.6\times 10^{52}\mathrm{photons}\mathrm{s}^1`$ at an age of $`1\times 10^7`$ yr for Model with IMF3 (i.e. where it intersects with Box A on Fig. 1). This corresponds to a total SFR of the hidden starburst of $`267M_{}`$ yr-1. In Arp 220, there are two concentrations of thermal dust emission with which the hidden starbursts are associated. Downes & Solomon (1998) evaluated the luminosity ratio between the western and the eastern nucleus as 1.5. We assume that the ratio of SFR between them is the same as this luminosity ratio. The SFR of the western nucleus is then $`160M_{}`$ yr-1 and that of the eastern nucleus is $`107M_{}`$ yr-1.

The results in this and the previous section depend intricately on our assumption of a star-formation rate that does not vary strongly with time (Point 1 in Section 2). An alternative way of explaining the deficiency is to have a star-formation rate that was extremely high at some time in the past but is low now. The highest mass stars, which produce most of the Lyman continuum, would since have died (these have very short lifetimes) and therefore the observed value of $`Q`$ is low. If the starburst was sufficiently powerful to generate many intermediate-mass stars, these could produce the observed $`L_{\mathrm{bol}}`$. Very large numbers of such stars would be required, since each star individually contributes only modestly to $`L_{\mathrm{bol}}`$. Such a scenario may then need to be finely tuned in terms of both the time-variation of the star-formation rate and the IMF so as to avoid overproducing the total mass in stars (this must be less than the dynamical mass), but this is still an important caveat to our results.

### 3.3. Visual Extinction

Once we have chosen a particular model, we can compute the visual extinction necessary to hide the hidden starbursts in the $`K`$-band $`HST`$ image of Scoville et al. (1998). We find that $`A_V`$ must be larger than 30 mag for the western nucleus and 40 mag for the eastern nucleus for either IMF1 where more than 0.4 of the ionizing photons are absorbed by dust or for IMF3 where none of the ionizing photons are absorbed by dust (see the discussions in Sections 3.1.1 and 3.1.3). This result is consistent with the recent MIR observations; i.e., $`A_V`$ of 45 mag (Genzel et al. 1998).

## 4. Conclusion

In all the cases that we consider, the hidden power source due to star formation at the 50% level or more in order to avoid an ionizing photon excess problem, and this starburst must be young ($`<7\times 10^7`$ yr). We derive a current star formation rate of $`267M_{}`$ yr-1, and an extinction $`A_V>30`$ mag to this hidden starburst.

YS thanks the Japan Society for Promotion of Science (JSPS) Research Fellowships for Young Scientist. NT thanks the PPARC for financial support. YT would like to thank R. -P. Kudritzki, Bob McLaren and Dave Sanders at Institute for Astronomy, University of Hawaii for their warm hospitality. This work was financially supported in part by Grant-in-Aids for the Scientific Research (Nos. 07044054, 10044052, and 10304013) of the Japanese Ministry of Education, Science, Sports and Culture.


format: text

text:
# Study of Charmonia near the deconfining transition on an anisotropic lattice with O(a) improved quark action

## 1 Introduction

It is generally believed that the Quantum Choromodynamics (QCD) exhibits a phase transition at some temperature $`T_c`$, and quarks and gluons confined in the low temperature phase are liberated to form the “quark gluon plasma”. At the beginning of 2000, CERN reported that the QGP state had been created in the heavy ion collision experiment. In these experiments, the $`J/\psi `$ suppression is regarded as the key signal of QGP formation. Since c-quark is heavy ( $`m_c=1.151.35`$ GeV) $`c\overline{c}`$ pair are hard to be generated except by the primary collisions of nucleons in the high energy heavy ion, i.e. the $`J/\psi `$ is not created with the thermal effects after the QGP formation. Therefor it is expected that the effects of deconfining clearly appear in this signal. Further investigation will be performed in RHIC project at BNL.

On the theoretical side, in spite of the various approaches, we are still far from the definite understanding of hadron properties near the transition and the fate of the hadronic states in the plasma phase. Since the phase transition changes the relevant degree of freedom of the system, the model approach which a priori assumes dynamical degrees of the system is difficult to treat physics near the phase transition. We investigate these problems using lattice QCD, which enables us to incorporate the nonperturbative effect of QCD from the first principle.

For a long time, in the lattice QCD simulations, hadron masses at finite temperature have been argued with the spatial correlation (screening mass). On the other hand, the study of the temporal correlation, which is related to the pole mass, has been started rather recently.

The previous work studied hadron properties near the $`T_c`$ based on light quarks. They caught a sign of chiral symmetry restoration above $`T_c`$, and a change of correlators in temporal and spatial direction near the $`T_c`$. From the discussion of wave function, however, they find the same strong spatial correlation at $`T>T_c`$ as that of below $`T_c`$, in other words hadronic mode at a long range survives in the deconfinement phase. At present, there is no well-established way of lattice simulations to attack the spectroscopy at $`T>0`$. However, our recent work seems one of the best approaches to the problem. It is interesting to apply this analysis to the heavy quarks. In this paper we focus on the heavy quarkonium state, which plays an important role as a signal of the quark gluon plasma formation.

Our goal of the investigation is the prediction about $`J/\psi `$ suppression or mass shift of charmonium as signals of QGP formation and the understanding of hadron properties and nature of QGP phase near the transition.

For the study of charmonium physics at $`T>0`$ on a lattice, there are several problems. Then we classify these problems into two category, and discuss them individually.

Precise calculation of temporal correlator of charmonium at $`T>0`$

Extraction of physical properties of charmonium from the correlators

Firstly we consider the former one. In the lattice QCD simulation at $`T>0`$, we set a temporal lattice extent to $`1/T`$. At high temperature, one needs the large lattice cutoff to work with the sufficient degrees of freedom in the temporal direction. In order to obtain the detailed information of temporal meson correlators at $`T>0`$, a high resolution in temporal direction is needed. The large lattice cutoff is also necessary to study a correlator of meson with the heavy quarks because of its rapid decreasing behavior. If one tries to overcome these difficulties with straightforward way, the tremendous large computational power is necessary. In order to get the sufficiently fine resolution with limited computer resources, we use the anisotropic lattice, which has a finer temporal lattice spacing $`a_\tau `$ than the spatial one $`a_\sigma `$.

In this work we adopt the same strategy as the previous work which was tractable to analyze light hadrons at $`T>0`$. This strategy is as follows. Firstly the mesonic operator are defined, then we observe its correlator and wave function (in the Coulomb gauge) at $`T=0`$. Next we investigate how they are affected by the temperature. In order to investigate the temperature effects for the state of interest, for example a ground and excited state of charmonium, we have to make the good mesonic operator which has the large overlap with the state. Because we should compare correlators at $`T0`$ within shorter temporal lattice extent. The wave function gives hints for existence of mesonic state at $`T>0`$. Especially we are interested in that of deconfinement phase. In the case of charmonium this wave function is an important quantity concerning the $`J/\psi `$ suppression.

We prepare two sets of gauge configuration whose lattice spacing are different. We control the temperature by changing temporal lattice extent $`N_t`$. Then above investigation is performed on these configurations.

This paper is organized as follows. In the next section we define the quark action on anisotropic lattice and discuss the dispersion relation of free quark and calibration of quark field. Sect. 3 is the preparation for study of charmonium correlator. Here various parameters of gauge configurations are determined and the calibration are performed. In Sect. 4 we report the charmonium spectroscopy and construction of optimized operator using variational analysis at $`T=0`$. Sect. 5 describes correlators at $`T>0`$ and the measurement of the wave function and compare these with results at $`T=0`$. The last section is the conclusion and discussion.

## 2 Quark action on the anisotropic lattice

### 2.1 Quark action

To treat the quark field on the lattice, we adopt the $`O(a)`$ improved Wilson quark formulation. To construct the quark action on the anisotropic lattice, we follow El-Khadra, Kronfeld and Mackenzie, for the following advantages. They expand the lattice Hamiltonian in the power of $`a`$, and determine the coefficient of each operator by matching the lattice Hamiltonian with the the continuum one except the redundant operators. The resultant action takes the same form as the clover quark action in the limit of $`m0`$. On the other hand, in the heavy quark mass region ($`m_q\mathrm{\Lambda }_{QCD}`$), the effective-theoretical treatment of the quark action enables us to use it for such a quark on a lattice of moderate cutoff. Since the full quark mass dependence is incorporated, the same form also covers the intermediate quark mass region, and then the small and the large mass regions are smoothly connected. Although our main target in this paper is the charm quark, it would be useful to take the other mass region into account for the future applications. In addition to these advantages, their argument is naturally in accord with the anisotropic lattice. They introduced the different hopping parameters for the spatial and the temporal directions, and proposed to tune them so that the rest and the kinetic mass take the same value. Such a treatment is inevitably required on the anisotropic lattice to assure that the anisotropy of the quark and the gauge fields coincide, especially if one employ the dispersion relation for the definition of the quark field anisotropy.

On the anisotropic lattice, the quark action takes almost the same form as in the Ref.  111 The notation in this paper is slightly different from the Ref.  . :

$$S_F=\underset{x,y}{}\overline{\psi }(x)K(x,y)\psi (y),$$
(1)
$`K(x,y)`$ $`=`$ $`\delta _{x,y}\kappa _\tau \left\{(1\gamma _4)U_4(x)\delta _{x+\widehat{4},y}+(1+\gamma _4)U_4^{}(x\widehat{4})\delta _{x\widehat{4},y}\right\}`$ (2)
$`\kappa _\sigma {\displaystyle \underset{i}{}}\left\{(r\gamma _i)U_i(x)\delta _{x+\widehat{i},y}+(r+\gamma _i)U_i^{}(x\widehat{i})\delta _{x\widehat{i},y}\right\}`$
$`\kappa _\sigma c_E{\displaystyle i\sigma _{4i}F_{4i}(x)\delta _{x,y}}r\kappa _\sigma c_B{\displaystyle \underset{ij}{}}{\displaystyle \frac{1}{2}}\sigma _{ij}F_{ij}(x)\delta _{x,y}.`$
The spatial and the temporal hopping parameters, $`\kappa _\sigma `$ and $`\kappa _\tau `$ respectively, are related to the bare quark mass $`m_0`$ and the bare anisotropy parameter $`\gamma _F`$ as follows.

$`\kappa _\sigma `$ $`=`$ $`{\displaystyle \frac{1}{2(m_0+\gamma _F+3r)}},`$ (3)
$`\kappa _\tau `$ $`=`$ $`\gamma _F\kappa _\sigma .`$ (4)
where the $`m_0`$ is in the spatial lattice unit. In the free quark case, the bare anisotropy $`\gamma _F`$ is taken to be the same value as the cutoff anisotropy $`\xi =a_\sigma /a_\tau `$. In practical simulation, the anisotropy parameter receives the quantum effect and should be tuned to give the same renormalized anisotropy for the fermion and the gauge fields. This “calibration” will be described later.

There have been used two choices of the value of the Wilson parameter $`r`$ for the anisotropic $`O(a)`$ improved quark action. In this work, we adopt the choice $`r=1/\xi `$. In this case, the temporal and the spatial directions are treated in the equal manner in the physical unit. As the result, the tree level dispersion relation holds the axis-interchange symmetry in the lowest order of $`\stackrel{}{p}^2`$. On the other hand, this choice decreases the masses of doublers which are introduced by the Wilson term to eliminate the unwanted poles at the edges of the Brillouin zone. The dispersion relation is examined in the later part of this section. Alternative choice, $`r=1`$, are adopted in the Ref. . In this case, the contribution of doublers would not cause any problem, in the cost of manifest axis-interchange symmetry. Since we aim to develop the form applicable to the whole quark mass region, $`r=1/\xi ^1`$ seems preferable especially in the light quark mass region.

Here it is useful to define $`\kappa `$ so that which has the same relation with the bare quark mass $`m_0`$ as the isotropic case:

$`{\displaystyle \frac{1}{\kappa }}={\displaystyle \frac{1}{\kappa _\sigma }}2(\gamma _F+3r4)(=2(m_0+4)).`$ (5)
For the light quark systems, the extrapolation to the chiral limit would be performed in $`1/\kappa `$. The coefficients of the clover terms, $`c_E`$ and $`c_B`$, depend on the Wilson parameter $`r`$. In our choice $`r=\xi ^1`$, $`c_E`$ and $`c_B`$ are unity at the tree level.

We apply the mean-field improvement proposed in the Ref. . On the anisotropic lattice, the mean-filed values of the spatial link variable $`u_\sigma `$ and the temporal one $`u_\tau `$ are different from each other. The improvement is achieved by rescaling the link variable as $`U_i(x)U_i(x)/u_\sigma `$ and $`U_4(x)U_4(x)/u_\tau `$. This replacement leads the following values for the coefficient of the clover terms.

$$c_E=\frac{1}{u_\sigma u_\tau ^2},c_B=\frac{1}{u_\sigma ^3}.$$
(6)
The determination of mean-field values of the link variable $`u_\sigma `$ and $`u_\tau `$ are described in the next section.

In this paper, our target mass region is around the charm quark mass. The temporal cutoffs in this work are $`4.5`$ and $`6.4`$ GeV, and well above the charm quark mass. For these quark mass and $`a_\tau ^1`$, the effective-theoretical treatment would not be necessary to be applied. Such consideration will be called for the calculations containing the $`b`$-quark on the same size of lattice. In the effective-theoretical treatment, the ratio of the spatial and the temporal hopping parameter is tuned so that they give correct dispersion relation of the nonrelativistic quark. On the anisotropic lattice, the calibration automatically incorporates this condition if one use the nonrelativistic dispersion relation as the anisotropy condition.

### 2.2 Dispersion relation of free quark

Now we consider how the dispersion relation of the free quark is changed by the introduction of anisotropy. Observing the action (2), one notices that the larger anisotropy $`\xi `$ causes the smaller spatial Wilson term. Then the question is how the contribution of the doubler eliminated by the Wilson term becomes significant. The action (2) leads the free quark propagator,

$$S(p)=\frac{1}{i\gamma _4\mathrm{sin}p_4+i\zeta \stackrel{}{\gamma }\stackrel{}{S}+m_t+(1\mathrm{cos}p_4)+\frac{1}{2}r\zeta \stackrel{}{\widehat{p}}^2},$$
(7)
where $`\zeta =\xi ^1`$, $`m_t=m_0\zeta `$, $`S_i=\mathrm{sin}p_i`$ and $`\widehat{p_i}=2\mathrm{sin}(p_i/2)`$. Then the dispersion relation of the free quark is

$$\mathrm{cosh}E(\stackrel{}{p})=1+\frac{\zeta ^2\stackrel{}{S}^2+(m_t+\frac{1}{2}r\zeta \stackrel{}{\widehat{p}}^2)^2}{2(1+m_t+\frac{1}{2}r\zeta \stackrel{}{\widehat{p}}^2)}.$$
(8)
Neglecting the higher order terms in $`\stackrel{}{p}`$ and in $`a`$, the relativistic dispersion relation $`E^2=m^2+\zeta ^2\stackrel{}{p}^2`$ holds for the small quark mass. ($`m^2=m_t^2(1m_t)`$, and $`m_t`$ is the bare parameter.)

Fig. 1 shows the dispersion relation (8) at $`\xi =5.3`$ and $`4.0`$ for various values of $`m_t`$. Now let us consider the practical cases that $`a_\tau ^1=4.5`$ GeV for $`\xi =5.3`$ (Set-I) and $`a_\tau ^1=6.4`$ GeV for $`\xi =4.0`$ (Set-II). These values are obtained in our numerical simulations, and described in the next section. In the heavy quarkonium, the typical energy and momentum exchanged inside the meson are in the order of $`mv^2`$ and $`mv`$ respectively. For the charmonium, $`v^20.3`$, then typical scale of the kinetic energy is around $`500`$ MeV. It is noted that $`m_t0.3`$ for Set-I and $`m_t0.2`$ for Set-II correspond to the charm quark mass. Let us consider two quarks inside meson with opposite momenta $`p=\pm \pi /a`$. Then $`2(E(p_z=a/\pi )E(0))0.5`$ GeV and $`1`$ GeV for Set-I and Set-II respectively. Although Set-I lattice may not be free from the systematic effect, Set-II would be successfully applicable to the low-lying charmonium system.

For comparison, we also examine the light quark mass region. For Set-I, $`m=0.02`$$`0.06`$ corresponds to 90-270 MeV, which is used in the Ref.  as the light quark mass region with the anisotropic Wilson quark action with $`r=1`$. $`E(\stackrel{}{p})E(0)`$ rapidly decrease at the edge of the Brillouin zone, and the height at $`z=a/\pi `$ is around 300 MeV. For two quarks with momenta $`p=\pm a/\pi `$, additional energy of doublers is $`600`$ MeV, and again seems not sufficiently large compared with the typical energy scale transfered inside mesons. In the case of Set-II, this value increases to 1.4 GeV, and seems to be applicable to the meson systems without large systematic effect.

### 2.3 Calibration of quark field

On anisotropic lattice, the anisotropy of quark configuration $`\xi _F`$ must be equal to that of gauge field $`\xi `$.

$$\xi (\gamma _G,\gamma _F)=\xi _F(\gamma _G,\gamma _F)$$
(9)
Since $`\xi `$ and $`\xi _F`$ are function of $`\gamma _G`$ and $`\gamma _F`$ in general, a nonperturbative determination of the combination of $`\gamma _G`$ and $`\gamma _F`$ which satisfy the condition (9) requires much effort. In the quenched case, however, these determination are rather easy to be performed, because $`\xi `$ can be determined independently of $`\gamma _F`$. After the determination of $`\xi `$, one can tune $`\gamma _F`$ so that the certain observable satisfies the condition (9). We call this procedure as “calibration”.

There are several determinations of $`\xi _F`$. In the Ref. , the ratio of the temporal to the spatial meson masses, $`\xi _F=m_\sigma /m_\tau `$, is used as such a observable. However this is not suitable for the present case, since the charm quark mass is larger than or comparable with the spatial lattice cutoffs. In order to determine $`\xi _F`$, we use the dispersion relation of the free meson. For a heavy quark, one may use the nonrelativistic dispersion relation, $`E=m+p^2/2m\xi _F^2`$. In this paper, we alternatively use the relativistic dispersion relation of meson for the calibration. This form is also available for the light quark mass region.

We assume that the meson is described by the following lattice Klein-Gordon action.

$$S=\underset{x}{}\frac{1}{2\xi _F}\varphi ^{}(x)\left[\xi _F^2D_4^{\mathrm{\hspace{0.17em}2}}\stackrel{}{D}^2+m_0^2\right]\varphi (x),$$
(10)
where $`m_0`$ is in the unit of $`a_\sigma `$. Then the free meson satisfies the dispersion relation

$$\mathrm{cosh}E(\stackrel{}{p})=1+\frac{1}{2\xi _F^2}(\stackrel{}{\widehat{p}}^2+m_0^2).$$
(11)
Using this relation, one can determine the anisotropy $`\xi _F`$ as

$$\xi _F^2=\frac{\stackrel{}{\widehat{p}}^2}{2(\mathrm{cosh}E(\stackrel{}{p})\mathrm{cosh}E(0))}$$
(12)
This condition forces the axis-interchange symmetry to the meson field.

## 3 Lattice setup

### 3.1 Gauge configuration

The numerical calculations performed on the two sets of anisotropic lattices, one is with the standard plaquette action and the other is with Symanzik type improved action at the tree level. These actions are represented with the following form;

$`S_G`$ $`=`$ $`{\displaystyle \frac{\beta }{\gamma _G}}{\displaystyle \underset{x,i<j3}{}}\left[c_{11}(1P_{ij}(x))+c_{12}(2R_{ij}(x)R_{ji}(x))\right]`$ (13)
$`+\beta \gamma _G{\displaystyle \underset{x,i3}{}}\left[c_{11}(1P_{i4}(x))+c_{12}(2R_{i4}(x)R_{4i}(x))\right],`$
where the plaquette $`P_{\mu \nu }(x)`$ and rectangular loop $`R_{\mu \nu }(x)`$ are defined as follows,

$$P_{\mu \nu }\frac{1}{3}\text{Re}\text{Tr}[U_\mu (x)U_\nu (x+\widehat{\mu })U_\mu ^{}(x+\widehat{\nu })U_\nu ^{}(x)]$$
(14)
$$R_{\mu \nu }\frac{1}{3}\text{Re}\text{Tr}[U_\mu (x)U_\mu (x+\widehat{\mu })U_\nu (x+2\widehat{\mu })U_\mu ^{}(x+\widehat{\nu }+\widehat{\mu })U_\mu ^{}(x+\widehat{\nu })U_\nu ^{}(x)].$$
(15)
The standard action is the case with $`c_{11}=1`$ and $`c_{12}=0`$, and the improved action is $`c_{11}=5/3`$ and $`c_{12}=1/12`$. The Set-I is the same configurations as used in the Ref. . The parameters for these configurations are summarized in Table 1. These parameters are adopted so that the spatial and temporal lattice extent are about 3 fm respectively.

Both numerical calculations are done on the quenched configurations which are generated by the pseudo-heat bath algorithm with 20,000 thermalization sweeps, the configurations being separated by 2,000 sweeps. These configurations are fixed to Coulomb gauge. The statistical errors are estimated using the jackknife method unless mentioned explicitly.

We determine the parameters for the gauge field configurations, renormalized anisotropy $`\xi (a_\sigma /a_\tau )`$ and spatial cutoff scale $`a_\sigma ^1`$. We define the ratios of the Wilson loop on spatial-spatial ($`\sigma \sigma `$) and spatial-temporal ($`\sigma \tau `$) plane.

$`R_\sigma (r,x)`$ $`=`$ $`W_{\sigma \sigma }(r,x)/W_{\sigma \sigma }(r+1,x),`$ (16)
$`R_\tau (r,t)`$ $`=`$ $`W_{\sigma \tau }(r,t)/W_{\sigma \tau }(r+1,t)`$ (17)
Then $`\xi `$ is determined so that the matching condition, $`R_\sigma (r,x)=R_\tau (r,t=\xi x)`$, is satisfied. From this analysis we conclude $`\xi =5.3(1)`$ (Set-I) and $`\xi =3.950(18)`$ (Set-II). Cutoff scales are determined from the static quark potential using the physical value of string tension $`\sqrt{\sigma _{\text{phys}}}=427\text{MeV}`$. The spatial cutoff scales are $`a_\sigma ^1=0.85(3)`$ GeV (Set-I) and $`1.610(14)`$ GeV (Set-II) respectively.

From the calculation of Polyakov loops and its susceptibilities, we find that the temperatures of lattices with $`N_t=18`$ (Set-I) and $`24`$ (Set-II) are both just above $`T_c`$. Then we estimate $`T_c=250260`$ MeV from the temporal cutoff scale on our lattices. These estimation are consistent with the other quenched lattice results. We prepare the finite temperature gauge configurations at just below $`T_c`$, just above $`T_c`$ and about $`1.5T_c`$ for each sets. These results are summarized in Table 2.

Here we notice that the hadronic correlator depends on which the sector Polyakov loop stays in $`Z_3`$ space at $`T>T_c`$. Since we treat the quenched QCD as an approximation of the full QCD, we choose the real sector for the value of Polyakov loop at $`T>T_c`$.

### 3.2 Mean-field value on the anisotropic lattice

We state the calculation of mean-filed value on the anisotropic lattice for the mean-field improvement of the clover coefficients.

There are two commonly used methods. One is determined from the expectation value of plaquette. This is widely used for its easiness to measure. On anisotropic lattice these are determined as follows.

$$u_\sigma =\frac{1}{3}\text{Re}\text{Tr}P_{ij}(x)^{\frac{1}{4}},u_\tau =\frac{1}{3u_\sigma }\text{Re}\text{Tr}P_{i4}(x)^{\frac{1}{2}}$$
(18)
The other is the trace of link calculated in the Landau gauge.

$$u_\sigma =\frac{1}{3}\text{Re}\text{Tr}U_i(x),u_\tau =\frac{1}{3}\text{Re}\text{Tr}U_4(x)$$
(19)
In the former case $`u_\tau `$ is greater than the unity on our lattice. Then the latter definition seems more reasonable. In this work we determined the mean-field values in the Landau gauge.

The Landau gauge fixing is realized by maximization of Eq.(20).

$$\underset{x}{}\text{Re}\text{Tr}\left(U_1(x)+U_2(x)+U_3(x)+c_{st}U_4(x)\right)$$
(20)
Here, the temporal coefficient $`c_{st}`$ appears on the anisotropic lattice. We adopt the self-consistent mean-field improvement of $`c_{st}`$. As the tree level we chose $`c_{st}=\xi ^2`$. Using the mean-field improved $`c_{st}=\xi ^2u_\sigma /u_\tau `$ we calculate the mean-field values recursively. We perform these calculation with 20 configurations. In our case, the result of 3rd measurement is consistent with the input $`u_\sigma /u_\tau `$ which is determined from the linear interpolation of the tree level and first mean-field improvement result. These results are summarized in Table 3 together with the mean-field improved $`c_E`$ and $`c_B`$. Here we mention that the mean-field improved $`\gamma _G`$ give the reasonable estimation for $`\xi `$ within 1% (Set-I) and 6% (Set-II) error.

### 3.3 Calibration result

We now turn to the calibration of quark field. It is performed by the Eq.(11) with the momentum $`\stackrel{}{p}=(0,0,0)`$ and $`(0,0,1)`$ for each mesonic channel, pseudoscalar (Ps) and vector (V). These results are summarized in Fig. 2. These calibrations are performed at the parameters which correspond to the pseudoscalar mass $`m_{Ps}2.6,3.0`$ and $`3.4`$ GeV respectively.

The detailed definition of meson correlator which is adopted here is shown in Sect. 4.1 In these calculations we use the smeared source function which is determined from the wave function at $`1/\kappa =9.868`$ (Set-I) and $`\kappa =8.797`$ (Set-II). Here we mention that the shape of wave function has very weak $`\kappa `$ dependence.

From the Fig. 2, we find that the $`\gamma _F`$ dependence of $`\xi _F`$ is almost linear and $`\kappa `$ dependence of its slope is small. The $`\gamma _F`$ which is determined from the pseudoscalar and the vector gives consistent value within the statistical error for each $`\kappa `$. Especially $`\xi _F`$ for Set-II are in good agreement with the pseudoscalar and the vector. We determine the $`\gamma _F`$ which satisfies the condition of Eq.(9) from the interpolation with the linear fitting. The error of $`\gamma _F`$ is estimated from the error of $`\xi `$ and $`\xi _F`$. These calibration results are summarized in Table 4. Here we mention that the mean-field improved $`\gamma _F`$ give the reasonable estimation for $`\xi _F`$ within 9% (Set-I) and 1% (Set-II) error.

## 4 Results at Zero temperature

### 4.1 Meson correlators

The mesonic operator which is used in this paper is in the following form,

$$𝒪_M^{(\omega )}(\stackrel{}{z},t)=\underset{\stackrel{}{y}}{}\omega (\stackrel{}{y})\overline{q}(\stackrel{}{z},t)\gamma _Mq(\stackrel{}{z}+\stackrel{}{y},t),$$
(21)
where $`\omega (\stackrel{}{y})`$ is a smearing function and $`\gamma _M=\gamma _5`$ and $`\gamma _1`$ for the pseudoscalar and the vector respectively. Using these mesonic operators we construct the meson correlators $`G_M^{\omega ^{}\omega }(t)`$ in the Coulomb gauge as the following,

$`G_M^{\omega ^{}\omega }(t)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{z}}{}}𝒪_M^{(\omega ^{})}(\stackrel{}{z},t)𝒪_M^{(\omega )}(0)`$ (22)
According to our strategy, we need to construct a good operator which has large overlap with the ground state ( or the excited state ) of heavy quarkonium. For this purpose we examine two types of correlators.

In the first type of correlator we use the smeared source and point sink ( $`\omega ^{}(\stackrel{}{y})\delta (\stackrel{}{y})`$ ), which is already used in the calibration described in Sect. 3.3 . This source function is defined from the measured wave function $`w(\stackrel{}{x})`$ so that the smearing reflects the actual distribution of quark and anti-quark. This smearing procedure works well for the suppressing higher excited state to the correlator. The wave function is well fitted to the form

$$w(\stackrel{}{y})\mathrm{exp}(a|\stackrel{}{y}|^p).$$
(23)
The determined parameters $`a`$ and $`p`$ are used to generate the smearing function. These correlators are used in the charmonium spectroscopy in the next subsection.

For the study at $`T>0`$ we need more systematic optimization of the correlators. We apply the variational analysis, and regard the diagonalized correlator as the optimized correlator. This analysis can extract not only the ground state but also the excited state. In this analysis we use the correlators with smeared source and sink, whose smearing functions are determined by solving the Schrödinger equation with the potential model. This technique is mentioned in Sect. 4.3 .

### 4.2 Charmonium spectroscopy

We show the results of spectroscopy for the heavy quarkonium at $`T=0`$, which is the basis of the study at $`T>0`$. In this section we also examine how our action works in the heavy quark system. The pseudoscalar and the vector meson channels are measured with the parameters listed in Table 4. Figure 3 show the effective mass plot at $`T=0`$. In this calculation we use the correlators with the smeared source introduced in the previous subsection. Here the parameters for the smearing function are $`a=0.7096,p=1.1267`$ for Set-I and $`a=0.2978,p=1.1918`$ for Set-II.

The plateaus of the effective mass appear beyond $`t10`$ (Set-I), $`20`$ (Set-II). Concerning the determination of smearing function our choice works well in heavy quark system. On the other hand in the light quark system the wider smearing function enhances the ground state contribution. In the small $`t`$ region, contributions from the excited states remain. However this small $`t`$ region is the main stage of the study at $`T>0`$. In order to estimate the temperature effects of the ground state contribution at the high temperature, we need further improvement of the mesonic operator. In the next subsection more systematic study with variational analysis is examined.

We fit the correlator to the single exponential form at $`t=30`$$`36`$ (Set-I), $`40`$$`48`$ (Set-II). These results are summarized in Table 5.

The parameters with $`1/\kappa =9.868`$ (Set-I) and $`1/\kappa =8.797`$ (Set-II) in Table 5 correspond approximately to the charm quark. Therefore we study the meson correlators with these parameters in the successive calculations.

Figure 4 shows the hyperfine splitting for charmonium ($`m_{J/\psi }m_{\eta _c}`$). For comparison, results of other group are shown simultaneously. Here we notice that these results largely depend on how to determine the lattice cutoff scale. Results of Fermilab action on the isotropic lattice, whose scale is determined from the physical value of the string tension, are roughly consistent with our results. Results of heavy relativistic action on anisotropic lattice and NRQCD show the similar tendency. These scales are determined from the Sommer scale and the 1P-1S splitting respectively. However, all the quenched results are roughly a half of experimental results $`117(2)`$MeV. The differences from the experimental value of the hyperfine splitting can be partly explained with the dynamical quark effects.

### 4.3 Variational analysis

We examine the variational analysis for constructing more optimized operators. With this analysis we can construct series of operators which have better matching to the states of interest.

Here we explain the principle of variational analysis briefly. Firstly, we suppose that the state generated by mesonic operator on the lattice is the linear combination of eigenstates for Hamiltonian. Practically we assume that the generated states consist of $`N`$ linearly independent states. Then we prepare $`N`$ mesonic operators, which have the same quantum numbers. These operators are constructed with different smearing function $`w(\stackrel{}{x})`$ in Eq.(21). Then we get the $`N\times N`$ correlator matrix, $`C_M^{ij}(t)`$, as the following,

$$C_M^{ij}(t)=\underset{\stackrel{}{z}}{}𝒪_M^{(\omega _i)}(\stackrel{}{z},t)𝒪_{M}^{(\omega _j)}{}_{}{}^{}(0).$$
(24)
Because of the symmetric nature of this matrix, we can get the diagonalized correlators which are optimized correlators for the state of interest.

In our study this variational analysis is applied in the minimal space including the ground, such as 1st excited and 2nd excited state. We have to prepare the smearing functions so that this analysis works well in this condition. Thus we construct the smearing functions $`\varphi _l(r)`$ using the Schrödinger equation with the potential model,

$`\left[{\displaystyle \frac{1}{2m_R}}{\displaystyle \frac{d^2}{dr^2}}+{\displaystyle \frac{l(l+1)}{2m_Rr^2}}+V(r)\right]y_l(r)`$ $`=`$ $`Ey_l(r),`$
$`y_l(r)`$ $`=`$ $`r\varphi _l(r),`$ (25)
where $`V(r)`$ is the static quark potential measured on our lattice and $`m_R=1.5/2`$ GeV. The spin interaction is neglected and we calculate only S-state ($`l=0`$). Figure 5 shows the lowest three solutions $`\varphi _0(r)`$ and the measured wave function at $`T=0`$ for Set-II.

We calculate the diagonalized correlator using the three types of smearing function in Fig. 5. The orthogonal matrices are obtained at each $`t`$, then we adopt the averaged one as the orthogonal matrix which is used for the calculation of the diagonalized correlator. With the orthogonal matrix we obtain the diagonalized correlator. Then its effective mass plots (Set-II) are shown in Fig. 6.

We fit the these data at $`t=10`$$`36`$ for Set-I and $`t=14`$$`48`$ for Set-II, and get the results summarized in Table 6.

The extracted mass of 1S state are consistent with previous results in Table 5. In the Table 6 results are presented in $`a_\tau `$ unit and physical unit for each sets, where the latter case includes the error of $`a_\tau `$.

The results of 2S-1S splitting are consistent with the experimental values within a statistical error. This is in contrast to the case of the hyperfine splitting in which we obtain the half of the experimental value, 2S-1S splitting is in good agreement with experimental value.

The variational analysis can directly extract one of correlators for the ground and excited states. Therefore this analysis is useful to investigate the excited state of hadron. In the next section we mainly use this analysis at $`T>0`$.

## 5 Results at Finite temperature

### 5.1 Temperature dependence of correlators

In accordance with our strategy the optimized operators which are constructed at $`T=0`$ are applied to the case at $`T>0`$. The optimized operators are defined from the variational analysis in the previous section. We construct the correlator matrix at $`T>0`$ with the same basis ( smearing function ) as $`T=0`$ ( Fig. 5). Then the orthogonal matrix defined at $`T=0`$ are operated to this correlator matrix, and its diagonal correlators are discussed in this subsection. Here we call these correlators as $`C_1(t)`$, $`C_2(t)`$ and $`C_3(t)`$ in the order of their effective masses. Figure 7 shows the effective masses of $`C_1(t)`$ and $`C_2(t)`$ for Set-II, although the effective mass of $`C_3(t)`$ is too noisy.

Below $`T_c`$ ( $`N_t=26`$ ) we can find the plateau of effective mass. The ground and the excited states seem to be observed below $`T_c`$. The masses of these states are almost the same as $`T=0`$ or slightly larger. Therefore we conclude that these correlators have little thermal effects at this temperature, and the spectral structure seems to keep the form at $`T=0`$. However it is difficult to identify the plateau precisely and determine the mass quantitatively with the present statistics. More detailed analysis with higher statistics may open a stage to discuss the potential mass shift of charmonium near to the $`T_c`$.

Above $`T_c`$ ( $`N_t=22`$ and $`16`$ ) effective masses have no clear plateau in whole $`t`$ region. These behaviors at least signal significant change of correlators when the system crosses $`T_c`$. This behavior appears noticeably in the effective mass of $`C_2(t)`$. In the case of the light quark system investigated in the Ref. , the effective masses increase as $`T`$ in the pseudoscalar and the vector channels. The observed behavior in present work, however, shows qualitatively different nature of the correlators.

As the comparison with above results, Fig. 8 shows the effective mass plots at $`T>0`$ for the correlators with smeared source by the wave function and point sink, described in Sect. 4.1 .

These results are consistent with the variational analysis. The remarkable change in the vector channel is apparently seen. Then the order of effective masses for the pseudoscalar and the vector mesons is reversed above $`T_c`$, the same as the case of free quark. This phenomenon is also brought about in the system of light quarks.

From the results of this subsection, we observed that the temporal mesonic correlators change drastically when the system goes through the deconfining transition. We can consider the several pictures which can explain (consistent with) the change for correlators. One of these pictures is that the mesonic bound state disappear in the deconfinement phase. This is the most interesting case as stated in our motivations. In order to discuss this picture we investigate the correlation between $`c`$ and $`\overline{c}`$ in the next subsection.

### 5.2 Wave function

In this subsection the $`c\overline{c}`$ bound state at $`T>0`$, especially in the deconfinement phase, is discussed in the light of “wave function”. The definition of the “wave function” in the Coulomb gauge is as follows.

$$w_M(\stackrel{}{r},t)=\underset{\stackrel{}{x}}{}\overline{q}(\stackrel{}{x}+\stackrel{}{r},t)\gamma _Mq(\stackrel{}{x},t)𝒪_M^{}(0)$$
(26)
Here this definition is the same form as the wave function at $`T=0`$. This wave function shows the spatial correlation between $`q`$ and $`\overline{q}`$, and gives us a hint of the mesonic bound state from its $`t`$ dependence. In the case of free quarks, $`q\overline{q}`$ has no bound state, then the wave function ought to broaden with $`t`$. On the other hand suppose quark and anti-quark form a bound state, the wave function holds the stable shape with $`t`$. We can discuss the existence of such the bound state by observing the $`t`$-dependence of the wave function. For this purpose we compare the correlation at spatial origin with another spatially separated point at each $`t`$. Therefore we define the wave function normalized at the spatial origin, $`\varphi _M(\stackrel{}{r},t)`$, as follows,

$$\varphi _M(\stackrel{}{r},t)=\frac{w_M(\stackrel{}{r},t)}{w_M(\stackrel{}{0},t)}.$$
(27)
From now on the wave function denotes this normalized definition.

Since the question is whether wave function has a stable shape or not, it is not necessary to use the optimized operator. Therefore the smeared source function with exponential form defined as with Eq.(23), is used for the analysis of $`t`$-dependence of the wave function.

Fig. 9 shows the results at $`T>0`$ with the smeared source function which is slightly wider than the observed wave function at $`T=0`$. The wave functions composed of free quark propagators are also shown together.

As is shown in the Fig. 9, the behaviors of the observed wave functions are clearly different from that of the free quark case at each temperature and in each mesonic channel. In the free quark case the wave functions are broadening as $`t`$ as expected. On the other hand, the observed wave functions are stable with the slightly narrower shape than source function. These behaviors are independent of the source function.

For a visible expression, we define the averaged orbital radius, $`r_0`$, as

$$r_0^3(t)=\underset{z}{}3z^2\varphi _M(z,t),$$
(28)
where we suppose a spherical symmetric wave function and the sum is over $`z`$ axis. These $`r_0(t)`$ of Set-I and Set-II are shown in Fig. 10.

Figure 10 shows the $`t`$ dependence of $`r_0`$ in the physical unit, where the error estimation takes into account the error of $`a_\sigma `$ and $`a_\tau `$. These results for Set-I and II are roughly consistent with each other. The behavior of observed $`r_0`$ are obviously different from free quark case, and shows the stable behavior with respect to $`t`$. This asymptotic value for the vector channel is larger than that of the pseudoscalar.

Below $`T_c`$, $`r_0`$ for each channel are almost same as the case of $`T=0`$. At $`T1.5T_c`$ the observed $`r_0`$ are slightly larger than that of $`T<T_c`$. However the wave function, even in the deconfinement phase, seem to have the stable $`r_0`$ independently of the source function. Therefore we conclude the same strong spatial correlation as below $`T_c`$ survives in the deconfinement phase at each mesonic channel.

The results of this subsection seem to suggest the presence of the hadronic state in the deconfinement phase, at least up to $`1.5T_c`$. This is the opposite situation to the picture mentioned in the previous subsection.

## 6 Conclusion

In this paper we explored the charmonium correlators in the Euclidean temporal direction at $`T>0`$ using the quenched lattice QCD simulation. The high resolution in the temporal direction is achieved by employing the anisotropic lattice. We examined the thermal effect on the correlators based on the following two quantities: the correlators between the optimized operators tuned at $`T=0`$, and the $`t`$-dependence of $`q\overline{q}`$ wave function.

In the calculation of the wave function in the Coulomb gauge, we found the strong spatial correlation even above $`T_c`$, up to at least $`1.5T_c`$. These results indicate the quark and the anti-quark tend to be close each other even in the deconfinement phase. The similar result is reported in the light quark mass region in the Ref. . In the case of charm quark mass these results suggest that the $`J/\psi `$ may not easily be resolved until $`1.5T_c`$. On the other hand, we also observed significant change of the nature of the correlators between the operators optimized at $`T=0`$. This was signaled by the drastic change of the behavior of effective mass. This situation is interesting, and at the same time puzzling. It is possible to consider several pictures which are able to explain our results. For example, the mesonic spectral function still have some peaks above $`T_c`$, and its width are broadened with thermal effects. Such a situation naturally explains the observed results in this work. The existence of hadronic modes just above $`T_c`$ were also suggested by the previous works. From the analysis of temporal meson correlators at $`T>0`$, the QCD vacuum still has non-perturbative nature above the phase transition, and far from the perturbative plasma state of quark and gluon. In spite of several systematic uncertainties our results in the present work have important and interesting implications on the fate of hadronic states above the critical temperature.

Our approach is directly applicable to the dynamical configuration. If the full QCD simulation on the anisotropic lattice is appropriately implemented. Then it is interesting to investigate what effect our results receive from the dynamical quarks. To achieve the high resolution in the temporal direction, we adopted the anisotropic lattice and employed the O(a) improved Wilson quark action on it. These implementation also useful to study the heavy particle such as the glueballs and scalar mesons. For the correlators of these states we are inevitably forced to extract the signals at the short time separation. The large temporal lattice cutoff enables us to simulate the heavy particles on lattices of moderate size with keeping the finite lattice cutoff effect small. The detailed information in the temporal direction is also significantly useful for the direct extraction of the spectral function from the lattice data. This approach may give us the further information on the spectral structure of mesons at $`T>0`$.

## 7 Acknowledgments

We thank the members of QCD-TARO Collaboration for interesting discussion. H.M. thanks T.Onogi, N.Nakajima and J.Harada for useful discussion. The simulation has been done on Intel Paragon XP/S and NEC HSP at the Institute for Nonlinear Science and Applied Mathematics, Hiroshima University, NEC SX-4 at Research Center for Nuclear Physics, Osaka University and Hitachi SR8000 at KEK. This work is supported by the Grant-in-Aide for Scientific Research by Monbusho, Japan ( No. 10640272, No. 11440080 ). H.M. is supported by the center-of-excellence (COE) program at RCNP, Osaka University.


text:
# Jahn-Teller, Charge and Magnetic Ordering in half-doped Manganese Oxides

## 1 Model and approximations

### 1.1 Hamiltonian

According to the arguments given in the introduction, we assume an infinite repulsion ($`U=\mathrm{}`$) between electrons on the same lattice site, and an infinite Hund coupling ($`J_H=\mathrm{}`$) between the localized $`t_{2g}`$ spins and the itinerant $`e_g`$ spins. One can therefore consider spinless electrons, their spin degree of freedom being unequivocally defined by the direction of the local $`t_{2g}`$ spins $`\stackrel{}{S}`$. Furthermore, we consider in this work a two-dimensional plane of the structure, with a half-filled band made of a single $`e_g`$ orbital. The effective model Hamiltonian is then:

$$H=H_{DE}+H_{Coul}+H_{ph}+H_{SE}+H_H$$
(1)
with

$`H_{DE}`$ $`=`$ $`{\displaystyle \underset{<ij>}{}}\stackrel{~}{t}_{ij}c_i^{}c_j`$
$`H_{Coul}`$ $`=`$ $`{\displaystyle \underset{<ij>}{}}V(n_in)(n_jn)`$
$`H_{ph}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}[K_bQ_{bi}^2+K_2Q_{2i}^2+K_sQ_{si}^2]`$
$`{\displaystyle \underset{i}{}}g_2Q_{2i}(n_in)+{\displaystyle \underset{i}{}}g_bQ_{bi}(n_in)L_s{\displaystyle \underset{<ij>}{}}Q_{si}Q_{2j}`$
$`H_{SE}`$ $`=`$ $`{\displaystyle \underset{<ij>}{}}[J_1J_2Q_s]\stackrel{}{S}_i\stackrel{}{S}_j`$
$`H_H`$ $`=`$ $`g\mu _b\stackrel{}{H}{\displaystyle \underset{i}{}}\stackrel{}{S}_i`$
The first term $`H_{DE}`$ represents the double exchange hopping of electrons on a square lattice. Here $`c_i^{}`$, $`c_i`$ are respectively the creation and annihilation operators for spinless electrons from a single band, and $`\stackrel{~}{t}_{ij}=t\mathrm{cos}(\theta _{ij}/2)`$ is the transfer integral between neighboring Mn sites whose ionic spins $`𝐒_i`$ and $`𝐒_j`$ make an angle $`\theta _{ij}`$ DeGen .

The second term $`H_{Coul}`$ describes the Coulomb repulsion between nearest neighbors ($`n_i=c_i^{}c_i`$ and $`n`$ is the average electron density, which is equal to $`1/2`$ in the present case).

The third term $`H_{ph}`$ is the elastic part, which includes the coupling of electrons to a Jahn-Teller (JT) mode $`Q_2`$ and of holes to a “breathing” mode $`Q_b`$ ($`g_2`$ and $`g_b`$ are the coupling strengths, $`K_2`$ and $`K_b`$ the spring constants). In the planar geometry considered here, the other Jahn-Teller mode $`Q_1`$ is not relevant. We have also introduced a shear mode $`Q_s`$, which is driven by $`Q_2`$. Such a shear deformation, which is experimentally observed at low temperatures, is essential to reconcile the alternating Mn4+ breathing and Mn3+ JT distortions wich develop in the ordered phases. A substantial shear deformation is indeed observed in La0.5Ca0.5MnO3 radaelli . It results in some $`MnOMn`$ bonds being shorter and other larger (”zig-zag” chains, see Fig. 1).

The term $`H_{SE}`$ represents the antiferromagnetic (AF) superexchange interaction $`J_1`$ between the ionic spins on neighboring sites, which are treated as classical. The additional term $`J_2Q_s`$ is a phenomenological implementation of the Goodenough rule: It can either enhance or reduce the AF coupling depending on the sign of the shear deformation, which accounts for the fact that longer (shorter) Mn-Mn bonds have a more (less) antiferromagnetic character Good . The last term $`H_H`$ takes into account the external magnetic field.

We shall study the Hamiltonian (1) in the mean-field approximation, describing the charge ordered (CO) phase as a charge density wave (CDW) with momentum $`(\pi ,\pi )`$. Let us call $`\overline{n}^A`$ and $`\overline{n}^B`$ the average electron densities in the two resulting sublattices, which correspond respectively to the Mn3+ and Mn4+ ions. We shall further assume that the JT coupling is only active on A sites, while the breathing deformations arise on B sites. With these approximations, the terms in the Hamiltonian which depend explicitely on $`(n_in)`$ reduce to

$$H_{MF}=\mathrm{\Delta }\underset{iAorB}{}(n_i^An_i^B)+const$$
(2)
where the order parameter $`\mathrm{\Delta }`$ is defined as

$$\mathrm{\Delta }=2V(\overline{n}^A\overline{n}^B)+(g_bQ_b+g_2Q_2)/2$$
(3)
and the chemical potential has been set to zero by adding a term $`\mathrm{\Delta }\mu =(g_2Q_2g_bQ_b)/2`$ to recover particle-hole symmetry (with these notations, the choice $`\overline{n}^A\overline{n}^B`$ corresponds to the $`Q`$’s being all positive).

The magnetic part is also treated in mean-field, according to de Gennes’ procedure DeGen , using a gaussian distribution for the angle of the classical spins with respect to the mean field direction. We consider the following magnetic phases: Ferromagnetic (F), paramagnetic (P), Néel anti-ferromagnetic (NAF), and CE-type ordering, with ferromagnetic zig-zag chains, coupled anti-ferromagnetically (CE). The most general unit cell which allows to describe all these phases is is made of 8 nonequivalent Mn sites in a plane (Fig. 1). In each of these magnetic configurations, the total free energy is minimized with respect to the following parameters: i) the magnetization on non-equivalent magnetic sites, ii) the average electron density $`\overline{n}^A`$ on sublattice A ($`\overline{n}^B`$ being just $`1\overline{n}^A`$) and iii) the lattice displacements.

### 1.2 Phenomenological treatment of screening effects

As mentioned in the Introduction, to have a realistic description of the phase diagram which includes both the metallic and charge ordered phases at half-filling, it is necessary to go beyond the simple mean-field approach described in the preceding section by including the effects of exchange and correlation. Such effects on the intersite Coulomb repulsion in the half-doping manganites have recently been analyzed within an RPA-like calculation sheng , which is known to be appropriate for interacting electron systems at metallic densities. Since a detailed study goes beyond the scope of this paper, we propose here a semi-phenomenological treatment of screening which allows a qualitative description of the transition between CO and metallic states, and which correctly reproduces the results of reference sheng . The method is further generalized to describe the screening of the electron-lattice interactions. Actually, the latter are due to the Coulomb repulsion between $`Mn`$ and $`O`$ ions and are therefore also screened in the metallic phase. This screening should be weaker than that of $`MnMn`$ interactions, since it involves $`MnO`$ rather than $`MnMn`$ charge fluctuations, but it should be sizeable.

The procedure will be carried out in two successive steps. The first step consists in writing a reasonable estimate for the exchange-correlation energy $`E_{xc}`$, which is defined as the correction to the ground-state energy beyond the Hartree mean-field result. In the second step, we shall define an effective hamiltonian $`H_{xc}`$ to be treated in the mean-field approximation, such that

$$H_{xc}=E_{xc}$$
(4)
This results in a modification of the atomic energy levels $`\pm \mathrm{\Delta }`$ (i.e. in a reduction of the CDW gap), and it yields a correction to the free energy which is precisely of the form $`E_{xc}`$.

#### 1.2.1 Exchange-correlation energy

Let us start by analyzing the simple ferromagnetic case at $`T=0`$, where the electron hopping is not renormalized by the mechanism of double exchange. In the metallic phase, which corresponds to a vanishing order parameter $`\mathrm{\Delta }`$, the leading correction comes from the exchange (Fock) terms. These terms are responsible for an increase of the carrier itinerancy, which can be viewed as a renormalization of the hopping parameter $`tt+Vc_i^{}c_{i+\delta }`$. Hence the kinetic energy is lowered by a quantity proportional to the interaction potential $`V`$, and one can write $`E_{xc}=aV`$ (the parameter $`a`$ is related to the dielectric constant of the system). On the other hand, in the CO phase, i.e. at strong $`\mathrm{\Delta }`$, the correlation energy corresponds to the interaction between density fluctuations on neighboring sites, each of them being proportional to $`\delta nt/\mathrm{\Delta }`$. Therefore, in this case the appropriate limiting formula is $`E_{xc}V(t/\mathrm{\Delta })^2`$. These results can be generalized to the screening of the electron-lattice interactions, by replacing $`Vg^2/K`$ and by introducing the corresponding order parameter $`\mathrm{\Delta }`$ as given by eq. (3).

A smooth interpolation between the weak and strong coupling behaviors can be obtained by writing the following expression for the exchange-correlation energy:

$$E_{xc}=\frac{aV+b(g_b^2/K_b+g_2^2/K_2)}{1+c\left(\frac{\mathrm{\Delta }}{t}\right)^2}$$
(5)
where $`a,b`$ and $`c`$ are phenomenological parameters (the ratio $`a/c=1.44`$ can be deduced from ref. sheng  and $`b/a=1/10`$ is chosen according to the ionic distances).

As was stated at the beginning of this section, this formula is only appropriate in the ferromagnetic case. It does not account for the fact that the mobility of the carriers taking part in the screening process is affected by the magnetic structure through the DE mechanism. We shall give here the arguments which allow a generalization of eq. (5) to the different kinds of magnetic orderings.

In the free-electron limit ($`\mathrm{\Delta }0`$), where screening is due to coherent band motion, one expects the correlation energy to be reduced by a factor $`\stackrel{~}{t}/t`$, where

$$\stackrel{~}{t}=t\mathrm{cos}(\theta _{ij}/2)$$
(6)
is the effective DE hopping parameter averaged in all space directions (this gives respectively $`1,8/9,1/2`$ in the $`F,P`$, $`CE`$ phases). The situation is slightly more complicated in the charge ordered phases, because $`E_{xc}`$ comes from incoherent hopping of the carriers to neighboring sites. According to Hund’s rule, such processes will be allowed only between sites with parallel spins, which defines an effective number of neighbors $`\stackrel{~}{z}z`$.

In the CE phase, for instance, the lattice can be divided into U (up) and D (down) sites, according to the spin direction. Since each site has 2 U and 2 D neighbors, a given electron can only hop to the 2 neighbors with the same spin direction, and consequently $`\stackrel{~}{z}/z=1/2`$. At finite temperatures, however, thermal fluctuations will reduce the absolute value of the local magnetization $`m`$. Accordingly, there will be a finite probability that a given U site has a $``$ spin, which is given by $`n_U^{}=(1m_U)/2`$ (an equivalent expression holds for D sites). The total probability for hopping away from a U site is therefore

$$2n_U^{}n_U^{}+2n_U^{}n_U^{}+2n_U^{}n_D^{}+2n_U^{}n_D^{}$$

where obviously $`n_U^{}=1n_U^{}`$. By adding the contributions for hopping processes starting from both U and D sites and dividing by 2, we obtain

$$\stackrel{~}{z}=\frac{1}{2}\left[4+(m_U+m_D)^2\right]$$
(7)
which correctly gives $`\stackrel{~}{z}/z=1,1/2,1/2`$ for the F,P,CE phases at $`T=0`$. Here the factors (6) and (7) introduce a feedback on the itinerancy of the electrons in the case of an applied magnetic field, which tends to align the spins ferromagnetically. This effect is essential in reducing the critical $`H`$ for the CE-FM transition at low temperatures, as we shall see below.

For each given magnetic configuration, instead of eq. (5), we shall use the following formula for the exchange-correlation energy:

$$\stackrel{~}{E}_{xc}=\frac{\left[\stackrel{~}{a}V+\stackrel{~}{b}(g_b^2/K_b+g_2^2/K_2)\right]}{1+\stackrel{~}{c}\left(\frac{\mathrm{\Delta }}{t}\right)^2}$$
(8)
where the screening parameters $`a,b,c`$ have been modified according to

$$\stackrel{~}{a}=a\frac{\stackrel{~}{t}}{t},\stackrel{~}{b}=b\frac{\stackrel{~}{t}}{t},\stackrel{~}{c}=c\frac{z}{\stackrel{~}{z}}\frac{\stackrel{~}{t}}{t}$$
(9)
We emphasize here that the terms in the numerator of Eq. (5) are rescaled by the $`\stackrel{~}{t}/t`$ factor since they arise from the coherent screening processes (mostly active when $`\mathrm{\Delta }0`$). On the other hand the local (incoherent) screening processes related to the term in the denominator of Eq. (5) are also rescaled by the effective number of accessible nearest neighbor sites.

#### 1.2.2 Mean-field potential from exchange and correlation

We wish to define an effective hamiltonian $`H_{xc}`$ to be treated in the mean-field approximation such that the correction to the free energy is equal to $`E_{xc}`$. To this purpose, we replace $`\mathrm{\Delta }`$ by an operator $`\widehat{\mathrm{\Delta }}`$ (e.g. the mean-field parameter $`n^A`$ is replaced by $`_{iA}n_i^A`$), and linearize the resulting expression. This gives

$`H_{xc}=B\stackrel{~}{c}{\displaystyle \frac{V}{t}}{\displaystyle \frac{\mathrm{\Delta }}{t}}{\displaystyle \underset{iAorB}{}}(n_i^An_i^B)+const`$ (10)
where we have defined

$$B=\left[\stackrel{~}{a}V+\stackrel{~}{b}(g_b^2/K_b+g_2^2/K_2)\right]/\left[1+\stackrel{~}{c}\left(\frac{\mathrm{\Delta }}{t}\right)^2\right]^2$$
(11)
The constant part in eq. (10) is

$$B\left\{1+\stackrel{~}{c}\frac{\mathrm{\Delta }}{t^2}\left[6V(\overline{n}^A\overline{n}^B)+(g_bQ_b+g_2Q_2)/2\right]\right\}$$
(12)
It is easy to verify that eq. (4) holds when $`\overline{n}^A=n^A`$ and $`\overline{n}^B=n^B`$.

One notices that a dielectric constant can be deduced from the screening of the gap, by writing

$$\mathrm{\Delta }\mathrm{\Delta }_{eff}=\mathrm{\Delta }B\stackrel{~}{c}\frac{V\mathrm{\Delta }}{t^2}$$

which gives

$$\epsilon =\frac{\mathrm{\Delta }}{\mathrm{\Delta }_{eff}}=\frac{1}{1c\frac{V}{t^2}\frac{\stackrel{~}{a}V+\stackrel{~}{b}(g_b^2/K_b+g_2^2/K_2)}{1+\stackrel{~}{c}\left(\frac{\mathrm{\Delta }}{t}\right)^2}}$$

## 2 Results

### 2.1 The phase diagram: Existence of a CE phase

The Hamiltonian in Eq. (1) is formed by several competing terms, and the corresponding phase diagram contains several phases, each one dominating in some region of the parameter space. To make the analysis simpler, we choose to vary together those parameters having similar physical effects. In particular, the electron-phonon couplings generically reduce the electron mobility and, at mean-field level, tend to give rise to a staggered charge ordering, acting similarly to the n.n. electron-electron repulsion $`V`$. Therefore in varying $`V`$ we keep constant the ratio $`V/(g^2/K)`$. For the sake of simplicity, we also keep a fixed $`J_2/J_1`$ ratio, although it is not the only possible choice.

Figure 2 reports typical phase diagrams at various temperatures as a function of the magnetic coupling $`J=J_1S^2/t`$ with $`J_1=J_2`$ ($`S=3/2`$) and of the repulsive e-e interaction $`V=0.5(g^2/K)`$. At low temperature (left panel) there is a metallic (i.e. without charge-ordering) ferromagnetic phase (FM) when the charge-ordering (CO) terms $`V`$ and $`g^2/K`$ are not too large. This FM phase is naturally suppressed by the increase of the antiferromagnetic (AF) superexchange coupling $`J_1`$. When the charge mobility is suppressed by the CO terms, one finds two distinct possible phases. At low values of the AF coupling the pure CO effects dominate and a ferromagnetic (F) CO phase occurs at sufficiently large values of $`V`$ (F-CO). The transition is first order, as found in Ref sheng , due to the exchange-correlation terms. On the other hand, by increasing the AF coupling, the CO ferromagnetism is destabilized and a CE phase takes place. This latter phase naturally realizes the best compromise between the electron mobility, favored by the ferromagnetic bonds, the CO, and the AF interactions increasing with $`J_1`$. The CE ordering arises due to competing lattice displacements. In particular a substantial shear mode is induced in the lattice to reconcile the (asymmetrical) JT deformations occurring in the Mn3+ ions with the (centrosymmetric) breathing deformations around the Mn4+ ions. The resulting lattice structure displays zig-zag chains formed by long bonds interlaced with zig-zag chains of short bonds. Then the peculiar CE magnetic structure naturally appears. In particular, according to Goodenough Good , orbital ordering makes the sign of the magnetic couplings to be correlated to the length of the bonds, with AF (F) magnetic couplings corresponding to short (long) bonds. Therefore, the lattice-driven chains with short bonds and with long bonds naturally translate into a lattice-driven CE magnetic structure.

The temperature evolution is represented in the three panels of figure 2. By increasing $`T`$ the weak ferromagnetism surviving in the CO phase at small $`J`$ is rapidly suppressed in favor of a a CO paramagnetic phase. Also the CE region, being due to a delicate balance between CO, FM, and AF interactions, is reduced rather rapidly. The FM phase at small values of $`J_1`$ is instead based on the double-exchange mechanism, which is more robust and, upon increasing $`T`$, is only weakly “invaded” by the CO paramagnetic phase. One also observes that compounds with a low-temperature CE magnetic structure, but laying near the CE-FM phase boundary, can undergo a first-order CE-to-FM transition upon increasing the temperature.

It is worth pointing out that, within our formal scheme, the CO paramagnetic phase is the only possible mean-field description of an insulating non-magnetic phase at moderate temperature values. The explored temperature range is indeed much too low to allow for a thermal disruption of the CO, which instead occurs at temperatures of the order of $`V`$. However, in strong coupling, one may speculate that a more refined description would allow for the disordering of the charge (possibly without spoiling the local lattice deformations around the charges) thus producing the disordered paramagnetic (and polaronic) phase which is observed in all manganites above a few hundreds of kelvins.

### 2.2 Sensitive role of the electronic bandwidth

A key issue is the role of the kinetic energy in the competition between the different phases. An extended experimental analysis of the phase diagrams of the various manganites tokura ; tokura2  suggests that the electronic bandwidth, among other parameters such as lattice disorder, plays a primary role in determining the stability and the competition between the FM and the insulating phases. In our model we investigated this relevant point by varying the bare hopping amplitude in front of the double-exchange term in the Hamiltonian (1). We also assume that the same mechanism inducing the variation of the nearest neighbour hopping of the itinerant electrons in the $`e_g`$ orbitals is responsible for variations of the hopping of the $`t_{2g}`$ electrons as well. This affects the superexchange couplings $`J_1`$ and $`J_2`$, in particular $`J_1`$ is expected to arise from second-order hopping processes of the $`t_{2g}`$ electrons $`J_1t^2/U_{t_{2g}}`$ (where $`U_{t_{2g}}`$ is some effective repulsion between electrons in the same doubly occupied $`t_{2g}`$ orbital). According to this assumption, when the hopping $`t`$ is increased without changing the intersite repulsion $`V`$ one moves downwards in the phase diagrams of Fig. 2, where the variable $`V/t`$ is reported on the $`y`$-axis. At the same time, however, the increasing $`t`$ produces an increase in the $`x`$-axis variable $`Jt`$. Therefore, by keeping all the Hamiltonian parameters fixed, but $`t`$ and $`J_1=J_2`$, one moves in the phase diagram along the dotted curves $`V/t=A/J`$ reported in Fig. 2. These curves correspond to similar physical systems, where the only nearest neighbour electronic hopping amplitudes have been varied.

It immediately appears that systems which slightly differ by the electronic hopping amplitude can have different magnetic structures at low temperature (Fig. 2, left panel). In particular, by increasing $`t`$ a first-order transition can occur at low temperature from a CE to a FM phase. Furthermore, the temperature evolution of systems with different bandwidth but laying near a phase boundary can be different. This is made apparent in Fig. 3, where we report the temperature-magnetic field phase diagrams for three systems with different hopping (and therefore also different magnetic couplings $`J_1`$) laying on the same dotted curve $`(V/t)=0.04/J`$. The three phase diagrams correspond to systems with slightly different magnetic couplings (differing at most by ten percent). Nevertheless, already at zero magnetic field, the three systems display a completely different evolution in temperature. The more insulating (i.e., with smaller $`t`$) system having $`J=0.037`$ (so that $`(V/t)=1.081`$) never becomes metallic at zero field, but only undergoes a first-order transition from a low-temperature CE phase to a paramagnetic insulating phase at $`T0.032t`$ (left panel, see also Figs.4,5). In the more metallic system (center panel) having $`J=0.038`$ (and $`V/t=1.053`$), the CE phase disappears at a lower temperature $`T0.02t`$ and it is replaced by an intermediate FM phase. The ferromagnetic order and the metallicity is then destroyed at a higher temperature $`T0.032t`$ where a paramagnetic CO phase takes place. We reiterate here that this latter phase is better to be intended as the mean-field representation of a disordered paramagnetic insulating phase. Finally, at even larger values of $`J=0.04`$, corresponding to $`V/t=1`$ the metallic phase is present already at low temperature and it survives up to a $`T0.05t`$.

The relevant role of the kinetic energy in stabilizing the uniform metallic FM phase at the expenses of the CO phases and particularly of the one with CE magnetic order is made even more apparent in the presence of a magnetic field. This is particularly clear in the first panel of Fig. 3, where the FM phase, which would be absent at zero field becomes the most stable solution at large enough $`H`$.

It is also worth mentioning that, due to the presence of screening, the metallic uniform solution is always a (local) minimum of the free energy. Therefore an (at least metastable) metallic solution is present even at zero field. The existence of a (local) minimum is a necessary condition for the occurrence of an hysteretic behavior at the transition. Of course the region of the hysteresis also depends on the height of the free energy barrier between the maxima, of domain walls and of other non-equilibrium properties. Nevertheless the region in T and H where two minima exist provides an (excess) estimate for the hysteresis region experimetally observed in half-doped Pr0.5Sr0.5MnO3, Nd0.5Sr0.5MnO3, (Nd1-ySmy)0.5Sr0.5MnO3 and Pr0.5Ca0.5MnO3 tokura ; tokura2 .

### 2.3 Discussion

The previous subsection illustrated the main results of the present work:

i) The CE phase does arise in the present one-orbital model and is crucially related to competition between the JT and breathing distortions involved in the charge-ordered state on $`Mn^{3+}`$ and $`Mn^{4+}`$ sites respectively. The shear lattice deformation results from this competition and couples to the magnetic degree of freedom, merely through orbital ordering.

ii) Exchange-correlation corrections are essential to stabilize a metallic ferromagnetic phase, due to substantial screening of both Coulomb and electron-lattice interactions.

iii) The kinetic energy is a most effective parameter in determining the relative stability of the various phases upon varying the temperature and the magnetic field.

As far as point i) is concerned, in the present work we show that the CE phase not only arises from electronic mechanisms based on the presence of (at least) two orbitals per Mn site. The existence of a CE phase in a model without orbital degrees of freedom is quite remarkable. It is indeed repeatedly stated in the literature solovyev ; jackeli ; okamoto  that the CE phase is stabilized by the kinetic energy gain arising from the orbital ordering forming ferromagnetic chains. Our results are not in contrast with this viewpoint, but underline that the above purely electronic mechanism is not the only possible one, and that the coupling with lattice degrees of freedom can be of primary importance. In this regard our results are related to previous Hartree-Fock mizokawa  and quantum Monte-Carlo yunoki  calculations, where the JT deformations were claimed to be relevant for the occurrence of a CE phase. Our low-temperature phase diagrams are qualitatively similar to the one reported in Fig. 2(c) in Ref.yunoki  once the distinction between order and disorder in the orbital degrees of freedom is discarded notascalaAF . Our contribution in this framework is to show that the lattice shear deformation is a relevant ingredient in its own even in the absence of cooperative mechanisms due to electronic or JT-induced orbital ordering.

Regarding point iii) we notice that a systematic analysis of the role of the hopping is relevant for the general understanding of the manganites. In the real materials of the general form R1-xAxMnO3 (where R and A are trivalent rare earth and divalent alkaline earth ions respectively) the bandwidth can be varied by changing the radius of the perovskite A site (where the R and A ions are located). Depending on the averaged ionic radius the bond angle of Mn-O-Mn deviates from 1800 in the orthorombic lattice. The smaller the radius of the A site is, the larger is this angle, which reduces the Mn-O overlap and the effective Mn-Mn hopping. A systematic experimental analysis of this effect is reported in Ref. tokura2 . The results summarized in Fig. 3 allow for a unified qualitative description of different half-doped materials. In particular the most insulating behavior in Fig. 3(a) is consistent with the generic features of Pr0.5Ca0.5MnO3. On the other hand, the center panel of Fig. 3 shows the same qualitative behavior of La0.5Ca0.5MnO3 or (Nd1-ySmy)0.5Sr0.5MnO3. Finally the most metallic system in the right panel is a good qualitative description of La0.5Sr0.5MnO3. Nevertheless, it has been pointed out note  that a rapid change in lattice constant $`K`$, rather than necessarily small changes of $`t`$, could be the clue for the very different behaviours of the systems (A,A)0.5MnO3. In our case this would correspond to an abrupt change along a vertical line in Figs. 2, and would enhance the first-order character of the transitions.

A semi-quantitative agreement can even be reached. In fact, the pure double-exchange ferromagnetic critical temperature ($`J/t=0`$) is, from our 2D mean-field calculation, $`T_c0.085t`$. A 3D estimate enhances this value by a factor $`3/2`$ owing to the number of nearest neighbours, yielding $`T_c^{DE}0.13t`$. For an average value $`t=0.3eV`$ one gets a transition temperature $`450K`$. It is reduced by the presence of the antiferromagnetic coupling, for instance in panel c of Figure 3, in zero field $`T_c0.05t`$ thus $`270K`$ in 3D, thus supporting De Gennes’ simple mean-field picture. Then one obtains in the center panel the value $`T_c^{CE}180K`$, and in the left panel $`T_c^{CE}170K`$. These values are reasonable, as compared with experimental ones, in particular one notices that $`T_c^{CE}`$ is strongly reduced compared to $`T_c^{DE}`$. This is due to the competition between the two order parameters (ferro and antiferro). Another way to understand it is to notice that in the CE phase the effective dimensionality is reduced by chain formation, together with charge localization this reduces the effective stength of double-exchange. On the other hand the effective antiferromagnetic exchange is close to that of stoechiometric CaMnO3 with a $`T_c`$ of $`120K`$.

We have also systematically investigated the role of the magnetic field in stabilizing the uniform FM phase. The typical energy differences involved in the first-order transitions are so small that accessible magnetic fields are sufficient to drive the transition from the insulating to the FM phases. Specifically, by taking a typical value of $`t0.3eV`$, one can see that $`H/t=0.015`$ (where $`H=g\mu _bSH`$) roughly corresponds to ten Teslas. This value agrees well with the typical values experimentally used to investigate the (T,H) dependence of the low-temperature CE insulating phase and the intermediate-temperature uniform FM phase tokura ; tokura2 .

## 3 Conclusion

Let us compare our approach with other models which have been proposed to describe the half-doping compounds. In our treatment, the main ingredient which is responsible for the CE-type magnetic ordering is the appearance of a shear deformation, with the consequent modification of the magnetic coupling along certain directions. In ref. yunoki , two different orbitals are retained for the $`e_g`$ electrons, but the shear deformation is absent. In their approach, the CE order arises because the orbitals prefer to have large overlaps along certain directions, thus favoring the kinetic energy along zig-zag chains where the spins are aligned ferromagnetically. In ref. jackeli , the electron-lattice interaction is absent, and it is again the anisotropic $`e_g`$ transfer amplitude of the two-orbital model which drives the CE state. In both cases, however, the AF coupling $`J0.1t`$ necessary to achieve the CE state is one order of magnitude higher than what is estimated from experiments, signaling that there must be some additional effect contributing to the CE ordering. Finally, we reiterate that self-consistent screening is necessary to explain that phases with marked charge order come in very close competition with metallic phases. We believe that this is a crucial feature of doped manganites, that further models addressing coexistence and texturing of those phases at small scales must take into account.


text:
# The Principle of Relativity and Modified Connection between the Mass of a Body and its Velocity: a Case of 𝑉>𝐶.

## Abstract

Following the basic idea expressed in , we assume that for any particle or body with mass $`M`$ its own time $`t`$ depends on therelative change $`\frac{\mathrm{\Delta }M}{M}`$ of that mass. Based on this assumption, one discusses possible existence of a hypothetical particles and bodies with a velocity $`V>C`$

The main idea of paper  is based on the assumption that the flow of time $`t`$ which is associated with a given particle or body with mass $`M`$ is a function of the relative change$`\left(\frac{\mathrm{\Delta }M}{M}\right)`$ of its mass:

$$dt=\mu ^1\frac{\mathrm{\Delta }M}{M}$$
(1)
The physical meaning of $`\mu `$is unknown: one does not know whether this parameter is a function of coordinates, or a velocity $`V`$, or electrical charge $`Z`$ etc. However, in

it was noted that in spite of quite rough estimations of values of$`\mu _n`$, $`\mu _w`$ and $`\mu _g`$ , a comparision of these parameters with dimensionless constants Cs , Cw and Cg for strong, weak and gravitational interactions shows that both curves on Fig.1 of the paper  are similar. That important observation was a serious reason to try-as much as possible-to clarify physical meaning of $`\mu `$-parameters. Besides that, below one discusses these conditions, which provide the velocity of a particle or body to be of $`V>C`$. To get more information about the $`\mu `$-parameters, let us assume that at the front of moving particle a layer of matter with higher density arises. Since one does not know neither the nature of assumed phenomenon nor any details of it, a simple model consisting of a potential barrier with a height of $`E_b`$ and a thickness of $`\mathrm{\Delta }x`$can be applied. This allows the value of$`E_b`$ to be roughly estimated. The probability $`D`$such barrier to be passed is given by well known equation:

$$D\mathrm{exp}\left(\frac{2}{h}\sqrt{2M(E_bE)}\mathrm{\Delta }x\right)$$
(2)
Here $`E`$ is the energy of moving particle, which was defined in  as follow:

$$E=\frac{\mu ^2}{2}(\mathrm{\Delta }S)^2\left(\frac{\mathrm{\Delta }M}{M}\right)^2$$
(3)
If $`E=E_b`$ one gets $`D=1`$ and the height $`E_b`$ of the barrier can be roughly

estimated. By using the values of $`\mathrm{\Delta }S10^{13}cm`$ and of $`\mu _n`$10$`{}_{}{}^{23}\mathrm{sec}_{}^{1}`$

 one finds that $`E_b10^{20}`$ $`M\left(\frac{\mathrm{\Delta }M}{M}\right)^2`$. Assuming that $`\left(\frac{\mathrm{\Delta }M}{M}\right)1`$,

one can get an estimation of the minimum heigt:$`E_b^{\mathrm{min}}`$ $`10^{20}`$ $`M`$. The same

approach can be used to estimate the value of $`E_b(\beta )`$ for $`\beta `$ -decay of a neutron. Using a value of$`\mu _w10^6\mathrm{sec}^1`$  and of $`\mathrm{\Delta }S10^{13}`$ $`cm`$ one obtains that $`E_b(\beta )10^{38}`$ $`M\left(\frac{\mathrm{\Delta }M}{M}\right)^2`$.

As a next step one needs to define a value of $`M`$ for the electron

antineutrino, which is unknown. To make our estimation more conservative one

can assume that this mass lays in a quite wide range of $`1eV<M<100eV`$. In addition, the value of $`E_b(\beta )`$ has to be comparable with the difference between the neutron and proton masses, which is of $`1,3MeV`$. This means that the value of $`\left(\frac{\mathrm{\Delta }M}{M}\right)^2`$. should be of $`10^{42}`$ $`10^{44}`$. Thus, the result leads us to unexpected conclusion: since $`\left(\frac{\mathrm{\Delta }M}{M}\right)_\nu 10^{21}`$ $`10^{22}`$, one should accept that the life-time of the electron antineutrino could be $``$ $`10^4`$times larger than the life-time of Universe!

It is interesting to note that if $`V<C`$, where $`C`$\- is the velocityof light, then the value of$`\mu `$ is limited :

$$\mu <\frac{\mathrm{\Delta }M}{M}C(\mathrm{\Delta }S)^1$$
(4)
On the other hand, since in  there are no any limitations the value of $`V`$to be larger than $`C`$, i.e.$`V>C`$, one can try to determine the conditions that provide this case. The$`\mu `$-parameter has to satisfy the following equation:

$$\mu >\frac{\mathrm{\Delta }M}{M}C(\mathrm{\Delta }S)^1$$
(5)
It is curios the values of $`\mu _n(V>C)`$ and $`\mu _w(V>C)`$ to be estimated.

One assumes that for strong interaction$`\mathrm{\Delta }S10^{13}cm`$ and $`\frac{\mathrm{\Delta }M}{M}1`$. The minimum estimated value of$`\mu _n(V>C)10^{23}\mathrm{sec}^1`$. Comparing both values of $`\mu _n(V<C)`$ and $`\mu _n(V>C)`$one can find that the difference between them is less than one order of magnitude.Similarly, one can get an estimation of the value of $`\mu _w(V>C)(1010^2)\mathrm{sec}^1`$ for weak interaction:. Thus, we find the following ratios:

$$\frac{\mu _n(V>C)}{\mu _n(V<C)}1$$
(6)
$$\frac{\mu _w(V>C)}{\mu _w(V<C)}10^8$$
(7)
Considering these quite curious ratios one can express some very interesting, but totally speculative thoughts. Indeed, since the values of both$`\mu _n`$-parameters are almost the same, is it possible that in some nuclear reactions an extremely small number of particles would penetrate trough such high potential barrier via a tunneling and would have $`V>C`$?

On the other hand, since the above assumed potential barrier for $`\beta `$ -decay is much lower compare to strong interacting particle reactions, one should expect higher yield of neutrinos with $`V>C`$.

References

. Elmir Dermendjiev, ” A Modified connection between the mass and time in the principle of relativity”, preprint nucl-th 0011028 (08. November 2000)