text:
# Superconducting Gap, Pseudogap and Their Relationship Investigated by Break-junction Tunneling Spectroscopy on 𝐵𝑖₂𝑆𝑟₂𝐶𝑎𝐶𝑢₂𝑂_{8+𝛿} Single Crystals

## Abstract

Temperature and doping dependent tunneling spectroscopy on $`Bi_2Sr_2CaCu_2O_{8+\delta }`$ single crystals has been investigated by using break-junction technique. The results provide direct evidence that the pseudogap persists into heavily overdoped regime, possibly in whole doping range of the superconducting phase. The gap value determined from the peak-to-peak of the conductance, $`\mathrm{\Delta }_{pp}`$, exhibits a jump-like increase as $`T`$ is raised across the critical temperature. The magnitude of the jump, the pseudogap opening temperature and superconducting gap decrease with raising doping. The $`\mathrm{\Delta }_{pp}`$ below $`T_c`$ can be naturally treated as the vector sum of two competing orders: a BCS-like superconducting gap and a residual pseudogap.

It is generally agreed that for underdoped cuprate high-$`T_c`$ superconductors a generic phenomenon is there is the pseudogap (PG) in their excitation spectrum for temperature ($`T`$) above the critical temperature $`T_c`$. The PG has been observed with various probes both spin and charge sensitive. Among the probes tunneling spectroscopy, one of the direct measures of quasiparticle excitations in a superconductor has attracted great deal of attention. However, thus far the results of tunneling, are still controversial. Different types of relation between the SG and the PG have been reported, such as smooth evolution, competing, coexistence of the two gaps, and distinct two gaps, with scanning tunneling microscope (STM), break-junction, and c-axis intrinsic junction. There is also the question that whether the PG exists in overdoped regime, especially in heavily overdoped region? Apparently, experimentally clarifying the above discrepancies is a pressing task and will set constraint to any theoretical models on the origin of the pseudogap and on the mechanism of high-$`T_c`$ superconductivity.

In this paper we present systematic investigation on the detailed evolution behavior of the gap $`\mathrm{\Delta }_{pp}`$, determined by a quarter of the span from the peak-to-peak in the tunneling conductance, on $`Bi_2Sr_2CaCu_2O_{8+\delta }`$ (Bi2212) single crystals with break-junction technique, specifically focusing on the overdoped regime. For all samples studied, including the highly overdoped ($`T_c`$ = 72 K), our low-noise tunneling conductance clearly demonstrates the existence of the pseudogap. Detailed $`T`$-dependent measurements show a jump-like increase of $`\mathrm{\Delta }_{pp}`$ as $`T`$ is raised across $`T_c`$ within a very narrow temperature range, indicating an unsmooth transformation from the superconducting state to the PG state. We show that the measured gap for $`T`$ below $`T_c`$ can be naturally treated as a vector sum of two competing components: the superconducting gap (SG) and the residual PG. The PG opening temperature, $`T^{}`$, the magnitude of the jump of $`\mathrm{\Delta }_{pp}`$ at $`T_c`$, $`\mathrm{\Lambda }`$, and the value of the SG, and PG decrease with increase of doping. The doping dependence of $`T^{}`$ implies that for Bi2212 the PG may exist in the whole overdoped regime provided there is the superconductivity. To our knowledge, this doping dependence of $`\mathrm{\Lambda }`$, as a characteristic parameter for the unsmooth transformation is newly observed by tunneling, and the present work is the first directly tunneling observation of the doping dependence of the $`T^{}`$.

We chose break-junction to carry out our tunneling investigation because it provides a clean and fresh interface easily, which avoids surface degradation and diminishes the barrier scattering. In addition, the break-junction is of the superconductor-insulator-superconductor (SIS) type junction, which enhances the gap feature in the tunneling conductance due to its unique convolution relationship from the electronic density of states (DOS) of two superconductors on both sides of tunnel barrier. Therefore, the peak position is less affected by thermal smearing and can be taken as a characteristic parameter. More importantly, previous studies at 4.2 K have provided reproducible and low-noise tunneling conductance consisting of sharp coherence peaks, a close to zero zero-bias conductance, and a significant peak-dip-hump structure. Those facts demonstrate that the break-junction technique is a reliable tool with high-energy resolution to investigate the tunneling spectroscopy.

The Bi2212 single crystals used in this study were of three doping levels. The near-optimally doped samples were the as-grown crystals grown by traveling-solvent-floating-zone method. The medium-overdoped samples were achieved by annealing the as-grown crystals ($`T_c`$ = 89 K), synthesized by directional solidification method, in flowing O2 at 500C for 120 hours. And the heavily overdoped samples, also from the as-grown Bi2212 single crystals grown by directional solidification, were obtained by annealing them under a pressure of 140 atm of O2 at 500C for 50 hours. The on-set $`T_c`$ of the three batches of samples was 90 K, 82 K and 72 K with transition width of 1 K, 1 K and 3 K, respectively, measured by AC susceptibility using Lakeshore AC Susceptometer-7000. Hereafter we label the samples as OP90, OV82 and OV72. All the characterizations were carried out on segments of the very single crystals used to measure tunneling spectroscopy. The tunneling conductance curve, $`dI/dVV`$, was measured by using a standard low-frequency ac lock-in technique with the four-probe method and a computer data acquisition system. To achieve low-resistive contacts ($`<`$ 2 $`\mathrm{\Omega }`$) between the crystals and the leads, gold film electrodes were sputtered onto each sample.

The break-junctions were formed in situ at 4.2 K. To realize it, thin (less than 0.1 mm) single crystal was glued onto phosphor-copper substrate with epoxy, and broken by applying a bending force perpendicular to the CuO2 plane with a precise differential-screw. After a break-junction was formed, a key problem left was to keep it from being unstable as temperature varied. Because there is mismatch between the thermal expansion coefficients among the Bi2212 crystal, epoxy and phosphor-copper substrate, which often causes the contact of two pieces of the broken off crystal become loosen or tighten, and leads to the junction conductance changed drastically when the temperature varies. It even breaks down the junction at high temperature. Through practice we found that by adjusting the thickness of epoxy, which depends on the dimensions of the substrate, we were able to overcome this problem and perform stable $`T`$-dependent investigation.

To diminish the complexity bringing about by Josephson current and to facilitate the analysis of experimental data, each time we carefully adjusted the differential screw to increase the junction resistance and depress the Cooper pair tunneling to undetectable small. Tunneling measurements were carried out at temperatures ranging from near 4.2 K to where the PG disappears.

Figure 1, 2 and 3 show the normalized tunneling conductance of different temperatures measured on OP90, OV82 and OV72, respectively. The normalization was done by divided each measured conductance by its conductance of high bias voltage. For clarity, the curves of $`TT_c`$ are sequentially shifted up with a unit value (in (a) panel of each figure), whilst those of $`TT_c`$ are plotted with same coordinate scale as they were measured (in (b) panel of each figure).

In Fig. 1 to 3 several points are notable. (1) At low $`T`$, the conductance for OP90 and OV82 exhibit the significant dip-hump structure, which has been observed in ARPES and tunneling measurements and is related to the ($`\pi `$, $`\pi `$) resonance in inelastic neutron scattering. The dip feature for OV72 is not as significant as of OP90 and OV82, but it is still visible in Fig. (3.c) by comparing the low-$`T`$ conductance to the ones of higher temperature. Those facts demonstrate the high quality of our junctions. (2) For $`TT_c`$, such as $`T`$ = 4.2 K, it is generally safe to take $`\mathrm{\Delta }_{pp}`$ as the SG, $`\mathrm{\Delta }`$(0). The low-T parameters, including $`\mathrm{\Delta }`$(0), and the dip position $`E_{dip}`$ in tunneling conductance for the three samples are listed in Table I. These values are consistent well with the results of Ref. , in which a detailed examination over a wide range of oxygen doping at 4.2 K was performed. On the other hand, the tunneling curve calculated from a weak-coupling d-wave DOS has significant deviation from the experimental date in the shape of conductance peak. This brings about difficulty to model fit the experimental curves even in the low $`T`$. (3) In Fig. 1, 2 and 3, the conductance curves at $`T>T_c`$ clearly demonstrate the existence of PG even in heavily overdoping level with the presence of peaks, and the low-noise background enables us to get the value of $`\mathrm{\Delta }_{pp}(T)`$, a quarter of the peak-to-peak span. (4) The opening temperature of the PG, $`T^{}`$, can be determined as the temperature at which the conductance peaks vanish, meanwhile the central dip in the conductance is filled-in. By so doing, the obtained $`T^{}`$ is 190, 150 and 120 K for sample OP90, OV82 and OV72, respectively, and we estimate the error bars are $`\pm `$15, $`\pm `$15, $`\pm `$10 K, respectively.

To gain more insights into the physics, we depict in Fig. 4 the $`\mathrm{\Delta }_{pp}(T)`$ taken from Fig. 1, 2 and 3, respectively, as a function of the normalized temperature $`T/T_c`$,. The filled squares, deltas and circles is corresponding to sample OP90, OV82 and OV72, respectively. The solid lines are the BCS $`T`$-dependent gaps corresponding to measured $`\mathrm{\Delta }`$(0) of OP90, OV82 and OV72, respectively. A general feature for the three junctions in Fig. 4 is that for $`2/3T_c<T<T_c`$, the $`\mathrm{\Delta }_{pp}(T)`$ decreases with increasing $`T`$ much slower than the BCS behavior. There are at least two factors that may cause this slow decrease. First, the life-time effect of the quasuiparticles, which smears the tunneling conductance and gives rise to a larger $`\mathrm{\Delta }_{pp}(T)`$ than the SG, $`\mathrm{\Delta }_s(T)`$. Second, the measured single gap is a vector sum of the SG and the residual PG for $`T`$ below $`T_c`$: $`\mathrm{\Delta }_{pp}(T)=(\mathrm{\Delta }_s^2(T)+\mathrm{\Delta }_p^2(T))^{1/2}`$. For simplicity, if we do not take into account the life-time broadening effect, and subtract the BCS gap from the $`\mathrm{\Delta }_{pp}(T)`$ according to the above formula we get the residual PG for $`T<T_c`$, which is shown in the inset of Fig. 4. The data for $`TT_c`$ in the inset are the same as those in the main part of Fig. 4. The vector sum scenario is suggested by a number of theoretical models stating that in superconducting state there exist two competing orders: the SG and the PG. As temperature is raised across $`T_c`$, $`\mathrm{\Delta }_{pp}`$ undergoes an abrupt jump within a very narrow temperature range ($``$2 K in OP90) and the $`T_c`$ can be regarded as where the maximum of $`d\mathrm{\Delta }_{pp}(T)/dT`$ is, or the middle of the jump. This behavior demonstrates that the maximum of DOS moves to higher energy unsmoothly as soon as long-range phase coherence disappears. To characterize this unsmooth behavior, we introduce $`\mathrm{\Lambda }`$, the jump in $`\mathrm{\Delta }_{pp}`$ within the measured temperature interval in the vicinity of $`T_c`$. It can be seen that $`\mathrm{\Lambda }`$ has the same trend with doping as $`T^{}`$ (see also Tab. I).

