text:
# Field-induced gap in the spin-1/2 antiferromagnetic Heisenberg chain: A density matrix renormalization group study

## I Introduction

The quasi-one-dimensional magnetic materials have attracted considerable experimental and theoretical interest since Haldane’s pioneering work  that pointed out the difference between the integer spin Heisenberg antiferromagnetic (HAFM) chains and the half-integer spin chains. By mapping the Heisenberg spin chains onto the $`O(3)`$ nonlinear sigma model, Haldane conjectured that the low-energy excitation spectrum displays a finite energy gap for the integer spin systems while for half-integer spin chains it is gapless. These conjectures have been verified by later detailed studies.

In the linear chain HAFM family, the spin-$`\frac{1}{2}`$ chain is of particular interest since most of its properties can be obtained exactly. These exact results serve as benchmarks in testing the validity of various approximation schemes. The low energy excitation spectrum in this system has no gap; the elementary excitations are spin-$`\frac{1}{2}`$ spinons; the ground state is quasi-long-range ordered, and the spin-spin correlations display power law decay.

The effect of an applied magnetic field on the spin-$`\frac{1}{2}`$ Heisenberg chain has also been studied to gain more insight into the physics of such a system. When the external magnetic field is present, the Hamiltonian of the system is written as

$$H=J\underset{i}{}(𝐒_i𝐒_{i+1}h_uS_i^z),$$
(1)
where $`J`$ is the coupling constant, $`𝐒_i`$ the spin-$`\frac{1}{2}`$ operator on site $`i`$, $`S_i^z`$ the $`z`$ component of $`𝐒_i`$, $`h_u=g\mu _BH/J`$ is the effective dimensionless uniform field, $`g`$ the average effective gyromagnetic ratio and $`H`$ the applied magnetic field. When $`h_u=0`$, the critical wave vector of the gapless excitation is located at 0 and $`\pi `$. The applied field will shift the critical wave vector of the gapless excitation away from 0 ( for transverse spin correlations) and $`\pi `$ ( for longitudinal spin correlations) to incommensurate values while the excitation remains gapless until the field is larger than its saturation value $`h_u=2.0`$. When $`h_u=2.0`$, the magnetization is saturated with all spins on the lattice sites oriented parallel to the applied uniform field. Further increase in the applied field leads to the opening of a gap in the low energy spectrum with its magnitude changing almost linearly with the applied field, corresponding to flipping one spin to its opposite direction.

The field-induced incommensurate state was first observed in the neutron scattering measurements on copper benzoate, Cu(C6D5COO)2$``$3D2O. Copper benzoate is a linear chain spin-$`\frac{1}{2}`$ AFM  with coupling constant $`J`$ 1.57 meV. In this material, the effective spin-$`\frac{1}{2}`$ Cu2+ ions form a linear chain structure, but the two neighboring copper sites are not totally equivalent. Because of the small value of the coupling constant $`J`$, it is possible to study the high field (large $`h_u`$) properties of Hamiltonian (1) and investigate how the induced incommensurate soft mode behaves with the changing magnetic field.

In the copper benzoate experiment,  in addition to the field-dependent incommensurate low energy modes, an unexpected non-zero energy gap induced by the magnetic field was detected. The value of the gap varies with the magnitude and the relative orientation of the applied magnetic field. It was first suggested that the unexpected gap is caused by the inequivalence of the Cu sites leading to an effective staggered $`g`$ tensor, which in turn gives rise to an effective staggered field in addition to the applied uniform magnetic field. Later, detailed analysis  shows that an additional staggered Dzyaloshinskii-Moriya interaction term provides similar contribution along with the staggered $`g`$ tensor, and both have the same order in magnitude. The Dzyaloshinskii-Moriya interaction together with the staggered $`g`$ factor can account for the observed non-zero energy gap in copper benzoate. After making some assumptions and ignoring the small exchange anisotropy, the effective Hamiltonian to describe the copper benzoate in the magnetic field can be written as:

$$H=J\underset{i}{}[𝐒_i𝐒_{i+1}h_uS_i^z(1)^ih_sS_i^x],$$
(2)
where $`h_s`$ is the induced effective dimensionless staggered field, which is the key term to account for the observed non-zero energy gap. The magnitude of $`h_s`$ depends on the magnitude and relative direction of the applied uniform field $`h_u`$ with respect to the sample.

Hamiltonian (2) has been studied using the bosonization approach, mapping on the Sine-Gordon model, and form-factor techniques. The gap and magnetization behavior, the dynamical magnetic susceptibility, the specific heat  and the electron spin resonance experiments have been analyzed using the sine-Gordon quantum field theory. There have also been numerical studies of the excitation energy and transition amplitudes in the staggered magnetic field, based on the Bethe Ansatz solutions. These studies focus mainly on the parameter range corresponding to the reported copper benzoate magnetic field experiment.  The field-induced gap and related magnetization, as well as spin correlations in generic spin-$`\frac{1}{2}`$ antiferromagnetic Heisenberg chains are of great interest but yet to be fully understood. Of particular interest is the study of a wider range of parameter and very high magnetic field conditions beyond those probed by the experiment and previous theoretical work, which may prove to be of significant importance in understanding the entire range of field-induced phenomena and the underlying physics in this interesting system.

In this paper, we report results of our numerical calculations of the ground state and the low energy excitations of Hamiltonian (2) using the density matrix renormalization group (DMRG) method. We study the most general case and take $`h_u`$ and $`h_s`$ as independent variables in our calculations. We find that the critical (saturation) uniform field $`h_u^c=2.0`$ serves as an important reference point in understanding the obtained results. When fixed $`h_u`$ is lower than or equal to $`h_u^c`$, the induced energy gap increases with $`h_s`$ as a power with exponent $`\frac{2}{3}`$, and for fixed $`h_u`$ larger than $`h_u^c`$, the small $`h_s`$ dependence is exponential. On the other hand, when $`h_s`$ is fixed, the $`h_u`$ dependence of the gap displays a minimum around $`h_u^c`$. When $`h_u`$ and $`h_s`$ increase simultaneously with a fixed ratio, the gap increases with the field when the ratio is small but develops a minimum around $`h_u^c`$ at larger ratios. The magnetization results are consistent with the intuitive expectation in general. The most interesting features are obtained for fixed $`h_u<h_u^c`$ and in small fixed $`h_s`$ cases. When $`h_u`$ is lower than $`h_u^c`$ and fixed, the existence of a small staggered field enhances the uniform magnetization instead of suppressing it; similar effects are also observed for small fixed $`h_s`$ case, where the staggered magnetization increases until $`h_u`$ approaches $`h_u^c`$. The uniform field induced incommensurate behavior is also studied, and the results show competing effects of the uniform and staggered fields, with the staggered field frustrating the incommensurate state.

In the following, the numerical DMRG results will be presented in Section II and a summary given in Section III.

## II DMRG Results

The DMRG method  is a powerful tool for the calculation of low lying states of quasi-one-dimensional systems and has been developed to calculate other properties of many strongly correlated systems. The accuracy of the DMRG calculations on spin chains is generally high. This has also been verified in the S=1 (Ref. ) and S=2 (Ref. ) Heisenberg chain calculations. In the present work, we employ the periodic boundary conditions (PBC) and use the infinite chain length algorithm of DMRG. We retain as many as 500 optimal states and compute up to chain length N=100 in each calculation. The largest truncation errors are of the order of 10-9 for zero uniform field calculations and 10-6 for non-zero uniform field calculations. To simplify the discussion, we set the coupling constant $`J`$ as the energy unit, $`h_u`$ as the effective uniform field and $`h_s`$ as the effective staggered field.

For the Hamiltonian considered here, the effect of the uniform field is to induce a uniform magnetization and shift the critical wave vector from $`\pi `$ (In this paper we will concentrate on this case for the pitch vector, where a peak in the static structure factor is expected), and the staggered field will induce a non-zero energy gap between the ground state and the lowest excited state. A non-zero staggered magnetization will also be expected for finite staggered fields. The nature of the ground state depends on the competition of the uniform and the staggered field. We will discuss the energy spectrum of the Hamiltonian first, and then show the ground-state magnetization behavior and investigate the incommensurate behavior of the spin correlation functions.

### A Energy Gap

When the staggered field is present alone in the spin-$`\frac{1}{2}`$ chain, that is, $`h_u`$=0 in Hamiltonian (2),

$$H=\underset{i}{}[𝐒_i𝐒_{i+1}(1)^ih_sS_i^x]$$
(3)
the x-component of the total spin $`S_{tot}^x`$ is conserved. This property can be used in DMRG calculations to reduce the dimension of the Hilbert space to be considered.

The spin-1 case of Hamiltonian (3) has been studied in detail. For the standard spin-1 Heisenberg chain, the low energy spectrum is gapful, and the lowest excited state is a spin triplet known as Haldane triplet. The presence of a staggered field will split the Haldane triplet into two branches, the transverse branch and the the longitudinal branch. Both branches are gapful and the gap increases with the staggered field.

