text:
# Numerical Study of Compressible Magnetohydrodynamic Turbulence in Two Dimensions

## 1 Introduction

The existence of cosmic magnetic fields in diffuse astrophysical plasmas (the interstellar media in galaxies, the intracluster media in clusters of galaxies and even the media associated with filaments and superclusters of galaxies) have been recognized for a while (for reviews, see e.g., Kronberg 1994; Beck et al. 1996; Zweibel & Heiles 1997; Carilli & Taylor 2002; also see Ryu et al. 1998; Clarke et al. 2001). Although the origin of such fields is not yet fully understood, turbulence is known to play an important role in the amplification and diffusion of existing magnetic fields. For instance, magnetic fields grow effectively by turbulent motion of conducting fluids (the $`\alpha `$ effect), but they are not further amplified when Maxwell stresses become strong enough to affect the turbulent motion itself. How strong the magnetic fields would be in order to influence turbulence is one of the intriguing topics in magnetohydrodynamic (MHD) turbulence. It has been argued through two-dimensional incompressible simulations that turbulent transport is reduced by weak magnetic fields whose energy is small compared to the kinetic energy of turbulent flows (Cattaneo & Vainshtein 1991; Cattaneo 1994). A similar suppression was also observed in three-dimensional simulations of incompressible flows (Tao et al. 1993). It is an important issue in astrophysics, since it is a part of the process for the generation and evolution of cosmic magnetic fields (see e.g., Ruzmaikin et al. 1988; Kulsrud et al. 1997; Kulsrud 1999).

In this work, we study compressible MHD turbulence by solving the ideal MHD equations with an isothermal equation of state. Since turbulent motion produces structures spanning a wide range of scales with accompanying energy transfer among different scales, high spatial resolution is required to cover a sufficient inertial range. Hence, in this paper we present high resolution two-dimensional simulations using up to $`1536^2`$ grid cells, leaving three-dimensional simulations for follow-up work. In addition, in order to achieve the highest possible magnetic Reynolds number and Reynolds number, no dissipation was included explicitly in our simulations. However, the resistivity and viscosity of numerical origin are still effective as small-scale dissipative channels. From the simulations, we examine the basic properties of the resulting MHD turbulence. We also address the dependence of the properties on resolution, including the amplification of magnetic field.

The rest of the paper is organized as follows. In the next section, we describe the problem setup as well as the code properties. In §3 and §4, we present the results of simulations of hydrodynamic and MHD turbulence. Finally in §5, the findings of this study are summarized.

## 2 Numerics

In simulations of compressible MHD turbulence, gas is heated by shocks and reconnection events, as well as numerical dissipation that mimics viscous and resistive influences (see below). Hence, in order to maintain the turbulence statistically in a steady-state, cooling should be applied to the internal energy. One simple way to handle it is to assume the isothermality of flows. The MHD equations of compressible, isothermal gas are

$$\frac{\rho }{t}+\stackrel{}{}(\rho \stackrel{}{v})=0,$$
(1)
$$\frac{\stackrel{}{v}}{t}+\stackrel{}{v}\stackrel{}{}\stackrel{}{v}+\frac{c_s^2}{\rho }\stackrel{}{}\rho \frac{1}{\rho }(\stackrel{}{}\times \stackrel{}{B})\times \stackrel{}{B}=\stackrel{}{f},$$
(2)
$$\frac{\stackrel{}{B}}{t}\stackrel{}{}\times (\stackrel{}{v}\times \stackrel{}{B})=0,$$
(3)
with an additional constraint

$$\stackrel{}{}\stackrel{}{B}=0,$$
(4)
for the absence of magnetic monopoles. Here, $`c_s`$ is the isothermal sound speed, and the units are chosen so that the factor of $`4\pi `$ does not appear in the equations. To enforce turbulence, a random force per mass, $`\stackrel{}{f}=f_x\widehat{x}+f_y\widehat{y}`$, was added. It has the following form:

$$f_{x,y}(x,y,t)=v_{\mathrm{amp}}\mathrm{cos}(\omega t+\delta _t)\mathrm{cos}(k_{\mathrm{𝑖𝑛},x}x+\delta _x)\mathrm{cos}(k_{\mathrm{𝑖𝑛},y}y+\delta _y),$$
(5)
where $`\delta _t`$, $`\delta _x`$, and $`\delta _y`$ are random phases in the range of $`0\delta _t,\delta _x,\delta _y\pi `$.

The above equations were solved using a multi-dimensional MHD code described in Kim et al. (1999), which is specifically designed for isothermal MHD. It is based on the explicit, finite-difference Total Variation Diminishing (TVD) scheme, which is a second-order accurate upwind scheme, and employs the minmod flux limiter. Simulations were performed in the computational domain of $`x=[0,L]`$, $`y=[0,L]`$, and $`L=1`$ with a periodic boundary condition using $`n_x^2=256^2`$, $`512^2`$, $`1024^2`$ and $`1536^2`$ grid cells. The total mass in the computational box is conserved, and the averaged density was set to be $`\rho =\rho _0=1`$. The values of other parameters in the simulations are the following: the isothermal sound speed $`c_s=1`$, the angular frequency of the random forcing $`\omega =2\pi `$, and the input wavenumbers $`k_{\mathrm{𝑖𝑛},x,y}`$ = $`8\pi `$. With this choice, we note that the input scale of the random forcing is one quarter of the computational box size and the period is one sound wave crossing time across the box. The amplitude of the random forcing, $`v_{amp}`$, was set so that without a magnetic field, hydrodynamic turbulence is saturated with the averaged Mach number $`M_s_{\mathrm{init}}v_{\mathrm{rms}}/c_s=1`$ (see §3). MHD turbulence was generated by introducing uniform magnetic fields of various strength into the saturated hydrodynamic turbulence (see §4).

Although our simulations do not contain explicit resistivity or viscosity, unavoidable numerical diffusion of magnetic field and momentum across cell produce effective numerical resistivity and viscosity, respectively. Kim et al. (1999) studied the character of numerical dissipation in the code used for this work through the decay of a two-dimensional Alfvén wave. They showed that for waves spanning different numbers of grid cells, $`n`$, the effective magnetic Reynolds number and Reynolds number are proportional to $`n^{1.66}`$, mimicking the “hyper-type” resistivity and viscosity ($`\eta ^4\stackrel{}{B}`$ and $`\mu ^4\stackrel{}{v}`$) rather than the “normal-type” resistivity and viscosity ($`\eta ^2\stackrel{}{B}`$ and $`\mu ^2\stackrel{}{v}`$) that characterize the physical dissipation of “collisional” fluids. They also estimated that the effective Reynolds numbers are larger than several hundreds if $`n8`$, so that the inertial range covers scales spanning 8 grid cells or more. We note that numerical dissipation estimates depend somewhat on the tests used. For instance, the test suggested by Zweibel et al. (2002) would have given somewhat smaller numerical dissipation. Hence, our estimation of the inertial range based on the Alfvén wave decay can be regarded to be rather conservative.

In most numerical studies of astrophysical turbulence based on explicit dissipation, the normal-type resistivity and viscosity have been included explicitly. However, dissipation processes in diffuse astrophysical plasmas are not well understood, and it is unlikely that the dissipation can be modeled by way of simple resistivity and viscosity coefficients. Hence, more careful consideration would be required in the studies of the properties of astrophysical turbulence when details on dissipative scales become important. In this paper we are concerned mostly with the properties of turbulence in the inertial range, which would be less sensitive to the form of dissipation, so we ignore dissipative complications outside this range.

## 3 Hydrodynamic Turbulence

Simulations of hydrodynamic turbulence were first performed for two reasons. One is to generate the initial states for MHD simulations. The other is to observe the change of turbulence properties from the the hydrodynamic case under the influence of magnetic fields. In the absence of magnetic fields, a trans-sonic turbulence with $`M_s_{\mathrm{init}}1`$ developed (see Table 2).

Figure 1 shows a typical resulting density distribution at an epoch $`t=15`$ and the time-averaged power spectrum of kinetic energy in the simulation with $`1536^2`$ grid cells. The existence of weak shocks is evident in the density image. Here, “weak” means the shock Mach number is small; typically $`M_{\mathrm{shock}}2`$. The density compression is $`\delta \rho /\rho _{\mathrm{rms}}0.45`$, when averaged over space and time in all the simulations of all resolutions (see Table 3). In a two-dimensional system, the forced turbulence of “incompressible” flows exhibits the dual energy spectrum known as the Kraichnan spectrum: in the range of $`k<k_{in}`$, the slope follows $`k^{5/3}`$ indicating the direct cascade of energy, while in the range of $`k>k_{in}`$, it follows $`k^3`$ meaning the direct cascade of enstrophy (see e.g., Lesieur 1997 and references therein). In the case of “high compressibility,” where strong shocks are common, the energy spectrum follows the $`k^2`$ slope in the inertial range, known as the Burgers spectrum (see e.g., Lesieur 1997 and references therein). In our results with weak compressibility, the energy spectrum is consistent with the Kraichnan spectrum as shown in the right panel of Figure 1, but the slope for $`k>k_{in}`$ is $`2.72`$; somewhat shallower than for incompressible turbulence. As listed in Table 4, the deviation of the slope from $`k^3`$ increases systematically as the resolution increases. This is because in higher resolution simulations, shocks are better resolved, and, hence, their effect is more evident in the energy power spectrum.

## 4 MHD Turbulence

Weak, uniform magnetic fields of various strength were added to the fully developed hydrodynamic turbulence. In the $`n_x^2=256^2`$, $`512^2`$ and $`1024^2`$ simulations, the flows of hydrodynamic turbulence at $`t=30`$ were taken as the initial states, while in the $`1536^3`$ simulations, the flow at $`t=15`$ (the right panel of Figure 1) was taken. The strength of the added magnetic fields was set so that the initial Alfvénic Mach number, $`M_A_{\mathrm{init}}v_{\mathrm{rms}}/c_{A,\mathrm{init}}`$, ranged from $`1000`$ to $`10`$ (it corresponds to the plasma $`\beta `$ of $`2\times 10^6`$ to $`200`$). The Alfvénic Mach number represents the ratio of Reynolds stresses to Maxwell stresses, so is a convenient way to characterize magnetic field strength in these flows.

The magnetic field strength grew rapidly, almost exponentially, during the initial transient period, and then saturated. Figure 2 shows the time evolution of rms magnetic field strength and flow velocity. The resulting flows of MHD turbulence after saturation can be naturally categorized into very weak field (VWF), weak field (WF) and strong field (SF) classes (Table 1). In the VWF cases, magnetic field influences are negligible, and the flow velocity does not change noticeably from purely hydrodynamic turbulence. In the WF cases, the magnetic field energy, $`E_{mag}`$, is still small compared to the flow kinetic energy, $`E_{kin}`$, but the magnetic field does locally affect small scale flow motions. Hence, the turbulence flow velocity is decreased from that of non-magnetic cases. Finally, the cases where $`E_{mag}`$ becomes comparable to $`E_{kin}`$ are classified as the SF class. We note that the change in the flow properties is gradual with the increase of the strength of the initial, uniform fields. Hence, boundaries between the classes are somewhat arbitrary. We set the criterion for the boundary between the VWF and WF classes as $`M_sv_{\mathrm{rms}}/c_s1`$; more specifically, $`M_s=0.9`$ was taken. On the other hand, the criterion for the boundary between the WF and SF classes was set as $`M_Av_{\mathrm{rms}}/c_A_{\mathrm{rms}}1`$, more specifically $`M_A=2`$ was taken. Among the SF cases, we observed examples where there a change in the magnetic field configuration takes place. In particular there are cases in which a dominant flux tube develops. These cases we categorized specifically as the field separation (FS) class. The growth of magnetic field strength and the characteristics of the four classes are further detailed in the following subsections.