Here one point we would like to emphasis is that the measured gap $`\mathrm{\Delta }_{pp}`$ in the normal state is not the true magnitude of the PG, because apparently the value of the PG should decrease or keep as constant when $`T`$ is raised. In general, high fluctuation or strong pair-breaking can bring out broadening conductance peak and widen separation between the peaks in the conductance. As to how to extract the true value of the PG and its temperature behavior from tunneling date is an open issue remaining to be address. The fact that for OP90 $`\mathrm{\Delta }_{pp}`$ keeps increasing as $`T`$ is raised reflects that in the normal state there is high fluctuation or strong pair-breaking which enhance as $`T`$ is raised. Moreover, the trend of doping dependence of $`\mathrm{\Delta }_{pp}(T)`$ demonstrates that the fluctuation or pair-breaking effect become weak in high doping region. Those features are consistent with the picture of precursor paring in electronic or spin channel in the normal state but we emphasize it exists in the overdoped regime, and directly show that the difference in electronic excitation between the superconducting state and the normal state near $`T_c`$ increases with lowering doping.

To further illustrate the present observation we plot the doping dependence of the $`\mathrm{\Delta }(0)`$ and $`T^{}`$ in Fig. 5, in which the dotted line is the BCS mean-field gap, $`\mathrm{\Delta }=2.14k_BT_c`$, the hole concentration $`p`$ is obtained from the empirical relation $`T_c/T_{c,max}=182.6(p0.16)^2`$, with $`T_{c,max}=92K`$. The filled circles are measured $`\mathrm{\Delta }`$(0) and the filled squares are the measured $`2.14k_BT^{}`$. The thick straight line is a guide to the eye. It is surprising that the extrapolation of the thick line just crosses to the abscissa at $`p`$ = 0.27, where the superconductivity disappears. That implies so long as it is a superconductor, there may be the PG state. A similar trend for Bi2212 can be found in Ref.. This phase diagram of Bi2212 is qualitatively different from that of YBa2Cu3O7-δ, to which it is shown the PG state only exists in the region where the hole concentration $`p`$ is less than a certain value, 0.19. Obviously this is an issue that needs to be further addressed.

To conclude, we have performed investigations on temperature and doping dependent tunneling spectroscopy of Bi2212 single crystals with break-junction technique, focusing on overdoped regime. The results provide direct evidence that for Bi2212 the PG persists into heavily overdoped regime, possibly in whole doping range accompanying the SG. Detailed measurements near $`T_c`$ have showed the measured $`\mathrm{\Delta }_{pp}`$ undergoes a jump-like increase within a very narrow temperature range as $`T`$ is raised across $`T_c`$, indicating an unsmooth transformation from the superconducting state to the PG state. The magnitude of the jump of $`\mathrm{\Delta }_{pp}`$, $`T^{}`$, the SG, and PG decrease with raising doping. We show the measure gap $`\mathrm{\Delta }_{pp}(T)`$ for $`T`$ below $`T_c`$ can be naturally treated as the vector sum of two competing orders: a BCS-like SG and a residual PG.

###### Acknowledgements.

This work was supported by a grant for State Key Program for Basic Research of China, and the National Center for Research and Development on Superconductivity.


format: text

text:
# REFERENCES

One of the most natural candidates for a symmetry of nature is parity (i.e. left-right) symmetry. While it is an established experimental fact that parity symmetry appears broken by the interactions of the known elementary particles, this however does not exclude the possible existence of exact unbroken parity symmetry in nature. This is because parity (and also time reversal) can be exactly conserved if a set of mirror particles exist. The idea is that for each ordinary particle, such as the photon, electron, proton and neutron, there is a corresponding mirror particle, of exactly the same mass as the ordinary particle. For example, the mirror proton and the ordinary proton have exactly the same mass. Furthermore the mirror proton is stable for the same reason that the ordinary proton is stable, and that is, the interactions of the mirror particles conserve a mirror baryon number. The mirror particles are not produced (significantly) in Laboratory experiments just because they couple very weakly to the ordinary particles. In the modern language of gauge theories, the mirror particles are all singlets under the standard $`GSU(3)SU(2)_LU(1)_Y`$ gauge interactions. Instead the mirror fermions interact with a set of mirror gauge particles, so that the gauge symmetry of the theory is doubled, i.e. $`GG`$ (the ordinary particles are, of course, singlets under the mirror gauge symmetry). Parity is conserved because the mirror fermions experience $`V+A`$ mirror weak interactions and the ordinary fermions experience the usual $`VA`$ weak interactions. Ordinary and mirror particles interact with each other predominately by gravity only.

At the present time there is a large range of experimental evidence supporting the existence of mirror matter (for a review see Ref.). Mirror matter is necessarily stable and dark and appears to provide a viable candidate for the inferred dark matter in the Universe as well as having important implications for early Universe cosmology. Mirror dark matter also has self interactions just like ordinary matter which may allow it to escape the fate of collisionless cold dark matter candidates such as hypothetical neutralinos which now appear to be ruled out by the observations. Moreover, mirror matter, like ordinary matter can form stars, planets and smaller bodies and there is interesting evidence for all these things. In particular mirror stars are a natural candidate for the observed MACHO gravitational microlensing events. Furthermore mirror planets would provide a simple explanation for the existence of close-in extrasolar planets which has been puzzling astronomers since their unexpected discovery in 1995. There is also evidence that the ‘dynamical mirror image’ system of an ordinary planet orbiting a mirror star has also been observed but interpreted as an ‘isolated’ planet because light from the mirror star was not detected.

The significance of mirror matter for astrophysics and cosmology is clear, perhaps of equal importance though is the implications of mirror matter for particle physics. While ordinary and mirror matter interacts with each other predominately by gravity, small non-gravitational interactions are actually possible. Due to constraints from gauge symmetry, renormalizability and parity symmetry it turns out that there are only 3 ways in which ordinary and mirror matter can interact with each other (besides gravity). This is via photon - mirror photon kinetic mixing, Higgs - mirror Higgs interactions and via ordinary neutrino - mirror neutrino mass mixing (if neutrinos have mass). While Higgs - mirror Higgs interactions will be tested if or when the Higgs particle is discovered, there is currently strong evidence for photon - mirror photon kinetic mixing and also ordinary neutrino - mirror neutrino mass mixing.

A simple consequence of the parity symmetry is that each of the ordinary neutrinos ($`\nu `$) will oscillate maximally into its mirror partner ($`\nu ^{}`$). This provides a very elegant explanation for the solar neutrino puzzle since the maximal $`\nu _e\nu _e^{}`$ oscillations imply an approximate 50% flux reduction for a large range of $`\delta m^2`$ which is in broad agreement with the solar neutrino data. Moreover this solution predicted the approximate energy independent recoil electron energy spectrum observed by super-Kamiokande as well as the $`50\%`$ flux reduction found in the Gallium experiments. In the case of the atmospheric neutrino anomaly the inferred 50% reduction of up-going $`\nu _\mu `$ is also nicely explained by maximal $`\nu _\mu \nu _\mu ^{}`$ oscillations. If the solar and atmospheric neutrino anomalies are due to oscillations into mirror neutrinos then oscillations between generations can be governed by small mixing angles which seems theoretically most natural. This reasoning is supported by the LSND experiment which has provided strong evidence for small angle $`\nu _e\nu _\mu `$ oscillations.

It is true, though, that the solution to the neutrino physics anomalies implied by the mirror matter theory does not give a perfect fit to every neutrino experiment. However this is probably a good thing, since it is unlikely that every experimental measurement is correct. In the case of solar neutrinos, the low Homestake result (1/3 c.f. 1/2 in the 5 other solar neutrino experiments) and also the recent SNO results do not favour the simplest mirror matter solution. In addition the atmospheric data slightly prefer $`\nu _\mu \nu _\tau `$ to $`\nu _\mu \nu _\mu ^{}`$ (although the extent to which $`\nu _\mu \nu _\mu ^{}`$ is disfavoured depends significantly on how the data is analysed). Because these disfavouring results are only at the 1.5-3.3 sigma level (and are largely dominated by systematics) they do not provide a strong case against the mirror matter theory. Importantly things will eventually become clear as more accurate measurements are done. The forthcoming $`NC/CC`$ SNO measurement should provide a solid result one way or the other.

Another important way that ordinary and mirror matter can interact with each other is via photon - mirror photon kinetic mixing. In field theory this is described by the interaction

$$=\frac{ϵ}{2}F^{\mu \nu }F_{\mu \nu }^{},$$
(1)
where $`F^{\mu \nu }`$ ($`F_{\mu \nu }^{}`$) is the field strength tensor for electromagnetism (mirror electromagnetism). This type of Lagrangian term is gauge invariant and renormalizable and can exist at tree level or maybe induced radiatively in models without $`U(1)`$ gauge symmetries (such as grand unified theories). One effect of ordinary photon - mirror photon kinetic mixing is to give the mirror charged particles a small electric charge. That is, they couple to ordinary photons with electric charge $`ϵe`$.

The most important experimental constraint on photon - mirror photon kinetic mixing is that it modifies the properties of orthopositronium. This effect arises due to radiative off-diagonal contributions to the orthopositronium, mirror orthopositronium mass matrix. This means that orthopositronium oscillates into its mirror partner. Decays of mirror orthopositronium are not detected experimentally which effectively increases the observed decay rate. Because collisions of orthopositronium destroy the quantum coherence, this mirror world effect is most important for experiments which are designed such that the collision rate of the orthopositronium is low. The only accurate experiment sensitive to the mirror world effect is the Ann Arbour vacuum cavity experiment. This experiment obtained a decay rate of $`\mathrm{\Gamma }_{oPs}=7.0482\pm 0.0016\mu s^1`$. Normalizing this measured value with the recent theoretical value of $`7.0399\mu s^1`$  gives

$$\frac{\mathrm{\Gamma }_{oPs}(exp)}{\mathrm{\Gamma }_{oPs}(theory)}=1.0012\pm 0.00023$$
(2)
which is a five sigma discrepancy with theory. It suggests a value $`ϵ10^6`$ for the photon - mirror photon kinetic mixing. Taken at face value this experiment is strong evidence for the existence of mirror matter and hence parity symmetry. It is ironic that the last time something important was discovered in high energy physics with a table top experiment was in 1957 where it was demonstrated that the ordinary particles by themselves appear to violate parity symmetry.

Of course this vacuum cavity experiment must be carefully checked by another experiment to make sure that mirror matter really exists. Actually this is quite easy to do. With the largest cavity used in the experiment of Ref. the orthopositronium typically collided with the cavity walls 3 times before decaying. If the experiment was repeated with a larger cavity then the mirror world effect would be larger because the decohering effect of collisions would be reduced. For example if a cavity 3 times larger could be used (which means that the orthopositronium would typically collide with the walls just once before decaying) then the mirror world would predict an effect 3 times larger.

There are several important implications of photon - mirror photon kinetic mixing with the relatively large value of $`ϵ10^6`$ suggested by the orthopositronium vacuum experiment. These include:

* Exploding mirror stars (mirror supernova) will emit a burst of (ordinary) gamma rays. This would occur because at the temperatures $`10\mathrm{MeV}`$ reached at the center of a typical supernova explosion the kinetic mixing will convert $`e^+e^{}e^+e^{}`$ which subsequently produces a relativistic fireball, which seems to qualitatively explain many of the features of the observed gamma ray bursts.
* Such a large value of $`ϵ10^6`$ will lead to the light mirror particles ($`e^\pm ,\gamma ^{},\nu ^{}`$) being brought into equilibrium with the ordinary particles above $`T=1\mathrm{MeV}`$ in the early Universe. While this is not a problem for the recent measurements of the Cosmic Microwave Background, it does suggest that standard BBN needs modification. For example, there might exist a large electron neutrino asymmetry which can compensate for the faster expansion rate leading to acceptable values of the light element abundances. Another possibility is that there might exist a large negative cosmological constant which will slow down the expansion rate at $`T1\mathrm{MeV}`$ .
* Mirror stars can become visible if they have some embedded ordinary matter. This is because the ordinary matter is heated by the mirror matter though photon - mirror photon kinetic mixing. Maybe the recently observed halo white dwarfs  (which are controversial) are really mirror stars or even mirror white dwarfs. Because of their age they may have accreted enough ordinary matter to be observable.

Perhaps the most remarkable possibility though is that there is some significant amount of mirror matter in our solar system. We don’t know enough about the formation of the solar system to be able to exclude the existence of a large number of Space Bodies (SB) made of mirror matter if they are small like comets and asteroids. The total mass of asteroids in the asteroid belt is estimated to be only about 0.05% of the mass of the Earth. A similar or even greater number of mirror bodies, perhaps orbiting in a different plane or even spherically distributed like the Oort cloud is a fascinating and potentially explosive possibility Large planetary sized bodies are also possible if they are in distant orbits. if they collide with the Earth. The possibility that such collisions occur and may be responsible for the 1908 Siberian explosion (Tunguska event) has been speculated in Ref.. The purpose of this paper is to study this possibility in detail and to point out the important ramifications of this idea which is that mirror matter should be present in the ground at the ‘impact’ sites and could be extracted as we will discuss.

If such small mirror bodies exist in our solar system and happen to collide with the Earth, what would be the consequences? If the only force connecting mirror matter with ordinary matter is gravity, then the consequences would be minimal. The mirror SB would simply pass through the Earth and nobody would know about it unless it was so heavy as to gravitationally affect the motion of the Earth. However if there is photon - mirror photon kinetic mixing as suggested by the orthopositronium vacuum cavity experiment, then the mirror nuclei (with $`Z^{}`$ mirror protons) will effectively have a small ordinary electric charge $`ϵZ^{}e`$. This means that the nuclei of the mirror atoms of the SB will undergo Rutherford scattering off the nuclei of the atmospheric nitrogen and oxygen atoms. In addition ionizing interactions can occur which can ionize both the mirror atoms of the space body and also the atmospheric atoms. The net effect is that the kinetic energy of the SB is transformed into light and heat (both ordinary an mirror varieties) and a component is also converted to the atmosphere in the form of a shockwave, as the forward momentum of the SB is transferred to the air which passes though or near the SB.

What happens to the mirror matter SB as it plummets towards the Earth’s surface depends on a number of factors such as its initial velocity, size, chemical composition and angle of trajectory. Of course all these uncertainties occur for an ordinary matter SB too. Interestingly it turns out that for the value of the kinetic mixing suggested by the Orthopositronium experiment, $`ϵ10^6`$, the air resistence of a mirror SB in the atmosphere is roughly the same as an ordinary SB assuming the same trajectory, velocity mass, size and shape (and that it remains intact). This occurs because the air molecules will lose their relative forward momentum (with respect to the SB) within the SB itself because of the Rutherford scattering of the ordinary and mirror nuclei as we will show in a moment. (Of course the atmospheric atoms still have random thermal motion). This will lead to a drag force of roughly the same size as that on an ordinary matter SB, implying an energy loss rate of

$$\frac{dE}{dx}=C_d\rho _{air}A\frac{v^2}{2},$$
(3)
where $`\rho _{air}`$ is the density of the air, $`v`$ is the velocity of the SB and $`A`$ is the cross sectional area. The drag coefficient, $`C_d`$ is of order unity - its precise value depending on the shape of the body. We will take $`C_d1`$. Eq.(3) is a standard result and quite easy to derive: The pressure of the atmosphere on the surface of the body increases linearly with the velocity of the body. Also the number of atoms striking the surface will increase linearly with the air density and also velocity (since the volume that the body sweeps out in a given time $`t`$ is just $`Avt`$). Eq.(3) implies that the bodies velocity decreases exponentially with distance ($`x`$),

$$v=v_ie^{x/D},$$
(4)
where $`v_i`$ is its initial velocity and

$`D={\displaystyle \frac{2R\rho _{SB}}{C_d\overline{\rho }_{air}}}10\left({\displaystyle \frac{R}{5\mathrm{meters}}}\right)\left({\displaystyle \frac{\rho _{SB}}{1g/cm^3}}\right)\mathrm{km}.`$ (5)
In this equation, $`\rho _{SB}`$ is the density of the SB and $`RV/A`$ is the ‘size’ of the body ($`V`$ is its volume). Note that we have used $`\overline{\rho }_{air}10^3g/cm^3`$ which is the air density at about 5 km altitude (the density at sea level is about twice this value) for a rough estimate of the mean density encountered as it travels through the atmosphere. The above calculation shows that the rate of energy loss of the SB in the atmosphere depends on its size and density. If we assume a density of $`\rho _{SB}1g/cm^3`$ which is approximately valid for a mirror SB made of cometary material (such as mirror ices of water, methane and/or ammonia) then the body will lose most of its kinetic energy in the atmosphere provided that it is less than roughly 5 meters in diameter. Of course things are complicated because the the SB will undergo mass loss (ablation) and also potentially fragment into smaller pieces and of course potentially melt & vaporize. Thus even a very large body (e.g. $`R100`$ meters as estimated for the Tunguska explosion) can lose its kinetic energy in the atmosphere if it fragments into small pieces.

An important difference between an ordinary and mirror SB is the rate and way in which it fragments, heats up and undergoes ablation because these properties depend very much on the interactions between the SB and the atmosphere. An ordinary matter SB undergoes huge pressure on its surface when it enters the atmosphere with cosmic velocity ($`30km/s`$) while in the case of a mirror matter body the effects of the pressure are distributed within the body to some extent, rather than just at the very surface. Let us now examine this in more detail.

Assume that the mirror matter SB is composed of atoms of mass $`M_A^{}`$ and the air is composed of atoms of mass $`M_A`$. The (mirror) electric charge in units of $`e`$ of the (mirror) nuclei, which we roughly assume to be half neutrons and half protons, will be $`Z=M_A/2M_P`$ ($`Z^{}=M_A^{}/2M_P`$), where $`M_P`$ is the proton mass. Let us assume that the trajectory of the SB is a straight line along the $`\widehat{z}`$ axis of our co-ordinate system. In the rest frame of the SB, the change in forward momentum of each of the on-coming atmospheric atoms is then The following equation is valid provided that $`M_A^{}M_A`$ but our conclusions will remain roughly the same for other cases of interest such as for $`M_A^{}M_A`$.

$$\frac{dP_z}{dt}=\mathrm{\Gamma }_{coll}M_A(v\mathrm{cos}\theta v)=2\mathrm{\Gamma }_{coll}M_Av\mathrm{sin}^2\frac{\theta }{2},$$
(6)
where $`\theta `$ is the scattering angle in the rest frame of the SB and $`\mathrm{\Gamma }_{coll}`$ is the collision rate of the atmospheric atom with the mirror atoms in the SB. Of course the collisions also generate transverse momentum (i.e. in the $`\widehat{x},\widehat{y}`$ directions) which is reduced by thermalization effects as the atoms in the atmosphere interact with themselves. For the present calculation we are only interested in the relative net momentum between the SB and the atmosphere and we can neglect this transverse motion in a rough approximation (which means that we can replace $`v`$ by $`v_z`$ below). The collision rate $`\mathrm{\Gamma }_{coll}`$ is given in terms of the cross section, relative velocity and number density in the usual way:

$$\mathrm{\Gamma }_{coll}=\sigma v_z\left(\frac{\rho _{SB}}{M_A^{}}\right).$$
(7)
Thus Eq.(6) becomes

$$\frac{dP_z}{dt}=2\left(\frac{M_A}{M_A^{}}\right)\frac{d\sigma }{d\mathrm{\Omega }}\rho _{SB}v_z^2\mathrm{sin}^2\frac{\theta }{2}d\mathrm{\Omega }.$$
(8)
There are various different processes which can contribute to the scattering cross section. For the velocities of interest, $`v\stackrel{<}{}70km/s`$, the cross section is dominated by Rutherford scattering§§§ Although the cross section is dominated by Rutherford scattering, ionizing collisions may also be important for generating light and perhaps may also allow the body to build up electric charge within. of the mirror nuclei of effective electric charge $`ϵZ^{}e`$ off the ordinary nuclei of electric charge $`Ze`$ modified for small angle scattering by the screening effects of the atomic electrons (at roughly the Bohr radius $`r_010^8\mathrm{cm}`$). It is given by (see e.g. ) We use standard particle physics units $`h/(2\pi )=c=1`$ unless otherwise stated. :

$$\frac{d\sigma }{d\mathrm{\Omega }}=\frac{4M_A^2ϵ^2e^4Z^2Z^2}{(4M_A^2v_z^2\mathrm{sin}^2\frac{\theta }{2}+1/r_0^2)^2}.$$
(9)
Thus we obtain from Eq.(8) and Eq.(9) the following differential equation for the distance travelled by each atmospheric atom ($`z`$) within the SB :

$$\frac{dP_z}{dt}=M_Av_z\frac{dv_z}{dz}Z^2Z^2\rho _{SB}\frac{ϵ^2e^44\pi }{M_A^{}M_Av_z^2}log_e\left(\frac{1}{M_Av_zr_0}\right),$$
(10)
which is valid for $`M_Avr_01`$. For $`M_A15M_P`$, $`M_Avr_050(v/30km/s)`$ which means that the above equation is approximately valid for the velocities of interest (the initial velocity, $`v_i`$, of a SB is typically between 15 and 60 km/s). Solving the above differential equation (neglecting the log factor which is of order 1) we find that the relative motion between the air molecules and SB is lost (upto random thermal motion) after travelling a distance within the SB of

$$z\frac{v^4M_A^2M_A^{}}{16\pi Z^2Z^2\rho _{SB}ϵ^2e^4}\left(\frac{10^6}{ϵ}\right)^2\left(\frac{v}{30km/s}\right)^4\mathrm{centimeters},$$
(11)
where we assumed $`\rho _{SB}1g/cm^3`$ and $`M_AM_A^{}15M_P`$ (with $`ZZ^{}7`$). For $`ϵ=10^6`$, Eq.(11) indicates that the atmospheric atoms lose essentially all of their relative momentum (of course they still have thermal motion) after penetrating a distance of the order of a few centimeters into the SB. (This distance may be somewhat greater for a body made of a heavy element such as mirror iron). If the SB remains intact then the above result implies that the air resistence of the mirror SB through the atmosphere is roughly the same as that of an ordinary matter SB, as we already assumed earlier and have now proved. This does not mean that only the outer regions of the mirror SB will be heated by the atmosphere. The atmospheric atoms still have rapid thermal motion which will penetrate deep into the mirror SB. This is of course completely unlike a SB made of ordinary matter which remains cool inside. This ‘internal heating’ of the mirror SB should make it easier for the body to fragment and/or possibly build up enough internal pressure to explode. However, because the huge pressure from the atmosphere is dissipated over some distance within the body rather than just at its surface, the rate of ablation of a mirror SB may be significantly less than that of an ordinary SB.