For the spin-$`\frac{1}{2}`$ case, the excitation spectrum for zero-field chain is gapless, and the presence of the staggered field may also open a finite gap between the ground state and the low-energy continuum. It is expected that the gapful excitations will also split into two branches, with the x-component of total spin $`S_{tot}^x`$=0 (longitudinal branch) and 1 (transverse branch).

The field dependence of the gap of the longitudinal and transverse branches for Hamiltonian (3) is shown in Fig. 1 (a). It is clear that the low energy spectrum becomes gapful as soon as the staggered field is present. The magnitude of both longitudinal and transverse modes increases when the staggered field becomes larger with the longitudinal mode goes up faster than the transverse one. This behavior is exactly the same as for the spin-1 chain. But for the spin-1 chain in the staggered field, the increase of the longitudinal gap is about three (two) times faster than the transverse one for small (large) staggered field. For the spin-$`\frac{1}{2}`$ chain case, this ratio is smaller. The two gaps are fitted using the equation

$$\mathrm{\Delta }=ah_s^b,$$
(4)
where $`a`$ and $`b`$ are fitting parameters. The least square fitting gives:

$`\mathrm{\Delta }_L`$ $`=`$ $`a_Lh_s^{b_L}=2.97h_s^{0.678},`$ (5)
$`\mathrm{\Delta }_T`$ $`=`$ $`a_Th_s^{b_T}=1.97h_s^{0.63},`$ (6)
where $`\mathrm{\Delta }_L`$ and $`\mathrm{\Delta }_T`$ denote the longitudinal and transverse gap, respectively. The ratio of the gap increase coefficient of the longitudinal mode to the transverse mode $`a_L/a_T`$ is about 1.5. The fitting curves Eq. (6) are also shown in Fig. 1(a). For larger staggered field, the fitting is almost perfect, while it is not very good for very small staggered field. (This deviation is not visible in Fig. 1 (a).) It should be noted here that the numerical error of DMRG results for small staggered field is much bigger than that for large field.

When the uniform field in Hamiltonian (2) is not vanishing, the competition of the uniform and the staggered field needs to be taken into account. The presence of the uniform field will affect the behavior of the energy gap. In the non-zero uniform field case, the calculation is more difficult, since now even $`S_{tot}^x`$ is not conserved any more. There is no good quantum number which can be used to reduce the relevant Hilbert space dimension. To investigate the gap behavior, the ground state and the lowest excited state must be calculated at the same time, and at least two states must be targeted in each calculation.

To study the effect of the uniform field on the gap behavior, we first calculate the staggered field dependence of the gap when uniform field $`h_u`$ is a non-zero constant. As described above, when the uniform field is present alone, there is a critical point $`h_u^c=2.0`$, where the magnetization is saturated. We need to perform calculations at different values of $`h_u`$, one in each part of the $`h_u`$ phase space. We choose $`h_u`$=0.5 for the small uniform field case and $`h_u`$=3.0 for the large field case. We also calculate the $`h_s`$ dependence of the gap at the uniform field critical point $`h_u^c`$=2.0.

The results for $`h_u`$=0.5, 2.0 and 3.0 are shown in Fig. 1 (b), (c), (d), respectively. For $`h_u`$=0.5 and 2.0, the spectrum is gapless when $`h_s`$=0, and an energy gap opens up when the staggered field is present. The gap increases with the staggered field $`h_s`$ following the same function Eq. (4) as in the $`h_u`$=0 case, although with different parameters. Least square fitting gives $`a`$=1.89, $`b`$=0.624 for $`h_u`$=0.5 and $`a`$=1.63, $`b`$=0.81 for $`h_u`$=2.0. Comparison with results in Fig. 1 (a) shows that both the coefficient and the exponent in the $`h_u=0.5`$ case differ only slightly from those of the transverse branch for $`h_u`$=0; for $`h_u`$=2.0, the difference becomes more pronounced. The calculated results show that in the uniform field gapless phase ($`h_u2.0`$), the induced gap is affected by the uniform field only when the uniform field is strong enough. For $`h_u`$=3.0, the gap dependence is different. The energy spectrum is gapful even when $`h_s`$=0. With the application of the staggered field, the gap increases. The increase of the gap is nearly exponential when the staggered field is not very strong, following $`\mathrm{\Delta }=\mathrm{exp}(ah_s^b)`$ with parameters a=1.186 and b=1.939. When the staggered field is large enough, the gap increase deviates from the exponential behavior.

In the experiment on Cu benzoate, the magnitude of the induced staggered field depends on the applied uniform field. Roughly speaking, when the relative orientation of the applied field is fixed, the induced staggered field increases linearly with the applied uniform field. To compare with experiment directly, we have also considered the case when $`h_s`$ and $`h_u`$ increase at the same time with the ratio $`h_s/h_u`$ being fixed. In Fig. 1 (e) and (f), we present the calculated results of the energy gap with the staggered field for the cases of $`h_u=h_s`$ and $`h_u=10h_s`$. In both cases, the field dependence of the gap for smaller staggered field can also be fitted with function (4). For $`h_u=h_s`$, the fitting holds at $`h_s=h_u<2.0`$ with $`a`$=1.786 and $`b`$=0.594, while for $`h_u=10h_s`$, it holds for $`h_s<0.1`$ with $`a`$=1.755 and $`b`$=0.613. The gap for $`h_u=h_s`$ case increases monotonically but that for $`h_u=10h_s`$ exhibits a minimum near $`h_s0.2`$ ($`h_u2.0`$), and then increases rapidly. The largest induced staggered field produced in experiment is about $`h_s`$ 0.05, and it is not big enough to detect the gap minimum shown in Fig. 1 (f). Experiments at higher magnetic field are needed to test this predicted phenomenon. Here we can see again that the critical point $`h_u^c=2.0`$ plays an important role in separating different regions where the scaling behavior of the field-induced gap shows qualitative difference.

The effect of the uniform field on the spin gap induced by the staggered field can be studied directly by calculating the uniform field dependence of the gap at fixed staggered field. From the above results, we have learned that the system stays in one single gapful phase when the staggered field is present alone, so we just choose two sets of $`h_s`$, $`h_s`$=0.05 as the small staggered field limit and $`h_s`$=1.0 as the large staggered field limit. In Fig. 1 (g) and (h), we show the gap behavior with the change of the uniform field for the two cases. It is clearly seen that in both cases the gap decreases when the uniform field increases from zero and reaches its minimum near $`h_u2.0`$, followed by a rapid increase with further increase in the uniform field. The increase of the gap after the minimum is approximately linear for both cases.

Detailed calculations around the gap minimum provide more information on the gap behavior. For $`h_s=0`$, the system is gapless until the uniform field reaches the saturation point $`h_u^c=2.0`$. When $`h_s`$ is not zero, the system has a finite energy gap, and the presence of a small uniform field may suppress the gap. However, for small uniform fields, the suppression of the gap is negligible; it becomes visible only when the uniform field is large enough. The gap minimum for non-zero $`h_s`$ occurs when the uniform field $`h_u`$ is slightly larger than 2.0. It occurs at $`h_u2.02`$ for $`h_s`$=0.05 and $`h_u2.1`$ for $`h_s`$=1.0. This means that the minimum gap uniform field value $`h_u^{min}`$ increases slowly with $`h_s`$. After the gap minimum, its behavior is dominated by the uniform field, and the gap increases almost linearly with $`h_u`$. We should emphasize here that the minimum gap for both cases is not zero, and for the $`h_s`$=1.0 case the decrease in magnitude is a little bigger than that for the $`h_s`$=0.05 case.

### B Magnetization

The existence of a non-zero uniform field will induce a magnetization in the system, and the existence of a non-zero staggered field will induce an additional staggered magnetization. The staggered magnetization $`M_s(N)`$ and the uniform magnetization $`M_u(N)`$ of the system with a finite chain length $`N`$ are defined as

$`M_s(N)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}(1)^iS_i^x,`$ (7)
$`M_u(N)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}S_i^z.`$ (8)
The results for the thermodynamic limit

$`M_s`$ $`=`$ $`\underset{N\mathrm{}}{lim}M_s(N),`$ (9)
$`M_u`$ $`=`$ $`\underset{N\mathrm{}}{lim}M_u(N).`$ (10)
can be obtained by studying different chain-length systems. The results of the magnetization and the staggered magnetization for the eight parameter sets ($`h_u`$, $`h_s`$) used in the previous subsection are shown in Fig. 2.

When the uniform magnetization is absent ($`h_u=0`$), the induced staggered magnetization increases with the staggered field and approaches the saturation value 0.5 when $`h_s\mathrm{}`$. From the field dependence of the staggered magnetization, we can extract the staggered magnetic susceptibility $`\chi ^{(s)}=\frac{M_s}{h_s}`$. In Fig. 2 (a), we can see clearly that the staggered magnetic susceptibility goes to infinity when the staggered field $`h_s0`$ instead of approaching a constant as in the case of the spin-1 chain. This is because the spin-$`\frac{1}{2}`$ chain is gapless for zero staggered field while it has a finite gap in the spin-1 case. The magnetization curves are fitted using the following function:

$$M_s=ah_s^b$$
(11)
with $`a`$=0.527 and $`b`$=0.277. The fitting line is also shown in Fig. 2 (a). In the $`h_s`$ range shown in Fig. 2 (a), the fitting is good, but it should be noted that the fitting Eq. (11) will not be valid for very large $`h_s`$, since it diverges when $`h_s\mathrm{}`$ instead of approaching the finite value 1/2 which is the strong staggered field limit of Hamiltonian (3).

When $`h_u`$ is finite and fixed, the staggered magnetization increases monotonically with the staggered field and approaches the saturation value $`M_s`$=0.5 when $`h_s`$ goes to infinity, the same as in the $`h_u=0`$ case. The magnetizations for $`h_u=`$0.5, 2.0, and 3.0 are shown in Fig. 2 (b), (c), and (d), respectively. Because the zero staggered field energy spectrum for $`h_u`$=0.5 and $`h_u`$=2.0 is gapless, the zero-field staggered magnetic susceptibility at $`h_s`$=0 is still divergent. For $`h_u`$=3.0, the $`h_s=0`$ system is gapful, so the zero-field staggered magnetic susceptibility has a finite value, $`\chi ^{(s)}(0)`$=0.299. For $`h_u`$=2.0 and 3.0, the uniform magnetization is saturated when the staggered field is absent; it decreases from the saturation value 0.5 when a staggered field is applied. When the staggered field is weak, the way the magnetization changes is different for the two cases. It decreases rapidly for $`h_u`$=2.0 but slowly for $`h_u`$=3.0. For $`h_u`$=0.5, at zero staggered field, the uniform magnetization has a finite value but is not saturated; when the staggered field increases, the uniform magnetization will also increase and reach a maximum when $`h_s`$ 0.5 before decreasing with further increasing staggered field.

In Fig. 2 (e) and (f), we show the field dependence of the magnetization for the $`h_u=h_s`$ and $`h_u=10h_s`$ cases, respectively. When the uniform and staggered fields increase simultaneously from zero, both staggered and uniform magnetizations increase. For both cases, the zero-field staggered susceptibility is divergent and the zero-field uniform susceptibility is zero, we should note here that the zero-field uniform susceptibility is not shown clearly in Fig. 2 (e) and (f), its zero value can only be obtained when the results for very small field is investigated. When $`h_u=h_s`$, the effect of the staggered field rises rapidly and is dominating at small fields. When the field increases, the effect of the uniform field becomes more important and must be taken into account. At $`h_s=h_u`$ 1.2, the staggered magnetization reaches its maximum, and it decreases with further increasing field. The uniform magnetization increases monotonically with the increasing fields. Both staggered and uniform magnetizations approach finite but nonsaturated values when the field goes to infinity. For $`h_u=10h_s`$, the effect of the staggered field is also dominating at small field. But at large field, the effect of the uniform field becomes dominant. The staggered magnetization decays fast to a small finite value $``$ 0.05 after reaching its maximum. The uniform magnetization increases monotonically and approaches a nearly saturated value 0.498 at high enough field.

When the staggered field is fixed, at $`h_u=0`$ the staggered magnetization is finite while the uniform magnetization is zero. The uniform field dependence of the magnetization for $`h_s`$=0.05 and $`h_s`$=1.0 is shown in Fig. 2(g) and (h). For $`h_s`$=0.05, the staggered magnetization increases with the uniform field first, reaches its maximum at $`h_u`$ 1.8, and then decays to zero rapidly. For $`h_s`$=1.0, the staggered magnetization decreases with the uniform field monotonically. But the decrease is not rapid for small $`h_u`$. In both cases, the uniform magnetization increases with the uniform field. At small uniform field, the uniform magnetization increases linearly with the uniform field. We have obtained zero-field uniform magnetic susceptibility $`\chi ^{(u)}(0)`$ = 0.1185 for $`h_s`$=0.05 and $`\chi ^{(u)}(0)`$ = 0.1283 for $`h_s`$=1.0.

For the $`h_u`$=0.5 ($`h_s`$=0.05) case, the existence of a small staggered (uniform) field enhances the corresponding uniform (staggered) magnetization instead of suppressing it. This phenomenon can be explained intuitively. In these cases, when a small uniform (staggered) field is applied, the coupling between neighboring spins is weakened. While this uniform (staggered) field is not strong enough to destroy the effect of the stronger staggered (uniform) field, it enhances the ratio between the effective uniform (staggered) field and the effective coupling constant and, consequently, the induced uniform (staggered) magnetization.

### C Correlation Function and Incommensurate Behavior

We define three ground-state spin correlation functions for chain length $`L`$, (i) $`𝒞_u`$ parallel to the uniform magnetic field, (ii) $`𝒞_s`$ parallel to the staggered field, and (iii) $`𝒞_y`$ along the remaining ($`y`$) axis, as

$`𝒞_u(ij)`$ $`=`$ $`S_i^zS_j^z,`$ (12)
$`𝒞_s(ij)`$ $`=`$ $`S_i^xS_j^x,`$ (13)
$`𝒞_y(ij)`$ $`=`$ $`S_i^yS_j^y.`$ (14)
The correlation function $`𝒞_y`$ is expected to display exponential decay because of the existence of the spin gap induced by the staggered field. $`𝒞_s`$ and $`𝒞_u`$ do not decay exponentially because of the effects of the non-zero staggered and uniform magnetization, respectively. $`𝒞_u`$ is also expected to show incommensurate behavior due to the existence of the uniform field. These correlation functions have the following form:

$`𝒞_u(l)`$ $`=`$ $`M_u^2+(1)^lA_1{\displaystyle \frac{e^{l/\xi }}{\sqrt{l}}}\mathrm{cos}(al+\theta _0),`$ (15)
$`𝒞_s(l)`$ $`=`$ $`(1)^lM_s^2+(1)^lA_2{\displaystyle \frac{e^{l/\xi }}{\sqrt{l}}},`$ (16)
$`𝒞_y(l)`$ $`=`$ $`(1)^lA_3{\displaystyle \frac{e^{l/\xi }}{\sqrt{l}}},`$ (17)
where $`M_u`$ and $`M_s`$ are the uniform and staggered magnetization, respectively, $`A_1`$, $`A_2`$, $`A_3`$, $`a`$, $`\theta _0`$ are $`l`$ independent constants, and $`\xi `$ is the correlation length. It should be noted here that the length $`\xi `$ may not be the same along three different directions.

From the spin correlation functions Eq. (17), we can extract the corresponding correlation length. When $`h_u`$=0, the obtained correlation length is a function of the staggered field. The product of the correlation length and the spin gap gives the spin wave velocity of the system. In Fig. 3, we show the relation between the inverse correlation length and the transverse gap. It is seen that the curve goes linearly which means that the spin wave velocity does not change with the magnitude of the staggered field. The linear fitting of the line gives the spin wave velocity $`v=\mathrm{\Delta }_T\xi _T=1/0.6364=1.5713`$, in good agreement with the exact spin wave velocity for the spin-$`\frac{1}{2}`$ chain $`\frac{\pi }{2}`$ 1.5708. As in the spin-1 chain case, this also serves as an independent check of the self-consistency of our calculations. From the figure, we can see for zero stagger field spin-$`\frac{1}{2}`$ chain, the spectrum is gapless, $`\mathrm{\Delta }=0`$, so the correlation length is infinite. For other cases in our calculations, the correlation length is difficult to obtain, because the chain length is short and the numerical error is bigger for those cases. But we can conclude from our results that in general the spin wave velocity is not a constant any more, and, instead, it changes with the applied uniform field.

The cosine function in $`𝒞_u`$ comes from the incommensurate behavior induced by the uniform field. Because of the existence of the staggered field, the net incommensurate behavior is a result of competition between $`h_u`$ and $`h_s`$. In Fig. 4, we present $`\varphi (l)=\sqrt{l}e^{l/\xi }(𝒞_u(l)M_u^2)`$ at chain length 100 as a function of $`l`$ for $`h_u`$=0.5 and $`h_s`$=0.05, where $`M_u`$ is obtained by the magnetization calculation discussed above. This clearly shows the existence of incommensurability in the system. In this case, the correlation length $`\xi `$ 4.05.

The incommensurate behavior in the correlation function leads directly to the peak shift from $`\pi `$ in the static structure factor $`𝒮(q)`$ which can be obtained from the correlation functions. For the uniform field ($`z`$) axis and the staggered field ($`x`$) axis, we can write:

$`𝒮_u(q)`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{l}{}}e^{iql}𝒞_u(l),`$ (18)
$`𝒮_s(q)`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{l}{}}e^{iql}𝒞_s(l).`$ (19)
Since we use the periodic boundary conditions in our calculations, the wave vector $`q`$ is well defined, $`q=\frac{2\pi n}{L},n=1,\mathrm{},L`$.