For quantitative discussions we computed, in addition to $`M_A`$ and $`M_s`$ (Table 2), the density compression, $`\delta \rho /\rho _{\mathrm{rms}}`$, and the intermittency (or the kurtosis of the field distribution) $`IB^4_{\mathrm{rms}}/B^2_{\mathrm{rms}}^2`$ (Table 3), as well as the slopes of the power spectra of kinetic and total energies (Table 4). Note that with a Gaussian distribution of magnetic field, $`\stackrel{}{B}\mathrm{exp}(x^2/2\sigma ^2)\widehat{y}`$, the value of $`I`$ is 2.39 when summed over the interval $`3\sigma x3\sigma `$, or 3.99 when $`5\sigma x5\sigma `$ is taken into account.

### 4.1 Growth of Magnetic Field Strength

There is no dynamo action in two-dimensions, and the total magnetic flux through given boundaries is conserved. Hence, the growth of magnetic field strength shown in Figure 2 is due to stretching and compression. The growth stops either when the magnetic energy reaches an equipartition with the kinetic energy, so that the back-reaction from magnetic field Maxwell stresses plays a dominant role, or when the separation between magnetic sheets is reduced to the diffusive scale, so that reconnection takes place (see e.g., Biskamp 1993 for details). If the initial, uniform magnetic field is weak enough, the growth is saturated by the latter cause before the magnetic energy reaches equipartition, as in the VWF and WF cases. Then, the resulting growth should depend on the diffusive scale and, hence, on the effective resistivity. Then the amplification of magnetic energy from the initial value is expected to be $`(k_{\mathrm{𝑖𝑛}}\delta )^2R_m`$ (see e.g., Biskamp 1993). Here, $`\delta `$ is the diffusive scale and $`R_m`$ is the magnetic Reynolds number.

Figure 3 shows the averaged amplification of magnetic energy in our simulations as a function of grid resolution. Only the VWF and WF cases were considered, because in the SF cases, the growth of magnetic field stopped as a consequence of Maxwell stresses, rather than reconnection. The amplification shown in the figure is proportional to $`n_x^{1.32}`$. Furthermore, with $`R_mn^{1.66}`$ for our code (see §2), the amplification is estimated to be approximately proportional to $`R_{m}^{}{}_{}{}^{0.8}`$, which is somewhat shallower than the expected dependence ($`R_m`$). However, this result agrees well with results of two-dimensional numerical simulations for incompressible MHD turbulence (Biskamp 1993).

### 4.2 VWF Cases

The VWF class includes those cases where the back-reaction of magnetic field into flow motions is insignificant. The criterion $`M_s0.9`$ corresponds to the rms velocity decrease by less than 10%. The characteristics of turbulence in these cases are summarized as follows: 1) The density compression is still large with $`\delta \rho /\rho _{\mathrm{rms}}0.38`$ (Table 3) and shocks still exist as shown clearly in the bottom-left panel of Figure 4. 2). The intermittency is large with $`I3.5`$ (Table 3), which indicates that the spatial contrast in the magnetic field strength distribution is high. In particular, the bottom-right panel of Figure 4 shows that the magnetic field is mostly thin tubes (or thin sheets in extension to the third dimension). 3) The Alfvén Mach number is $`M_A4.5`$ (Table 2), and hence the kinetic energy is, at least, an order of magnitude larger than the magnetic energy. Also, the power of kinetic energy, $`P_{k}^{}{}_{}{}^{kin}`$, is larger than that of magnetic energy, $`P_{k}^{}{}_{}{}^{mag}`$, over most wavenumbers, as shown in the bottom panels of Figure 5. However, $`P_{k}^{}{}_{}{}^{mag}P_{k}^{}{}_{}{}^{kin}`$ in a small range of large wavenumbers. This is because the magnetic field first built up on small scales. However, it saturated due to reconnection before the power could extend to larger scales. 4) The slope of $`P_{k}^{}{}_{}{}^{kin}`$ was not affected much by the magnetic field. The change was less than a few percent as shown in Table 4. But the slope of the power of the total energy, $`P_{k}^{}{}_{}{}^{tot}`$, changed significantly compared to the hydrodynamic case, especially in the MHD cases close to the boundary between the VWF and WF classes. This is because the small scale power of magnetic energy pushed $`P_{k}^{}{}_{}{}^{tot}`$ up, and hence $`P_{k}^{}{}_{}{}^{tot}`$ does not have a single, well-defined slope any more.

### 4.3 WF Cases

The WF class was categorized by the property that globally the magnetic energy is smaller than the kinetic energy, but yet the back-reaction of magnetic field is not negligible. As the criterion to represent this property, $`M_s0.9`$ and $`M_A2`$ were adopted (Table 2). The characteristics of turbulence in these cases are the following: 1) As shown in the middle-left panel of Figure 4, the occurrence of shocks is reduced and their strength is weakened compared to the hydrodynamic case. The density compression is $`0.28\delta \rho /\rho _{\mathrm{rms}}0.38`$ (Table 3), which is smaller than that in the VWF cases. 2) The contours of magnetic field lines plotted in the middle-right panel of Figure 4 show the appearance of almost circular magnetic islands. At the same time, the magnetic field is less concentrated than in the VWF cases, with smaller intermittency $`2.2I3.5`$ (Table 3). 3) Although in the small wavenumber region the magnetic energy power is still smaller than the kinetic energy power, $`P_k^{mag}P_{k}^{}{}_{}{}^{kin}`$, at large wavenumbers, $`P_{k}^{}{}_{}{}^{mag}P_{k}^{}{}_{}{}^{kin}`$, as shown in the middle panels of Figure 5. Hence, turbulent transport can be suppressed by magnetic fields on the small scales where $`P_{k}^{}{}_{}{}^{mag}P_{k}^{}{}_{}{}^{kin}`$. Consequently, the flow velocity decreases noticeably from that of hydrodynamic turbulence (see Figure 2). 4) The power of magnetic energy, $`P_{k}^{}{}_{}{}^{mag}`$, changed the shape of $`P_{k}^{}{}_{}{}^{tot}`$ as well as that of $`P_{k}^{}{}_{}{}^{kin}`$ in power spectra, and the slopes are not simply defined over the entire inertial range of wavenumber any more (see Figure 5). We calculated the slopes of $`P_{k}^{}{}_{}{}^{kin}`$ and $`P_{k}^{}{}_{}{}^{tot}`$ over $`2\times k_{in}kk_{24}`$, where the slopes are relatively well defined, and listed them in Table 4. Here, $`k_{24}`$ is the wavenumber corresponding to 24 grid cells. The slope of $`P_{k}^{}{}_{}{}^{kin}`$ is significantly smaller than those of the hydrodynamic and VWF cases. The slope of $`P_{k}^{}{}_{}{}^{tot}`$ is quite small. But this is because $`P_{k}^{}{}_{}{}^{mag}`$ peaks in the wavenumber range over which the slope was calculated (see Figure 5).

### 4.4 SF Cases

In the SF cases, the magnetic energy grows to become comparable to the kinetic energy, with $`M_A2`$ at saturation (Table 2). Hence, the flow velocity is influenced significantly by magnetic field (Figure 2), with $`M_s0.75`$ at saturation (Table 2). The characteristics of SF turbulence are the following: 1) With the significantly reduced Mach number, shocks are rare, as shown in the top-left panel of Figure 4. The density compression is accordingly reduced to $`\delta \rho /\rho _{\mathrm{rms}}0.28`$ (Table 3). 2) As shown in the top-right panel of Figure 4, circular magnetic flux islands, or loops, are common. This trend of changing topology in magnetic field lines, from tubes to loops, with increase of magnetic field strength was observed also in the two-dimensional simulations of incompressible MHD turbulence (see e.g., Biskamp 1993). At the same time, the spatial contrast of magnetic field strength distribution becomes low with smaller intermittency $`I2.2`$ (Table 3). 3) The power of magnetic energy exceeds that of kinetic energy with $`P_k^{mag}P_k^{kin}`$ over all $`kk_{in}`$ as shown in the top panels of Figure 5, although still $`P_k^{mag}<P_k^{kin}`$ in $`kk_{in}`$. 4) The slopes of $`P_{k}^{}{}_{}{}^{kin}`$ and $`P_{k}^{}{}_{}{}^{tot}`$, calculated over $`2\times k_{in}kk_{24}`$ as in the WF cases, are listed in Table 4. Across the above wavenumber range, $`P_{k}^{}{}_{}{}^{kin}`$, $`P_{k}^{}{}_{}{}^{tot}`$ and $`P_{k}^{}{}_{}{}^{mag}`$ all exhibit a single slope close to $`1.5`$; namely, the slope of the Iroshnikov-Kraichnan spectrum (the top panels of Figure 5).

We note that the emergence of the slope of $`1.5`$ is achieved only in the very high resolution simulations with $`1024^2`$ and $`1536^2`$ grid cells, not in the simulations with $`256^2`$ and $`512^2`$ grid cells. In fact, there is no case categorized as the SF class in the simulations with $`256^2`$ and $`512^2`$ grid cells at all (see Table 1). That is, our two-dimensional simulations of MHD turbulence start to show converged behaviors with $`1024^2`$ grid cells or more.

### 4.5 FS Cases

The cases with $`M_A_{\mathrm{init}}30`$ show FS class behavior, wherein a dominant flux tube forms, and, as a consequence, the magnetic field separates from the background flow (see Fig. 6). In such FS cases, flows become anisotropic and their characteristics are different from the other cases. The appearance of this behavior can be understood as follows. Recall that there are three ideal 2D MHD invariants; the total energy, $`E_{tot}`$, the mean square of the magnetic potential, $`A`$, and the cross helicity,

$$K=\frac{1}{2}\stackrel{}{v}\stackrel{}{B}d^2\stackrel{}{x}.$$
(6)
Among them, the power spectrum of $`A`$ exhibits an inverse cascade, while the power spectra of $`E_{tot}`$ and $`K`$ exhibit a normal cascade (see e.g., Biskamp 1993). Hence, just as in three dimensions, large scale magnetic field power can be built up in the regime $`kk_{in}`$ in two-dimensional MHD turbulence. In our FS cases this power build-up proceeds through the formation of a magnetic flux tube with a coherent length larger than the scale associated with $`k_{in}`$. However, as pointed out in §4.1, there is no dynamo action in our simulations. Hence, the formation of such a flux tube should involve reconnection followed by the expulsion of gas out of the flux tube, reducing the gas density (so, also the gas pressure) inside the tube, while maintaining approximate pressure equilibrium.

We note that the critical value of $`M_A_{\mathrm{init}}`$ for the FS class is independent of numerical resolution in our simulations (see Table 1). This is because the existence of the flux tube depends on its ability to resist ram pressure bending, which depends on large scale flows, so is independent of numerical resolution. On the other hand, since the formation of the dominant flux tube involves reconnection, we expect the time to reach the state where the magnetic field separates from the background flow would depend on resistivity. The flux tube was developed in a few tens eddy turnover times in our simulations. However, in astrophysical environments such as interstellar media or intracluster media, the resistivity is much smaller than that of our simulations. The classical magnetic Reynolds number in the interstellar medium can be as large as $`10^{50}`$ or so (see e.g., Spitzer 1979). Hence, we expect that the emergence of the flux tube would take too long in astrophysical environments, and instead, turbulence of the SF class is more likely to develop. So the failure of SF cases in our $`256^2`$ and $`512^2`$ simulations is probably a limitation of numerical simulations; in particular from very large numerical resistivity compared to astrophysical environments and perhaps also from the periodic boundary we employed.