Incidentally, if $`ϵ\stackrel{<}{}10^8`$ instead of the value $`10^6`$ indicated by the orthopositronium vacuum cavity experiment, a small or moderate sized SB would not lose significant energy in the atmosphere because the atmospheric atoms would pass through the body without losing much of their relative momentum. In this case the SB would release most of its energy underground in the Earth’s crust. The distance over which this would occur would simply be given roughly by Eq.(11) with the replacement $`\rho _{SB}\rho _E`$ ($`\rho _E`$ is the density of the Earth) and $`M_AM_A^{}`$, which is

$`L{\displaystyle \frac{v_i^4M_A^{}^2M_A}{16\pi \rho _EZ^2Z^2ϵ^2e^4}}\left({\displaystyle \frac{v_i}{30km/s}}\right)^4\left({\displaystyle \frac{10^9}{ϵ}}\right)^2km,`$ (12)
which was advertized earlier in Ref..

Returning to the most interesting case of large photon \- mirror photon kinetic mixing, $`ϵ10^6`$ which is indicated by the orthopositronium experiment, our earlier calculation suggests that most of the kinetic energy of a mirror matter SB is released in the atmosphere like an ordinary matter SB if it is not too big ($`\stackrel{<}{}5`$ meters) or fragments into small objects. It seems to be an interesting candidate to explain the 1908 Tunguska explosion (as well as smaller similar events as we will discuss in a moment). The Tunguska explosion toppled approximately 2,100 square kilometers of trees in a radial pattern (i.e. like spokes on a wheel) with an atmospheric release of energy estimated to be the TNT equivalent of roughly 1000 atomic bombs. There was also evidence that the inner 300 square kilometers of trees was burned from above. The broad features of the event suggest a huge explosion in the atmosphere at an altitude of between about 2.5 and 9 km which produced a downward going spherical shockwave. The spherical shockwave toppled the trees in the radial pattern and the heat from the explosion caused the flash burn of the trees. An interesting feature of this event is the lack of any extraterrestrial fragments or any (ordinary) crater(s). The estimated mass of the SB is of the order of 100 thousand tons. That is no typo. It is a remarkable result that such a large amount of extraterrestrial material apparently vanished without leaving behind significant remnants. Over the last 75 years about 35 scientific expeditions to the Tunguska site have been made with many types of search techniques, but all coming back empty handed. There have also been searches for microparticles in tree resin with some success. However their tiny abundance is hardly consistent with what might have been expected. It seems therefore to be a real possibility that the Tunguska event was due to a mirror matter SB which would not leave any ordinary fragments (the observed microparticles, if there are indeed of extraterrestrial origin, may simply be due to a small proportion of ordinary matter accreted within the mirror matter SB). Furthermore, the internal heating of the mirror SB by the interactions of the atmospheric atoms within the SB may actually cause the required atmospheric explosion.

It is also interesting to note that there is evidence that smaller ‘Tunguska-like’ events are actually quite common, occurring on a yearly basis. Such events have been catalogued by Ol’khovatov with the most recent such event occurring only a few months ago in Jordan. There are many events (see e.g. Ref.) where low altitude ‘fireballs’ are observed, yet such fireballs (if they are due to an ordinary matter SB) should originate from huge and enormously bright fireballs higher up in the atmosphere because of ablation and fragmentation. These bright parents of low altitude fireballs are inexplicably not observed. Even more remarkable is that these ‘fireballs’ have been observed in some cases to actually hit the ground (we will discuss an explicit example of this in a moment), yet no meteorite fragments were recovered. The strange properties of these events has lead to purely geophysical explanations. For example, it has been proposed that they are due to some poorly understood coupling between tectonic and atmospheric process rather than to some type of SB. Mirror matter represents an exciting and fun alternative possibility which can be tested in a number of ways as we will now briefly discuss.

First, it requires large photon - mirror photon kinetic mixing of the order given by the orthopositronium experiment for the mirror SB to release its energy in the atmosphere. Thus, we could simply repeat the orthopositronium experiment to make sure that mirror matter exists with the required kinetic mixing. More work could be done in trying to understand the detailed properties of mirror matter space bodies interacting with the Earth’s atmosphere which might allow the idea to be more rigorously compared with observations. For example, the 1997 Greenland event was observed with satellites and a ground based video camera. This event has been estimated to be due to a 36,000 Kg SB which fragmented and exploded over Greenland. No fragments or even meteoritic dust in the snow was found by search teams. The study also found that the SB had an anomalous ablation coefficient which might be something which could be used to possibly test the mirror matter hypothesis for these space bodies.

Perhaps the most spectacular way to test the idea though is to actually find it! Mirror matter could be searched for in the ground at the various impact sites. Any mirror matter fragments would have melted when they hit the ground and reformed becoming mixed with ordinary matter at some distance underground. The small effective ordinary electric charges of the mirror electrons ($`ϵe`$) which is given to them by the photon - mirror photon kinetic mixing should easily lead to enough electrostatic repulsion (which is linear in $`ϵ`$) to resist gravity which means that the mirror matter will eventually stop (if it solidifies). There may be some amount close to the surface which could potentially be extracted and purified. Importantly, many of these sites are very localized and very accessible. For example, in the recent Tunguska-like event which occurred in Jordan (about 50 kilometers from the capital Amman) only a few months ago  the fireball was observed (by a crowd of about 100 people in a funeral procession) to break up into two pieces and observed to actually hit the ground! The two sites where the ‘objects’ landed featured a half burnt tree and a half burnt rock (see Ref. for the remarkable pictures) but no ordinary crater and no ordinary matter fragments.  Potentially a mirror matter SB could leave a type of impact crater depending on the chemical composition of the SB and also on the nature of the Earth’s surface at the impact site. One could take samples of earth below the burnt tree (or the parts of the burnt tree itself) and try to extract mirror atoms. This might be possible by taking samples and putting them into a centrifuge which should allow the mirror matter to be separated from the ordinary matter (or at least greatly purified). It would be a very exciting experiment and lots of fun too!

Finally, mirror matter should have all sorts of useful industrial applications. Of course it is premature to speculate too much along these lines until it is actually discovered, but the point is that its possible existence is not merely of interest to people who want to understand the fundamental laws of nature or find out what the Universe is made of. Unlike Higgs particles or top quarks it may actually be a very useful new material with all sorts of practical applications. This provides another important motivation to search for it, either by repeating the orthopositronium experiment in vacuum or by digging it out of the ground. Of course I love Higgs particles and top quarks too but it is also important to remember that pure research in particle and astrophysics can sometimes lead to discoveries with widespread implications for society, in addition to the intrinsic merits and long term importance of such pure science itself.

Acknowledgement

The author is an Australian Research Fellow. The author would like to thank Z. Ceplecha for patiently answering some of my questions.


text:
# Antiferromagnetic integer-spin chains in a staggered magnetic field: approaching the thermodynamic limit through the infinite-size DMRG.

## I Introduction

The behavior of antiferromagnetic (AFM) Heisenberg spin chains is dominated by quantum fluctuations which suppress magnetic long-range order. However, after Haldane’s proposal , it has become clear that the value of the spin $`S`$ plays a crucial role. Integer and half-integer spin cases are indeed completely different. The spectrum of the half-integer spin chain is gapless in the thermodynamic limit, and the system is critical, in the sense that the linear response to an infinitesimal staggered magnetic field, coupled to the would-be order parameter, is divergent (at zero temperature). Conversely, the integer-spin chain has a spectrum which stays gapped in the thermodynamic limit. The first excited state has a finite distance from the groundstate (the so-called Haldane gap), and the system is not critical, i.e., the linear response to an infinitesimal staggered magnetic field is finite. The groundstate, which is characterized by a finite correlation length, is a spin liquid.

The two behaviors are only reconciled as $`S\mathrm{}`$. This limiting case, which corresponds to the suppression of quantum fluctuations, leads to the closure of the Haldane gap and to the divergence of the staggered susceptibility even in the integer-spin case, which becomes thus indistinguishable from the half-integer-spin case, as it is naively expected at large $`S`$.

All theoretical approaches are well controlled only in the limit of large $`S`$, while the increasing relevance of quantum fluctuations makes the predictions less accurate as $`S`$ is reduced. On the other hand, the introduction of the staggered magnetic field gradually freezes quantum fluctuations, and the problem arises of the description of the spin-liquid state in a staggered magnetic field and its evolution to the frozen state. A non-linear $`\sigma `$-model approach has been recently developed to investigate this evolution . This approach relies on the Haldane ansatz, and becomes exact in the limit $`S\mathrm{}`$, in the absence of staggered magnetic field. However, the validity of the ansatz must be limited by some finite field, at which the spin-liquid description breaks down. Moreover, at small $`S`$, the starting ansatz is (quantitatively) not appropriate even at zero field, and corrections must be considered, when describing the field-driven crossover from the spin liquid to the polarized AFM state.

This paper is therefore devoted to the numerical analysis of the response of the AFM integer-spin chains to a staggered magnetic field, by means of the density-matrix renormalization group (DMRG) , which gives basically exact results for one-dimensional systems, that are not biased by any a priori assumption, and represent an almost ideal numerical test for analytical results.

The case $`S=1`$ has been most extensively studied within DMRG, but the case $`S=2`$  represents quite a hard numerical task even in the absence of a magnetic field, and lacks any analysis in the presence of the staggered field, that will be the main object of our numerical investigation.

Our starting Hamiltonian reads

$$=J\underset{i=1}{\overset{L1}{}}𝐒_i𝐒_{i+1}H\underset{i=1}{\overset{L}{}}(1)^iS_i^z,$$
(1)
where $`𝐒_i`$ is the spin-$`S`$ operator on site $`i`$, $`S=1,2,\mathrm{}`$ is integer, and $`L`$ is the number of sites in the chain. $`J>0`$ is the AFM Heisenberg coupling, $`H`$ is the amplitude of the staggered magnetic field along the $`z`$ axis, which is coupled with the staggered magnetization $`M^z=_i(1)^iS_i^z`$. In view of the forthcoming DMRG study, we explicitly assumed open boundary conditions (OBC) in writing the first term in (1).

The Hamiltonian (1) is manifestly not invariant under reflection with respect to the mid point of the chain when $`L`$ is even. However, in the standard implementation of the DMRG algorithm , $`L`$ is even, and the reflection symmetry is used to reduce the numerical effort by identifying the right block with the reflection of the left one (for more details, see the discussion in Sec. II). Although there is no problem in implementing a non-symmetric algorithm , the Hamiltonian (1) can be easily made invariant under reflection by means of a simple gauge transformation (i.e. a local rotation of the reference frame), given by $`(1)^iS_i^z\stackrel{~}{S}_i^z`$, and, e.g., $`(1)^iS_i^x\stackrel{~}{S}_i^x`$, while $`\stackrel{~}{S}_i^y`$ coincides with $`S_i^y`$. Once the transformation is performed, the Hamiltonian (1) is re-cast in the form