In Fig. 5, we present the static structure factor $`𝒮_u`$ at $`h_u=0.5`$ and $`h_s`$=0.05 for the even chain length from 60 to 100. The results for different chain length systems fall onto the same curve. This success is due to the small correlation length of the system considered ($`\xi 4.05`$). The chain lengths used here are much larger than the correlation length, and the finite size effect is very weak. Fig. 5 shows a two peak structure symmetric about $`\pi `$ which is obtained by Eq. (18). For this case, we have obtained the critical wave vector shift $`\delta q=|q\pi |`$ 0.224 $``$ 0.07 $`\pi `$. The accuracy of $`\delta q`$ obtained from the peak deviation of $`𝒮_u(q)`$ is good; it is mainly limited by the finite chain lengths used in the calculations and the error is estimated to be less than 1 %.

For all parameters studied in these calculations, using the peak position of the static structure factor $`𝒮_u(q)`$ to determine the existence of the incommensurate state, we found the critical wave vector shift in three sets of parameters, $`h_u=0.5`$, $`h_u=10h_s`$, and $`h_s=0.05`$. The critical wave vector shift $`\delta q`$ versus field in these cases is shown in Fig. 6. Since the largest total chain length in our calculations is N=100, we cannot detect the incommensurate state if $`\delta q<\frac{2\pi }{100}=0.02\pi `$. For the $`h_u=0.5`$ case, the critical wave vector shift becomes smaller when the staggered field $`h_s`$ increases from zero. It is expected to go to zero when the staggered field is large enough. In other words, the existence of the staggered field will frustrate the incommensurate state. For the other two cases $`h_u=10h_s`$ and $`h_s=0.05`$, the critical wave vector shifts away from $`\pi `$ with increasing uniform field but eventually returns to $`\pi `$ around $`h_u^c`$.

For comparison, we examine our calculated results with parameters corresponding to the reported experimental case on Cu benzoate.  The highest magnetic field reached in the experiment is about 7T which corresponds to $`h_u`$ 0.52 . The energy gap observed at this field is $`\mathrm{\Delta }`$ 0.4 meV, considering the coupling constant for the material $`J`$=1.57 meV, $`\mathrm{\Delta }/J`$ 0.255. From our calculated results for $`h_u`$ = 0.5, we estimate that the induced staggered field $`h_s`$ 0.04. The critical wave vector shift $`\delta q`$ for $`h_u`$=0.5, and $`h_s`$=0.04 is $`\delta q`$ 0.274 $``$ 0.09 $`\pi `$. The critical wave vector shift for $`h_u`$=0.52, $`h_s`$=0.04 will be slightly larger than 0.09 $`\pi `$. In the experiment, the largest wave shift for H=7T is about 0.12 $`\pi `$. This comparison shows that the calculated results are in good agreement with the experiment.

## III summary

We have carried out systematic calculations using the density matrix renormalization group method to study the behavior of the energy gap, magnetization, and spin correlation functions of spin-$`\frac{1}{2}`$ antiferromagnetic Heisenberg chain in the presence of a uniform and a perpendicular staggered magnetic field. An extensive examination of the parameter space has revealed many interesting features beyond those reported in previous studies. In particular, results at very high magnetic field show quantitatively and even qualitatively different behaviors in the energy gap and magnetization, from those found for lower fields. For the $`h_u=10h_s`$ case, which is close to the real parameters in the experiments, we can see from Fig. 1 and Fig. 2 that the competition of the staggered and uniform field is visible only when $`h_s>0.15`$, that is $`h_u>1.5`$. In Cu benzoate, $`h_u=1.5`$ corresponds to a applied field $`H`$ 21 Tesla. Further experimental investigation at magnetic fields higher than 21 Tesla is needed to test these predictions.

The field-induced energy gap is dominated by the staggered field when the uniform field is smaller than the standard spin-$`\frac{1}{2}`$ chain saturation field $`h_u^c`$=2.0. When the uniform field is comparable or larger than $`h_u^c`$, the effect of the uniform field becomes important and must be taken into account. The uniform field introduces frustration effects and creates a local minimum in the energy gap near $`h_u^c`$ in several cases.

The magnetization results clearly reveal the competition of the uniform and staggered fields. However, for some uniform (staggered) field, the existence of a small staggered (uniform) field enhances the uniform (staggered) magnetization instead of suppressing it. The competition of the two kinds of field can also be seen from the incommensurate behavior with the staggered field suppressing the incommensurate state and moving the critical wave vector closer to the zero-field value $`\pi `$.

The specific heat experiment on Cu benzoate shows that the field-induced gap scales with approximately 2/3 power of the applied magnetic field. The analytic results yield the same power law dependence for Hamiltonian (2). In our numerical results, for $`h_u=0`$ case, the scaling power $`b_L`$ =0.678 for the longitudinal gap and $`b_T`$=0.63 for the transverse gap. They are in good agreement with the experiment and the analytic value. The existence of the non-zero uniform field modifies the power law relation. When the uniform field is small, the modification is almost negligible. Further high field experiments are needed to examine the predicted effect of the applied uniform field on the scaling behavior of the induced gap.

###### Acknowledgements.

J. Lou would like to thank Prof. M. Oshikawa, Prof. T. K. Ng, Dr. Xiaoqun Wang, and Dr. Tao Xiang for useful discussion. This work was supported by the Department of Energy at the University of Nevada, Las Vegas and by the Chinese Natural Science Foundation.


format: text

text:
# Nonrelativistic ionization energy for the helium ground state.

## Abstract

The helium ground state nonrelativistic energy with 24 significant digits is presented. The calculations are based on variational expansion with randomly chosen exponents. This data can be used as a benchmark for other approaches for many electron and/or three-body systems.

From the early days of the quantum mechanics the ground state ionization energy of the helium atom was a benchmark for approximate methods of solving nonrelativistic Schrödinger equation for a few–body system. One of the earliest variational calculations has been performed by Hylleraas Hylleraas  in 1929 and it yields an about 5 significant digits for the nonrelativistic ionization potential. In 1957, Kinoshita Kinoshita  presented a 7 digit number obtained with a 39 parameter function, which along with higher order corrections including the Lamb shift calculations confirmed a very good agreement with the best experimental value. Since that time with the development of computer power the accuracy grows very rapidly. We would like to mention here the two most recent calculations. The first Sims  is aimed to elaborate an efficient variational method for the many electron atoms. The second Drake  is to find an effective and economical way for studying the helium and helium-like two electron atoms.

In this short report we want to present a new very accurate value for the nonrelativistic energy of the helium ground state. In our calculations we strictly follows a method described in var00 . The two modifications are to be stated. First, a sextuple precision arithmetics (about 48 decimal digits) implemented as a Fortran-90 module has been used instead of a multiprecision package written by Bailey. The use of this module gives an about 5-fold improvement in computational time and allows for to increase significantly a length of a basis set. The module is based on a representation of a sextuple precision number by a set of 3 double precision numbers. It is assumed that an exact sum of these double precision numbers is some sextuple precision number sextuple . Second, we have taken a multilayer variational wave function with 5 independent sets of variational parameters (instead of 4 as in var00 ), which consecutively approximates one after another smaller and smaller distances of electrons with respect to a nucleus.

A variational wave function is expanded in a form var00

$$\begin{array}{c}\psi _0=\underset{i=1}{\overset{N/2}{}}\left\{U_i\text{Re}\left[\mathrm{exp}(\alpha _ir_1\beta _ir_2\gamma _ir_{12})\right]+W_i\text{Im}\left[\mathrm{exp}(\alpha _ir_1\beta _ir_2\gamma _ir_{12})\right]\right\}\hfill \end{array}$$
(1)
where $`\alpha _i`$, $`\beta _i`$ and $`\gamma _i`$ are complex parameters generated in a quasi-random manner:

$$\alpha _i=\left[\frac{1}{2}i(i+1)\sqrt{p_\alpha }(A_2A_1)+A_1\right]+i\left[\frac{1}{2}i(i+1)\sqrt{q_\alpha }(A_2^{}A_1^{})+A_1^{}\right],$$
(2)
$`x`$ designates the fractional part of $`x`$, $`p_\alpha `$ and $`q_\alpha `$ are some prime numbers, $`[A_1,A_2]`$ and $`[A_1^{},A_2^{}]`$ are real variational intervals which need to be optimized. Parameters $`\beta _i`$ and $`\gamma _i`$ are obtained in a similar way. The actual values of these parameters for the calculation with the largest basis set of 5200 functions is presented in Table I. As is seen from the Table fine tuning of variational parameters is not required that greatly facilitates calculations.

Table II demonstrates a convergence of the variational expansion with the number of basis functions. Extrapolated value has been obtained by means of the simple extrapolation formula,

$$E(\text{extrap})=E(N)C\times N^\nu ,$$
(3)
where parameters $`C`$ and $`\nu `$ are taken from the best fit of the last 4 or 5 calculations.