## 5 Summary & Conclusion

We performed high resolution two-dimensional simulations of isothermal MHD turbulence using up to $`1536^2`$ grid cells. Compressibility was taken into account, with Mach numbers $`M_s1`$. In order to maximize the magnetic Reynolds number and the Reynolds number, no explicit dissipation was included. Instead, resistivity and viscosity of numerical origin were utilized. Our findings are summarized as follows:
1) The growth of magnetic energy from the initial value depends on effective resistivity and, hence, on resolution. We measured a maximum magnetic energy amplification proportional to $`n_x^{1.32}`$. After the growth saturated, the turbulence could be categorized into three classes, which we labeled very weak field (VWF), weak field (WF), and strong field (SF) classes, depending on the strength of the uniform component of magnetic fields, or equivalently the initial Alfvén Mach number, $`M_A_{\mathrm{init}}`$. Each class is characterized by different turbulence properties.
2) In the WF cases, although $`E_{mag}<E_{kin}`$ globally, the rms flow velocity decreases due to the magnetic field, which can become locally important. That is, even a weak magnetic field reduces turbulent transport (Cattaneo & Vainshtein 1991; Cattaneo 1994).
3) Although the inertial range power spectra of MHD turbulence are not represented by power-laws any more, part of the inertial range still develops a slope close to $`1.5`$; that of the Iroshnikov-Kraichnan spectrum, in the SF cases. Such a slope developed only in highest resolution simulations with $`1024^2`$ and $`1536^2`$ grid cells.
4) In the lower resolution simulations with $`256^2`$ and $`512^2`$ grid cells, there are no cases categorized as the SF class. Instead, with $`M_A_{\mathrm{init}}30`$, those simulations developed into FS cases, where a large flux tube dominates the flow structures and separates the flow from the magnetic field. This FS property is the result of large numerical resistivity along with the periodic box used.

The conclusions of our work are the following:
1) For the VWF, WF, and SF cases, we observed properties of turbulence which are consistent with those of incompressible counterparts, whenever comparisons are made. Hence, we conclude that weak compressibility of $`\delta \rho /\rho _{\mathrm{rms}}`$ up to $`0.45`$ would not be important in characterizing MHD turbulence.
2) Converged behavior in simulated turbulence, such as that of the SF class, starts to appear only in simulations with very high resolution, $`1024^2`$ grid cells or more, in our simulations. The fact that simulations with $`1024^2`$ grid cells or more are necessary for our two-dimensional study of MHD turbulence points to the need of very high resolution in three-dimensional MHD turbulence studies as well. For hydrodynamic turbulence, Porter et al. reported that converged results emerged with $`512^3`$ or more grid cells in their simulations. However, we note that numerical dissipation would differ in different codes based on different schemes. As a result, the convergence behavior would be different, too. These findings emphasize the importance in turbulence studies of high resolution grids and codes with the smallest possible numerical dissipation. Any explicit treatment of dissipation will require even higher resolution simulations, since numerical dissipation and diffusion would have to be smaller than the explicit terms.

The work by HL and DR was supported by an international collaboration grant from Korea (KOSEF F01-2000-000-10009-0) and Germany (DFG). The work by JK was supported by Strategic National R&D Program (M10222000005-02B0600-00400) from MOST, Korea. The work by TWJ was supported by the National Science Foundation (AST96-19438; AST00-71167) and the University of Minnesota Supercomputing Institute. Simulations were made through the support by “The Supercomputing Application Support Program of KISTI”.


format: text

text:
# High-resolution x-ray study of the nematic - smectic-A and smectic-A - smectic-C transitions in 8̄𝑆5-aerosil gels.

## I Introduction

The phase transition behavior of many thermotropic liquid crystals is quite well characterized deGennes\_Prost  making liquid crystals ideal materials for studying the effects of weak perturbations, including especially quenched random fields and interactions. In a recent series of papers Park ; Leheny ; Germano2 ; Clegg , it has been shown that the nematic to smectic-A (N - SmA) transition in liquid crystal-aerosil dispersions is well described as a continuous symmetry (XY) transition with quenched random fields. This facilitates studying this important class of transitions. Early calorimetric studies have shown a systematic variation in the critical fluctuations as the aerosil concentration increases. These experiments and accompanying theory have focused on the isotropic to nematic (I - N) Germano1  and the N - SmA transitions Park ; Leheny ; Germano2 ; Clegg ; Germano1 . Here we present results on a substantially different material at the N - SmA transition and we present the first x-ray data for the smectic-A to smectic-C (SmA - SmC) transition for liquid crystal-aerosil dispersions. The SmA - SmC transition for heptyloxybenzylidene(7O.4)-aerosil dispersions has previously been investigated using ac-calorimetry H\_and\_G .

4-$`n`$-pentylphenylthiol-4’-$`n`$-octyloxybenzoate ($`\overline{8}S5`$) is a non-polar liquid crystal with a succession of ordered phases (Fig. 1). At $`T_{NA}^0=336.6`$ K a transition from the orientationally ordered N phase to the layered SmA phase occurs. For this particular liquid crystal the pretransitional correlation volume is highly anisotropic Safinya81 . As we shall demonstrate in the present work, this gives rise to a strong asymmetry in the powder-averaged line shape observed with x-ray diffraction. This effect is enhanced by the relatively small fourth-order correction to the correlation function in the transverse direction. Thus, the spherically averaged smectic line shape is very different from that observed for materials such as octylcyanobiphenyl (8CB) and octyloxycyanobiphenyl (8OCB). Consequently, this $`\overline{8}S5`$ study provides a test of the generality of the line-shape analysis used in Refs. Leheny  and Clegg  versus that used in Ref. Aerogel\_Science . As the temperature is lowered further, $`\overline{8}S5`$ exhibits a SmA - SmC transition ($`T_{AC}^0=329.3`$ K) where the molecules develop a tilt with respect to the layer normal. The tilt angle is the order parameter for this phase transition. The transitional behavior has mean-field rather than critical characteristics albeit with an anomalously large sixth order term in the free energy Birgeneau83 , which is typical for the Landau tricritical to mean-field crossovers observed for most SmA - SmC transitions M\_and\_G .

A liquid crystal is perturbed with weak-to-intermediate strength quenched disorder when aerosil particles are introduced. These particles are $`7`$nm diameter silica spheres with hydrophilic surfaces. When aerosil is dispersed in a liquid crystal, hydrogen bonds form between the aerosil particles creating a low density gel. On very long length scales the gel structure exhibits fractal correlations Germano1 ; however, the gel has no particular correlations commensurate with the wave vector of the SmA structure and hence is effectively random for our purposes. Here we present high-resolution x-ray diffraction measurements for $`\overline{8}S5`$ samples containing a range of concentrations of aerosil particles. Both the SmA and SmC type liquid crystal correlations in these samples have been characterized as a function of temperature.

Previous x-ray studies of the N - SmA transition in liquid crystal-aerosil dispersions demonstrated that the transition responded to the aerosil as though its principal effect was to apply a random field that pinned the phase of the density wave Park ; Leheny ; Clegg . The static random-field fluctuations are more effective than thermal fluctuations at destroying long-range smectic order, and this was evident in the development of the structure factor as a function of temperature and aerosil concentration. Here we wish to investigate whether the non-polar $`\overline{8}S5`$ is less strongly perturbed than the polar cyano materials 8CB and 8OCB were. In addition, the tilt order parameter at the mean-field SmA - SmC transition couples to the aerosil in a different manner than the SmA order parameter.

Section II describes the experimental techniques used to make and study the samples. In Section III the analysis protocol is outlined and the data and model parameters are presented. Finally in Section IV we draw conclusions from our measurements.

## II Experimental techniques

The $`\overline{8}S5`$ was synthesized at Kent State University and the transition temperatures were observed to be close to those previously recorded for the bulk material Schantz ; Neubert81 . The aerosil particles (hydrophilic type-300) were provided by Degussa and were dried prior to sample preparation. The aerosil concentration was characterized by $`\rho _S`$, the mass of aerosil divided by the volume of liquid crystal. Samples with $`\rho _S=0.0,0.025,0.039,0.073,0.096,0.220,`$ and $`0.350`$ g cm-3 were prepared (we drop the units hereafter). The corresponding pore volume fractions and mean pore sizes are given in Table 1. As with previous studies of 8OCB-aerosil dispersions Clegg , it was necessary to avoid crystallization of the $`\overline{8}S5`$ samples prior to measurement since this results in phase separation. The aerosil particles were dispersed in a solution of $`\overline{8}S5`$ in spectroscopic grade acetone; the system was then sonicated and kept well above the crystallization temperature until the measurements were performed. The x-ray studies were undertaken once the acetone solvent had been evaporated. The sample environment and experimental geometry have been described in detail previously Clegg . The studies were made using 8 keV x-rays at the X20A beam line at the National Synchrotron Light Source. Of great interest here was the asymmetry of the smectic scattering peak and its evolution with temperature and disorder strength. To insure that the observed effects were due to the sample and not the instrumental resolution, measurements were made at both positive and negative scattering angles. As with previous studies, the number of scans was kept limited to avoid undue radiation damage to the sample. The results obtained in the N, SmA, and SmC regimes are reported below.

## III Results and analysis

Figure 2 shows a succession of scans of $`\overline{8}S5`$ with $`\rho _S=0.096`$ as the temperature is decreased. At high temperatures in the nematic phase (Fig. 2(a)) there is a broad peak corresponding to short range pretransitional SmA fluctuations. The observed powder average has a characteristic asymmetry. On reducing the temperature the correlated regions grow in size as indicated by the narrower line profile (Fig. 2(b)). Unlike bulk $`\overline{8}S5`$, the liquid crystal enters a SmA state with a finite correlation volume. Further cooling takes the sample toward the SmC state. Prior to the onset of the molecular tilt, the peak begins to broaden (Fig. 2(c)). As the sample enters the SmC regime the peak moves to higher wave vectors and begins to narrow (Fig. 2(d)). At the lowest temperatures studied, the angle of molecular tilt exceeds $`20^{}`$ while the line-shape asymmetry is reversed compared to that in the SmA phase (Fig. 2(e)). The solid lines in Fig. 2 are the results of fits to a model to be described below.

In keeping with previous studies of the N-SmA transition in liquid crystal - aerosil dispersions Leheny ; Clegg , the line shape has been taken to be composed of two terms: the first is based on the line shape for the SmA thermal fluctuations in the pure liquid crystal and the second is the same line shape squared, which corresponds to the contribution of static random-field fluctuations. The latter are caused by the randomly positioned aerosil surfaces pinning the phase of the density wave. The results of fits in Fig. 2(a-c) and in Fig. 3 show that this line shape describes the data in the SmA region very well. The full expression is