$$\stackrel{~}{}=J\underset{i=1}{\overset{L1}{}}\left[\stackrel{~}{S}_i^z\stackrel{~}{S}_{i+1}^z+\frac{1}{2}\left(\stackrel{~}{S}_i^+\stackrel{~}{S}_{i+1}^++\stackrel{~}{S}_i^{}\stackrel{~}{S}_{i+1}^{}\right)\right]H\underset{i=1}{\overset{L}{}}\stackrel{~}{S}_i^z,$$
(2)
where $`\stackrel{~}{S}_i^\pm =\stackrel{~}{S}_i^x\pm \mathrm{i}\stackrel{~}{S}_i^y`$, and the external field is now coupled to the uniform magnetization $`\stackrel{~}{S}^z=_i\stackrel{~}{S}_i^z`$. The token we pay to recover reflection symmetry is the appearance of terms of the form $`S_i^+S_{i+1}^+`$ and $`S_i^{}S_{i+1}^{}`$, in which the spin is simultaneously raised or lowered on neighboring sites.

The original Hamiltonian (1) commutes with the $`z`$-component of the total spin $`S^z=_iS_i^z`$ and, therefore, the gauge-transformed Hamiltonian (2) commutes with the staggered magnetization $`\stackrel{~}{M}^z=_i(1)^i\stackrel{~}{S}_i^z`$. It is easily realized that the groundstate of the model lies within the $`\stackrel{~}{M}^z=0`$ subspace at large values of the field $`H`$. At zero field, the groundstate is instead degenerate (in the presence of OBC), but, since the total spin is obviously integer, a state with $`\stackrel{~}{M}^z=0`$ is certainly present in the groundstate multiplet. We explicitly checked that (at least a component of) the groundstate always lies in the $`\stackrel{~}{M}^z=0`$ subsector for every value of the amplitude of the external field $`H`$. Henceforth, for the sake of definiteness, we assume $`H>0`$, the case of negative $`H`$ being trivially recovered by changing $`\stackrel{~}{S}^z\stackrel{~}{S}^z`$ and $`\stackrel{~}{M}^z\stackrel{~}{M}^z`$ in the following discussion. In particular the groundstate energy and magnetization are even and odd functions of the external field respectively. In the rest of the paper we refer to the properties of the model described by (2), and in particular to the uniform magnetization along the external field, $`\stackrel{~}{S}^z`$. The correspondence with the original model (1) is straightforward.

The plan of the paper is the following: in Sec. II we discuss the general problem of extracting the thermodynamic properties within the DMRG approach, in connection with the interplay of finite-size effects and truncation of the Hilbert space, and we specialize the discussion to the case of the Haldane spin-liquid phase at small magnetic field; in Sec. III we present our numerical results, mostly on the $`S=2`$ case, and discuss the various regimes which characterize the model at small (Sec. III A) and large (Sec. III B) magnetic field, showing how these can be also characterized by the peculiar form of the finite-size corrections to the magnetization; concluding remarks are found in Sec. IV.

We point out that throughout the paper, whenever we discuss finite-size corrections, we specifically refer to corrections in the presence of OBC, which are those affecting our DMRG results.

## II Extracting the Thermodynamic properties

After the introduction, by S. R. White , of the DMRG, a wide majority of the DMRG studies focused on the so-called finite-size algorithm, in which the system size is first increased until a chosen size is reached, and then stopped. Subsequently, a few further “sweeps” are performed to improve the basis for that system size. On the other hand, in the infinite-size algorithm the size is always increased at each step. It is clear that, for a given lattice size, the finite-size algorithm gives the best estimate (at fixed Hilbert space). Conversely, if one is really interested in the thermodynamic limit, the information obtained with the infinite-size algorithm are of extreme value , and become crucial in the integer-spin case at issue in this paper, as we show below. We point out that extreme care has to be taken in checking the convergence of the procedure with respect to the truncation of the Hilbert space, before any sensible statement on the physical meaning of the numerical results is made.

The algorithm we used is the standard implementation of the infinite-size DMRG, which is described in  (see also ). Here we limit ourselves to a quick summary which serves as an introduction to the notations which we use in the paper. As it is customary we start with a 4-site chain which is divided into a left and a right part (blocks). The density matrix of, say, the right block in the groundstate of the entire system (superblock) is diagonalized and $`N_k`$ eigenstates, corresponding to the largest eigenvalues, are kept as the most representative of the block, while the other eigenstates are truncated away, thus reducing the size of the Hilbert space. In the presence of reflection symmetry with respect to the mid point of the chain, the left block is simply the reflection of the right block. Once the truncated left and right blocks are obtained, two more sites are added in the middle of the chain, and the new left and right blocks are defined, each as the “old truncated block + newly added site” system. The procedure is then iterated. After $`N_{RG}`$ iterations the chain contains $`L=2N_{RG}+2`$ sites. At each step, a measure of the reliability of the procedure is provided by the truncation error $`=1_{i=1}^{N_k}w_i`$, where $`w_i`$ are the eigenvalues of the density matrix, sorted in decreasing order. It follows from the definition that $`01`$, and that, in the most general case, $`=1`$ only when all the basis states are kept at each iteration, unless, under special circumstances, a simplification occurs in the spectrum of the density matrix (see below, Sec. III B). The success of DMRG when dealing with one-dimensional systems is due to their peculiar topology in the presence of short-range interactions. In most cases, indeed, the truncation error can be made small by keeping a number of states $`N_k`$ which is much smaller than the full size of the block Hilbert space \[$`(2S+1)^{L/2}`$ in the present case\], and can be considered, for all practical purposes, independent of the system size $`L`$ at fixed $``$. We point out, however, that a small value of $``$ does not by itself imply a good convergence with respect to the truncation of the Hilbert space, first of all because the convergence must be directly checked on the observed physical quantity, which may have a peculiar dependence on $``$, and also because different observables have different rates of convergence .

The observables are extracted by means of the standard DMRG procedure . In this paper we mainly discuss the groundstate energy

$$E_0(J,H;L,N_k)=\left[\stackrel{~}{}\right]^T,$$

and magnetization

$$M(J,H;L,N_k)=\left[\stackrel{~}{S}^z\right]^T,$$

where the symbol $``$ denotes the expectation value over the DMRG approximate groundstate wavefunction, for the system of size $`L`$, with $`N_k`$ states retained in the density-matrix truncation procedure, and the superscript $`T`$ means that the corresponding operator has been properly truncated onto the superblock basis.

In the following we define the groundstate energy per site in units of $`J`$, $`e_0=E_0/(JL)`$, and the magnetization per site in units of $`S`$, $`m=M/(SL)`$, such that $`0m1`$. We also introduce the notation $`hH/J`$, since $`e_0`$ and $`m`$ depend on $`J`$ and $`H`$ only through their ratio $`h`$.

A glance at the magnetization curves as a function of the system size $`L`$, for various magnetic fields, and a fixed number of states (e.g., $`N_k=50`$, see Fig. 1), shows that there are two different behaviors at small and large field respectively. The convergence of the magnetization to the thermodynamic limit is non-monotonic at small $`h`$, and $`m`$ is characterized by a well-pronounced bump at a characteristic length $`L_b`$ which increases with increasing $`S`$ and decreases with increasing $`h`$. This property strongly suggests that system sizes smaller than $`L_b`$ are definitely not representative of the thermodynamic limit, and that only for $`LL_b`$, where a monotonic convergence to the thermodynamic limit is recovered, the finite system is a reasonable representation of the infinite-length chain. For $`S=1`$, at the smallest field we investigated, $`h=0.002`$, we found $`L_b25`$ (see lowest curve in the left panel of Fig. 1). This effect is more dramatic as $`S`$ increases from 1 to 2. Fig. 1 clearly indicates that, for $`S=2`$, and for the smallest applied field $`h=0.0001`$, $`L_b70`$, and that, even for $`L200`$, the effect of the bump is not negligible, in the sense that the system appears to be far from an asymptotic convergence (see lowest curve in the right panel of Fig. 1). We shall address the physics underlying this behavior in Sec. III.

It should be clear that the presence of the bump strongly limits the possibility to draw reliable conclusions from finite-size DMRG calculations, unless the latter are pushed to exceedingly large sizes, with huge numerical effort. On the other hand, only a careful extrapolation procedure can give insights on the actual thermodynamic limit. The infinite-size algorithm is certainly suitable for this kind of study.

Within our procedure we extract the thermodynamic limit for any given number of states $`N_k`$, by fitting the large-$`L`$ behavior of the energy and magnetization as $`e_0(h;L,N_k)=e_0^{\mathrm{}}(h;N_k)+A/L`$ and $`m(h;L,N_k)=m^{\mathrm{}}(h;N_k)+B/L`$ respectively. The long-living corrections are $`O(1/L)`$ due to the OBC. We explicitly checked that the linear fit in $`1/L`$ at small field is accurate only if $`LL_b`$. For the smallest field (and $`S=2`$), we had to consider $`L1600`$ to obtain a reliable result, and the situation gets worse at larger $`S`$ .

Once $`e_0^{\mathrm{}}(h;N_k)`$ and $`m^{\mathrm{}}(h;N_k)`$ are obtained, we carefully check the convergence with respect to the Hilbert space, for both $`m`$ and $`e_0`$, by taking the extrapolation to the limit $`N_k\mathrm{}`$ in the region where the observables display a linear dependence on $`1/N_k`$ . Up to $`N_k=150`$ states were needed in the case $`S=2`$ for the smallest field $`h=0.0001`$, to obtain the linear scaling. We emphasize that this procedure explicitly deals with the convergence to the limit of infinite Hilbert space of the expectation values of the different observables, and is not equivalent to fixing a small truncation error a priori.

## III Results and discussion

In this section we mostly devote our analysis to the case $`S=2`$, for which no results in the presence of the magnetic field can be found in the literature.

The groundstate energy per site $`e_0`$ extrapolated to the infinite-size limit ($`L\mathrm{}`$) and to the full Hilbert space ($`N_k\mathrm{}`$), is a smooth monotonically decreasing function of the magnetic field, correctly interpolating between the zero-field value $`e_0(h=0)=4.76126(1)`$, and the large-field asymptotic behavior $`e_0(h1)Sh`$. The small-field region is plotted in Fig. 2, and the large-field region is shown in Fig. 3. Our result for $`e_0(h=0)`$ is obtained by an extrapolation as a function of $`h`$, and the (relatively) large error bar is due to this procedure. The value compares reasonably well with Ref.  ($`e_0=4.761244(1)`$). Our value is slightly lower, suggesting that our extrapolation procedure is slightly more efficent that the one in Ref. .

The magnetization along the magnetic field can be either computed as outlined in Sec. II, by directly evaluating the expectation value of the corresponding operator $`\stackrel{~}{S}^z`$ on the (approximate) DMRG groundstate wavefunction, or exploiting the Hellmann-Feynman theorem. The latter implies $`\stackrel{~}{S}_z=E_0/H`$ when $`E_0`$ is an exact eigenvalue. $`m`$ may therefore be estimated by a numerical derivative of the groundstate energy per site $`e_0`$ with respect to the field amplitude $`h`$. It must be pointed out that when the wavefunction is not an exact eigenstate, as in the case of DMRG, the equality does not (necessarily) hold. However, if the state obtained within DMRG is close to an exact eigenstate, the Hellmann-Feynman relation is expected to hold within good accuracy.