In Table III a comparison with the most recent and most accurate values is presented. Our result extends accuracy of the previous calculations by more than three orders of magnitude.

The author would like to thank J.S. Sims and G.W.F. Drake for stimulating to publish this work.


text:
# 1 Introduction

## 1 Introduction

There is mounting evidence showing the existence of a light scalar isoscalar meson . However, its nature as a genuine meson state or as a resonance dynamically generated in the $`\pi \pi `$ scattering is still under debate. A review of the current status of the discussion can be found in ref. .

In any case, it is such a broad resonance that its effects are hardly visible in any phase shifts or decay plots. The situation could be quite different at finite densities and/or temperatures where the $`\sigma `$ could become lighter and narrow . Apart from its implications related to chiral symmetry restoration, this topic is of special interest due to the relevance of the $`I=J=0`$ channel in the nucleon-nucleon interaction. Any substantial change of this meson properties could alter the in medium nucleon-nucleon interaction and therefore our current understanding of the nuclear matter in heavy ion collision, neutron stars or even in normal nuclei if the changes occur at low densities.

The density dependence of the $`\sigma `$ properties has been studied in several models. Bernard et al.  have found a decrease of the mass as a function of the density using a generalized Nambu Jona-Lasinio model, until it converges to the mass of the pion, its chiral partner. It is noticeable, that in the same model the $`\rho ,\omega `$ and $`\pi `$ mesons have almost constant masses.

In ref. , Hatsuda et al. studied the $`\sigma `$ propagator in the linear $`\sigma `$ model and found an enhanced and narrow spectral function near the $`2\pi `$ threshold caused by the partial restoration of the chiral symmetry, where the $`m_\sigma `$ mass would approach $`m_\pi `$. The same conclusions were reached using the nonlinear chiral lagrangians in ref. .

Similar results, with large enhancements in the $`\pi \pi `$ amplitude around the $`2\pi `$ threshold, have been found in a quite different approach by studying the $`swave`$, $`I=0`$ $`\pi \pi `$ correlations in nuclear matter . In these cases the modifications of the $`\sigma `$ channel are induced by the strong $`pwave`$ coupling of the pions to the particle-hole ($`ph`$) and $`\mathrm{\Delta }`$-hole ($`\mathrm{\Delta }h`$) nuclear excitations. It was pointed out in  that this attractive $`\sigma `$ selfenergy induced by the $`\pi `$ renormalization in the nuclear medium could be complementary to additional $`swave`$ renormalizations of the kind discussed in  calling for even larger effects.

On the experimental side, there are also several results showing strong medium effects in the $`\sigma `$ channel at low invariant masses in the $`A(\pi ,2\pi )`$  and $`A(\gamma ,2\pi )`$  reactions. At the moment, the cleanest signal probably corresponds to the $`A(\gamma ,2\pi ^0)`$ reaction, which shows large density effects that had been predicted in both shape and size in ref. , using a model for the $`\pi \pi `$ final state interaction along the lines of the present work.

Our aim in this paper is to study the $`\sigma `$ mass and width at finite densities in the context of the model developed in . These works, which provide a very economical and successful description of a wide range of hadronic phenomenology, use as input the lowest orders lagrangian of chiral perturbation theory  and calculate meson-meson scattering in a coupled channels unitary way. The $`\sigma `$ meson is not a basic ingredient on this approach, it is dynamically generated and appears as a pole of the $`\pi \pi `$ scattering amplitude in the second Riemann sheet.

The nuclear medium effects on the scalar isoscalar channel were implemented in this framework in refs. . As in other approaches, large medium effects were found. Namely, the imaginary part of the $`\pi \pi `$ scattering amplitude showed a clear shift of strength towards low energies as the density increases.

In the next section we present, for the sake of completeness, a brief description of the model used for the $`\pi \pi `$ interaction in both vacuum and dense medium, which is already published elsewhere. The following section describes the method used to calculate the $`\pi \pi `$ scattering amplitude in the complex energy plane and to search for the $`\sigma `$ pole position at finite densities.

## 2 $`\pi \pi `$ interaction

In this work we consider only the scalar isoscalar ($`\sigma `$) channel and follow the simple method of ref.  for $`\pi \pi `$ interaction in vacuum and refs.  for the nuclear medium effects. Additional information on this and related approaches for different spin isospin channels can be found in refs. .

### 2.1 Vacuum

The basic idea is to solve a Bethe Salpeter (BS) equation, which guarantees unitarity, matching the low energy results to the chiral perturbation theory ($`\chi PT`$) predictions. We consider two coupled channels, $`\pi \pi `$ and $`K\overline{K}`$ and neglect the $`\eta \eta `$ channel which is not relevant at the low energies we are interested in.

The BS equation is given by

$$T=V+VGT.$$
(1)
Eq. 1 is a matrix integral equation which involves the two mesons one loop divergent integral, (see Fig. 1), where $`V`$ and $`T`$ appear off shell. However, for this channel both functions can be factorized on shell out of the integral. The remaining off shell part can be absorbed by a renormalization of the coupling constants as it was shown in refs. . Thus, the BS equation becomes purely algebraic and the $`VGT`$ originally inside the loop integral becomes then the product of $`V`$, $`G`$ and $`T`$, with $`V`$ and $`T`$ the on shell amplitudes independent of the integration variables, and $`G`$ given by the expression

$$G_{ii}(P)=i\frac{d^4q}{(2\pi )^4}\frac{1}{q^2m_{1i}^2+iϵ}\frac{1}{(Pq)^2m_{2i}^2+iϵ}$$
(2)
where $`P`$ is the momentum of the meson-meson system. This integral is regularized with a cut-off adjusted to optimize the fit to the $`\pi \pi `$ phase shifts ($`\mathrm{\Lambda }=1.03`$ GeV).

The potential $`V`$ appearing in the BS equation is taken from the lowest order chiral Lagrangian

$$L_2=\frac{1}{12f^2}<(_\mu \mathrm{\Phi }\mathrm{\Phi }\mathrm{\Phi }_\mu \mathrm{\Phi })^2+M\mathrm{\Phi }^4>$$
(3)
where the symbol $`<>`$ indicates the trace in flavour space, $`f`$ is the pion decay constant and $`\mathrm{\Phi }`$, $`M`$ are the pseudoscalar meson and mass $`SU(3)`$ matrices. This model reproduces well phase shifts and inelasticities up to about 1.2 GeV. The $`\sigma `$ and $`f_0(980)`$ resonances appear as poles of the scattering amplitude. The coupling of channels is essential to produce the $`f_0(980)`$ resonance, while the $`\sigma `$ pole is little affected by the coupling of the pions to $`K\overline{K}`$ .

### 2.2 The nuclear medium.

As we are mainly interested in the low energy region, which is not very sensitive to the kaon channels, we will only consider the nuclear medium effects on the pions. The main changes of the pion propagation in the nuclear medium come from the p-wave selfenergy, produced basically by the the coupling of pions to particle-hole ($`ph`$) and Delta-hole ($`\mathrm{\Delta }h`$) excitations. For a pion of momentum $`q`$ it is given by

$$\mathrm{\Pi }(q)=\frac{\left(\frac{D+F}{2f}\right)^2\stackrel{}{q}^2U(q)}{1\left(\frac{D+F}{2f}\right)^2g^{}U(q)}$$
(4)
with $`g^{}=0.7`$ the Landau-Migdal parameter, $`U(q)`$ the Lindhard function and $`(D+F)=1.257`$. The expressions for the Lindhard functions are taken from ref. .

Thus, the in medium BS equation will include the diagrams of Fig. 2 where the solid line bubbles represent the $`ph`$ and $`\mathrm{\Delta }h`$ excitations.

In fact, as it was shown in , the contact terms with the $`ph`$ ($`\mathrm{\Delta }h`$) excitations of diagrams (b)(c)(d) cancel the off shell contribution from the meson meson vertices in the term of Fig. 2(a). Hence, we just need to calculate the diagrams of the free type (Fig. 1) and those of Fig. 2(a) with the amplitudes factorized on shell. Therefore, at first order in baryon density, we are left with simple meson propagator corrections which can be readily incorporated by changing the meson vacuum propagators by the in medium ones.

There are some other medium corrections, like the one depicted in Fig. 3

which has been considered in several papers like  using the linear sigma model.

In our approach, we start with only pseudoscalar meson fields. As it was discussed previously, the sigma is generated dynamically through the rescattering of the pions. A possible analog to the diagram of Fig. 3 is given by the diagrams of Fig. 4.

We evaluate the contribution of these diagrams by using the chiral Lagrangians  involving the octet of baryons and the octet of pseudoscalar mesons. We find these terms to be proportional to $`\overline{p}\gamma ^\mu p\overline{n}\gamma ^\mu n`$ and thus they vanish in symmetric nuclear matter.

Also we would need to consider diagrams of the type shown in Fig. 5