$`S_{NA}(𝐪)=`$ $`{\displaystyle \frac{\sigma _1}{1+\xi _{}^2(q_{}q_{}^0)^2+\xi _{}^2q_{}^2+c\xi _{}^4q_{}^4}}`$
$`+`$ $`{\displaystyle \frac{\sigma _2}{[1+\xi _{}^2(q_{}q_{}^0)^2+\xi _{}^2q_{}^2+c\xi _{}^4q_{}^4]^2}}`$ (1)
which is powder averaged and convoluted with the measured resolution function. In Eq. (1), $`\sigma _1`$ and $`\sigma _2`$ are the amplitudes of the thermal and random-field terms respectively; $`\xi _{}`$ and $`\xi _{}`$ are the correlation lengths parallel and perpendicular to the layer normal; likewise $`q_{}`$ and $`q_{}`$ are the components of the wave vector $`𝐪`$; $`q_{}^0`$ is the peak position and $`c`$ is the coefficient of the fourth order correction term. The amplitude of the random-field term increases as the temperature is reduced (Fig. 3).

As with previous studies Park ; Leheny ; Clegg , it is assumed that the dependence of $`\xi _{}`$ and $`c`$ on $`\xi _{}`$ remains the same in the aerosil samples as it is in the pure materials. It should be noted that the asymmetric line shape is a consequence of powder averaging Eq. (1) for a highly anisotropic correlation volume together with a small value for $`c`$ in $`\overline{8}S5`$ which, in turn, means that there are long tails in the transverse direction. This $`\overline{8}S5`$ line shape cannot be reproduced using the approximation scheme described in Ref. Aerogel\_Science  which does not incorporate explicitly the bulk liquid crystal structure factor. Figure 4 shows the parallel correlation length, $`\xi _{}`$, values extracted from the fits. At high temperatures the short-ranged smectic correlations develop in the nematic phase in a similar manner to that in pure $`\overline{8}S5`$ and to other liquid crystal-aerosil gels. At temperatures below the bulk N - SmA transition temperature the behavior differs from that of the pure material. Estimates of the correlation lengths parallel and perpendicular to the layer normal close to the N - SmA pseudotransition res\_note  are listed for each sample in Table 1. The parallel correlation lengths are very much greater than those for 8CB-aerosil Park ; Leheny  and 8OCB-aerosil Clegg  at the same $`\rho _S`$ values. These correlation lengths are very long compared to the size of the pores, more than ten times longer for most samples. The highly anisotropic correlation volume is evident from the perpendicular correlation length values. These values are comparatively short, roughly half the mean pore size (Table 1). At lower temperatures the behavior in the SmA regime is qualitatively very different from that observed for 8CB-aerosil and 8OCB-aerosil systems Leheny ; Clegg . The size of the correlated SmA regions in a given $`\overline{8}S5`$-aerosil sample decrease markedly as the temperature is lowered toward the SmA - SmC transition whereas $`\xi _{}`$ is roughly independent of temperature in the SmA region for the materials with no SmC phase Leheny ; Clegg . Close to the SmA-SmC transition there is a tilt instability which can be induced by compression and this may make the liquid crystal less resistant to the disordered environment Ribotta74 .

Figure 5 illustrates the low temperature line shape for the finite range SmC ordered state when $`\rho _S=0.073`$. The asymmetry is pronounced and is different from that observed when the structure factor has SmA character. The solid line in Fig. 5(a) is the result of a fit of the data to the model given in Eq. (1). This model does not describe the line shape at all well. The solid line in Fig. 5(b) is the result of a fit using the SmC structure factor as the basis for the two component line shape. The SmC structure factor is taken from studies of the pretransitional N-SmC fluctuations in $`\overline{7}S5`$-$`\overline{8}S5`$ mixtures, the so called Chen-Lubensky structure factor MM87 :

$`S_{AC}(𝐪)={\displaystyle \frac{\sigma _1}{1+\xi _{}^2(q_{}q_{}^0)^2U_pq_{}^2+V_pq_{}^4}}`$
$`+{\displaystyle \frac{\sigma _2}{[1+\xi _{}^2(q_{}q_{}^0)^2U_pq_{}^2+V_pq_{}^4]^2}}`$ (2)
with $`U_p,V_p>0`$. This is identical to Eq. (1) except that the coefficient of the $`q_{}^2`$ term is assumed to be negative. Since the line shape and the peak position have quite different responses to the porous environment it is more convenient to restate this structure factor displaying explicitly the component of the peak position perpendicular to the long-axis of the molecules, $`q_{}^0`$:

$`S_{AC}(𝐪)={\displaystyle \frac{\sigma _1}{1+\xi _{}^2(q_{}q_{}^0)^2+D_p^2(q_{}^2(q_{}^0)^2)^2}}`$
$`+{\displaystyle \frac{\sigma _2}{[1+\xi _{}^2(q_{}q_{}^0)^2+D_p^2(q_{}^2(q_{}^0)^2)^2]^2}}`$ (3)
This is seen to give a good account of the data in Fig. 5 and Fig. 2(d & e). The best fit between this structure factor and the data was found by varying the thermal amplitude $`\sigma _1`$, the random-field amplitude $`\sigma _2`$, the parallel correlation length $`\xi _{}`$, the perpendicular parameter $`D_p`$, and the peak position $`q_c=\sqrt{((q_{}^0)^2+(q_{}^0)^2)}=2\pi /d_c`$. The two components, $`q_{}^0=q_ccos\varphi `$ and $`q_{}^0=q_csin\varphi `$ were determined using the previously determined relationship Safinya80  between the tilt angle, $`\varphi `$, and the layer spacing in the SmA phase, $`d_a`$, that is $`\varphi =1.2\sqrt{(2(1d_c/d_a))}`$.

The SmA-SmC transition in $`\overline{8}S5`$ is a mean-field transition with the tilt angle as the order parameter Safinya80 . Hence the pretransitional fluctuations are weak in the pure material. Figure 6 shows the tilt angle, as a function of temperature, for the full range of disorder strength studied. The tilt angle was determined from the SmC peak position, $`q_c`$. The effective transition temperature, $`T_{AC}^{}`$, was found by fitting the Landau-Ginzburg tricritical model Birgeneau83  to the temperature dependence of the tilt angle. When plotted versus $`TT_{AC}^{}`$, the tilt angles for all samples collapse onto a single curve (Fig. 6). The robustness of $`\varphi `$ in gels suggests that, while the disorder may be coupling to the position and presumably the orientation of the smectic layers, it is not strongly affecting the molecular orientations parms\_note . Figure 7 shows parameters extracted from the fits of Eq. (3) to the data acquired below $`T_{AC}^{}`$. Only results for low $`\rho _S`$ are presented. In Fig. 7(a) the correlation volume, $`\xi _{}D_p`$, is shown as a function of temperature in the SmC region. Increasing the aerosil concentration progressively decreases the extent of the SmC correlations. For higher $`\rho _S`$ the line shape becomes distorted as will be described below. In Fig. 7(b) the ratio of the random field to thermal amplitudes is presented in the form $`a_2/\sigma _1`$, where $`a_2=\sigma _2/\xi _{}D_p`$. The ratio $`a_2/\sigma _1`$ is independent of normalization and is predicted to have a power law dependence on $`\rho _S`$ Germano2 . Increasing aerosil concentration gives rise to a larger random-field contribution to the structure factor. Note that both Figs. 7(a) and 7(b) show $`\rho _S`$ dependence for temperatures just below $`T_{AC}^{}`$ that diminishes at lower temperatures. The low temperature behavior could be an early symptom of the distortion that becomes clear for the higher density samples.

Table 1 lists the SmA-SmC transition temperatures as a function of $`\rho _S`$. There is an initial sharp drop followed by a recovery. This effect is not understood, but qualitatively similar variations have been observed at N-SmA and N-I transitions in other liquid crystal-aerosil gels Germano1 .

For the highest densities of aerosils and the lowest temperatures the peak shape corresponding to SmC fluctuations becomes distorted. A shoulder develops on the low tilt angle side of the reflection peak, as shown in Fig. 8(a). The shoulder becomes more pronounced as the temperature is reduced and as $`\rho _S`$ is increased. Figure 8(b) shows the result of subtracting the result of a fit (solid line Fig. 8(a)) from the data. A peak is observed which falls between the SmA peak position $`q_a=q_{}^0`$ and the SmC peak position $`q_c`$. This behavior indicates that a minority of regions become fixed at lower tilt angles.

## IV Conclusions

X-ray diffraction results have been presented for $`\overline{8}S5`$ aerosil dispersions over the temperature range for which the bulk liquid crystal has N-SmA and SmA-SmC transitions. The gel of aerosil particles is observed to destroy the SmA and SmC phases in the sense that the ordered states becomes short ranged. Both the short-range order SmA and SmC structure factors are well described by a model assuming random fields. The correct thermal fluctuation structure factor has to be used as the thermal term and its square as the random-field term. Estimates at low $`\rho _S`$ and more accurate values for high $`\rho _S`$ res\_note  indicate that $`\xi _{}`$, the correlation lengths parallel to the layer normal, are very long compared to those of 8CB and 8OCB for the same pore sizes. This suggests that the non-polar molecule experiences a weaker pinning field than the polar materials do. The highly anisotropic correlation volume for this material makes the fitting procedure simpler. It also makes clear that the two-Lorentzian model used by Bellini and coworkers Aerogel\_Science  cannot be applied in this case. That the two component line shape gives a good account of data for the SmA - SmC transition in $`\overline{8}S5`$-aerosil gels is intriguing.

For other liquid crystal-aerosil gels exhibiting N - SmA transitions, it has been possible to make detailed comparisons with theoretical models. For 8CB-aerosil Park ; Leheny ; Germano2  and 8OCB-aerosil Clegg , the parameters extracted from the line shape compared favorably with the predictions of Ref. A\_and\_P , giving the appearance that the gel applies a random field which pins the phase of the smectic density wave. The low temperature behavior of the correlation length, due to enhanced anomalous elasticity, predicted by Ref. R\_and\_T  has not been observed for LC-aerosil gels. In $`\overline{8}S5`$, the SmA regime appears to be strongly affected by the impending SmA - SmC transition. The correlation length is never constant, as for a random-fields model, or, increasing on cooling as predicted for anomalous elasticity. Currently there are no theoretical predictions for the effects of a porous random environment on the SmA - SmC transition.

The SmA-SmC tilt (order parameter) behavior is preserved almost unmodified in the aerosil samples. This suggests that any effect that the aerosils may have on the molecular orientations has little consequence for the transition and that this is true even for quite high aerosil concentrations. The bare correlation volume for SmA materials exhibiting a SmC phase can be estimated from the elastic properties of the material and is typically (70 Å)3 Safinya80 . For the aspect ratio of the $`\overline{8}`$S5 molecule this is likely to imply a bare correlation length parallel to the molecule of close to ten molecular lengths. This long-range of interaction is highly likely to minimize the influence of random disorder. If some small fraction of the molecules are forced to point in uncorrelated random directions they will tend to have no net effect: the large bare correlation length implies an averaging over the disorder. For high $`\rho _S`$, the appearance of correlations peaked at a wave vector mid-way between those of the SmA and SmC phases indicates that small regions of the sample become stuck at a lower tilt angle.

## V Acknowledgments

We are grateful to S. LaMarra for technical assistance, to N.A. Clark for discussing unpublished research, and to M. Ramazanoglu for helpful comments. Funding in Toronto was provided by the Natural Science and Engineering Research Council and support was provided at WPI by the NSF under award DMR-0092786 and at JHU under award DMR-0134377. Dr Mary E. Neubert and David G. Abdallah Jr. were supported by NSF grant DMR 89-20147.