The magnetization values obtained following the two procedures are shown in Fig. 4, and fall on a single regular curve. This remarkable agreement represents then a strong evidence of the convergence of our method with respect to the truncation of the Hilbert space, that may be used as a further check, besides the control of the truncation error $``$, which is standard in DMRG, and the extrapolation to $`N_k\mathrm{}`$ explained in Sec. II.

The magnetization is a smooth function of $`h`$ for all the values of the magnetic field, and correctly reproduces both the zero-field limit $`m=0`$, and the $`h1`$ limit (shown in Fig. 5), where the curve saturates to $`1`$. Nonetheless, a relatively sharp crossover at $`h0.005`$ is signaled in Fig. 4 by a rapid change in the slope of the curve. The small-$`h`$ region is characterized by a spin-liquid behavior, and the crossover can be naturally associated with a change in the nature of the groundstate from the spin-liquid state to a polarized Néel state, induced by a gradual freezing of quantum fluctuations. A further support of this interpretation comes from the analysis of the curves of $`m`$ as a function of $`L`$ shown in Fig. 1. The small-$`h`$ region coincides indeed with the range in which the bump is observed in the curves. This finding sheds light on the physical origin of the bump. The bump can in fact be associated with a spin-liquid behavior, since it is observed only in this region.

If we now turn to the magnetization curves of Fig. 1, we can interpret the non-monotonicity in a simple way. When the system has a spin-liquid groundstate, it is characterized by a spin gap, and the correlation functions decay with a characteristic lengthscale $`\xi `$. If $`L\xi `$, the system is too small to display the spin-liquid features, and behaves as if $`\xi \mathrm{}`$, i.e., similarly to a half-integer spin system, which responds more strongly to a small staggered external field. As a result, this leads to an “overshooting” of the magnetization, that assumes values larger than the thermodynamic limit. Only when $`L\xi `$, the finite-system is representative of a spin liquid: the screening of the field becomes effective and the magnetization relaxes to its bulk value. This leads to an interpretation of $`L_b`$ as an intrinsic lengthscale related to some correlation length $`\xi `$. A more quantitative study of this relationship requires a model for the evolution of the Haldane spin-liquid phase in the small-field region, and is currently underway .

### A The small-field region

As mentioned above, it is well known that the integer-spin AFM Heisenberg chain has a spin-liquid groundstate, characterized by a spin gap, the so-called Haldane gap. It is quite natural that the spin-liquid behavior extends to small finite fields. As a consequence, at the smallest field $`h0.0001`$, the magnetization is linear in the field, as opposed to the half-integer spin case, which is critical (i.e. with divergent linear susceptibility) at zero field. Our (lower-bound) estimate of the linear susceptibility is $`\chi (h=0)=\frac{m}{h}|_{h=0}=1070`$, in extremely good agreement with the calculations in Ref. . In this region the non-linear $`\sigma `$-model is therefore the most correct field theory. Unfortunately, the Haldane gap decreases by increasing the spin $`S`$, so that the field-driven crossover from the spin liquid to the antiferromagnet is expected to occur at an external field which decreases with increasing $`S`$ ($`h_c0.08`$ for $`S=1`$ and $`h_c0.005`$, for $`S=2`$, as stated above). For $`h>h_c`$, the discrepancy with the non-linear $`\sigma `$-model results increases, even if a broad intermediate region exists, where the system is no longer a spin liquid, but the groundstate is not yet well described as a perturbation of the fully polarized antiferromagnet. In this intermediate region, the convergence of the magnetization to the thermodynamic limit is monotonic, and the asymptotic value is approached from below, contrary to the spin-liquid regime. Finally, a much smoother crossover takes place at larger fields, signaling the onset of a large-field “perturbative” region, where $`m`$ monotonically approaches its limiting value from above, as we discuss in Sec. III B. We point out that the field-driven freezing of the quantum fluctuations is seen in the gradual simplification of the spectrum of the density matrix, i.e., a continuous reduction of the truncation error $``$ by increasing $`h`$ at fixed $`N_k`$.

### B The large-field region

In the limit $`h\mathrm{}`$ the groundstate is trivially polarized. In this limit the spectrum of the density matrix is therefore dominated by a single state, with eigenvalue $`w_1=1`$, and the truncation error $``$ vanishes as soon as $`N_k>1`$. It is therefore reasonable to expect that at large magnetic fields, $`h1`$, the spectrum of the density matrix will be much simpler than in the spin-liquid phase, making the convergence to the full Hilbert space relatively easy. As a consequence, the DMRG data are exceedingly accurate at $`h=1`$, even by keeping $`N_k50`$ states (giving $`10^{13}`$), and at $`h=10`$, keeping $`N_k30`$ states (giving $`10^{16}`$). Moreover, in this limit we can check our numerical data by comparison with the analytical results of standard perturbation theory around the $`J=0`$ fully polarized groundstate. The groundstate wavevector for the Hamiltonian (2) at $`J=0`$ is $`|0=|S_1\mathrm{}|S_L`$, with a groundstate energy $`E_0^{(0)}=HSL`$. As the perturbation Hamiltonian is

$$\stackrel{~}{𝒱}=J\underset{i=1}{\overset{L1}{}}[\stackrel{~}{S}_i^z\stackrel{~}{S}_{i+1}^z+\frac{1}{2}(\stackrel{~}{S}_i^+\stackrel{~}{S}_{i+1}^++\stackrel{~}{S}_i^{}\stackrel{~}{S}_{i+1}^{})],$$

the groundstate energy up to second order in $`J`$ is given by

$$E_0E_0^{(0)}+\stackrel{~}{𝒱}_{00}+\underset{k>0}{}\frac{|\stackrel{~}{𝒱}_{k0}|^2}{E_0^{(0)}E_k^{(0)}}.$$

The excited states which contribute to the sum over $`k`$ are all degenerate, and correspond to wavevectors of the form $`|k=|S_1\mathrm{}|S1_k|S1_{k+1}\mathrm{}|S_L`$, $`k=1,\mathrm{},L1`$, with unperturbed eigenvalues $`E_k^{(0)}=HSL+2H`$. As a result, for $`h1`$ we find

$$e_0hS\left[1+S\frac{1}{h}+\frac{S}{2}\left(\frac{1}{h}\right)^2\right]+\frac{S^2}{L}\left(1+\frac{1}{2h}\right),$$
(3)
where the infinite-size value (first term in the r.h.s.) and the leading finite-size corrections $`O(1/L)`$ (second term in the r.h.s.), are explicitly shown. The agreement with the DMRG data for $`h>1`$ is perfect, as shown in Fig. 3.

Analogously one can obtain the magnetization along the field, by calculating the groundstate wavefunction up to second order in perturbation theory, and then evaluating the expectation value of the operator $`\stackrel{~}{S}^z`$. It is however much simpler to make use of the Hellmann-Feynman theorem once more and calculate $`\stackrel{~}{S}^z=E_0/HSLS^2(L1)/(2h^2)`$, obtaining

$$m\left(1\frac{S}{2h^2}\right)+\frac{S}{2h^2L}.$$
(4)
We have thus found that the convergence to the thermodynamic limit \[the first term in the r.h.s. of Eq. (4)\] is monotonic and the finite-size corrections \[the second term in the r.h.s. of Eq. (4)) are positive, in perfect agreement with the DMRG data (see Fig. 5\]. The deviation from the saturation value in the thermodynamic limit vanishes as $`1/h^2`$, in disagreement with the $`1/\sqrt{h}`$ behavior found within the nonlinear $`\sigma `$-model approach . Since the calculation of Ref.  relies on the Haldane ansatz, the disagreement is expected at large magnetic fields, where the groundstate is Néel-like, rather than spin-liquid-like.

## IV Conclusions

We have presented extensive DMRG calculations for integer-spin antiferromagnetic Heisenberg chains in an external staggered magnetic field. The infinite-size DMRG algorithm has proven much more suitable than the finite-size version, to extract the thermodynamic properties of these systems, due to the presence of sizeable characteristic lengthscales, which induce anomalous finite-size effects.

We have been able to determine three regions in the phase diagram of the $`S=2`$ chain, as a function of the applied staggered field. The first region, for $`0<h0.005`$ is the spin-liquid region, in which the spectrum is characterized by the Haldane gap, the spin susceptibility is linear in the field, with a coefficient $`\chi (0)=1070`$, that compares extremely well with a non-linear $`\sigma `$-model calculation . In this region the convergence of the magnetization to the thermodynamic limit is not monotonic as the system size $`L`$ is increased, due to the existence of a characteristic lengthscale $`\xi `$. Thus, quite large systems must be considered to correctly describe the spin-liquid phase. Indeed, for $`L\xi `$ the system is expected to behave as if $`\xi \mathrm{}`$, as it is the case for half-integer spin, with a strong increase of the magnetization. It is only for $`L\xi `$ that the screening of the external field in the spin-liquid phase becomes effective and the magnetization can relax to its bulk value.

The second region ($`0.005h1`$) is an intermediate crossover region, where the system is no longer a spin liquid, due to the field-driven freezing of quantum fluctuations, but the field is not large enough to give rise to a linear antiferromagnet. As a result, both the non-linear $`\sigma `$-model (valid in the Haldane phase), and perturbation theory around the linear antiferromagnet (obviously valid for $`h\mathrm{}`$), are not accurate. In this region the magnetization curve monotonically approaches the thermodynamic limit, and the limit is reached from below (contrary to the Haldane phase).

Finally, for $`h>1`$, the magnetic field is large enough to completely freeze quantum fluctuations and give rise to an almost saturated Néel antiferromagnet, that can be well described by perturbation theory in $`1/h`$. As predicted by perturbation theory, the magnetization decreases monotonically with increasing system size and approaches its thermodynamic limit from above.

A more detailed and quantitative comparison between our numerical results and the field-theoretical results of the nonlinear $`\sigma `$-model is presently in progress, aiming to clarify the identification of the characteristic length $`L_b`$ observed in the DMRG data with some correlation length, and the precise region of validity of the analytical approaches .

###### Acknowledgements.

It is a pleasure to thank F. Becca for a substantial contribution to the improvement of the DMRG code. We also thank L. Capriotti, E. Ercolessi, G. Morandi, F. Ortolani, P. Pieri, M. Roncaglia, for stimulating discussions.


text:
# The evolution of 𝛀_HI and the epoch of formation of damped Lyman-𝛼 absorbers.

## 1 Introduction

One of the fundamental phenomena still poorly understood in cosmology is the detailed process relating to the origin of structure formation after the epoch of recombination. The basic dilemma is, that while the directly observable baryonic content of galaxies at the present epoch is concentrated in stars, in the past, logically, this must have been in the form of gas. Therefore the only way to obtain a self consistent and complete picture of the galaxy formation process is to combine studies of the star light and the star formation rate with studies of the gas content of the Universe to learn about the underlying metal production and gas consumption rates.