We get a negligible contribution to the the $`\pi \pi `$ scattering in the nuclear medium from this diagram  and for simplicity we do not include it in this work.

The $`\pi \pi `$ scattering amplitude obtained using this model exhibits a strong shift towards low energies. In Fig. 6, we show the imaginary part of this amplitude for several densities. Quite similar results have been found using different models  and it has been suggested that this accumulation of strength, close to the pion threshold, could reflect a shift of the $`\sigma `$ pole which would approach the mass of the pion.

Other pion selfenergy contributions related to $`2ph`$ excitations, and thus proportional to $`\rho ^2`$, can be incorporated in the pion propagator. As we are most interested in the region of low energies we can take as estimation the corresponding piece of the optical potentials obtained from pionic atoms data, following the procedure of ref.  and substituting in Eq. 4

$$\left(\frac{D+F}{2f}\right)^2U(q)\left(\frac{D+F}{2f}\right)^2U(q)4\pi C_0^{}\rho ^2$$
(5)
with $`\rho `$ the nuclear density and $`C_0^{}=(0.105+i0.096)m_\pi ^6`$. Its effects are small except at large densities as can be appreciated by comparing Fig. 6, with Fig. 7 of ref.  where this piece is included.

## 3 $`\sigma `$ meson mass and width as a function of the baryon density

In this section we will consider the lightest pole position of the $`\pi \pi `$ scattering amplitude in the S=I=0 channel in the complex energy plane, which corresponds to the $`\sigma `$ resonance. In our model, the vacuum pole position occurs at a mass around 500 MeV and the width is around 400 MeV. A compilation of values for these magnitudes obtained in several modern analyses can be found in ref. . See also Fig. 1 of ref.  and .

The analytical structure of the meson meson scattering amplitude is driven by the meson loops, which are given by the formula

$$G_{\pi \pi }=i_0^{\mathrm{}}\frac{d^4q}{(2\pi )^4}D_\pi (q^0,\stackrel{}{q})D_\pi (\sqrt{s}q^0,\stackrel{}{q})$$
(6)
where $`D_\pi `$ is the $`\pi `$ propagator which, in the nuclei, incorporates the p-wave coupling of the pions to particle-hole and $`\mathrm{\Delta }`$-hole excitations. We can simplify the calculation by using the Lehmann representation for the meson propagators

$$D_\pi (q^0,\stackrel{}{q})=\frac{2}{\pi }_0^{\mathrm{}}𝑑x\frac{xImD_\pi (x,\stackrel{}{q})}{(q^0)^2x^2+iϵ}$$
(7)
After some algebraic manipulation we obtain

$$G_{\pi \pi }=_0^{\mathrm{}}\frac{dE}{2\pi }(\frac{1}{P^0E+iϵ}\frac{1}{P^0+Eiϵ})F(E)$$
(8)
where $`F(E)`$ is a real function given by

$$F(E)=\frac{d^3q}{(2\pi )^3}_E^E\frac{dx}{\pi }ImD_\pi ((E+x)/2,\stackrel{}{q})ImD_\pi ((Ex)/2,\stackrel{}{q})$$
(9)
which in vacuum takes the simple form

$$F_{vac}(E)=\frac{p_\pi }{4\pi E}=\frac{1}{4\pi }\sqrt{\frac{1}{4}\frac{m_\pi ^2}{E^2}}$$
(10)
As we are interested in the positive energy region, we can concentrate our attention on the first term of Eq. 8 which is the only one that has a cut. We can rewrite that term, neglecting the $`ϵ`$, like

$$\frac{F(E)}{P^0E}\frac{F(Real(P^0))}{P^0E}+\frac{F(E)F(Real(P^0))}{P^0E}$$
(11)
The $`E`$ integration of the first term it is the only piece that depends on the Riemann sheet and can be readily evaluated. In the second sheet the integral takes the value

$$_{0,(II)}^{\mathrm{}}𝑑E(\frac{F(Real(P^0))}{P^0E})=_{0,(I)}^{\mathrm{}}𝑑E(\frac{F(Real(P^0))}{P^0E})+2\pi iF(Real(P^0))$$
(12)
Substituting this result in Eq. 8 we get the following expression for the two meson propagator in the second sheet

$$G_{\pi \pi }^{(II)}=iF(Real(P^0))+_0^{\mathrm{}}\frac{dE}{2\pi }(\frac{1}{P^0E}\frac{1}{P^0+E})F(E)$$
(13)
The rest of the pieces entering the BS equation are analytical single valued functions and therefore we can proceed to look for poles in the scattering amplitude.

## 4 Results

The function $`F(E)`$ includes the phenomenological information on the pion selfenergy in the nuclear medium. The results for several densities are shown in Fig. 7.

The $`F(E)`$ strength is related to the imaginary part of the $`\sigma `$ propagator and therefore reflects the energy and density dependence of the different $`\sigma `$ decay channels. We can start classifying them according to their density dependence (See. Fig.8).

At low densities, the $`\sigma `$ meson can only decay into two pions, $`\sigma \pi +\pi `$. Therefore we have a threshold at $`E=2m_\pi `$. As the density grows the $`\sigma \pi +(Nh)`$ decay channel becomes relevant and we have a new threshold at $`E=m_\pi `$, which is clearly visible in the curve corresponding to $`\rho =\rho _0/8`$. Finally, at larger densities, mechanisms such as $`\sigma 2(Nh)`$ which go like $`\rho ^2`$ become important. They are possible even at very low energies.

Using the two meson propagator of Eq. 8 we can solve the BS equation, obtain the meson-meson scattering amplitude and look for the poles in the complex energy plane. The results for the $`\sigma `$ pole position are shown in Fig. 9 for densities up to 1.5 $`\rho _0`$. Note, however, that the calculation is more reliable at low densities because some contributions of order $`\rho ^2`$ or higher are missing.

We find that both mass and width decrease as the density increases, reaching a mass around 250 MeV and a similar width at 1.5 times the nuclear density. At large densities the $`2ph`$ pieces of the pion selfenergy become relevant decreasing further both mass and width. These results are little sensitive to the kaon channel. Its omission affects less than 1 percent the mass and increases the width around 5 percent. The results of Fig. 9 can be cast in terms of an effective potential which can be approximated by

$$V_\rho =a\frac{\rho }{\rho _0}+b\left(\frac{\rho }{\rho _0}\right)^2$$
(14)
with $`a=358i\mathrm{\hspace{0.17em}108}`$ MeV and $`b=140+i\mathrm{\hspace{0.17em}23.6}`$ MeV.

Qualitatively similar results for the mass are found in other models. See for instance Fig. 10 of ref. . However, we should stress that there the changes respond to a reduction of the pion decay constant ($`f`$) value which we have kept constant. The basic ingredient that drives the mass decrease in our calculation is the p-wave interaction of the pion with the baryons in the medium.

Although the $`\sigma `$ mass drops significantly, the width stays relatively large, even when the mass is close or below the $`2\pi `$ threshold. This width comes mostly from medium decay channels, namely, the decay into a single pion and a $`ph`$ excitation or $`2ph`$ excitations. Therefore we should not expect any signal of narrow $`\sigma `$ mesons in the medium as expected in  or in models that include only purely mesonic decays. Nonetheless the low mass could modify the long and medium range nucleon nucleon interaction in nuclei or nuclear matter.

## 5 Conclusions

In this work, we have searched for the lightest pole of the $`\pi \pi `$ scattering amplitude, the $`\sigma `$ meson, in the complex energy plane as a function of the baryon density.

The growing consensus that the $`\sigma `$ is rather a $`\pi \pi `$ resonance of dynamical origin than a genuine QCD state built from $`q\overline{q}`$ pairs, finds its support in recent chiral unitary formulations of the $`\pi \pi `$ interaction which we have followed in the present work. This approach has been extended to include the interaction of the pions with a nuclear medium, which allows us to trace the density dependence of the $`\sigma `$ pole.

A quite general agreement has been reached about these medium modifications of the $`\pi \pi `$ interaction between different groups, either working in the present approach or using other $`\pi \pi `$ scattering formulations, as long as they satisfy some minimal chiral constraints. All these models are, however, involved numerically to the point that extrapolating analytically the results to the complex plane is a difficult task, particularly when some pieces of the $`\pi `$ selfenergy, like the $`2p2h`$ contributions, are taken from a comparison with pionic atoms data and not from a theoretical model.

In this work, the analytical extrapolation was made possible by the use of the Lehmann representation for the $`\pi `$ propagators. This allowed us to express the scattering matrix in terms of the pion propagators evaluated on the real energy axis.

The results obtained show a dropping of the $`\sigma `$ mass as a function of the density, down to values close to two pion masses at normal nuclear density. The width decreases moderately and reaches values around 300 MeV at $`\rho =\rho _0`$. Yet, this decrease is surprising in view of the fact that in the nuclear medium there are more channels for the $`\sigma `$ decay, like pion-$`ph`$ or 2$`ph`$. The presence of these channels partly makes up for the strong reduction of the $`2\pi `$ channel due to the smaller phase space available when the $`\sigma `$ mass drops.