text:
# Remarkable long-range-systematic in the binding energies of 𝛼-nuclei II

## 1 Introduction

Twenty years ago we published a remarkable systematic of nuclear binding energies . The two-nucleon separation energies of $`\alpha `$-nuclei where found to be approximately constant and equal to the 2-proton, 2-neutron and PN-separation energies of the $`{}_{}{}^{4}He`$-nucleus. At these days this systematics could only be followed up to $`{}_{}{}^{72}Kr`$. In the past 20 years the table of nuclear binding energies was considerably extended111Nuclear Structure and Decay Data, National Nuclear Data Center, Brookhaven National Laboratory, Upton, NY 11973-5000; http://www.nndc.bnl.gov/nndc/nndcnsdd.html.: It is now interesting how far this systematic can be followed also in these new data. This is indeed well possible as will be shown in this publication.

The nuclear separation energies are defined by the difference of the binding energies of initial nucleus minus the sum of that of the final nuclei. I correct it in the most simple way for the different Coulomb energies according to the trivial formula 2 :

$`S_{2N}^{nuclear}`$ $`=`$ $`(B_i+U_i^C)(B_f+U_f^C)U_{2N}^C`$ (1)
$`U^C`$ $`=`$ $`{\displaystyle \frac{3}{5}}(e^2/r_0)Z(Z1)/A^{1/3},r_0=1.25\text{ fm},`$ (2)
where $`B_i`$ and $`B_f`$ are the experimental binding energies of the initial and final nucleus respectively, $`2N`$ corresponds to any pair of nucleons \[proton-proton ($`\pi ,\pi `$), neutron-proton ($`\nu ,\pi `$), and neutron-neutron ($`\nu ,\nu `$) pair\]. As can be clearly seen in figure 1, the removal of a 2N-pair costs about the same energy as for a bare $`\alpha `$-particle $`30`$MeV. In contrast, the removal of an additional $`\pi \pi `$, $`\nu \nu `$ or $`\pi \nu `$ pair costs much less energy. Moreover, this strongly depends on the size the mother nucleus, resp. the number of $`\alpha `$-particles surrounding this pair.

Figure 2 show the separation energies of an $`\alpha `$-particle out of a $`n_\alpha `$-nucleus. The peaks at $`n_\alpha =4,7,14`$, and the disappearance of peaks above $`n_\alpha =14`$, remind of the special stable $`3`$-dim clusters of closed packed spheres at $`n=4,7,13,55`$. The latter ones with icosaeder structure are observed also for atomic clusters of the heavier noble gases like Xen. For the lighter ones like Ar clusters the shift from $`13`$ to $`14`$ is well known.

## 2 Acknowledgement

I am grateful to W.v.Oertzen for pointing at the new and far more extended table of nuclear binding energies.


text:
# Diamagnetic Persistent Currents and Spontaneous Time-Reversal Symmetry Breaking in Mesoscopic Structures

\[

## Abstract

Recently, new strongly interacting phases have been uncovered in mesoscopic systems with chaotic scattering at the boundaries by two of the present authors and R. Shankar. This analysis is reliable when the dimensionless conductance of the system is large, and is nonperturbative in both disorder and interactions. The new phases are the mesoscopic analogue of spontaneous distortions of the Fermi surface induced by interactions in bulk systems and can occur in any Fermi liquid channel with angular momentum $`m`$. Here we show that the phase with $`m`$ even has a diamagnetic persistent current (seen experimentally but mysterious theoretically), while that with $`m`$ odd can be driven through a transition which spontaneously breaks time-reversal symmetry by increasing the coupling to dissipative leads.

\]

The interplay of disorder and interactions is a rich source of unexplained phenomena in the bulk, especially in two dimensions, despite three decades of theoretical effort. In mesoscopic systems one is confronted with phenomena not seen in the bulk, such as Coulomb Blockade oscillations of the zero-bias conductance, or persistent currents in mesoscopic rings in a small external magnetic field(for reviews see). It has been realized in the last several years that one needs to take interactions seriously in order to understand the experimental Coulomb Blockade peak spacing statistics. This understanding has led to Universal Hamiltonian treatments for weak-coupling (electron gas parameter $`r_s=1/a_0\sqrt{\pi \rho }1`$, where $`a_0`$ is the Bohr radius and $`\rho `$ is the electron density), in which a constant charging interaction and a constant exchange interaction are kept in addition to the single-particle energy.

Recently, two of us investigated the stability of the Universal Hamiltonian to other interactions of the Landau Fermi-liquid type, which are expected to be present in ballistic quantum dots with chaotic boundary scattering, but are not in dots deep in the diffusive limit. The Landau interaction $`u(\theta (𝐤)\theta (𝐤^{}))`$ depends only on the angle between the momenta $`𝐤,𝐤^{}`$ of the interacting particles, and can be parameterized in two dimensions by a set of Landau parameters $`u_m`$:

$$u(\theta \theta ^{})=u_0+\underset{m=1}{\overset{\mathrm{}}{}}u_m\mathrm{cos}m(\theta \theta ^{})$$
(1)
The values of the Landau parameters depend on the strength and range of interactions, and can only be accessed numerically for large $`r_s`$ (strong interactions). It was found that while the Universal Hamiltonian was stable in the regime $`u_m>u_m^{}=1/2\mathrm{ln}2`$, it became unstable to a mesoscopic Pomeranchuk phase for $`u_m<u_m^{}`$. Soon afterwards one of us and R. Shankar showed how to construct the large-$`N`$ theory of the strong-coupling regime (the dimensionless conductance $`g`$ plays the role of $`N`$). Details of this treatment will appear soon.

In this Letter, treating Coulomb Blockade and persistent currents within the same approach, we show that the mesoscopic Pomeranchuk phases display unexpected signatures in the persistent current, including a diamagnetic persistent current (seen experimentally but so far unexplained) in a model without superconductivity for $`m`$ even, and spontaneous time-reversal symmetry breaking for $`m`$ odd.

We want the effective Hamiltonian in an energy window of width the Thouless energy $`E_T`$ around the Fermi energy. This is the regime of validity of Random Matrix Theory (RMT). For ballistic structures $`E_T=\mathrm{}v_F/L`$, where $`v_F`$ is the Fermi velocity and $`L`$ is the linear system size. The dimensionless conductance is defined as the number of single-particle energy states (of mean spacing $`\mathrm{\Delta }`$) in this window $`g=E_T/\mathrm{\Delta }`$. Our effective Hamiltonian has a noninteracting part representing the chaotic scattering at the walls, and a Fermi-liquid-like interacting part which conserves momentum. The order parameter $`𝝈`$ of the Pomeranchuk phase is a two-dimensional vector whose magnitude $`\sigma `$ and direction $`\chi `$ represent the size and direction of the maximum Fermi surface distortion. The shape of the deformed Fermi surface is given by $`\sigma \mathrm{cos}(m\theta \chi )`$, where $`m`$ is the angular momentum channel in which the instability occurs. To determine the behaviour of $`𝝈`$ the fermions are integrated out and an effective action is obtained, the dominant part of which is self-averaging for large $`g`$.

$`\overline{S}_{eff}=`$ $`g^2{\displaystyle 𝑑t\left(\frac{|𝝈|^2}{2\mathrm{\Delta }}\left(\frac{1}{|u_m|}\frac{1}{|u_m^{}|}\right)+\lambda (𝝈^2)^2\right)}+`$ (2)
$`g{\displaystyle \frac{d\omega }{2\pi }|𝝈(\omega )|^2f(\omega )}`$ (3)
The “kinetic” term $`f(\omega )`$ behaves like $`\omega ^2`$ for $`\omega \mathrm{\Delta }`$ and like $`|\omega |`$ for $`\omega \mathrm{\Delta }`$, indicating the Landau damping of $`𝝈`$. For large $`g`$ the saddle-point of this action dominates, since fluctuations are down by $`1/g`$. A strong enough attractive Landau parameter $`u_mu_m^{}`$ leads to symmetry-breaking.

We have numerically evaluated the effective potential for the collective variable $`𝝈`$ and analyzed its dependence on external magnetic flux. This effective potential can be obtained as the ground state energy $`(𝝈)`$ of a noninteracting fermion Hamiltonian where $`𝝈=\sigma (\mathrm{cos}\chi \widehat{i}+\mathrm{sin}\chi \widehat{j})`$ appears as a parameter:

$$(H_\sigma (\chi ))_{\alpha \beta }=\epsilon _\alpha \delta _{\alpha \beta }g\sigma M_{\alpha \beta }(\chi )$$
(4)
where the first term encodes the chaotic scattering (with the eigenvalues $`\epsilon _\alpha `$ controlled by RMT), and $`M_{\alpha \beta }(\chi )`$ represents the coupling between the collective mode $`𝝈`$ and particle-hole excitations of the fermions. In the large-$`g`$ limit, the ground state energy of the system in the strong-coupling regime is simply the value of the effective potential at its global minimum. This energy automatically contains the contributions which are specific to the particular disorder realization in addition to the self-averaging contributions determnined earlier. The averaged energy landscape of $`𝝈`$ in the strong-coupling phase would be a “Mexican Hat”, with rotational symmetry. The sample-specific contributions break this symmetry completely for $`m`$ even, leading to a single minimum, and to a two-fold symmetry for $`m`$ odd, as can be seen in Figure 1.

To see the relation of the above to the persistent current we note that

$$I(\varphi )=\frac{F(\varphi )}{\varphi }.$$
(5)
where $`\varphi `$ is the external flux piercing the sample. At zero temperature the Free energy $`F`$ is just the ground state energy $``$, so we desire to obtain $``$ as a function of $`\varphi `$. We obtain this by taking the noninteracting part of the Hamiltonian from the RMT ensemble of crossover Hamiltonians parametrized as

$$H_{cross}=\sqrt{\frac{1}{1+C^2\varphi ^2}}\left(H_S+C\varphi H_A\right)$$
(6)
where $`H_{S,A}`$ represent symmetric ($`𝒯`$ preserving) and antisymmetric ($`𝒯`$ breaking) random matrices drawn from their respective normalized ensembles, and $`C`$ is a factor of order unity which depends on the shape of the sample and the precise nature of the chaotic scattering at the boundary.

Figure 2 shows the dependence of $``$ on the crossover parameter $`C\varphi `$ for $`m`$ even and odd, and clearly shows the diamagnetic behaviour for $`m`$ even. To see why this is special, consider what is known about persistent currents. In a mesoscopic ring penetrated by a flux, the ground state energy has to be periodic in the flux, since an integer number of flux quanta can be gauged away.

$$I_{pers}(\varphi )=\frac{F}{\varphi }=I_1\mathrm{sin}(2\pi \varphi /\varphi _0)+I_2\mathrm{sin}(4\pi \varphi /\varphi _0)+\mathrm{}$$
(7)
where $`\varphi _0=h/e`$ is the flux quantum. Only the even moments $`I_{2n}`$ survive disorder-averaging.

In the noninteracting case the typical, fluctuating values of the Fourier coefficients are (for small $`n`$) $`I_{n,typ}E_T/\varphi _0`$ while the average is $`I_{2n}\mathrm{\Delta }/\varphi _0`$. Experiments typically measure $`I_{2n}`$ and a few other low-order harmonics.