Quasar absorbers provide a powerful observational means to study the early stages of galaxy formation independently of their intrinsic luminosity. The hydrogen absorbers are normally divided into three classes according to their column density: damped Lyman-$`\alpha `$ (DLA) systems with N(HI) $`2\times 10^{20}`$ atom cm-2; Lyman-limit Systems (LLS) with N(HI) $`1.6\times 10^{17}`$ atom cm-2; and the Lyman-$`\alpha `$ Forest with N(HI) ranging from $`10^{12}`$ to $`1.6\times 10^{17}`$ atom cm-2. The damped Lyman-$`\alpha `$ absorption lines (DLAs) are of particular importance since they contain the bulk of the neutral gas in the Universe at high redshift and are the major directly observable baryonic component at these redshifts.

For historical reasons based on the observed HI column density distribution in local galactic discs, Wolfe et al. (1986) introduced a defining limit for DLAs of N(HI) $`2\times 10^{20}`$ atom cm-2. This definition arises from the fact that 21 cm observations of local spirals show the column density drops sharply beyond this threshold (Bosma 1981). These systems exhibit strong damping wings, however, technically any absorption system with a doppler parameter $`b<`$100 km/s and N(HI) $`>10^{19}`$cm-2 will exhibit damping wings. As we will show, this ’low redshift’ DLA definition needs to be extended at high redshift to include systems down to $`10^{19}\mathrm{atom}\mathrm{cm}^2`$.

DLAs were originally thought to be the precursors of present day $`L_{}`$ disk galaxies (Wolfe et al. 1986), a scenario which is supported by some current models (Prochaska $`\&`$ Wolfe 1998). However, Haehnelt, Steinmetz & Rauch (1998) and Ledoux et al. (1998) have shown that the rotating disk interpretation for DLA systems is not unique and that the velocity structure observed by Wolfe and co-workers can also be explained by infalling sub-galactic clumps in collapsing dark matter halos with virial velocities of $``$100 $`\mathrm{km}\mathrm{sec}^1`$.

LLS are a lower column density superset of DLAs, which at $`z<1`$ are probably associated with galactic halos (Steidel, Dickinson & Persson 1994). At high redshift LLS, irrespective of their physical nature, are an important contributor to the UV opacity of the Universe since they essentially block all radiation shortward of 912 Å in the rest frame.