Altogether, the changes in the mass and width of the $`\sigma `$ are quite large, compared to other typical meson or baryon properties in a dense medium. Therefore, we expect drastic signals in the invariant mass distributions of two pions (or photons) in 2$`\pi `$ (2$`\gamma `$) production experiments if those pions (photons) can couple strongly to the scalar isoscalar channel. Some present data seem to support these claims, although further experiments and calculations should be performed to test these findings.

## 6 Acknowledgements

This work is partly supported by DGICYT contract no. BFM2000-1326.


text:
# Passive Sliders on Growing Surfaces and (anti-)Advection in Burger’s Flows

## Abstract

We study the fluctuations of particles sliding on a stochastically growing surface. This problem can be mapped to motion of passive scalars in a randomly stirred Burger’s flow. Renormalization group studies, simulations, and scaling arguments in one dimension, suggest a rich set of phenomena: If particles slide with the avalanche of growth sites (advection with the fluid), they tend to cluster and follow the surface dynamics. However, for particles sliding against the avalanche (anti-advection), we find slower diffusion dynamics, and density fluctuations with no simple relation to the underlying fluid, possibly with continuously varying exponents.

The advection of a passive scalar, such as dye, a temperature field, or dust particles, in a turbulent flow is a scientific problem of great interest, since the fluctuations of the scalar field appear to have no simple relation to the statistics of the flow siggiareview . This holds not only for fully developped turbulence, as emerging from the Navier-Stokes equation, but even for the much simpler case of a random Gaussian velocity field which has power-law correlations in space and is delta-correlated in time kraichnan .

In this article, we introduce and study a related problem of particles passively sliding on a fluctuating landscape, such as a growing surface. To exploit its mathematical similarity to the better known problem of advection of a passive scalar, we start by considering a velocity field $`\stackrel{}{v}(\stackrel{}{x},t)`$ which follows the Burger’s equation, given by

$$\frac{\stackrel{}{v}}{t}+\lambda (\stackrel{}{v})\stackrel{}{v}=\nu ^2\stackrel{}{v}+\zeta _h(\stackrel{}{x},t).$$
(1)
We also impose the additional condition

$$\stackrel{}{v}=h,$$

which ensures a gradient flow with no vortices. This condition is preserved by the random stirring with the gradient of a Gaussian white noise $`\zeta _h`$, with correlations

$$\zeta _h(\stackrel{}{x},t)\zeta _h(\stackrel{}{x}^{},t^{})=2D_h\delta ^d(\stackrel{}{x}\stackrel{}{x}^{})\delta (tt^{}).$$

The noise term also prevents the formation of shock waves. Shock waves occur in the absence of a stirring force, or with stirring correlated over large spatial scales mezard . In contrast to the Navier-Stokes equation, Burger’s equation has no pressure term, and the fluid is therefore compressible. Our use of Burger’s equation is not so much motivated by describing a particular physical situation, but because its flows are simpler than Navier-Stokes turbulence, which is a yet unsolved problem. When expressed in terms of the scalar variable $`h(\stackrel{}{x},t)`$, Eq. (1) becomes the Kardar-Parisi-Zhang (KPZ) equation KPZ  for interface growth,

$$\frac{h}{t}=\nu ^2h+\frac{\lambda }{2}(h)^2+\zeta _h(\stackrel{}{x},t),$$
(2)
where $`h(\stackrel{}{x},t)`$ is the height of interface at position $`\stackrel{}{x}`$ at time $`t`$. Fluctuations of the surface height obey the scaling relation FV85

$$\left[h(\stackrel{}{x},t)h(\stackrel{}{x}^{},t^{})\right]^2=\left|\stackrel{}{x}\stackrel{}{x}^{}\right|^{2\chi }g\left[\frac{tt^{}}{\left|\stackrel{}{x}\stackrel{}{x}^{}\right|^z}\right].$$
(3)
In one dimension, the only case that we shall consider here, the roughness exponent is $`\chi =1/2`$, and the dynamic exponent is $`z=3/2`$. These values are exact due to a fluctuation-dissipation condition, and the Galilei invariance of Eq. (1) fns .

A particle moved by the flow is characterized by the Langevin-equation

$$\frac{d\stackrel{}{x}}{dt}=a\stackrel{}{v}+\zeta _x(t),$$
(4)
where the white noise $`\zeta _x`$ describes the effects of the temperature and has a correlation function

$$\zeta _x(t)\zeta _x(t^{})=2\kappa \delta (tt^{}).$$

For $`a=\lambda `$, the particle is advected in the same way as an element of the fluid. For $`0<a<\lambda `$, the particle is slower than the neighboring fluid elements, while a negative $`a`$ describes a particle moving against the flow (anti-advection). In the language of the fluctuating interface (with $`\stackrel{}{v}=h`$), Eq. (4) describes a particle sitting on the surface and sliding with a velocity proportional to the local slope, for instance under the influence of gravitation. (This equation also describes the dynamics of a single domain wall separating two phases dro00 .) The corresponding equation for the particle density $`\rho (\stackrel{}{x},t)`$ is

$$\frac{\rho }{t}=\kappa ^2\rho +a(\rho h)+\zeta _\rho (\stackrel{}{x},t),$$
(5)
where we have again added a white noise, with

$$\zeta _\rho (\stackrel{}{x},t)\zeta _\rho (\stackrel{}{x}^{},t^{})=2D_\rho ^2\delta ^d(\stackrel{}{x}\stackrel{}{x}^{})\delta (tt^{}).$$

In this paper, we consider only noninteracting particles that move independently of each other. The case of many particles with hard-core repulsion is treated in Ref. das01 .

The difficulty in analysis of the stochastic Eqs. (2,5) arises from the non-linear couplings. A perturbative treatment of the non-linearities results in divergences that have to be controlled, e.g. by a renormalization group (RG) treatment. Such a treatment is well-known for the KPZ equation (2), and can be extended to the advection equation (5). In particular, at the one loop order, the reduced parameters $`\alpha =a/\lambda `$ and $`k=\kappa /\nu `$, satisfy the RG flows

$`{\displaystyle \frac{d\alpha }{dl}}`$ $`=`$ $`\alpha \left[z_\rho {\displaystyle \frac{3}{2}}+{\displaystyle \frac{2(\alpha \alpha ^2)}{(1+k)^2}}\right],`$
$`{\displaystyle \frac{dk}{dl}}`$ $`=`$ $`k\left[z_\rho 2+{\displaystyle \frac{\alpha ^2(3k)}{k(1+k)^2}}\right],`$
$`{\displaystyle \frac{dD_\rho }{dl}}`$ $`=`$ $`D_\rho \left[z_\rho 2\chi _\rho 1+2{\displaystyle \frac{\alpha ^2}{k(1+k)}}\right].`$ (6)
Similar to Eq. (3), the exponents $`\chi _\rho `$ and $`z_\rho `$ describe the scaling of particle density fluctuations via

$$\left[\rho (x,t)\rho (x^{},t^{})\right]^2\left|xx^{}\right|^{2\chi _\rho 2}g\left[\frac{tt^{}}{\left|\stackrel{}{x}\stackrel{}{x}^{}\right|^{z_\rho }}\right].$$
(7)
Standard dynamic RG calls for a single exponent $`z`$ for the relative scaling of time and space. The reason for including a distinct $`z_\rho `$ for the dynamics of particles shall be explained below. Other than this difference, the above RG equations are a special case of the flows derived in Ref. deniz  for a polymer drifting through a random medium.

Equations (6) show quite distinct behavior depending on the sign of $`\alpha =a/\lambda `$. While $`\alpha >0`$ is the natural choice for advection by a fluid, both signs are possible for particles sliding on a fluctuating surface. Indeed, we first encountered $`\alpha <0`$ (anti-advection) in connection with domain walls in a model of growing surfaces dro00 , inspiring the further simulations described below. For $`\alpha >0`$, the RG equations have an attractive fixed point at finite parameter values, with the standard dynamic exponent $`z=z_\rho =3/2`$. Using $`z_\rho =3/2`$, it is easy to show that the first two RG equations (6) yield the fixed point values $`\alpha ^{}=1`$ and the $`k^{}=1`$, while the third sets $`\chi _\rho =3/4`$. Note that for $`\alpha =k=1`$, the density fluctuations share the Galilei-invariance of the fluid flows, and since symmetries are preserved under renormalization, this fixed point could have been predicted without any calculations. However, this symmetry does not fix the scaling dimension $`\chi _\rho `$ (in contrast to the values of $`\chi `$ and $`z`$), and we shall in fact show below that the observed value of $`\chi _\rho =1/2`$ is different from the above one-loop prediction. The case $`a/\lambda >0`$ is the one that is usually considered in related work, as for instance in Ref. fns , which contains a paragraph on passive scalar advection in randomly stirred flows; in Ref. DE99  which discusses (among others) a situation that is equivalent to the advection of a single particle by Burger’s flow; in Ref. CSC  which focusses on the world lines of particles, and in Ref. das01 .