Interactions, when included in renormalized first-order perturbation theory, produce $`I_2\mu ^{}\frac{E_T}{\varphi _0}`$ where $`\mu ^{}`$ (of order $`1`$) is the dimensionless Cooper-channel interaction at low energies. Thus, interactions enhance the average persistent current, but if $`\mu ^{}>0`$ $`I_2`$ should be paramagnetic, while if $`\mu ^{}<0`$ it should be diamagnetic. This prediction, while of the same order-of-magnitude as the experiments, has the wrong sign. Materials that show no sign of superconductivity (implying that $`\mu ^{}>0`$) show a diamagnetic $`I_2`$. Many explanations have been proposed to account for this puzzle (a recent one being Ref.), but the question remains open, as summarized in Ref.. In this context a diamagnetic persistent current of order $`E_T/\varphi _0`$ in a model without superconductivity is striking. The fact that our treatment is nonperturbative in the interactions enables us to evade the usual sign. Our approach, while suggestive, is not directly applicable to the experiments on $`Au`$ and $`Ag`$ rings since the samples are not likely to be in the strong-coupling regime, and are not fully in the ballistic limit (the elastic mean free path is of the same order as the system size). On the other hand, our theory would apply directly to ensembles of ballistic $`GaAs`$ rings of the type used in Refs. , but at stronger coupling.

Let us now turn to $`m`$ odd. The exact degeneracy of the two global minima separated by $`\pi `$ in the angle $`\chi `$ can be seen from Figure 1, and can be proved analytically using the relation $`H_\sigma ^{}(\chi )=H_\sigma (\chi +\pi )`$. The two degenerate minima are related by the time-reversal transformation $`𝒯`$. A particular value of $`\chi `$ leads to a distortion of the Fermi surface along the direction specified by $`\chi `$. Under $`𝒯`$, $`𝐤𝐤`$ and a distortion of the Fermi surface for odd $`m`$ maps to an inequivalent state at $`\chi +\pi `$ with the same energy, since the underlying Hamiltonian is $`𝒯`$-invariant. The ground state of a Hamiltonian quantum system with a two-fold degenerate potential is the symmetric combination of the two minima. This applies to the isolated mesoscopic structure, whose dynamics is Hamiltonian at energy scales smaller than $`\mathrm{\Delta }`$ (for energy scales in the range $`\mathrm{\Delta }\omega E_T`$ the dynamics is dissipative with ohmic dissipation, see Eq. (3)). The splitting between the symmetric and antisymmetric combinations is the tunneling amplitude between the two minima, here $`\mathrm{\Delta }e^g`$. The two minima correspond to states carrying opposite persistent currents, and are macroscopically distinguishable.

The coupling of the mesoscopic structure to the leads produces ohmic dissipation at arbitrarily low energies. This is precisely the case of the Caldeira-Leggett model considered and solved by Chakravarty, and Bray and Moore. The effective action of our model at low energies ($`\omega \mathrm{\Delta }`$) is

$`g{\displaystyle 𝑑t(V(\chi (t))+\frac{1}{2}(d\chi /dt)^2)}+`$ (8)
$`{\displaystyle \frac{2g\mathrm{\Gamma }^2\sigma ^2}{\pi \mathrm{\Delta }^2}}{\displaystyle \frac{dtdt^{}}{(tt^{})^2}\mathrm{sin}^2\left(\frac{\chi \chi ^{}}{2}\right)}`$ (9)
where $`\chi `$ is the angle of $`𝝈`$ in the Mexican Hat, $`V(\chi )`$ is the doubly-degenerate realization-specific potential of Figure 1, and $`\mathrm{\Gamma }`$ is the level width induced by coupling to the leads. The long-range interaction in imaginary time comes from a $`|\omega |`$ “kinetic”term, which in turn arises from the Landau damping of $`𝝈`$ due to decay into particle-hole pairs at arbitrarily low energies, possible because each formerly sharp level $`\alpha `$ is broadened by coupling to the leads.

The model has a weak-dissipation phase in which the ground state is still the symmetric superposition of the two minima, and a strong-dissipation phase in which the particle is localized in one minimum. The transition between the two phases occurs for $`g(\mathrm{\Gamma }/\mathrm{\Delta })^2𝝈^21`$. For large enough $`g`$ even a weak coupling to the leads ($`\mathrm{\Gamma }\mathrm{\Delta }/\sqrt{g}`$, since $`𝝈^21`$) is sufficient to meet this criterion, and leads to localization in one minimum of the twofold degenerate effective potential even at zero temperature, corresponding to a spontaneous breaking of $`𝒯`$. The $`𝒯`$-breaking transition could be monitored by measuring the peak-height statistics.

If one turns on an external flux, one minimum moves up in energy as the flux increases while the other moves down. The ground state (which moves down) displays a paramagnetic persistent current of order $`E_T/\varphi _0`$. In the isolated dot, or in the case with weak dissipation, one starts with a symmetric superposition of the two minima as the zero-flux ground state. As $`\varphi `$ increases the system crosses over to fully $`𝒯`$-broken dynamics when the energy difference of the two minima is greater than their splitting in zero field, which is $`\mathrm{\Delta }e^g`$. Thus, the crossover will occur for an external flux $`\varphi _X\varphi _0e^g`$, as compared to the noninteracting crossover flux $`\varphi \varphi _0/\sqrt{g}`$. For strong enough dissipation, the ground state already breaks $`𝒯`$, and the variation of $``$ contains a term first-order in $`\varphi `$, which implies a spontaneous persistent current at zero flux.

Figure 3 shows the ensemble-averaged second Fourier coefficient $`I_2`$ for $`m`$ even and odd. Ideally one would find a periodic behaviour of $``$ with $`\varphi `$. This cannot be captured by the crossover Hamiltonian, but must be put in by hand. Since we do not know the number $`C`$ connecting the crossover parameter to the flux $`\varphi `$, there is an inherent ambiguity in this procedure. We show $`I_2`$ for three different choices for $`C`$. As can be seen, the qualitative results are unaffected.

Finite temperature produces additional interesting effects for $`m`$ odd. For temperatures exceeding the tunnel splitting $`\mathrm{\Delta }e^g`$ the superposition states are irrelevant, and one can think of the two minima as separately thermally occupied. Because of their (near) degeneracy, there will be a Curie-like susceptibility of $`1/T`$ of the persistent current (and therefore the magnetization) to external flux.

In summary, we have shown that some surprising signatures of the mesoscopic Pomeranchuk regimes show up in the persistent current. There is a diamagnetic persistent current (for $`m`$ even) without any superconductivity. The $`m`$ odd case undergoes a spontaneous time-reversal symmetry-breaking transition as the coupling to the leads is increased, and displays a spontaneous persistent current (at zero flux) in the $`𝒯`$-broken phase. It would be very interesting to explore the behaviour of persistent currents in the quantum critical regime, where large fluctuations of the order parameter $`𝝈`$ and finite quasiparticle lifetime at low energies are expected.

We thank R. Shankar for illuminating discussions and E. Mucciolo for pointing out Ref. , and the National Science Foundation for partial support under grants DMR 98-04983 (DH and HM) and DMR 0071611 (GM).


text:
# Elimination of spiral chaos by periodic force for the Aliev-Panfilov model

## Abstract

Spiral chaos appears in the two dimensional Aliev-Panfilov model. The generation mechanism of the spiral chaos is related to the breathing instability of pulse trains. The spiral chaos can be eliminated by applying periodic force uniformly. The elimination of spiral chaos is most effective, when the frequency of the periodic force is close to that of the breathing motion.
PACS numbers:05.45.a,05.45.Xt, 87.19.Hh

Some types of cardiac arrhythmia are characterized by rotating waves, which are similar to spiral waves found in excitable media . Control and elimination of arrhythmia and the spiral waves are medically important. The control of regular spiral patterns in excitable media has been studied with several methods. Meandering of the spiral core can be controlled by a periodic parameter modulation , impulses and periodic force have been applied to suppress spiral waves , and local and global feedback have been applied to excitable systems to eliminate spiral waves . A serious cardiac arrhythmia such as ventricular fibrillation is related to spiral chaos where many spiral cores are spontaneously generated. To eliminate spiral chaos is more important.

Periodic forcing to spatio-temporal chaos was used to restore regular waves for diffusively coupled chemical oscillators and the complex Ginzburg-Landau equation . In this paper, we attempt to control and eliminate numerically spiral chaos by applying a periodic force and find an effective frequency for the elimination. We use the Aliev-Panfilov model for the cardiac cell :

$`{\displaystyle \frac{e}{t}}`$ $`=`$ $`K(ea)(e1)er+^2e,`$
$`{\displaystyle \frac{r}{t}}`$ $`=`$ $`[ϵ+\mu _1r/(\mu _2+e)][rke(eb1)].`$ (1)
Here $`e`$ stands for the membrane potential and $`r`$ stands for the conductance of the inward current. This model is a phenomenological model which represents certain feature of impulse propagation in cardic tissue. The parameter values of $`K,a,ϵ,\mu _1,\mu _2`$ and $`b`$ are evaluated based on the real experiment. The model equation exhibits spiral breakup and spiral chaos in a certain parameter range .

We perform numerical simulations of Eq. (1). The parameters except for $`a`$ are fixed as $`K=8,ϵ=0.01,\mu _1=0.11,\mu _2=0.3`$ and $`b=0.1`$. Firstly, we show numerical results for a one-dimensional system, that is, $`^2e`$ in Eq. (1) is replaced by $`^2e/x^2`$. The numerical simulation was performed using the finite difference method with $`\mathrm{\Delta }t=0.005`$ and $`\mathrm{\Delta }x=0.5`$. The system size $`L`$ is 500, and the periodic boundary conditions are used. We have found that a pulse train propagates with a constant velocity for large $`a`$, however, the pulse train exhibits a breathing instability by decreasing $`a`$. Figure 1(a) displays the time evolution of $`e(x)`$ at $`a=0.063`$ for the pulse train with wavenumber $`k=2\pi 10/L0.126`$. The pulse width is temporally oscillating. The oscillatory instability of a single pulse in a one-dimensional ring of excitable media is called the alternans instability, since the pulse width almost alternates every time when the pulse turns around the ring . The oscillatory instability of a pulse train with spatial period $`2\pi /k`$ almost corresponds to the alternans instability in a ring of size $`2\pi /k`$. The alternans instability in a ring of size $`2\pi /k=L/10`$ occurs as a supercritical Hopf bifurcation at $`a=a_{c0}0.0647`$. However, the oscillatory instability of the pulse train occurs at slightly larger $`a=a_c0.0664`$ for $`k=2\pi 10/L`$, accompanying spatial modulation. As seen in Fig. 1(a), the phase of the breathing motion is not synchronized for all pulses in the pulse train, as the result of the spatial modulation. As $`a`$ is decreased further, the breathing amplitude becomes larger, and the spatial modulation grows. Finally, one pulse is annihilated and the wavenumber is decreased to $`k=2\pi 9/L`$. The time evolution of the wavenumber decrease process is shown in Fig. 1(b) for $`a=0.044`$. It is characteristic of the cardiac tissue that the pulse width is easily varied with spatial periods of pulse trains. This characteristics is involved in the Aliev-Panfilov model. The critical parameters for the breathing instability and the wavenumber changing bifurcation strongly depend on the pulse interval. Figure 2 displays the two bifurcation curves (the breathing instability and the wavenumber decreasing bifurcation) as a function of wavenumber $`k`$. As the pulse interval $`2\pi /k`$ is narrower, the instabilities occur more easily, and the two bifurcation curves approach each other. The wavenumber changing process occurs just after the breathing instability for $`k0.2`$.