In this paper we present an analysis based on a new sample of high column density absorbers (DLA candidates and LLS) from a recently completed high-redshift survey (Péroux et al. 2001) combined with data from the literature (Storrie-Lombardi et al. 1996c; Storrie-Lombardi $`\&`$ Wolfe 2000 and refs therein). The combined high redshift dataset is based on observations of $``$100 quasars with z$`>`$4 and includes complete samples of 29 DLA candidates and 37 LLS with z$`>`$3.5. The layout of this paper is as follows: §2 shows how the LLS can be used to constrain the cumulative number of absorbers above log $`N(HI)=17.2`$ atom cm-2 and examines the redshift evolution of the column density distribution. The cosmological neutral gas evolution is presented in §3 and implications of our results on theories of structure formation are detailed in §4. Unless otherwise stated, this paper assumes $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_M=0.3`$ and $`H_0=65`$ km s-1 Mpc-1.

## 2 Column Density Distribution

Our new sample of high redshift Lyman-limits systems (Péroux et al. 2001) is combined with data collected from previous surveys (Storrie-Lombardi et al. 1994) in order to determine the number of LLS per unit redshift (Figure 1). We have chosen to use LLS with $`z2.4`$ in our analysis because of the very small number of systems known at low redshift ($`0.4z2.4`$), all of which are below $`z=1.1`$. Moreover, the paper is mainly concerned with the evolution of the high column density absorber population at z$`>`$2.

In the case of the colour selected quasars observed by Storrie-Lombardi et al. (1996c), Péroux et al. (2001) objects with $`z<4.2`$ are excluded to minimise colour selection bias. Only LLS with optical depth $`\tau >1`$ are used resulting in 67 systems in 124 quasars without broad absorption features, not within $`3000`$ km s-1 of the quasar emission redshift. The latter systems are excluded as they might be associated with the quasar itself (Storrie-Lombardi et al. 1994). Details on the LLS selection process and redshift estimates can be found in Péroux et al. (2001). Our new sample adds 26 LLS systems at z$`>`$3.5 to the 11 systems used by Storrie-Lombardi et al. (1994). The analysis follows that of Storrie-Lombardi et al. (1994) and the LLS number density is fitted using a power law of the form $`N(z)=N_o(1+z)^\gamma `$. The parameter values for our fit, $`N_o=0.07_{0.04}^{+0.13}`$ and $`\gamma =2.45_{0.65}^{+0.75}`$, are determined using a maximum likelihood analysis. Interestingly the number density distribution in the Lyman-$`\alpha `$ forest shows a break at $`z1.5`$ below which the distribution is flat ($`\gamma 0.2`$, Weymann et al. 1998) and above which the distribution is steep ($`\gamma =2.19\pm 0.27`$, Kim, Cristiani & D’Odorico 2001), comparable to the LLS results. Our new determination agrees very well with the previous results of Storrie-Lombardi et al. (1994) but there is a significant difference at high redshift when compared with the extrapolated results of Stengler-Larrea et al. (1995). This is not surprising since their analysis contained no data above z=3.5

We now analyse the combined samples of Lyman-limit systems and confirmed/candidate damped Lyman-$`\alpha `$ systems. The cumulative number of absorbers per unit distance interval for two redshift ranges is presented in Figure 2. The absorption distance interval, $`dX`$, (Bahcall & Peebles, 1969; Tytler, 1987) is used to correct to co-moving coordinates and thus depends on the geometry of the Universe since

$$X(z)=_0^z(1+z)^2\left[(1+z)^2(1+z\mathrm{\Omega }_M)z(2+z)\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/2}𝑑z$$
(1)
This equation differs from equation (3) in Storrie-Lombardi et al. (1996a) in that we have included the terms for a non-zero Lambda Universe. The data with N(HI) $`2\times 10^{20}`$ atom cm-2 are DLA candidates taken from our recent high-redshift survey (our observations more than double the redshift path surveyed at $`z`$$`\genfrac{}{}{0pt}{}{_>}{^{}}`$$`3.5`$ – see Péroux et al. 2001) and DLAs published by Storrie-Lombardi & Wolfe 2000. The data used for the analysis are tabulated in Appendix A & B.

The power law fit to the observed number of LLS per unit redshift is used to calculate the expected number of LLS systems:

$$LLS_{expected}=\underset{i=1}{\overset{n}{}}_{z_{min}}^{z_{max}}N_o(1+z)^\gamma 𝑑z$$
(2)
where $`z_{min}`$ and $`z_{max}`$ define the redshift path along which quasar absorbers were searched for. The LLS line profiles cannot be used to directly measure their column densities because in the range $`1.6\times 10^{17}`$ to $`2\times 10^{20}`$ atom cm-2 the curve of growth is degenerate. Nevertheless, the expected number of LLS provides a further constraint on the cumulative number of quasar absorbers and clear evidence that a simple power law is not a good representation of the observations. We thus choose to fit the data with a $`\mathrm{\Gamma }`$-distribution (cf. Schechter (1976) function used in studies of the galaxy luminosity function) as introduced by Pei $`\&`$ Fall (1995) and Storrie-Lombardi, Irwin & McMahon (1996b):

$$f(N,z)=(f_{}/N_{})(N/N_{})^\beta e^{N/N_{}}$$
(3)
where $`N`$ is the column density, $`N_{}`$ a characteristic column density and $`f_{}`$ a normalising constant. A maximum likelihood analysis is used to derive the parameters in various redshift ranges (see Table 1). The $`1`$, $`2`$ and $`3`$-$`\sigma `$ confidence contours are shown in Figure 3 for $`z<3.5`$ and $`z>3.5`$. The distributions are clearly different at the $``$3$`\sigma `$ level and indicate that there are less high column density systems (N(HI) $`>10^{21}`$ atom cm-2) at high-redshift, z$`>`$3.5, compared with 2$`<`$z$`<`$3.5 confirming the results from Storrie-Lombardi, McMahon & Irwin (1996a) and Storrie-Lombardi & Wolfe (2000). For comparison with previous works, Figure 4 shows the column density distribution, $`f(N,z)`$, at $`z_{\mathrm{abs}}>3.5`$ together with the double power law of Petitjean et al. (1993).

## 3 Cosmological Mass Density Evolution

The mass density of absorbers can be expressed in units of the current critical mass density, $`\rho _{crit}`$, as:

$$\mathrm{\Omega }_{DLA}(z)=\frac{H_o\mu m_H}{c\rho _{crit}}_{N_{min}}^{\mathrm{}}Nf(N,z)𝑑N$$
(4)
where $`\mu `$ is the mean molecular weight and $`m_H`$ is the hydrogen mass. The total HI may be estimated as:

$$_{N_{min}}^{\mathrm{}}Nf(N,z)𝑑N=\frac{N_i(HI)}{\mathrm{\Delta }X}$$
(5)
If a power law is used to fit $`f(N,z)`$, up to $`90\%`$ of the neutral gas is in DLAs (Lanzetta, Wolfe & Turnshek 1995), although an artificial cut-off needs to be introduced at the high column density end because of the divergence of the integral. If instead a $`\mathrm{\Gamma }`$-distribution is fitted to $`f(N,z)`$ this removes the need to artificially truncate the high end column distribution and can be used to probe in more detail the neutral gas fraction as a function of column density and how this changes with redshift. Figure 5 shows the cumulative mass fraction as a function of column density for $`z>`$ and $`<3.5`$. At $`z`$ $`\genfrac{}{}{0pt}{}{_>}{^{}}`$$`3.5`$, up to $`45\%`$ of the neutral gas is in systems with $`10^{19}<N(HI)<10^{20.3}`$. We refer to these systems as sub-DLAs. These observations supports a picture where at z$`>`$3.5, we are directly observing the formation of high column density neutral hydrogen DLA systems from lower column density units. The value of $`N_{}`$ decreases by a factor of $``$5 from $`1.6\times 10^{21}\mathrm{atom}\mathrm{cm}^2`$ at z$`<`$3.5 to $`2.9\times 10^{20}\mathrm{atom}\mathrm{cm}^2`$ at z$`>`$3.5.

Figure 6 displays $`\mathrm{\Omega }_{DLA}`$ contained in DLAs (filled circles) and the total amount of neutral gas (DLAs plus sub-DLAs) for a non-zero $`\mathrm{\Lambda }`$ Universe (grey stars). Our observations are consistent with no evolution in the redshift range z= 2 to 5. Empty Universe models produce similar results.

## 4 Discussion and conclusions

Our study is mainly based on optically selected quasars, so it is obvious that quasars that lie behind dusty DLA systems will be underrepresented if they exist. Pei & Fall (1995) have used self-consistent closed-box/inflow-outflow galactic chemical evolution models to show that the fraction of missing DLAs at $`z=3`$ ranges from 23% to 38%. However, dust is unlikely to be important at high-redshift because of the short time available for its production. Dust will have a larger effect at z=2 compared with z$`>`$3 and thus it could not cause the form of evolution in f(N,z) that we observe. But a recent work by Ellison et al. (2001) has shown that measurements of $`\mathrm{\Omega }_{DLA}`$ based on optically selected quasars survey are not significantly affected by dust at z$``$2.

It is possible that the overall cosmological evolution of $`\mathrm{\Omega }_{DLA}`$ is dominated by feedback processes rather than by gas consumption due to star formation. The observed mass of HI at any redshift may relate more to the recombination timescale for the ionized HII and the timescale for cooling and collapse into molecular hydrogen ($`\mathrm{H}_2`$) and thereafter into stars. Models from Efstathiou (2000) indicate that DLA systems might predominantly be due to the outer parts of galaxies which do not even participate in star formation.

It can be seen in Figure 6 that $`\mathrm{\Omega }_{DLA}`$ is significantly below the current estimates of $`\mathrm{\Omega }_M`$ in stars in the nearby Universe. This is a significant change in the situation compared with Storrie-Lombardi et al. 1996c (Figure 1b). The main reason for this change is that for the currently favoured Lambda-dominated cosmology, the mass in HI at high redshift drops by a factor of $``$50% compared with an $`\mathrm{\Omega }_M`$=1 Universe. This is purely a geometric effect and mainly affects the relative normalisation between z=2 and z=0.

In any case, the observations in the redshift range 2 to 5 show no evolution in the total amount of neutral gas in contrast to the earlier results of Lanzetta, Wolfe and Turnshek (1995), who found that $`\mathrm{\Omega }_{DLA}`$(z$``$3) was twice $`\mathrm{\Omega }_{DLA}`$(z$``$2). Under simple assumptions of closed box evolution, this could be interpreted as indicating there is little gas consumption due to star formation in DLA systems in this redshift range. Similarly, at z $`>`$ 2, Prochaska, Gawiser & Wolfe (2001) conclude that there is no evolution in the metallicity of DLA systems from column density-weighted Fe abundance measurements in DLAs. Since metallicity studies focus on the higher column density systems they may be giving a biased or incomplete view of global galactic chemical evolution at z$`>`$3. These metallicity evolution results are still very much open to debate as another study by Savaglio (2000) shows that the metallicity content of DLAs and sub-DLAs does decrease with redshift when one excludes the highest column density systems (N(HI)$`>6\times 10^{20}`$) from the analysis. It is important to note that at z$`>`$3.5, 90% of the HI lies below this limit. Moreover, the current practice of using column density weighted metallicities neglects the fact that the metallicity observations are biased towards high HI column density systems and hence do not necessarily trace the global metallicity evolution.

The detailed comparison of high redshift (z$`>`$2) observations with lower redshifts is currently quite difficult. In the local Universe 21-cm emission measurements are used to determine the neutral gas mass density. However, direct detections of DLA systems are paradoxically more difficult for several reasons: the observed DLA wavelengths are shifted to the UV requiring HST observations, the redshift path available per quasar is lower and the number of DLA systems per unit redshift is lower. In fact the most recent determination of $`\mathrm{\Omega }_{DLA}`$ at z$``$1 by Rao & Turnshek (2000) is higher than the previous estimate by Lanzetta, Wolfe and Turnshek (1991).

Nevertheless, since DLAs are the most directly observable baryonic mass systems at high redshift, their properties are an important constraint for theories of galaxy formation (Kauffmann $`\&`$ Charlot 1994, Kauffmann $`\&`$ Haehnelt 2000). For example, Somerville, Primack $`\&`$ Faber (2000) have used the integrated total neutral gas (HI+HeII) to constrain models (Figure 6) which vary in their recipes for star formation (due to collisional starburst, constant efficiency quiescent star formation or accelerated quiescent star formation). The models take into account the “cold gas” which includes neutral as well as molecular gas.

To summarise, we find that at z$`>`$3.5 the fraction of neutral gas mass in DLAs above the commonly used Wolfe et al. (1996) DLA definition of N(HI)$`>2\times 10^{20}`$ atom cm-2 is only 55$`\%`$ compared with $``$90% at z$`<`$3.5 and that the remaining fraction of the neutral gas mass lies in systems below this limit, in “sub-DLAs” with column density $`10^{19}<`$ N(HI) $`<2\times 10^{20}`$ atom cm-2. We also find that the total mass is conserved over the redshift range 2 to 5. This is strongly indicative that we are observing the assembly of high column density systems from lower column density units and independently of the precise physical nature of damped Lyman-$`\alpha `$ shows that we are observing the epoch of their formation or initial collapse. A systematic study of the kinematic and metallic properties of DLA systems with z$`>`$3.5 and N(HI) above $`10^{19}`$ atom cm-2 is an obvious program for the new generation of echelle spectographs on 8m class telescopes. It will also be important to directly establish the N(HI) column density distribution function below 2 $`\times 10^{20}`$ atom cm-2. We are currently undertaking such a program using VLT UVES archival echelle data of high redshift quasars (Dessauges-Zavadsky et al. 2003, Péroux et al. 2003). The data used for the analysis in this paper are available from http://www.ast.cam.ac.uk/~quasars.

## 5 Acknowledgments

We thank Bob Carswell, Mike Fall, Sandro D’Odorico, Patrick Petitjean, Max Pettini, Martin Rees, and Joop Schaye for helpful discussions. CP is grateful to PPARC and the Isaac Newton Trust for support and RGM thanks the Royal Society and the Raymond and Beverly Sackler Foundation for support.

## Appendix A The sample of quasars with confirmed or candidate DLA absorbers

Table 2 lists the 114 confirmed or candidate DLA absorbers used in our study, their redshift and HI column densities together with the redshift path surveyed for each quasars. The following corrections to the low redshift sample derived by Lanzetta, Turnshek and Wolfe (1995) have been made: the Q1329+4117 DLA at $`z_{abs}`$=0.5193, and the Q2112+059 DLA at $`z_{abs}`$=0.2039 have been ruled out (Jannuzi et al. 1998, Fynbo, Moller & Thomsen, 2001), while Q0302$``$223 has a DLA at $`z_{abs}`$=0.1014 (Jannuzi et al. 1998).

## Appendix B The sample of quasars without DLA absorbers

Table 2 lists the quasars from our study where no DLAs where discovered. The table also provides the redshift path surveyed for each quasars.


text:
# Right orderable residually finite 𝑝-groups and a Kourovka notebook problem

## Abstract.

A. H. Rhemtulla proved that if a group is a residually finite $`p`$-group for infinitely many primes $`p`$, then it is two-sided orderable. In problem 10.30 of the Kourovka notebook 14 th. edition, N. Ya. Medvedev asked if there is a non-right-orderable group which is a residually finite $`p`$-group for at least two different primes $`p`$. Using a result of Dave Witte, we will show that many subgroups of finite index in $`\mathrm{GL}_3()`$ give examples of such groups. On the other hand we will show that no such example can exist among solvable by finite groups.

###### Key words and phrases:

right orderable group, residually finite $`p`$-group

A group $`G`$ is right orderable if it has a total ordering $``$ such that $`xyxgyg`$ whenever $`g,x,yG`$. It is *two-sided orderable* if in addition $`xygxgy`$ whenever $`g,x,yG`$. For much information on right ordered groups, see the books . It was proved in  that a group which is a residually finite $`p`$-group for infinitely many primes $`p`$ is two-sided orderable. In problem 10.30 of the Kourovka notebook 14 th. edition , N. Ya. Medvedev asks if there is a non-right-orderable group which is a residually finite $`p`$-group for at least two different primes $`p`$. We shall prove

###### Theorem 1.

Let $`𝒫`$ be a finite set of primes. Then there exists a non-right-orderable group which is a residually finite $`p`$-group for all $`p𝒫`$.

The proof depends on a theorem of Witte , and we will see that the groups in Theorem 1 can be taken to be subgroups of finite index in $`\mathrm{GL}_3()`$, the group of 3 by 3 invertible matrices with integer entries.

Going in the opposite direction, we shall prove the following result.

###### Theorem 2.

Let $`p,q`$ be distinct primes and let $`G`$ be a solvable by finite group. If $`G`$ is a residually finite $`p`$-group and a residually finite $`q`$-group, then $`G`$ has a series

$$1=G_0\mathrm{}G_n=G$$

with $`G_iG`$ and $`G_{i+1}/G_i`$ torsion free abelian for all $`i`$. In particular $`G`$ is right orderable.

###### Acknowledgement.

This work was carried out while I was at the Sonderforschungsbereich in Münster. I would like to thank Wolfgang Lück for organizing my visit to Münster, and the Sonderforschungsbereich for financial support.

###### Proof of Theorem 1.

For each prime $`p`$, let $`G_p=\{A\mathrm{GL}_3()AImodp\}`$, the congruence subgroup of level $`p`$ (here $`I`$ denotes the identity matrix of $`\mathrm{GL}_3()`$). Then $`G_p`$ is a residually finite $`p`$-group. Set $`G=_{p𝒫}G_p`$. Then $`G`$ is a residually finite $`p`$-group for all $`p𝒫`$. Furthermore $`G`$ has finite index in $`\mathrm{GL}_3()`$, so by \[5, Proposition 3.3\] we see that $`G`$ is not right orderable. This completes the proof of Theorem 1. ∎

###### Proof of Theorem 2.

Observe that $`G`$ is torsion free. Let $`H`$ be a normal solvable subgroup of finite index in $`G`$ with minimal derived length. We shall prove the result by induction on the derived length of $`H`$, the result being obvious when this is zero because then $`H=1`$ and thus $`G`$ will be finite. We now assume that the derived length of $`H`$ is at least one.

Let $`A_0`$ denote the penultimate term of the derived series of $`H`$. Then $`A_0`$ is a normal abelian subgroup of $`G`$ and $`H/A_0`$ has strictly smaller derived length than that of $`H`$. Let $`A`$ be a maximal normal abelian subgroup of $`G`$ containing $`A_0`$.

We shall let $`\widehat{K}`$ denote the pro-$`p`$ completion of a group $`K`$. Then the exact sequence $`1AGG/A1`$ yields an exact sequence

$$\widehat{A}\widehat{G}\widehat{G/A}1.$$

Let $`B`$ denote the image of $`\widehat{A}`$ in $`\widehat{G}`$. Since $`A`$ is abelian, we see that $`\widehat{A}`$ is abelian and we deduce that $`B`$ is abelian. Also we may view $`G`$ as a subgroup of $`\widehat{G}`$ because $`G`$ is a residually finite $`p`$-group. Therefore $`BG`$ is an abelian normal subgroup of $`G`$ containing $`A`$ and we conclude that $`BG=A`$. But $`BG`$ can also be described as the kernel of the natural map $`G\widehat{G/A}`$. Therefore $`G/A`$ is isomorphic to a subgroup of $`\widehat{G/A}`$ and we deduce that $`G/A`$ is a residually finite $`p`$-group. Similarly $`G/A`$ is a residually finite $`q`$-group. Since the derived length of $`HA/A`$ is strictly less than the derived length of $`H`$, induction shows that $`G`$ has a series

$$A=G_1\mathrm{}G_n=G$$

with $`G_iG`$ and $`G_{i+1}/G_i`$ torsion free abelian for all $`i`$. By setting $`G_0=1`$, we obtain the required series for the first part of Theorem 2.

The assertion that $`G`$ is right orderable now follows from \[3, Theorem 7.3.2\] and the fact that torsion free abelian groups are right orderable. ∎