For $`a/\lambda <0`$ and $`z_\rho =z`$, the RG flows drive the parameters $`\alpha `$ and $`k`$ towards zero, such that to leading order $`\alpha ^2/k=1/6`$ (and $`\chi _\rho =5/12`$). The vanishing of the parameter $`a`$ under the coarse-graining action of RG suggests that on larger scales the particles are less able to perform their sliding motion. A potential interpretation of the lack of a standard fixed point is that the particles do not follow the fluctuating surface, and are instead governed by a slower timescale. (A similar vanishing of parameters under RG occurs for decoupled equations with distinct time scales.) Allowing for a distinct value of $`z_\rho `$ in the physically meaningful interval $`(3/2,2)`$, does indeed lead to nonzero fixed points for $`\alpha `$ and $`k`$. Interestingly, the resulting exponent is non-universal, varying continuously along the fixed line $`\alpha ^{}(k)`$. This conclusion appears to be supported by the simulations described below. However, we shall make no quantitative claims based on this RG study, both because of the unconventional use of two timescales, and the “strong-coupling” nature of the fixed point, which precludes controlled perturbation expansion in any dimension.

Next, we present the results of computer simulations. The growing interface is obtained from a restricted solid-on solid model, by adding “bricks” to a growing “wall”. Starting from a flat surface and imposing periodic boundary conditions, bricks are added to randomly chosen sites such that no overhangs are formed, and with the center of each brick atop the intersection of two bricks in the layer below. Such a growth rule leads to surfaces that are in the same universality class as the KPZ equation. The parameter $`\lambda `$ is negative in this model as the average growth velocity is less for a sloped surface. The advected particles are intially placed at random surface positions. Whenever a brick is added to the site of a particle, it is moved to the lower one of its two nearest-neighbor sites with probability $`(1+K)/2`$, and to the higher nearest-neighbor site with probability $`(1K)/2`$. The parameter $`K`$ is adjustable, but if not explicitly mentioned in the following, it is set to $`K=1`$. When both nearest-neighbor sites are at the same height, the particle is moved to the right or left with equal probability. The downhill character of these moves corresponds to $`a>0`$ in Eq. (4). The case $`a<0`$ is implemented similarly, moving particles uphill instead of downhill.

After a time period long enough for establishing a stationary state, the diffusion of the particles and the particle densities were quantified. Figure 1 shows the mean square distance $`x^2(t)t^{2/z_\rho }`$ travelled by a particle as function of time, for different values of the simulation parameters.

As suggested by the RG analysis, the particles can follow the surface fluctuations for $`a/\lambda >0`$, and their diffusion is characterized by the same dynamical exponent $`z_\rho =3/2`$ as the surface. For $`a/\lambda <0`$ (particles moving downhill), the diffusion exponent $`z_\rho `$ is larger, and it seems to be nonuniversal, i.e., it depends on the value of the parameter $`K`$.

This qualitative difference between the cases $`a/\lambda >0`$ and $`a/\lambda <0`$ is also visible in the correlations of particle density. To examine this quantity, we divided the system into $`L/l`$ bins of size $`l`$ and counted the number of particles in each bin. Let us choose a particle at random, and denote by $`N(l,L)`$ the number of particles that are in the same bin. This quantity is plotted in Fig. 2 as function of $`l`$, for $`a/\lambda >0`$. Using the scaling form of Eq. (7), we have

$$N(l,L)_{l/2}^{l/2}𝑑x\rho (x,t)\rho (0,t)𝑑xl^{2\chi _\rho 1}L^{22\chi _\rho },$$
(8)
where the dependence on $`L`$ follows from $`N(L,L)`$ being the total number of particles in the system, which was chosen proportional to the system size $`L`$. The curves in Fig. 2 are in fact proportional to $`L`$, and become constant for larger $`l`$, both results consistent with $`\chi _\rho =1/2`$ in Eq. (8). This value is different from the 1-loop RG result obtained above.

The corresponding results for $`a/\lambda <0`$ are shown in Fig. 3 (for $`K=1`$). Since particles cannot fully follow the surface fluctuations, the proportion of particles within a given interval around a randomly chosen particle increases with increasing system size, and we measured a nontrivial scaling exponent $`\chi _\rho 0.85`$. For smaller values of $`K`$, the scaling behavior is not so clear. The asymptotic scaling of the factor that maps curves for different $`L`$ onto each other could not be determined reliably. It appears, however, that $`\chi _\rho `$ is nonuniversal, becoming larger with decreasing $`K`$ (and increasing particle diffusion exponent $`z_\rho `$).

Particle clustering, which is the most prominent feature of our system, does not occur in Navier-Stokes turbulence, where the fluid is incompressible. Instead, fluid elements are streched and folded, leading ultimately to a homogeneous particle distribution, unless new particles are continuously added.

In order to better understand the nonuniversal and nontrivial scaling for $`a/\lambda <0`$, we measured directly the velocity of particles as a function of an external slope. A particular slope was imposed on the system by choosing periodic boundary conditions with a fixed height difference between the two ends. The average speed of particles was then evaluated as function of slope and system size. As expected, for $`a/\lambda >0`$ the mean velocity is simply proportional to the slope. For $`a/\lambda <0`$, however, the velocity shows a nontrivial scaling behavior of the form

$$v=L^yh𝒞(h\sqrt{L}),$$
(9)
with $`y0.14`$, as shown in Figure 4. The scaling function $`𝒞(x)`$ is constant for small $`x`$ and increases as $`x^{2y}`$ for large $`x`$, such that $`v(h)^{1+2y}`$ is independent of $`L`$.

(Note that the data for $`L=1024`$ show large fluctuations and seem sensitive to the choice of the random number generator, and therefore not all data points collapse with those for smaller $`L`$.) Just as for the density correlations, the data for smaller values of $`K`$ did not show nice scaling behavior for the system sizes in our simulations. Equation (9) implies that the factor $`a`$ relating velocity to slope \[see Eq. (4)\] decreases as a function of system size, consistent with the RG prediction. However, while the 1-loop result suggests an inverse logarithmic dependence of $`a`$ on the length scale (implying $`y=0`$), the value of $`y`$ found in the simulations is nonzero.

The form of Eq. (9) can be justified by the following argument: The growth rules for the surface imply that sequences of brick addition usually move from local minima in the uphill direction. The reason is that the addition of a brick generates a potential growth site (where a brick can be added without generating overhangs) at the nearest uphill position. The particles that try to slide downhill are therefore faced with an upward avalanche of growth mounds of different sizes that hamper their downhill motion. (This is analogous to a moving fluid density fluctuation faced by a particle advecting against the flow.) The horizontal size of the largest mounds, $`l`$, scales as $`L\stackrel{~}{𝒞}(h\sqrt{L})`$, with $`\stackrel{~}{𝒞}(x)`$ being constant for small $`x`$ and $`\stackrel{~}{𝒞}(x)x^2`$ for large $`x`$. If we assume that the particle velocity is slowed down by some power of the size of the largest mounds, we obtain Eq. (9). The exponent $`y`$ may then have a non-universal dependence on the manner in which the particle goes over these barriers.

Finally, we shall try to compute the exponent $`y`$ by relating characteristic time scales of surface and particle fluctuations. Consider a portion of the surface of size $`l`$, with a typical slope of order $`|h|l^{1/2}`$. For this slope to change by surface fluctuations, we should wait a time of the order $`tl^zl^{3/2}`$. During this time the advecting particle will move against the slope by an amount

$$xa(l)l^{1/2}tt^{2(1y)/3},$$

leading to the scaling relation

$$z_\rho =\frac{3}{2(1y)}.$$
(10)
For $`y0.14`$, this gives $`z_\rho 1.74`$, which is the value measured in simulations of particle diffusion for $`K=1`$.

In conclusion, we have found surprisingly intricate scaling phenomena for passive sliders on stochastically growing surfaces, and the equivalent problem of passive scalar motion in Burger’s turbulence. In particular, there is a qualitative difference between the cases $`a/\lambda >0`$ and $`a/\lambda <0`$, describing respectively passive advection with or against the flow. The former leads to clustering, while the latter yields non-trivial exponents. This conclusion is bolstered by RG analysis, simulations, and scaling arguments. While the simulations provide convincing indication of a nonuniversal value of the exponent $`z_\rho `$ characterizing particle diffusion(for $`a/\lambda <0`$), the evidence is less compelling for as the density correlations. The absence of nice scaling for most values of simulation parameters may indicate slow crossovers, but then also leaves open the possibility that the exponent $`z_\rho `$ is universal even in the case $`a/\lambda <0`$.

###### Acknowledgements.

We thank Kay Wiese for helpful discussions. B.D. was supported by the Deutsche Forschungsgemeinschaft (DFG) under Contract No Dr300-2/1. M.K. is supported by the NSF through grant number DMR-01-18213.