Next, we consider the instability of the spiral pattern in two dimensions. A stationary rotating spiral becomes unstable in a certain parameter range as shown by Panfilov . We have performed a numerical simulation of Eq. (1) in two dimensions. The system size is $`300\times 300`$ and no flux boundary conditions are used. A spiral pattern is stable for $`a>0.13`$ for the fixed parameters $`K=8,ϵ=0.01,\mu _1=0.11,\mu _2=0.3`$ and $`b=0.1`$. The initial condition is a regular spiral pattern obtained numerically for $`a=0.13`$. Figure 3(a) displays a time evolution of $`e(x,y)`$ on the line $`y=Lx`$ at $`a=0.115`$. At this parameter, pulse trains with wavenumber $`k0.185`$ are emitted from the spiral core. A one-dimensional pulse train with wavenumber $`k0.185`$ is unstable for the breathing motion at $`a=0.115`$ as is shown in Fig. 2. The pulse trains indeed exhibit breathing motion, and the wavenumber decreasing processes occur at $`t520`$. In one-dimensional system, the wavenumber decreasing process leads to a more stable structure with smaller wavenumber, however, in two dimensions, the wavenumber decreasing process leads to formation of additional topological defects, the spiral breaks up, and then spiral chaos appears. This is another route of spiral breakup, although it is similar to the spiral breakup via a wavenumber changing process by the Eckhaus instability . Figure 4(b) displays a snapshot of $`e`$ at $`t=800`$, where $`e`$ takes a larger value than 0.4 in the shaded region. Two main spirals and several small spirals are generated as a result of the spiral breakup.

To eliminate the pulse train and the spiral chaos, we apply an external periodic force. The model equation is written as

$`{\displaystyle \frac{e}{t}}`$ $`=`$ $`k(ea)(e1)er+^2e+F\mathrm{sin}(\omega t),`$
$`{\displaystyle \frac{r}{t}}`$ $`=`$ $`[ϵ+\mu _1r/(\mu _2+e)][rke(eb1)],`$ (2)
where $`F`$ and $`\omega `$ are the amplitude and frequency of the external periodic force. We first show a numerical result of a one-dimensional system. Figure 4(a) displays a time evolution for $`F=0.009`$, $`a=0.06`$, $`k=2\pi 10/L`$ and $`\omega =2\pi /60`$. The initial condition is a breathing state for Eq. (1) without the external force. The breathing motion of the pulse width is entrained to the external force. The amplitude of the breathing motion grows, and finally the pulse train structure collapses completely for $`t>230`$. Thus, the pulse train could be eliminated by applying the external periodic force. We have investigated a critical value of the amplitude $`F`$ for the complete collapse. Figure 4(b) displays the critical values $`F_c`$ as a function of the period $`T=2\pi /\omega `$ for $`a=0.06`$. The phase diagram has a shape like the Arnold tongue for the forced entrainment, and the critical value $`F_c`$ takes the smallest value at $`T57.5`$. The period is close to the period $`T^{}61`$ of the natural breathing motion without the periodic forcing. These results are interpreted as a kind of resonance. That is, if the period of the external force is close to the period of the natural breathing motion, the effect of the external force is enhanced and the pulse train collapses easily.

We have applied a periodic force to eliminate the spiral chaos in two dimensions. The system size is $`200\times 200`$ and the parameter $`a`$ is 0.1. The spiral chaos appears for this parameter. We have used a snapshot of the spiral chaos as an initial condition for the forced system. In the simulation of the two dimensional system, we have applied the periodic force $`F\mathrm{sin}\omega t`$ to the spiral chaos only for $`0<t<6\pi /\omega `$ (three periods), and observed whether the spiral chaos collapses or the spiral chaos is maintained after the periods of the external forcing. The three periods are sufficient to see the effect of the external forcing. Figure 5 displays three snapshots of pattern of $`e`$ at $`t=0,100`$ and $`300`$ for $`a=0.1,F=0.036`$ and $`\omega =2\pi /60`$. The spiral chaos collapses completely and the uniform state $`e(x,y)=0`$ and $`r(x,y)=0`$ has been obtained after the application of external force. We have numerically obtained the critical value $`F_c`$ for the collapse of the spiral chaos. The critical strength $`F_c`$ is shown in Fig. 6 as a function of the period $`T=2\pi /\omega `$. The critical curve takes minimum at $`T50`$. It is close to the period $`55`$ of the natural breathing motion at this parameter.

To summarize, we have performed numerical simulations of the Aliev-Panfilov model. We have found a breathing instability for the pulse trains in one dimension. The breathing instability leads to spiral breakup in two dimensions. The periodic force is applied to eliminate the pulse trains and the spiral chaos. We have found that the most effective frequency to eliminate the wave patterns is close to the natural frequency of the breathing motion. This is interpreted as a kind of the resonance effect. The periodic forcing with the most effective frequency may be effective as a method of mild defibrillation.


text:
# 1 The quark probability distribution at finite temperature for μ=0.

Quark Probability Distribution at Finite Temperature and Density
Lianyi He, Guang Bian, Jinfeng Liao and Pengfei Zhuang

Physics Department, Tsinghua University, Beijing 100084, China

## Abstract

Quark deconfinement phase transition at finite temperature and density is investigated in the frame of quantum mechanics. By solving the Schrödinger equation for a heavy quark in a thermal mean field, we calculate the quark probability distribution as a function of temperature and density. The confined wave function in the vacuum expands outward rapidly when the temperature and density are high enough. The obtained phase transition line agrees qualitatively with the result of lattice QCD.

It is widely accepted that the collective effect of a multiparton system can change the vacuum structure of Quantum Chromodynamics (QCD), and the partons confined in a hadron bag can move out when the temperature and density are high enough and form a new state of matter, the so called quark-gluon plasma (QGP). This kind of new matter may be produced in the early stage of a high-energy heavy-ion collision. The dense partons, produced in the vacuum at high temperature or compressed at high baryon density, make the quark wave functions overlap in a space with dimension of the colliding nuclei. An important question in relativistic heavy ion collisions is how to identify the new state of matter if it is created in the early stage. According to the numerical calculation of lattice QCD, the critical temperature of quark deconfinement is about 150 MeV when the density effect is excluded.

The idea that the overlap of wave functions, or the space extension of wave functions at finite temperature and density is the original reason of quark deconfinement phase transition arises naturally an interesting question to study the behavior of quark distribution function in the frame of quantum mechanics. In this letter, we will study a quark moving in a thermal mean field describing multiparton interactions, and obtain the quark probability distribution by solving the Schrödinger equation at finite temperature and density. The temperature and density at which the quark wave function expands outward rapidly can be considered as the critical point of quark deconfinement. Because of the limitation of non-relativistic quantum mechanics, we focus on the motion of a heavy quark at finite temperature and density, corresponding to, for example, the dissociation of $`J/\psi `$ as a signature of QGP in relativistic heavy ion collisions.

In the string-like models that describe the quark confinement well, one uses a linear potential $`V=kr`$ to express the interaction between two quarks. To simplify the numerical calculation in the following, we use a three dimensional square well

$`V(r)=\{\begin{array}{cc}0\hfill & 0<r<a\hfill \\ V_0\hfill & a<r<b\hfill \\ \mathrm{}\hfill & r>b\hfill \end{array}`$ (4)
instead of the linear potential. Here $`a`$ is the boundary of the region where the motion of the quark is asymptotically free, and can be taken as the radius of the hadron constructed by the considered quarks, $`b`$ is the space scale of the system, and $`V_0`$ should be so large that the quark wave function can be confined in the asymptotically free region $`r,a`$. Since the region of collective interaction in high-energy heavy-ion collisions is of the order of the colliding nuclei, we treat $`b`$ as the radius of a heavy nucleus. In our numerical calculations, we take the parameters of the confinement potential as $`a=1`$ fm, $`b=5`$ fm, and $`V_0=10`$ GeV.

At finite temperature and density, the quark will be affected by many-body interactions in the region $`r>a`$ in addition to the confinement potential. In mean field approximation, the thermodynamic potential of the quasi-particle system is

$$\mathrm{\Omega }(T,\mu )=g\frac{d^3𝐩}{(2\pi )^3}\frac{p^2}{3E_q}\left(f_q(T,\mu )+f_{\overline{q}}(T,\mu )\right)$$
(5)
with the quark and antiquark densities

$`f_q(T,\mu )`$ $`=`$ $`{\displaystyle \frac{1}{e^{\frac{E_q\mu }{T}}+1}},`$
$`f_{\overline{q}}(T,\mu )`$ $`=`$ $`{\displaystyle \frac{1}{e^{\frac{E_q+\mu }{T}}+1}},`$ (6)
where $`g`$ is the quark degree of freedoms of flavors, colors and spins, $`T`$ and $`\mu `$ are respectively the temperature and baryon chemical potential, $`E_q=\sqrt{m_{eff}^2+p^2}`$ is the quark energy with effective mass $`m_{eff}`$ obtained in the mean field. Because the contribution from heavy quarks to the thermodynamic potential is very small compared with light quarks, we will only take $`u`$ and $`d`$ quarks into account, i.e., $`g=12`$.

From the rules of thermodynamics, the pressure of the system is

$$P(T,\mu )=\mathrm{\Omega }(T,\mu ).$$
(7)
Considering the relation between the potential $`V_{th}`$ and the force $`\stackrel{}{F}_{th}`$ that the heavy quark we study feels in the region $`r>a`$,

$$\stackrel{}{F}_{th}(\stackrel{}{r})=V_{th}(\stackrel{}{r}),$$
(8)
i.e.,

$$\frac{dV_{th}(r)}{dr}=4\pi r^2P(T,\mu ),$$
(9)
and assuming that the potential falls down with increasing $`r`$ and vanishes on the boundary of the system, namely $`V_{th}(r=b)=0`$, the potential $`V_{th}`$ can be expressed as

$$V_{th}(r|T,\mu )=_r^b𝑑r^{}4\pi r^2P(T,\mu )=\frac{4}{3}\pi (b^3r^3)\mathrm{\Omega }(T,\mu ).$$
(10)
With the obtained potential, the motion of the heavy quark at finite temperature and density is described by the Schrödinger equation

$$i\frac{}{t}\psi (\stackrel{}{r},t)=\left(\frac{1}{2m}^2+V_{eff}(r|T,\mu )\right)\psi (\stackrel{}{r},t),$$
(11)
where $`m`$ is the mass of the heavy quark, and $`V_{eff}`$ is the effective potential defined as

$`V_{eff}(r|T,\mu )=\{\begin{array}{cc}0\hfill & 0<r<a\hfill \\ V_0+V_{th}(r|T,\mu )\hfill & a<r<b\hfill \\ \mathrm{}\hfill & r>b\hfill \end{array}`$ (15)
By separating the variables of the wave function $`\psi (\stackrel{}{r},t)=R(r)Y(\theta ,\phi )e^{iEt}`$, and letting $`\varphi (r)=rR(r)`$, we obtain the radial equation for a stationary state

$$\frac{1}{2m}\frac{d^2\varphi (r)}{dr^2}+\left(V_{eff}(r|T,\mu )+\frac{l(l+1)}{2mr^2}\right)\varphi (r)=E\varphi (r),$$
(16)
where $`E`$ is the eigenvalue of the quark energy. Since we focus in this letter on quark deconfinement at finite temperature and density, we discuss in the following the ground state of the Schrödinger equation only. The corresponding radial equation is then simplified as

$$\frac{1}{2m}\frac{d^2\varphi _0(r)}{dr^2}+V_{eff}(r|T,\mu )\varphi _0(r)=E_0\varphi _0(r).$$
(17)
We first discuss qualitatively the deconfinement phase transition through the behavior of the effective potential as a function of $`T`$ and $`\mu `$. From the $`T`$\- and $`\mu `$-dependence of the thermodynamic potential, we can immediately obtain the conclusion that the confinement potential in the vacuum at $`T=\mu =0`$ will be suppressed by the thermal mean field $`V_{th}<0`$, and the motion of the quark will extend outside when the temperature and density increase. Since the mean field decreases monotonously from the maximum on the boundary $`r=a`$ to the minimum on the other boundary $`r=b`$, the suppression at $`r=a`$ is most strong, and the potential at $`r=b`$ does not change at finite temperature and density.

In order to illustrate quantitatively the above estimated quark deconfinement phase transition, we should solve the Schrödinger equation (17) at finite temperature and density to get the quark probability distribution

$$\rho (r|T,\mu )r^2R^2(r|T,\mu ).$$
(18)
We take the considered heavy quark in the Schrödinger equation as a $`c`$ quark which forms $`J/\psi `$ with $`m=1.6`$ GeV, and the effective mass of the light $`u`$ and $`d`$ quarks in $`V_{eff}`$ as their constituent quark mass $`m_{eff}=0.35`$ GeV.

In the limit of $`\mu =0`$, the probability distribution as a function of temperature obtained by solving numerically the Schrödinger equation (17) is shown in Fig.(1). In the vacuum with $`T=\mu =0`$, the chosen parameters for the square well guarantee that almost the whole quark wave function is restricted in the region $`r<a`$, that means quark confinement. As temperature increases the quark motion extends outside very slowly in the beginning, then the extension is accelerated rapidly in the neighborhood of $`T=120`$ MeV, which indicates that the two $`c`$ quarks forming $`J/\psi `$ begin to dissociate in the hot mean field, and finally the quark moves in the whole region of the system $`0<r<b`$. The case $`T=0`$ is similar to that of $`\mu =0`$, the quark probability distribution as a function of $`r`$ for different chemical potentials is indicated in Fig.(2). The probability expands outward rapidly in the neighborhood of $`\mu =445`$ MeV.

It is well known that for a particle moving in a central field the most probable radius describes the bound degree of the stationary state. The most probable radius $`r_0(T,\mu )`$ of the $`c`$ quark motion at finite temperature and density defined by

$$\frac{\rho (r|T,\mu )}{r}|_{r=r_0}=0$$
(19)
is indicated as a function of $`T`$ for $`\mu =0`$ in Fig.(3). It is clear to see that around $`T=120`$ MeV $`r_0`$ goes up from $`a/2`$ to $`b/2`$ very fast. We can determine the critical temperature $`T_c`$ and chemical potential $`\mu _c`$ of the quark deconfinement phase transition as the point at which the most probable radius $`r_0(T_c,\mu _c)`$ reaches $`a`$, the boundary of the asymptotically free region in the vacuum,

$$r_0(T_c,\mu _c)=a.$$
(20)
The critical line which separates the deconfinement phase from the confinement phase is shown in Fig.(4). While the most probable radius $`r_0`$ jumps from $`a/2`$ to $`b/2`$ in a very narrow region around $`T_c`$, it is still continuous at the critical point, and the phase transition we study is a continuous one. When the temperature and chemical potential exceed the critical values, the region where the quark is asymptotically free will expand, i.e., the position of $`a`$ will move outward.

In summary, by describing the collective effect as a thermal mean field and combining it with the confinement potential in the vacuum, we have demonstrated the physical picture of quark deconfinement in the frame of quantum mechanics through solving the Schrödinger equation and calculating the quark probability distribution at finite temperature and density. When the temperature and density are high enough, the quark wave function extends outward from inside the hadron. The critical temperature and density around which the wave function extension happens rapidly separates the confinement from deconfinement. The obtained critical temperature $`T_c=120`$ MeV at $`\mu =0`$ agrees qualitatively with the lattice QCD calculation.

Acknowledgments: The work was supported by Grant Nos. NSFC19925519, 10135030, and G2000077407.


text:
# References

ON THE ANOMALOUS STRUCTURES OF THE VECTOR LEPTONIC CURRENTS

Rasulkhozha S. Sharafiddinov

Institute of Nuclear Physics, Uzbekistan Academy of Sciences, Tashkent, 702132 Ulugbek, Uzbekistan

Each of existing types of the electric charges come forwards in the system as the source of a kind of the dipole moment. Therefore, to investigate these regularities we have established the compound structures of Dirac and Pauli form factors. They state that the electron possesses as well as the anomalous electric charge.

Owing to the vector nature of virtual photon, the elastic scattering of electrons by spinless nuclei depends on the Dirac $`F_{1e}(q^2)`$ and Pauli $`F_{2e}(q^2)`$ form factors of light leptons . They are of course the functions of the square of four - dimentional momentum transfer. However, in spite of large number of works dedicated to the interaction between the electron and field of emission, thus far remains many uncertainties both in the structures and in the behavior of these currents. Usually it is accepted that $`F_{2e}(0)`$ is equal to the electron anomalous magnetic moment , and its full magnetic moment is defined by the combination  of form factors

$$\mu _e^{full}=\frac{F_{1e}\left(0\right)}{2m_e}+F_{2e}\left(0\right).$$
(1)
It appears that here $`F_{1e}(0)`$ gives the electric charge leading to the appearance of the electron normal magnetic moment.

Such a procedure were based actually on the assumption  of that the functions $`F_{1e}(q^2)`$ and $`F_{2e}(q^2)`$ are not the Fourier transforms of the spatial distributions of the electric charge and magnetic moment of a particle. This is explained by some consequences of the classical model of an extensive electron .

According to the classical theory of electromagnetic mass , the availability of the eigenenergy $`E_0`$ of the electron electrostatic field implies the existence of the electric part of the electron rest mass:

$$m_e^{em}=\frac{E_0}{c^2}.$$

The opinion has been speaked out that all the mass of the electron is equal to its electromagnetic mass. Such an idea called simply a hypothesis of field mass and testifies in favor of the unsteadity of charge distribution of the electron.

We start from the duality of matter that the mass and charge of a particle correspond to the most diverse form of the same regularity of the nature of this field . It states that each of all possible types of charges arises as a consequence of the availability of a kind of the inertial mass . Thereby such a mechanism leads to the appearance of the intraelectron interratio between the forces of the electric and unelectric nature. Therefore, the charge distribution of the electron must be steady.

The purpose of the present work is to discuss some consequences and implications implied from the above - mentioned regularities of the nature of matter. They give of course the justification of that in the same presentation as the form factors $`F_{1e}(q^2)`$ and $`F_{2e}(q^2)`$ was used are not in the states to explain the observed vector picture of the electron. For understanding the mechanism of the anomalous interaction of Pauli at the fundamental level, one must elucidate the compound structures of these functions.

From such a purpose, we not only must write the form factors $`F_{ie}(q^2)`$ in the form

$$F_{ie}\left(q^2\right)=f_{ie}\left(0\right)+A_{ie}\left(\stackrel{}{q^2}\right)+\mathrm{}$$
(2)
but also need conclude that each of existing types of the electric charges come forwards in the system as the source of a kind of the dipole moment. Herewith the independent components $`f_{ie}(0)`$ coincide with the normal size of the electric charge and magnetic moment of the electron:

$$f_{1e}\left(0\right)=e_e^{norm},f_{2e}\left(0\right)=\mu _e^{norm}=\frac{e_e^{norm}}{2m_e^{norm}},$$
(3)
where and further it is necessary to keep in mind that $`e_e^{norm}`$ for a particle (antiparticle) has the negative (positive) sign.

The second terms $`A_{ie}(\stackrel{}{q^2})`$ characterize the dependence of form factors on the square of three - dimensional momentum transfer $`\stackrel{}{q^2}`$ and at the emission of a real photon ($`\stackrel{}{q^2}`$=0) are reduced to the values

$$A_{1e}\left(0\right)=e_e^{anom},A_{2e}\left(0\right)=\mu _e^{anom}=\frac{e_e^{anom}}{2m_e^{anom}}.$$
(4)
Here $`m_e^{norm}`$ and $`m_e^{anom}`$ are the Coulomb normal and anomalous masses. Insofar as the full electric mass is concerned, we will start from the fact that the Coulomb mass and charge of a particle correspond to two form of the same regularity of its electric nature. Therefore, we conclude  that

$$m_\nu ^E=m_\nu ^{norm}+m_\nu ^{anom}+\mathrm{}.$$
(5)
So, it is seen that any of form factors $`F_{1e}(q^2)`$ and $`F_{2e}(q^2)`$ includes in self both normal and anomalous interactions between the electron and field of emission. In other words, they must be Fourier transforms of the spatial density of charge and moment. Their value at zero four - dimensional momentum transfer $`(q^2=0)`$ defines the full static size of the electric charge and magnetic moment of a particle:

$$F_{1e}\left(0\right)=e_e^{full}=e_e^{norm}+e_e^{anom}+\mathrm{},$$
(6)
$$F_{2e}\left(0\right)=\mu _e^{full}=\mu _e^{norm}+\mu _e^{anom}+\mathrm{}.$$
(7)
By following the compound structures of form factors (3), we get

$$f_{2e}\left(0\right)=\frac{f_{1e}\left(0\right)}{2m_e^{norm}}.$$
(8)
Exactly the same one can found from (4) that

$$A_{2e}\left(0\right)=\frac{A_{1e}\left(0\right)}{2m_e^{anom}}.$$
(9)
According to the presented here point of view, the electron possesses as well as the anomalous electric charge which has an estimate of

$$e_e^{anom}=\frac{\alpha }{2\pi }\left(\frac{m_e^{anom}}{m_e^{norm}}\right)e_e^{norm}$$
(10)
in assuming that the size of $`A_{2e}(0)`$ is equal to the electron Schwinger magnetic moment:

$$A_{2e}\left(0\right)=\mu _e^{anom}=\frac{\alpha }{2\pi }\frac{e_e^{norm}}{2m_e^{norm}}.$$
(11)
To brighter reveal our ideas one must apply to the process of elastic scattering of electrons and their neutrinos by spinless nuclei as to the source of unique information about structures of leptonic currents. It is already clear from (2) that in the case of one - photon exchange only the independent components of form factors are responsible for the interaction with matter. Therefore, a study of the behavior of light leptons $`(l=e,\nu _e)`$ in the nucleus charge field leads us to the equation

$$2m_l^{norm}\frac{f_{2l}\left(0\right)}{f_{1l}\left(0\right)}=\pm 1.$$
(12)
Comparison of (12) with (8) say in favor of correspondence principle which states that each terms of the expansions (2) correspond to the definite approximations . Under such circumstances the possibility of the inclusion of the anomalous phenomena $`A_{ie}(\stackrel{}{q^2})`$ in the discussion is realized only in the second Born approximation. Nevertheless, without loss of generality, we must have in view of that any non - zero component of the interaction of Pauli implies the availability of a kind of the Dirac interaction. Of course, the above - noted regularities of vector picture of the electron and its neutrino open up new possibilities for developments of our sights at the nature of matter.