text:
# VLA 8.4-GHz monitoring observations of the CLASS gravitational lens B1933+503

## 1 Introduction

B1933+503 is the most complex arcsec-scale gravitational lens system known and the most noteworthy system discovered in the Cosmic Lens All-Sky Survey (CLASS; Browne et al. 1998), a survey of approximately 15,000 flat-spectrum radio sources observed with the VLA at 8.4 GHz with a resolution of 200 milliarcsec (mas). Sykes et al. (1998) report the discovery of this lens and present VLA and MERLIN maps that reveal up to ten components, four of which are compact and have flat spectra while the rest are more extended and have steep spectra. The ten components are believed to be the multiple images of a background source that consists of a flat-spectrum core (quadruply imaged) and two compact lobes symmetrically situated on opposite sides of the core (one quadruply imaged and the other doubly imaged). A Hubble Space Telescope (HST)/WFPC2 image ($`I`$-band) of B1933+503 (Sykes et al. 1998) shows a faint galaxy with a compact core (the lensing galaxy) but none of the images of the background object are detected. Spectra taken with the Keck telescope have given a redshift for the lensing galaxy of 0.755 (Sykes et al. 1998). Sub-mm observations of B1933+503 (Chapman et al. 1998) pointed to the fact that the source of the lens system should be lying above a redshift of 2 and that has been confirmed recently by observations made with the United Kingdom Infra-red Telescope (UKIRT) that have yielded a redshift of 2.62 for the background source (Norbury et al., in preparation).

HST NICMOS F160W-band observations (Marlow et al. 1998) of B1933+503 have uncovered the infra-red counterparts to two of the four compact components. These same two components were also detected in a VLBA 1.7-GHz map with a resolution of 6 mas. Marlow et al. suggest that the failure to detect the other two flat-spectrum components in the infra-red can be due to either rapid variability on a time-scale less than the time delay or due to extinction in the interstellar medium (ISM) of the lensing galaxy. The absence of these same two components in the radio could arise due to scatter-broadening of the images when propagating through the ISM of the lens galaxy as well as from time variability.

B1933+503 offers a unique opportunity to construct an accurate mass model due to the wealth of observational constraints that are available. As well as the large number of images, yielding both positional and flux density ratio constraints, some of these are extended. Perhaps more crucial is that three epochs of VLA 15-GHz data (Sykes et al. 1998) suggest that the flat-spectrum components are variable, by as much as 33 per cent over the period July to September 1995. This means that it should be possible to determine the relative time delays from flux density monitoring. Taken together this makes B1933+503 a promising system for a determination of the Hubble constant, $`H_0`$ (Refsdal 1964). The observational constraints given by Sykes et al. (1998) have been used as a starting point for the model presented by Nair (1998). The time delays of the flat-spectrum components with respect to component 1 derived from this model are approximately 8 d for component 3, 7 d for component 4 and 9 d for component 6.

In this paper we present the results of a four-month monitoring campaign using the VLA in its ‘A’ configuration at a frequency of 8.4 GHz. This frequency was chosen over 15 GHz (where the magnitude of flux density variations are likely to be greater) for a number of reasons. As well as the sensitivity of the VLA being greater at 8.4 GHz, the phase stability of the array is much poorer at 15 GHz. In addition, VLA 15-GHz monitoring of the lens system B0218+357 (Biggs et al. 1999) showed that at this frequency the data are subject to calibration errors of a few per cent due to the gains of the antennas varying with elevation. In Section 2 we describe the observations and observing strategy as well as the reduction methods that we used on the data. In Section 3 we show the light curves that were obtained from the data and we follow this with a discussion of the results.

## 2 Observations and Data Reduction

B1933+503 was observed between February 16 1998 and June 1 1998, during which time the VLA was in its ‘A’ configuration. The 37 epochs were separated, on average, by 2.8 d. The observations were made simultaneously in two frequency bands (IFs) each with a bandwidth of 50 MHz and centre frequencies of 8.4351 and 8.4851 GHz respectively. This gives an angular resolution of approximately 200 mas.

The observing strategy was kept as consistent as possible over the period of the monitoring. At each epoch the B1933+503 observations were split into two or three scans, each of approximately 6 min duration. Before and after each of the B1933+503 scans were separate ($``$one min duration) scans of the steep-spectrum source B1943+546 (Patnaik et al. 1992) which was used as the phase and amplitude calibrator. B1943+546 has all the requirements of a good primary calibrator being close to the target, having a strong correlated flux density and being relatively unresolved with the VLA at this frequency. The duration of the combined observations of both sources at a typical epoch was 30 min, but several epochs consisted of 20 min observing time.

Calibration was performed using the NRAO Astronomical Image Processing Software package aips and each IF was calibrated separately. The first step was to make the best possible map of 1943+546 as this calibrator is slightly resolved at the observing frequency. This was done by performing a first order calibration of several best epochs of 1943+546 data assuming the flux density of 610 mJy measured by Patnaik et al. (1992). This source has a steep spectrum and is found from VLBI maps (Polatidis et al. 1995; Xu et al. 1995) to be dominated by emission from two extended lobes on scales of $``$50 pc. Therefore it is unlikely that the source’s flux density will have changed significantly since this measurement or that it will have changed over the period of the monitoring. These epochs were then combined and iteratively mapped and self-calibrated. The structure found from this map agrees well with that seen in the above mentioned VLBI maps. The CLEAN components from the map were then used as a model input to the aips task calib to derive telescope gain solutions at each epoch of observation. These gain solutions were then applied to the target B1933+503.

The positions and approximate flux densities of each component of B1933+503 were first established by making a map of the source. In order to make a high dynamic-range map of the source we combined ten epochs of observations with the dbcon task. As the images of the radio core are seen in higher resolution MERLIN maps to be compact, a traditional CLEAN was not used to make the map. Instead, a single delta component (point source) was placed at the approximate position of each of the compact components in the dirty map and the positions and flux densities of these were then optimised by model-fitting. Once the best fit had been found, the remaining extended flux was then CLEANed in the usual way. The resulting map is shown in Fig. 1 with eight of the ten components labelled. The remaining two components cannot be easily distinguished due to the low resolution of the image. Higher resolution observations and lens modelling have shown that component 2 actually consists of two merging images, denoted 2a and 2b by Nair (1998). The remaining component (1a) can be barely discerned in Fig. 1 as a slight extension to the north-east of component 1. This extension is much more obvious in a MERLIN 1.7-GHz map (Sykes et al. 1998).

The flux densities of the four compact components of B1933+503 were determined by model-fitting directly to the $`(u,v)`$ data, and to each IF individually, using the difmap (Pearson et al. 1994; Shepherd 1997) model fitter. Each epoch of monitoring data was fit using the model derived from the map with all flux densities and positions held constant, except for the flux densities of the four compact components which were allowed to float. The model fit was interspersed with two iterations of phase self-calibration followed by a final overall amplitude correction i.e. one correction for each antenna for the entire observing period, normalising against the extended structure whose brightness should not vary over the monitoring observations. These amplitude corrections are never more than a few per cent per telescope. The reduced $`\chi ^2`$ indicated a reasonably good fit to each epoch, with a mean value of 1.35 and an rms scatter of 0.04 for each IF’s data. At this point the two flux density measurements at each epoch from each IF are averaged together.

## 3 Results

The averaged total flux density light curves for components 1, 4, 3 and 6 are shown in Fig. 2 and are plotted from the top of the figure downwards in the order in which they are expected to vary. The range encompassed by the y-axis is, for each component, equal to the same fraction ($`20\%`$) of the average flux density calculated over the period of the monitoring. The error bars attached to each data point have been estimated from a consideration of the difference between the flux densities derived from each IF at each epoch. For each component, a frequency histogram of the flux density difference at each epoch is approximately Gaussian with a mean that does not differ from zero by more than 0.02 mJy. The 1$`\sigma `$ error derived using this method for each component is very similar and has been fixed at 0.07 mJy in Figure 2. This is in fact only a lower limit to the uncertainty in the flux densities as it does not take into account errors that are correlated in each IF. These include errors arising from overall amplitude miscalibration and those that are a function of the $`(u,v)`$ coverage of the data which is essentially the same for each IF of a given epoch.

From the variations seen in Figure 2 it does not appear that the source has varied significantly during the course of our monitoring campaign. This conclusion is supported in a number of ways. First, the mean flux density and rms percentage scatter about this mean for components 1, 4, 3 and 6 are: 4.35 mJy, 15.73 mJy, 3.25 mJy and 3.43 mJy, and 1.6, 0.6, 2.6, and 2.2 per cent respectively. What is apparent is that the rms scatter in the light curves is inversely proportional to the mean flux density of the component. If there were actually a dominant signal of intrinsic variability then it would be expected that the fractional scatter in the data would be equal for each component. This is not the case. This can be seen most clearly by comparing the ‘variability’ seen in the brightest component (4) which is much less than that seen in the other three fainter components. Also, a preliminary time delay analysis using both a chi-squared minimisation and the Dispersion method of Pelt et al. (1996) result in best-fit time delays that are totally at odds with those predicted by lens models of this system. The rms scatters quoted above are in general larger than expected from thermal noise considerations alone (by factors of $``$2–3) and indicate that the data are corrupted by errors stemming from deconvolution and model fitting. Errors introduced by incorrect amplitude calibration from epoch to epoch must be present at a low level ($`0.6\%`$) as these would produce the same fractional scatter in each component.

Table 1 lists the flux densities of the four flat-spectrum components measured from the map in Fig. 1 together with previous values obtained at 8.4 GHz with the VLA in ‘A’ configuration on July 6, 1995 (Sykes et al. 1998) as well as those determined from the original VLA snapshot observation taken in 1994. We see that although there were no changes between the two most recent epochs there are signs that the overall flux density of these compact components in 1994 was higher than that in the last two epochs. As flux density variability has been reported at 15 GHz by Sykes et al. (1998) and because a comparison of the MERLIN and VLA maps by the same authors support variations in the relative flux densities of the same components we think that this source has been, and will probably again be, variable.

## 4 Concluding remarks

We have presented the results of a four month program of monitoring the ten-component lens system B1933+503 with the VLA which have not allowed us to measure the independent time delays in this otherwise promising system for determining $`H_0`$. The component light curves show fluctuations of 0.7–2.5 per cent which we do not believe represent source variability. Instead they are more likely to result from a combination of thermal noise, deconvolution and model fitting errors. Since the system has already shown indications of variability (Sykes et al. 1998), we think that we have caught this system in a quiescent phase. For this reason, we have embarked on a second monitoring campaign that commenced in June 1999, also with the VLA at 8.4 GHz. An advantage that we have now over the previous attempt is that the flux density ratios between the components are known to high accuracy. Since we did not detect any variability during the 4 month monitoring campaign, which is a period significantly longer than the predicted time delay of 7–9 d, we have determined the intrinsic flux density ratios (Table 1). Therefore, should this second set of observations again not show reliable features that can be identified in each light curve, but do show a monotonic increase or decrease in flux density, a time delay may still be determined. This has previously been demonstrated by Koopmans et al. (1999) for the CLASS lens B1600+434.

## Acknowledgments

This research was supported by European Commission, TMR Programme, Research Network Contract ERBFMRXCT96-0034 “CERES”. The VLA is operated by the National Radio Astronomy Observatory which is supported by the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. ADB acknowledges the receipt of a PPARC studentship.


format: text

text:
# Magneto-optical control of bright atomic solitons

## I Introduction

The experimental demonstration of Bose-Einstein condensation in atomic vapors has rapidly lead to spectacular new advances in atom optics. In particular, it has enabled its extension from the linear to the nonlinear regime , very much like the laser led to the development on nonlinear optics in the 1960s. It is now well established that two-body collisions play for matter waves a role analogous to that of a Kerr nonlinear crystal in optics. The first experimental verification of this analogy was four-wave mixing . In addition, it is also possible to nonlinearly mix optical and matter-waves, as demonstrated recently in matter-wave superradiance  and in the first realization of a phase-coherent matter-wave amplifier . In addition, it is known that the nonlinear Schrödinger equation which describes the condensate in the Hartree approximation support soliton solutions . For the case of repulsive interactions normally encountered in BEC experiments, the simplest solutions are dark solitons, that is, “dips” in the density profile of the condensate . These dark solitons have been recently demonstrated in two experiments  which appear to be in good agreement with the predictions of the Gross-Pitaevskii equation.

While very interesting from a fundamental physics point-of-view, dark solitons would appear to be of limited interest for applications such as atom interferometry, since what would be desirable is to achieve the dispersionless transport of a spatially localized ensemble of atoms, rather than a ”hole.” In that case, bright solitons are much more interesting. However, the problem is that large condensates are necessarily associated with repulsive interactions, for which bright solitons might appear to be impossible since the nonlinearity cannot compensate for the kinetic energy part (diffraction) in the atomic dynamics. While this is true for atoms in free-space, this is not the case for atoms in suitable potentials, eg. optical lattices. This is because in that case, it is possible to tailor the dispersion relation of the atoms in such a way that their effective mass becomes negative. For such negative masses, a repulsive interaction is precisely what is required to achieve soliton solutions. This result is known from nonlinear optics, where such soliton solutions, called gap solitons, have been predicted and demonstrated .

We showed in previous work  that such bright solitons are also possible in matter-wave optics using spinor condensates, but did not discuss explicitly how to excite and control them. The present paper addresses these questions, and shows that a combination of optical and magnetic fields can be used to generate solitons of various velocities, and subsequently, to control them, split and recombine them, etc. This magneto-optical control results from the use of spinor condensates, and can be achieved at minimal cost in terms of atomic loss. In particular, we illustrate how to realize a Mach-Zehnder interferometer for bright atomic solitons. This opens up the way to intriguing new ways to manipulate and transport coherent matter in ways which complement those offered by standard optical tweezers or by atomic wave guides .

Amongst Dan Walls many contributions to nonlinear and atom optics he was one of the first to recognize the emerging area of nonlinear atom optics, and he introduced the idea of atomic solitons traveling in laser beams already in $`1994`$ along with Zhang and Sanders . The current wave of interest in atomic solitons therefore finds its origin in Dan’s seminal work, and we all greatly miss the originality, insight, and energy he brought to the field.

The paper is organized as follows: Section II briefly reviews the analysis leading to the predictions of gap solitons , while section III discusses some of their most important characteristics. Section IV exploits these properties to develop tools to excite gap solitons and control their dynamical behavior, leading to the demonstration of an atomic Mach–Zehnder interferometer in Section V. Finally, section VI is a summary and outlook. While the main text uses a simple model of light-matter interaction for notational clarity, a more realistic coupling scheme, using the full hyperfine structure of the Sodium $`3S_{1/2}`$$`3P_{3/2}`$ transition, is presented in Appendix A.

## II Physical Model

To set the stage for our analysis, we first briefly review the main ingredients of the theory of gap solitons of Ref. . The system we consider consists of a Bose-Einstein condensate interacting with two counterpropagating, focused Gaussian laser beams of equal frequency $`\omega _l`$ but opposite circular polarizations, see Fig. 1. The optical dipole potential associated with the applied laser beams is assumed to provide tight transverse confinement for the BEC in the $`(X,Y)`$ plane, thereby forming a cigar shaped condensate of transverse cross-sectional area $`A_T`$. In the following, we confine our discussion to the one-dimensional dynamics of the BEC along the $`Z`$-axis for simplicity.

In addition to supplying a transverse optical potential, the laser beams can drive two-photon transitions between different Zeeman sublevels of the atomic ground state. For illustrative purposes we consider the case of Sodium and the two-photon coupling of the Zeeman sublevels $`|1=|F_g=1,M_g=1`$ and $`|1=|F_g=1,M_g=1`$. For example, starting in the $`|1`$ state this process involves the absorption of a $`\sigma _+`$ photon from the right propagating laser beam followed by emission of a $`\sigma _{}`$ photon into the left propagating laser beam.

We must of course assume that the excited states involved in the atom-field interaction are far-detuned from the applied laser frequency, a necessary requirement to avoid the detrimental effects of spontaneous emission. This however raises a serious issue, since for the alkali atoms that we have in mind the two-photon coupling strength vanishes in the limit when the detuning is large compared to the excited–state hyperfine splitting . But this difficulty can be circumvented: We show in the Appendix that by using a four-photon scheme involving an additional $`\pi `$-polarized field incident perpendicular to the $`Z`$-axis, one can achieve an effective two-photon coupling between the states $`|\pm 1`$ which survives in the limit of large detunings. We therefore proceed with our two-photon coupling model and refer the interested reader to the Appendix for details.

By restricting our attention to the coupled states $`|\pm 1`$ the effective single–particle Hamiltonian for our model system can be written as

$$H_{\mathrm{𝑒𝑓𝑓}}=\frac{P_Z^2}{2m}+g\mathrm{}\delta ^{}\left[|11\left|e^{2iK_lZ}+\right|11|e^{2iK_lZ}\right],$$
(1)
where we have omitted constant light-shift terms. Here $`P_Z`$ is the atomic center–of–mass momentum operator along $`Z`$, $`m`$ the mass of the atom, $`K_l=\omega _l/c`$ is the magnitude of the field wave vector along $`Z`$, $`g`$ is a coupling constant between the ground and excited states characteristic of the atom and transition involved, $`\delta ^{}=𝒟^2^2/\mathrm{}^2\delta `$, with the detuning $`\delta =\omega _l\omega _a`$, $``$ is the laser field amplitude, and $`𝒟`$ is the reduced dipole moment for the $`3S_{1/2}`$$`3P_{3/2}`$ transition.

The first term of the effective Hamiltonian (1) describes the quantized atomic center-of-mass motion, and the remaining terms give the effective coupling of the two Zeeman sublevels via the applied laser fields. The exponential terms $`\mathrm{exp}(\pm 2iK_lZ)`$ arise from the fact that the two-photon transitions involve the absorption of a photon from one light field and reemission into the other. Finally, introducing a spinor macroscopic condensate wave function $`\mathrm{\Psi }(Z,t)=[\mathrm{\Psi }_1(Z,t),\mathrm{\Psi }_1(Z,t)]^T`$ normalized to the number of atoms $`N`$, and including the many–body effects via a mean-field nonlinearity, we obtain the coupled Gross–Pitaevskii equations

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=H_{\mathrm{𝑒𝑓𝑓}}\mathrm{\Psi }+U|\mathrm{\Psi }|^2\mathrm{\Psi },$$
(2)
where $`U=4\pi \mathrm{}^2a_{sc}/m`$, $`a_{sc}`$ is the s-wave scattering length, $`|\mathrm{\Psi }|^2=|\mathrm{\Psi }_1|^2+|\mathrm{\Psi }_1|^2`$, and we have assumed that the magnitude of the self- and cross-nonlinearities are equal for simplicity.

It is convenient to re-express Eq. (2) in dimensionless form by introducing the scaled variables $`\tau =t/t_c`$, $`z=Z/l_c`$ and $`\psi _j=\mathrm{\Psi }_j/\sqrt{\rho _c}`$ where

$$t_c=\frac{1}{g\delta ^{}},l_c=\frac{t_c\mathrm{}K_l}{m},\rho _c=\left|\frac{g\mathrm{}\delta ^{}}{U}\right|.$$
(3)
Equations (2) then become

$`i{\displaystyle \frac{}{\tau }}\left(\begin{array}{c}\psi _1\\ \psi _1\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}M^2& e^{2ik_lz}\\ e^{2ik_lz}& M^2\end{array}\right)\left(\begin{array}{c}\psi _1\\ \psi _1\end{array}\right)`$ (10)
$`+`$ $`\text{sgn}\left(g\delta ^{}/U\right)\left|\psi \right|^2\left(\begin{array}{c}\psi _1\\ \psi _1\end{array}\right),`$ (13)
where $`M=g\delta ^{}m/2\mathrm{}K_l^2`$ is a mass-related parameter such that $`k_l=K_ll_c=1/2M`$. Throughout this paper we use a characteristic laser intensity of $`I=50\text{ W/cm}^2`$ for each of the two counter-propagating laser beams and a wavelength of $`\lambda _l=985\text{ nm}`$. In addition we use $`g=1/4`$, which is simply the Clebsch-Gordan coefficient for the $`F_g=1F_e=1`$ transition. This results in the characteristic scale values $`t_c=685\text{ }\mu \text{s},l_c=12.1\text{ }\mu \text{m},\rho _c=8.59\text{ }\mu \text{m}^3`$ and $`M=0.0065`$.

## III Gap Solitons

The spatially modulated coupling between the optical fields and the condensate induces a single-particle band structure with regions of negative effective mass. As mentioned in the introduction, this leads to the possibility of bright atomic solitons even for repulsive interactions . Their energy lies in forbidden gaps of the linear band structure, hence the name gap solitons.

### A Analytical Solutions

Approximate analytic expressions for the gap solitons can be obtained by expressing the spinor condensate components as

$$\psi _{\pm 1}(z,\tau )=e^{\pm ik_lz}e^{i\tau /4M}\varphi _{\pm 1}(z,\tau ),$$
(14)
where the field envelopes $`\varphi _{\pm 1}`$ are assumed to be slowly varying in space compared to $`1/k_l`$. Neglecting then the second-order spatial derivatives yields the coupled partial differential equations

$`i\left({\displaystyle \frac{}{\tau }}\pm {\displaystyle \frac{}{z}}\right)`$ $`\left(\begin{array}{c}\varphi _1\\ \varphi _1\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}\varphi _1\\ \varphi _1\end{array}\right)`$ (24)
$`\pm (|\varphi _1|^2+|\varphi _1|^2)\left(\begin{array}{c}\varphi _1\\ \varphi _1\end{array}\right),`$
where the choice $`\pm 1=\text{sgn}(g\delta ^{}/U)`$. For a red-detuned laser and our choice of $`g`$ this becomes $`\pm 1=\text{sgn}(U)`$. Aceves and Wabnitz  have shown that these dimensionless equations have the explicit two-parameter gap soliton solutions (see also Ref. )

$`\varphi _1`$ $`=`$ $`\pm {\displaystyle \frac{\mathrm{sin}(\eta )}{\beta \gamma \sqrt{2}}}\left({\displaystyle \frac{e^{2\theta }+e^{i\eta }}{e^{2\theta }+e^{\pm i\eta }}}\right)^vsech\left(\theta {\displaystyle \frac{i\eta }{2}}\right)e^{\pm i\sigma },`$ (25)
$`\varphi _1`$ $`=`$ $`{\displaystyle \frac{\beta \mathrm{sin}(\eta )}{\gamma \sqrt{2}}}\left({\displaystyle \frac{e^{2\theta }+e^{i\eta }}{e^{2\theta }+e^{\pm i\eta }}}\right)^vsech\left(\theta \pm {\displaystyle \frac{i\eta }{2}}\right)e^{\pm i\sigma },`$ (26)
with $`1<v<1`$ is a parameter which controls the soliton velocity, $`0<\eta <\pi `$ is a shape parameter, and

$$\beta =\left(\frac{1v}{1+v}\right)^{\frac{1}{4}},\gamma =\frac{1}{\sqrt{1v^2}},$$
(27)
$$\theta =\gamma \mathrm{sin}(\eta )(zv\tau ),\sigma =\gamma \mathrm{cos}(\eta )(vz\tau ).$$
(28)
Since we are interested in creating bright solitons in the presence of repulsive interactions we restrict ourselves to $`\text{sgn}(U)=+1`$, corresponding to the choice of the upper sign in the analytic solutions.

The characteristic length scale associated with the solitons is $`l_c`$, so that the approximate solitons (26) are valid for $`K_ll_c=1/2M>>1`$. This inequality is well satisfied for our choice of parameters, which gives $`1/2M=76.9`$.

### B Characteristic Properties

From the dependence of the hyperbolic–secant on $`\theta =\gamma \mathrm{sin}(\eta )(zv\tau )`$ in Eqs. (26), we identify the gap soliton parameter $`v=V_g/V_R`$ as the group velocity $`V_g`$ of the soliton in units of the recoil velocity $`V_R=l_c/t_c=\mathrm{}K_l/m`$. Since $`1<v<1`$, the magnitude of the group velocity is bounded by the recoil velocity, which is $`V_R=1.77\text{ cm/s}`$ for the present case of Sodium. From Eqs. (26), one can extract further important soliton properties, such as the number $`N_s`$ of atoms in the soliton shown in Fig. 2, and the soliton width $`W_s=w_sl_c`$ shown in Fig. 3. Specifically, the number of atoms in a particular gap soliton is given by

$$N_s=A_T𝑑Z[|\mathrm{\Psi }_1(Z,t)|^2+|\mathrm{\Psi }_1(Z,t)|^2],$$
(29)
where $`A_T`$ is the effective transverse area. Fig. 2 illustrates two general soliton properties, namely that increasing the shape parameter $`\eta `$ increases the atom number, but that faster solitons have lower atom number. Similarly, Fig. 3 shows that increasing the shape parameter $`\eta `$ decreases the soliton width and that faster solitons have narrower widths.

Given that the analytic gap soliton solutions (26) hold for broad envelopes, we confine our attention to the case $`\eta <1`$, with atom numbers in the range $`N_s10^410^5`$, and soliton widths $`W_s40100`$ $`\mu `$m for the parameters at hand.

We can gain further insight into the structure of the gap solitons by taking the more extreme limit $`\eta <<1`$ of Eqs. (26). Using the definitions (14), and returning to dimensional units we have then

$`\mathrm{\Psi }_1(Z,0)`$ $`=`$ $`{\displaystyle \frac{\eta }{\beta \gamma }}\sqrt{{\displaystyle \frac{\rho _c}{2}}}sech(Z/W_0)(1)^ve^{i(K+K_l)Z},`$ (30)
$`\mathrm{\Psi }_1(Z,0)`$ $`=`$ $`{\displaystyle \frac{\beta \eta }{\gamma }}\sqrt{{\displaystyle \frac{\rho _c}{2}}}sech(Z/W_0)(1)^ve^{i[(KK_l)Z+\pi ]},`$ (31)
where

$$W_s=3.44W_0=3.44\left(\frac{l_c\sqrt{1v^2}}{\eta }\right),K=\frac{\gamma v}{l_c}.$$
(32)
Here $`W_s=3.44W_0`$ is the soliton width, the factor $`3.44`$ being the numerical conversion from the width of the hyperbolic secant to the $`1/e^2`$ width of the distribution, and $`K=k/l_c`$ a velocity-dependent wave vector shift. This expression agrees well with the features displayed in Fig. 3, in that the width decreases with increasing $`\eta `$ and $`v`$. The soliton atom number obtained by combining Eqs. (29) and (31),

$$N_s=2(A_Tl_c\rho _c)\eta (1v^2),$$
(33)
correctly predicts the scaling properties of Fig. 2.

An essential point to keep in mind is that the gap solitons are coherent superpositions of the two Zeeman sublevels, and the approximate solutions (31) contain important information on the phase and amplitude relations that need to be created between them to successfully excite and manipulate gap solitons. In particular, they show that there is always a spatially homogeneous $`\pi `$ phase difference between the two states. In addition, the two components have the spatial wave vectors

$$K_{\pm 1}=K\pm K_l,$$
(34)
with $`K=\gamma v/l_c`$. For our parameters, $`K_l=6.38`$ $`\mu `$m-1 and $`|K|<1/l_c=0.08`$ $`\mu `$m-1. Despite the fact that it is so small, $`|K|`$ is an important factor since it controls the soliton velocity. Finally, it follows from dividing the amplitudes of the two components that

$$\left|\frac{\mathrm{\Psi }_1(Z,t)}{\mathrm{\Psi }_1(Z,t)}\right|^2=\frac{1}{\beta ^4}=\left(\frac{1+v}{1v}\right),$$
(35)
which shows that their relative occupation depends on the soliton velocity parameter $`v`$. For $`v=0`$ the sublevels are equally populated, but as $`v1`$ the $`|1>`$ sublevel has a larger population, and vice versa for $`v1`$.

The characteristic time scale for the evolution of the gap solitons can be determined from the plane-wave exponential factors in Eqs. (26). Converting to dimensional form the soliton period $`t_s`$ is defined as the time to accumulate a $`2\pi `$ phase, or in the limit $`\eta 0`$

$$t_s=2\pi t_c\sqrt{1v^2}.$$
(36)
Physically, $`t_s`$ corresponds to the internal time scale for the gap soliton. In order to observe a soliton-like behavior, it is therefore necessary to investigate the atomic propagation over several periods.

We conclude this section by noting that Aceves and Wabnitz  have shown that the gap soliton solutions (26) are stable solutions of Eqs. (24) in that they remain intact during propagation, even when perturbed away from the exact solutions. However, one should remember Eqs. (24) are only an approximation to the exact system of Eqs. (10) so that in general the gap solitons are solitary wave solutions only. As such, they are not guaranteed to be absolutely stable.

## IV Gap soliton control

### A State manipulation

Summarizing the previous section, gap solitons require the right population in each Zeeman sublevel, a phase difference of $`\pi `$ between these sublevels, and appropriate plane-wave factors $`e^{iK_{\pm 1}Z}`$. The shapes of the Hartree wave functions corresponding to the two Zeeman sublevels are hyperbolic secant, which we approximate by a Gaussian in the following. They could for example be initialized in an optical dipole trap . Manipulating the gap solitons is therefore reduced to the problem of controlling the populations and phases throughout the spinor condensate. This is achieved via a magneto-optical control scheme involving a combination of pulsed coherent optical coupling and of phase-imprinting using spatially inhomogeneous magnetic fields.

The coherent optical coupling can be achieved e.g. by a laser pulse of frequency $`\omega _l`$ propagating perpendicularly to the $`Z`$-axis and with linear polarization perpendicular to that axis. For sufficiently short pulses, one can neglect changes in the center-of-mass motion of the atoms during its duration, leading to a very simple description. We assume for simplicity a plane-wave rectangular pulse of duration $`t_p`$ and of spatial extent large compared to the soliton. The Hamiltonian describing the coupling between this pulse and the condensate is then the same as in Eq. (1) without the linear momentum exchange terms $`\mathrm{exp}(\pm 2iK_lZ)`$ and the kinetic energy term, and with $`\delta ^{}\delta _p^{}`$. The state of the system after the pulse is then easily found to be

$`\left(\begin{array}{c}\mathrm{\Psi }_1(t_p)\\ \mathrm{\Psi }_1(t_p)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{cos}\chi & i\mathrm{sin}\chi \\ i\mathrm{sin}\chi & \mathrm{cos}\chi \end{array}\right)\left(\begin{array}{c}\mathrm{\Psi }_1(0)\\ \mathrm{\Psi }_1(0)\end{array}\right)`$ (43)
$``$ $`M_L(\chi )\left(\begin{array}{c}\mathrm{\Psi }_1(0)\\ \mathrm{\Psi }_1(0)\end{array}\right).`$ (46)
where $`\chi =g\delta _p^{}t_p`$ is the excitation pulse area and the operator $`M_L`$ can be used to control the population transfer by an appropriate choice of $`\chi `$.

The required phase relation between the two states can be achieved via Zeeman splitting. Considering for concreteness a spatially inhomogeneous rectangular magnetic field pulse of duration $`t_B`$ we have, neglecting again all other effects,

$$i\mathrm{}\frac{}{t}\mathrm{\Psi }_{\pm 1}(Z,t)=\pm \mu _Bg_F\left(B_0+B^{}Z\right)\mathrm{\Psi }_{\pm 1}(Z,t),$$
(47)
where $`g_F`$ is the Landé g–factor of hyperfine ground state, $`\mu _B`$ is the Bohr magneton, $`B_0`$ the spatially homogeneous component of the magnetic field, and $`B^{}`$ its gradient, the direction of the magnetic field being along the $`Z`$-axis. The application of this field results in the state

$`\left(\begin{array}{c}\mathrm{\Psi }_1(t_B)\\ \mathrm{\Psi }_1(t_B)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}e^{i(\vartheta +K_BZ)}& 0\\ 0& e^{i(\vartheta +K_BZ)}\end{array}\right)\left(\begin{array}{c}\mathrm{\Psi }_1(0)\\ \mathrm{\Psi }_1(0)\end{array}\right)`$ (54)
$``$ $`M_B(\vartheta ,K_B)\left(\begin{array}{c}\mathrm{\Psi }_1(0)\\ \mathrm{\Psi }_1(0)\end{array}\right),`$ (57)
where $`\vartheta =(\mu _Bg_F/\mathrm{})B_0t_B`$ and $`K_B=(\mu _Bg_F/\mathrm{})B^{}t_B`$ are the imprinted phase shift and phase gradient (or wave vector), respectively. That is, the application of the magnetic pulse results in a phase difference of $`2\vartheta `$ between the two Zeeman sublevels, and in addition it imparts them wave vectors $`\pm K_B`$.

We note that although we have treated the pulsed excitations above in an impulsive manner to illustrate their action, our simulations describe the actions of the pulses correctly to the equations of motion (10). The numerics confirm the accuracy of the impulsive approximation for the parameters at consequence of the fact that we consider pulse durations significantly shorter than the soliton period (36) of $`t_s=4.3`$ ms.

### B Gap soliton excitation

To illustrate how stationary and moving solitons can be excited using the proposed magneto-optical state scheme, we start from a scalar condensate in the $`|1`$ state, $`\mathrm{\Psi }(Z,0)=[0,\mathrm{\Psi }_0(z)]^T`$, with spatial mode

$$\mathrm{\Psi }_0(Z)=\frac{N_s}{\sqrt{A_T}}\left(\frac{2}{\pi W_s^2}\right)^{1/4}e^{Z^2/W_s^2},$$
(58)
with $`W_s`$ and $`N_s`$ the width and atom number of the gap soliton desired, see Figs. 2 and 3. This Gaussian is chosen to approximate the hyperbolic-secant structure of the analytic gap soliton solution in Eq. (31). For a stationary solution we need to prepare the Zeeman sublevels with equal population and with a $`\pi `$ phase difference. We further need to impose wave vectors which are equal in magnitude but opposite in sign $`K_{\pm 1}=\pm K_l`$. This can be achieved by applying a laser pulse with area $`\chi =\pi /4`$, followed by a magnetic pulse with $`\vartheta =\pi /4`$ and $`K_B=K_l`$. The state then transforms as $`(t=t_p+t_B)`$

$`\mathrm{\Psi }(Z,t)`$ $`=`$ $`M_B(\pi /4,K_l)M_L(\pi /4)\mathrm{\Psi }(Z,0)`$ (59)
$`=`$ $`{\displaystyle \frac{e^{\frac{3i\pi }{4}}}{\sqrt{2}}}\left(\begin{array}{c}e^{iK_lZ}\\ e^{i\pi }e^{iK_lZ}\end{array}\right)\mathrm{\Psi }_0(Z).`$ (62)
Using the atomic parameters of Sodium, this situation can be realized for a $`10\text{ }\mu \text{s}`$ light pulse with $`I=2.69\text{ KW/cm}^2`$, and a $`200\text{ }\mu \text{s}`$ magnetic field pulse with $`B_0=0.89\text{ mG}`$ and $`B=72.5\text{ G/cm}`$ . Fig. 4 shows the resulting stable evolution of the total density $`\left|\mathrm{\Psi }(Z,t)\right|^2`$. As a result of the Gaussian approximation to the exact solution there are some slow oscillations imposed on the motion, but the solution remains centered at $`Z=0`$ and stationary over a time $`t=200`$ ms, much longer than the soliton period $`t_s=4.3`$ ms.

The excitation of a moving soliton is slightly more complicated since the velocity dependent wave vector $`K=\gamma v/l_c`$ in Eq. (31) is no longer zero. To deal with this situation, we first impart the wave vector $`K`$ to the initial $`|1`$ state in Eq. (58) using a magnetic pulse with $`\vartheta =0,K_B=K`$, which for Sodium and a $`200`$ $`\mu `$s magnetic pulse can be realized using $`B_0=0,B^{}=+0.094`$ G/cm. Here we specifically take $`v=0.1`$, so that $`K=8306`$ m-1. We further recall that moving solitons must have unequal populations of the two Zeeman sublevels. For $`v=0.1`$, we have from Eq. (35) that the ratio between the $`|1`$ and $`|1`$ populations should be $`1.22`$. This is achieved by a coherent optical coupling with $`\chi =0.835`$, and $`I=2.69\text{ KW/cm}^2`$, corresponding to a pulse duration of $`10.6`$ $`\mu `$s. Finally, we impart the wave vectors $`\pm K_l`$ and the $`\pi `$ phase difference between the states the Zeeman sublevels with a magnetic pulse, as for the stationary soliton above. Summarizing, the full magneto-optical control sequence is described at $`(t=2t_B+t_p)`$ by

$`\mathrm{\Psi }(Z,t)`$ $`=`$ $`M_B(\pi /4,K_l)M_L(\chi )M_B(0,K)\mathrm{\Psi }(Z,0)`$ (63)
$`=`$ $`e^{\frac{3i\pi }{4}}\left(\begin{array}{c}\mathrm{sin}\chi e^{i(K+K_lZ)}\\ \mathrm{cos}\chi e^{i\pi }e^{i(KK_lZ)}\end{array}\right)\mathrm{\Psi }_0(Z).`$ (66)
Fig. 5 shows the resulting numerical simulation of a gap soliton with $`v=0.1`$ ($`V_g=0.18`$ cm/s), illustrating its stable propagation is exhibited over many soliton periods. We have used the same scheme to launch gap solitons over the full range of velocities.

### C Soliton splitting

The next application of magneto-optical control that we consider is soliton splitting. To achieve this goal, we take advantage of the fact that for fast solitons almost all of the population is in one Zeeman sublevel, and the other state can be viewed as a small perturbation. For example, for $`v=0.5`$ the ratio of the populations is already $`6`$. This implies that the relative phase between the two states is no longer of importance, so we need only concentrate on getting the plane-wave factors right.

Assume for concreteness that we start from an initial condition with $`N_0`$ atoms in the $`|1`$ state, and apply a laser pulse of area $`\pi /4`$ to transfer half the population to the $`|+`$ state. We can then apply a magnetic pulse to impose the wave vectors $`\pm (K_l+K)`$, to the $`|\pm 1`$ states, with $`K=\gamma v/l_c`$ corresponding to a given velocity $`|v|`$. Now if we choose $`N_0=2N_s`$, with $`N_s`$ the atom number for that velocity, then for $`|v|`$ large enough we may expect to see oppositely moving solitons emerge from the initial state. This is illustrated in Fig. 6, which clearly illustrates the emergence of two solitons with opposite velocities. Fig. 7, which shows the density profiles for the individual Zeeman sublevels for a time $`t=82`$ ms, confirms that that each gap soliton indeed comprises a dominant Zeeman sublevel plus a small component of the other state. During the early stages of propagation, the solitons rearrange their phase and shape before settling down. This is accompanied by some slowing down and the familiar shedding of “radiation.” The emerging solitons are therefore slower than their “design velocity.” For the specific example of Fig. 6, the actual velocity is found to be $`1.08\text{ cm/s}`$.

### D Soliton reversal

In addition to offering the possibility of exciting moving gap solitons, magneto-optical control can also reverse their direction of propagation along $`Z`$. This can again be achieved by making use of the fact that for fast solitons almost all the population is in one Zeeman sublevel, so that we need only concentrate on getting the plane-wave factors right. In the numerical simulation of Fig. 8 we create a gap soliton with $`v=0.59`$ ($`V_g=1`$ cm/s) and let it propagate for $`100`$ ms. This soliton consists essentially of a plane-wave factor $`e^{i(K+K_l)z}`$ multiplying state $`|1`$. At $`t=100`$ ms we apply a laser pulse of area of $`\pi /2`$ to transfer all the population to state $`|1`$ and change the phase to $`e^{i(K+K_l)z}`$ using a magnetic field pulse. This results in the same soliton as before the control sequence, but with opposite velocity. Note however the loss of some atomic population to “radiation” in the process.

## V Atomic Mach–Zehnder Interferometer

As an illustration of the potential use of gap solitons employing magneto-optical control here we consider a nonlinear atomic Mach-Zehnder interferometer. Solitons present some advantages for atom interferometry in that they are many-atom wavepackets which are immune to the effects of spreading, hence allowing longer path lengths, and also reduced signal-to-noise for large atom numbers. Typically many-body effects limit the utility of high-density wavepackets due to spatially varying mean-field phase shifts, but solitons have the cardinal virtue that they have fixed spatial phase variations (for our situation this applies for faster solitons). Thus solitons may provide a key to making maximal use of high density sources for atom interferometry. Indeed they have long been advocated for all-optical switching applications due to these very properties.

Our specific demonstration of a nonlinear atomic interferometer based on solitons involves an initial scalar condensate that is split into two oppositely moving solitons along $`Z`$, see Fig. 9 for $`t<60`$ ms. At $`t=60`$ ms laser and magnetic pulses are applied which act to reverse the direction of the two solitons. The process of reversing causes some loss of atoms in both solitons as before. The two reversed soliton components come together again at $`t=120`$ ms. Since the colliding solitons are predominantly in opposite orthogonal Zeeman sublevels, the interference pattern appearing during the collision is due to the contamination of each soliton by the other state, see e.g. Fig. 7). Fig. 10 shows the interference at the soliton collision in the total density (solid line) and also the individual Zeeman sublevels, with a fringe contrast of around 30%. These results demonstrate the potential use of gap solitons for realizing nonlinear atom interferometers with high brightness sources.

## VI Summary and outlook

Employing a spinor rather than a scalar condensate gives the opportunity to externally manipulate bright atomic solitons by conceptually simple magneto-optical methods without losing the stable solitonic behavior. Specifically, the characteristic properties of the gap soliton solutions lead to realistic manipulation schemes that have been demonstrated explicitly in a one-dimensional situation. This scheme has proven successful in exciting solitons with different velocities at a minimal atom loss rate as well as splitting, and reversing their direction of propagation. The combination of these techniques resulted in the demonstration of an atomic Mach–Zehnder interferometer, which might be of interest in atom optical sensors. We remark that although this was not explicitly discussed here, we have also found numerical evidence for stable three–dimensional gap solitons, when adding an optical potential for confinement in the transverse direction.

It is now becoming amply evident that much of the future of atom optics lies in integrated systems, or ”atom optics on a chip” . Optical and magnetic waveguides and beam splitters have recently been demonstrated, but the coupling of a condensate into these guides remains a major experimental challenge. The use of bright solitons, with their potential to transport bright matter waves in a controlled fashion, might offer a solution to this problem, and merits further investigations. Further topics of interest include the manipulation and control of bright solitons using phase imprinting methods that employ the AC–Stark effect, as have already been realized for dark solitons.

###### Acknowledgements.

We have benefited from numerous discussions with M. G. Moore and P. S. Jessen. The authors would also like to thank J. V. Moloney for CPU time. This work is supported in part by the U.S. Office of Naval Research under Contracts No. 14-91-J1205 and No. 14-99-1-0806, by the National Science Foundation under Grant No. PHY98-01099, by the U.S. Army Research Office, and by the Joint Services Optics Program. S.P. was partly supported by the German National Merit Foundation.

## A Four–Photon Coupling

As mentioned in the Introduction, the Raman coupling between the $`|F_g=1,m_g=1`$ and $`|F_g=1,m_g=1`$ Zeeman sublevels vanishes in the case of far-off resonance light, a result of destructive interference with other $`3P_{3/2}`$ hyperfine excited states . In this limit the system reduces to an effective $`|J_g=1/2|J_e=3/2`$ transition where no $`\mathrm{\Delta }m=\pm 2`$ transitions are allowed. This difficulty can be overcome by adding a $`\pi `$–polarized laser beam from the side (e.g. along the $`X`$-axis), which will not transfer any momentum along the transverse direction. A proper choice of detunings allows one to use another hyperfine ground state to mediate the $`\mathrm{\Delta }m=\pm 2`$ transition, without populating it significantly.

One disadvantage of this approach is that since this is now a four-photon process, the intensity of the lasers must be increased while still avoiding spontaneous emission. Choosing a frequency difference between the $`\pi `$–light and the $`\sigma `$-light a few MHz larger than the splitting between the hyperfine ground states as shown in Fig. 11, three time scales govern the dynamics of the system. They are determined by (a) the coupling of the ground state to the excited states, with detunings in the THz range; (b) the effective coupling between the hyperfine ground states, with detunings in the GHz resp. MHz range; (c) the desired effective $`\mathrm{\Delta }m=\pm 2`$ coupling. As it turns out, coupling the $`|F_g=2,m_g=1`$ and $`|F_g=2,m_g=1`$ levels is a better choice than staying in the $`F_g=1`$ manifold, as this minimizes the coupling to the outer states $`|F_g=2,m_g=\pm 2`$. Employing the intermediate $`|F_g=1,m_g=0`$ state makes the system aware of the hyperfine structure and leads to a non–vanishing coupling.

Separating the time scales leads to an effective coherent evolution between the $`|F_g=2,m_g=\pm 1`$ states. It takes the same form as in Eq. (1), but with

$`g`$ $``$ $`g_{\mathrm{𝑒𝑓𝑓}}=1/192,`$ (A1)
$`\delta ^{}`$ $``$ $`\delta _{\mathrm{𝑒𝑓𝑓}}^{}={\displaystyle \frac{𝒟^4_\sigma ^2_\pi ^2}{\mathrm{}^4\delta ^2\mathrm{\Delta }}},`$ (A2)
where $`_\sigma `$ is the amplitude for the $`\sigma `$–light resp. $`_\pi `$ for the $`\pi `$–light and the detunings $`\delta `$ and $`\mathrm{\Delta }`$ are defined in Fig. 11.

The value of $`g_{\mathrm{𝑒𝑓𝑓}}`$ involves Clebsch–Gordan and Wigner–6J coefficients for all hyperfine transitions and $`\delta _{\mathrm{𝑒𝑓𝑓}}^{}`$ is clearly a four–photon term due to the product of four electric field envelopes. Choosing the wavelength of the $`\sigma `$–light $`\lambda _\delta =985\text{ nm}`$ as before, $`\mathrm{\Delta }=\omega _\pi \omega _\sigma =1.85\text{ GHz}`$ and laser intensities $`I_\sigma =100\text{ KW/cm}^2`$ and $`I_\pi =158.3\text{ KW/cm}^2`$ we obtain an effective coupling $`1/(g_{\mathrm{𝑒𝑓𝑓}}\delta _{\mathrm{𝑒𝑓𝑓}}^{})=685.0\text{ }\mu \text{s}`$, as in Sec. II.

Numerical simulations of the full level scheme of Sodium is in excellent agreement with the derived coupling in Eq. (A1) and shows that the intermediate level $`|F_g=1,m_g=0`$ only contains around 0.1 percent of the atomic population. All other levels are negligibly populated, they are indeed far–off resonance. Since $`|F_g=1,m_g=0`$ is a ground state, spontaneous emission is not an issue.

The intensities might be reduced if $`\delta `$ were chosen smaller, but that would of course alter the properties of the soliton, since they depend on the wavelength of the optical fields. Note that we cannot apply permanent magnetic fields to remove the degeneracy of the hyperfine ground states and thus push non–participating states out of resonance since this would destroy the solitons. The light shifts induced by the lasers are the same for each manifold and will not remove this degeneracy.


text:
# 1 Introduction

## 1 Introduction

Exact inflationary solutions  are of interest because they allow one to study physical effects in a simple manner than would otherwise be allowed from numerical studies. In addition, they allow one to write down exact forms for decaying cosmological constant models  as we shall see below. It is thus always interesting if one can think of new methods for developing exact solutions to inflation models. The results of the present work are threefold. Firstly a new method for generating exact inflationary solutions is given. Secondly a new inflation solution is presented. Thirdly, from this exact solution, a new form of decaying cosmological constant is presented.

Exact inflationary solutions have been studied before  and we would like to point out the difference between previous work and the present work. It is well known that inflationary potentials of the form $`V(\varphi )e^{\pm \varphi }`$ can be solved exactly and these types of models were studied in . This type of potential is considered in the present work (Section 3.1), but only as a check on the more general solution method that we develop. The works by Lidsey , Ellis and Madsen  and Lucchin and Matarrese  consider more complicated potentials and are all similar in that they specify either $`H\frac{\dot{a}}{a}`$  or $`a(t)`$  initially and then go on to derive the corresponding $`V(\varphi )`$. Ideally however, one would wish to do this the other way around, i.e. specify a form for $`V(\varphi )`$ initially and then find the corresponding expressions for $`a(t)`$ and $`\varphi (t)`$. Our present work does not quite achieve this, but we come closer to this goal than previous work . We specify a function $`F(\varphi )`$ (see later) which gives $`V(\varphi )`$. Then we are able to solve for $`a(t)`$ and $`\varphi (t)`$. As mentioned above the other works  do this the other way around, which is less satisfactory.

The inflationary scenario holds that at a very early time the universe underwent a rapid period of quasi-exponential expansion. This idea was proposed as a way of eliminating the necessity of fine-tuning of initial conditions in the standard hot big bang model, as this rapid expansion would account for the flatness of the observed universe and also explain why regions of the universe that appear causally disconnected are in thermal equilibrium . The model that is studied here assumes the interaction of a classical scalar (inflaton) field. It also assumes that on large scales the universe has a uniform energy density and that the structure is the same in every direction. These assumptions of homogeneity and isotropy and the equations of General Relativity give rise to the so called Friedmann models which seem to describe the large-scale behavior of our universe. These models allow for the inclusion of a cosmological constant term $`\mathrm{\Lambda }`$, however in this study it is assumed that this term is identically zero and that the expansion is driven by interaction with the scalar field (i.e. that the energy density of the scalar field can act as an effective cosmological term).

The size of the universe is given in terms of the scale factor $`a`$. In the Friedmann-Robertson-Walker cosmological model the evolution of the scale factor $`a`$ is specified by

$$H^2\left(\frac{\dot{a}}{a}\right)^2=\frac{8\pi G}{3}\rho \frac{k}{a^2}+\frac{\mathrm{\Lambda }}{3},$$
(1)
where $`H`$ is the Hubble parameter and overdots indicate time derivatives and $`\rho `$ is the energy density. The curvature term $`\frac{k}{a^2}`$ will be neglected in this model for the sake of simplicity and also because this term will quickly approach zero for any inflationary solution. In models whose only matter source is a scalar field $`\varphi `$ the energy density has the form

$$\rho _\varphi =\frac{1}{2}\dot{\varphi }^2+V,$$
(2)
where $`VV(\varphi )`$ is the interaction potential. During a period where the energy density of the scalar field dominates, equation (1) takes the form

$$H^2\left(\frac{\dot{a}}{a}\right)^2=\frac{8\pi G}{3}\left(\frac{1}{2}\dot{\varphi }^2+V\right)$$
(3)
giving the differential equation determining $`a(t)`$ as

$$H\frac{\dot{a}}{a}=\pm \frac{g}{3}\sqrt{\dot{\varphi }^2+2V}$$
(4)
where

$$g\sqrt{12\pi G}$$
(5)
The $`\pm `$ signs in (4) indicate an expanding or contracting solution for $`a(t)`$. The equation which specifies the time evolution of the field density is

$$\ddot{\varphi }+3H\dot{\varphi }+V^{}=0$$
(6)
with $`V^{}\frac{dV}{d\varphi }`$. Equations (4) and (6) form a set of coupled equations which specify the expansion in this model. Eliminating $`H\frac{\dot{a}}{a}`$ yields

$$\ddot{\varphi }+g\dot{\varphi }\sqrt{\dot{\varphi ^2}+2V}+V^{}=0$$
(7)
where we have included only the expanding solution from (4). The contracting solution would involve a $``$ sign in front of the square root in the above equation. In principle one could specify a potential $`VV(\varphi )`$ and solve this equation for $`\varphi (t)`$. One could then find the functional form of all other dynamical quantities, i.e. $`a(t)`$, $`\rho (a)`$, etc. . In practice this is difficult and few exact solutions have been found.

## 2 General Solution

In this paper a new closed-form solution is presented and the method of solution is described in detail as it may provide opportunities to find additional exact solutions. The solution is obtained by finding equivalent coupled equations that are written in terms of an arbitrary function $`F`$. One can specify $`F(\varphi )`$ to be any convenient function. This will determine the functional form of the interaction potential, which together with $`F`$, determine the scalar field density $`\varphi (t)`$.

Defining the function $`FF(\varphi )`$ via

$$\dot{\varphi }\pm \sqrt{(F1)V}$$
(8)
gives equation (7) as

$$\ddot{\varphi }\pm gV\sqrt{F^21}+V^{}=0$$
(9)
Equations (8) and (9) form a set of coupled equations which are equivalent to equation (7). They can be uncoupled by differentiating (8) with respect to time and substituting for $`\ddot{\varphi }`$ in (9). Solving for $`\ddot{\varphi }=\frac{1}{2}\frac{d\dot{\varphi }^2}{d\varphi }`$ from (8) we have

$$\ddot{\varphi }=\frac{1}{2}\left[(F1)V^{}+VF^{}\right]$$
(10)
(with $`F^{}\frac{dF}{d\varphi }`$ ) and inserting this expression for $`\ddot{\varphi }`$ into (9) leads to

$$(F+1)V^{}+VF^{}\pm 2gV\sqrt{F^21}=0$$
(11)
It would be nice if we could just specify the interaction potential $`V(\varphi )`$ and solve this first-order differential equation for $`F(\varphi )`$. However in the absence of such a solution, we instead choose the less convenient method of doing it the other way around, i.e. by specifying $`F(\varphi )`$ and then solving for $`V(\varphi )`$. It can be seen that if one chooses $`FF(\varphi )`$, then equation (11) is always separable and the potential is given by

$$V=\beta \mathrm{exp}\left(\frac{F^{}\pm 2g\sqrt{F^21}}{F+1}𝑑\varphi \right)$$
(12)
with $`FF(\varphi )`$ and $`\beta `$ is a constant. This may be simplified to

$$V=B\mathrm{exp}\left(2g\sqrt{\frac{F1}{F+1}}𝑑\varphi \right)$$
(13)
where $`B`$ is a constant. It is important to note that equations (8), (9), (11), (12), (13) are consistent only if one takes either all the upper signs or all the lower signs. One cannot take upper signs in one equation and lower signs in another.

The choice of function $`F`$ is arbitrary. When one specifies $`F(\varphi )`$ then one has the potential $`V(\varphi )`$ as given in the above equation and also the solution for $`\varphi (t)`$ and $`a(t)`$ as given in equations (8) and (4). In other words one has a general method of finding exact inflationary solutions. This expression (13) for $`V(\varphi )`$ and the procedure for getting $`\varphi (t)`$ by solving (8) is the first main result of the present paper.

## 3 Special Cases

### 3.1 $`F=constant`$

It should be noted, for example, that if one chooses $`F=constant`$, then one is led to the solution for $`\varphi (t)`$ given by Barrow , although the constants could be chosen differently. Equation(13) gives

$$V=B\mathrm{exp}\left(2g\sqrt{\frac{F1}{F+1}}\varphi \right)$$
(14)
with solutions from (8) of the form

$$\varphi (t)=\pm \frac{1}{g}\sqrt{\frac{F+1}{F1}}\mathrm{ln}\left[\pm \frac{g(F1)\sqrt{B}}{\sqrt{F+1}}(tC)\right]$$
(15)
with $`C`$ being a constant and the upper or lower signs in (15) in front of the term $`\frac{1}{g}`$ correspond to the upper or lower signs in (14) or (8). However the $`\pm `$ signs occurring in the square brackets occur regardless of what sign is taken in (14) or (8). That is, (14) or (8) have two solutions for each choice of sign, the two solutions corresponding to the two signs in the square brackets. Equation (15) is simply the Barrow solution . One can now easily verify that (15) and (14) satisfy (7).

Substituting (14) and (15) into (4) yields

$$a(t)t^{\frac{1}{3}\left(\frac{F+1}{F1}\right)}$$
(16)
and eliminating $`\varphi (t)`$ in favor of $`t`$ using (15) gives, from (2),

$$\rho (a)\frac{1}{a^{6\left(\frac{F1}{F+1}\right)}}$$
(17)
which is of the form of a decaying cosmological constant . The direct relation between $`V(\varphi )`$ and $`\rho (a)`$ is perhaps made clearer if we define a constant $`\lambda `$ via

$$V(\varphi )B\mathrm{exp}(\lambda \varphi )$$
(18)
which, upon comparing to (14) is

$$\lambda 2g\sqrt{\frac{F1}{F+1}}$$
(19)
so that

$$a(t)t^{\frac{4g^2}{3\lambda ^2}}$$
(20)
giving

$$\rho (a)\frac{1}{a^{\frac{3\lambda ^2}{2g^2}}}$$
(21)
which is consistent with equation (37) of reference . If $`\lambda `$ is determined by say some other particle physics theory then $`\rho (a)`$ is specified. In the next section we pick $`\lambda =g`$ for simplicity. Doing that here would give

$$a(t)t^{\frac{4}{3}}$$
(22)
and

$$\rho (a)\frac{1}{a^{\frac{3}{2}}}$$
(23)
thus showing that the Barrow solution results in a power law decaying cosmological constant. Such a power law decay is the assumption in many decaying cosmological constant models.

### 3.2 $`F=\mathrm{cosh}(\lambda \varphi )`$

One of the other main results of the present paper is the exact solution for $`\varphi (t)`$ resulting from the choice $`F\mathrm{cosh}(\lambda \varphi )`$. Then equation (13) gives

$$V(\varphi )=C\left(1+\mathrm{cosh}\lambda \varphi \right)^{2g^{}1}$$
(24)
where $`C`$ is a constant and with the $`\pm `$ sign corresponding to the $`\pm `$ sign in (8) and the expanding solution (7) and

$$g^{}\frac{g}{\lambda }$$
(25)
Thus equation (8) becomes

$$\dot{\varphi }=\pm \sqrt{C}\frac{\mathrm{sinh}\lambda \varphi }{(1+\mathrm{cosh}\lambda \varphi )^{1\pm g^{}}}$$
(26)
This equation is only consistent if all upper or all lower signs are taken, i.e. one should not mix upper and lower signs. This equation is most easily solved if the denominator disappears, i.e. the constant $`\lambda `$ is chosen so that $`\lambda =g\sqrt{12\pi G}`$ and with the lower signs, corresponding to the lower sign in definition (8). Thus for the rest of this paper, whenever $`\lambda `$ appears it is taken to be $`\lambda =g`$.  The upper $``$ sign in equation(24), corresponding to the upper $`+`$ sign in definition (8) gives

$$V(\varphi )=\frac{C}{(1+\mathrm{cosh}\lambda \varphi )^3}$$
(27)
and the lower $`+`$ sign in (24) corresponding to the lower $``$ sign in definition (8) gives

$$V(\varphi )=C(1+\mathrm{cosh}\lambda \varphi ).$$
(28)
We will study this latter solution (28) whose graph of $`V(\varphi )`$ versus $`\varphi `$ looks similar to a parabola. Using this result (28) in equation (8) leads to

$$\dot{\varphi }=\sqrt{C}\mathrm{sinh}\lambda \varphi $$
(29)
Solving for $`\varphi (t)`$ yields

$`\varphi (t)={\displaystyle \frac{2}{\lambda }}\mathrm{coth}^1\left\{\mathrm{exp}\left[\lambda \sqrt{C}(tD)\right]\right\}`$
$`={\displaystyle \frac{1}{\lambda }}\mathrm{ln}\left\{{\displaystyle \frac{\mathrm{exp}\left[\lambda \sqrt{C}(tD)\right]+1}{\mathrm{exp}\left[\lambda \sqrt{C}(tD)\right]1}}\right\}`$ (30)
where $`D`$ is a constant. The form of the potential suggests that $`\varphi (t)`$ be a function that decreases from an initial maximum value similar to the chaotic inflation model. We can choose $`D=0`$ and so (30) can be written as

$$\varphi (t)=\frac{1}{\lambda }\mathrm{ln}\left[\frac{\mathrm{exp}(\lambda \sqrt{C}t)+1}{\mathrm{exp}(\lambda \sqrt{C}t)1}\right].$$
(31)
This expression for $`\varphi (t)`$ is just a monotonically decreasing function of $`t`$. Again one can easily verify that (31) and (28) satisfy (7).

The form of the potential (28) is similar to the chaotic inflation potential, $`V(\varphi )=\frac{1}{2}m^2\varphi ^2`$, in that the initial value of the scalar field determines how the expansion proceeds and also that the true vacuum state (global potential minimum) is located at $`\varphi =0`$ . For the chaotic inflation potential the value of the scalar field oscillates around this minimum with $`\varphi (t)`$ taking negative values. The hyperbolic potential given here does not lead to this behavior, instead the scalar field (31) smoothly approaches zero at late times.

This solution can be used to determine the evolution of the scale factor with time by solving equation (4) for the expanding ($`+`$ square root) solution. Substituting (28) and (29) into (4) yields

$$a(t)=\left[\mathrm{exp}(2\lambda \sqrt{C}t)1\right]^{\frac{1}{3}},$$
(32)
for $`a(0)=0`$.

Equations (28), (31) and (32) constitute the second main result of the present paper, in that a new exact inflationary solution has been given for the potential specified by $`V(\varphi )=C(1+\mathrm{cosh}\lambda \varphi )`$.

Many authors have studied models where the cosmological term decays as the scale factor increases and have focused on examples where the functional relationship is one of power-law decay with one or more terms . Here the effective cosmological term arises from the energy density of the scalar field. From equations (2) and (28,31,32) we obtain a functional form for $`\rho (a)`$ namely

$$\rho (a)=2C\left(1+\frac{2}{a^3}+\frac{1}{a^6}\right)$$
(33)
(Note that $`\lambda `$ does not appear in this expression.) This is the third main result of the present paper presenting a new form of decaying cosmological constant. This differs from previous models  in that a constant is added onto two power law decay terms.

## 4 Summary

In summary, we have introduced a general formula (13) for finding exact inflationary solutions. This formula has been checked by showing that it reproduces Barrow’s solution . Secondly, a new exact solution for a hyperbolic interaction potential has been presented. Other exact solutions can be found by choosing different functional forms for $`F`$ appearing in (13), and as such should prove useful in studying other inflationary models. Finally a new decaying cosmological constant model, which includes a constant term plus two power law decay terms, results from this exact solution.

Acknowledgements
ATK is funded by the Wisconsin Space Grant Consortium and NSF grant PHY-9507740.


text:
# The Density Matrix Renormalization Group technique with periodic boundary conditions

## I Introduction

In 1992 the density matrix renormalization group (DMRG) technique in real space was invented by S. R. White  and has been mostly applied to the diagonalization of one-dimensional (1D) quantum chain spin Hamiltonians. Three years later the DMRG method was redesigned by T. Nishino  and applied to classical spin 2D models. The DMRG method for the classical models is based on renormalization of the transfer matrix. It is a variational method maximizing the partition function using a limited number of degrees of freedom and its variational state is written as a product of local matrices .

DMRG has been used for many various quantum models. It provides results with remarkable accuracy for larger systems than it is possible to study using standard diagonalization methods. The 2D classical systems treated by the DMRG method exceeds the classical Monte Carlo approach in accuracy, speed, and size of the systems . A further DMRG improvement of the classical systems is based on Baxter’s corner transfer matrix , the CTMRG , and its generalization to any dimension .

Applications of the DMRG technique for calculation of the thermodynamic properties in the 2D classical systems has been done by . Treating of non-symmetric transfer matrices or non-hermitian quantum Hamiltonians has also been studied by the DMRG technique .

It was shown that DRMG method yields very accurate estimations of ground state energy of finite quantum chains and free energy of classical strips of finite width with open boundary conditions. We have developed the DMRG method with periodic boundary conditions for strips of classical spins and shown that, similarly as for quantum chains, it gives these quantities with much less degree of accuracy. Nevertheless, the DMRG method is mostly used for prediction of physical quantities and critical properties of infinite systems in connection with finite size scaling or extrapolation of the results from finite-size systems to infinite ones. The objective of this paper is to study the DMRG and exact methods with two different boundary conditions for finite strips of various widths and compare their results with known exact results for infinite 2D system. It is shown that while for the exact diagonalization of finite-strip transfer matrices scaling properties of the system improve, for DMRG approach there exists an optimum width for each degree of approximation. The developed approach is tested on 2D Ising model.

Paper is organized as follows: in Sec. II we mention briefly the DMRG for the open boundary conditions; Sec. III contains the modification of the DMRG method for the periodic boundary conditions; in Sec. IV we present the results obtained by the DMRG method with periodic and open boundary conditions, the exact diagonalization method and show how to determine the optimum strip width for finite-size scaling, and in Sec. V the results will be summarized.

## II DMRG with open boundary conditions

The transfer matrix approach is a powerful method for exact numerical calculation of thermodynamical properties of lattice spin models defined on finite-width strips. If the width of the strip is too large and the capacity of the computer is exceeded, the DMRG method is found to be useful for an effective reduction of the transfer matrix size. It can be used for calculation of global quantities such as free energy as well as of a spatial dependence across the strip of local quantities, e.g. spin correlation functions.

The properties of an infinite strip of finite width $`L`$ are given by the solution of ‘left’ eigenvectors and corresponding eigenvalues of the transfer matrix equation

$$\underset{\{\sigma \}}{}\mathrm{\Psi }_l\left(\{\sigma \}\right)T\left(\{\sigma \}|\{\sigma ^{}\}\right)=\lambda \mathrm{\Psi }_l\left(\{\sigma ^{}\}\right),$$
(1)
where $`\{\sigma \}`$ is a set of $`L`$ spins $`\{\sigma _1,\sigma _2,\mathrm{},\sigma _L\}`$ defined on a row and $`\{\sigma ^{}\}`$ is a set of $`L`$ spins on the adjacent row. The transfer matrix is a product of Boltzmann weights given by the lattice Hamiltonian. For non-symmetric transfer matrices besides the left eigenvectors $`\mathrm{\Psi }_l^i`$, the right eigenvectors $`\mathrm{\Psi }_r^i`$ should be calculated, as well.

Reducing the size of the transfer matrix the standard DMRG technique proceeds in two regimes:
(1) In the process of iterations, the infinite system method (ISM) pushes both ends of the transfer matrix further so that each step of the ISM enlarges the lattice size by two sites. The transfer matrix (superblock) is constructed from three blocks: left $`T_l`$ and right $`T_r`$ transfer matrices (blocks) and a Boltzmann weight $`W_B`$, in particular

$$T_{[2j+2]}=T_l^{(j)}W_BT_r^{(j)},$$
(2)
where the index on the left hand side denotes the number of sites in one row of the whole superblock $`T`$ at the $`j`$th step of iteration. The Boltzmann weight usually is a function of several spins interacting among each other, e.g. for the Ising model with nearest-neighbor interactions the Boltzmann weight has the form

$$W_B(\sigma _1\sigma _2|\sigma _1^{}\sigma _2^{})=\mathrm{exp}\left\{\frac{J}{k_BT}(\sigma _1\sigma _2+\sigma _1^{}\sigma _2^{}+\sigma _1\sigma _1^{}+\sigma _2\sigma _2^{})\right\}.$$
(3)
In the first step of the ISM (for details, see ) $`T_l^{(1)}=T_r^{(1)}=W_B`$ is put. The whole procedure has $`L/21`$ steps

$$T_l^{(1)}W_BT_r^{(1)}T_l^{(2)}W_BT_r^{(2)}\mathrm{}T_l^{(L/21)}W_BT_r^{(L/21)}.$$
(4)
and stops when the desired strip width of $`L`$ sites is reached.

The first steps of the iteration scheme (4) are exact but if the superblock matrix $`T`$ becomes too large, a reduction procedure, to keep the size of superblock constant, should be introduced.

The first step of (4) introduces open conditions at the strip boundaries. If the temperature of the system is lower than the critical one and the strip width is wide enough, the symmetry of the system is spontaneously broken (order parameter becomes non-zero), and after reaching the fixed point of the iteration procedure, the system does not depend on the boundary conditions any more. The calculations with periodic boundary conditions described in the next Section give in this regime the same result as with the free ones.
(2) the finite system method (FSM) improves numerical accuracy of ISM result by left and right moves (sweeps) according to the following prescription:

$$T_l^{(L/21)}W_BT_r^{(L/21)}T_l^{(L/2)}W_BT_r^{(L/22)}\mathrm{}T_l^{(L3)}W_BT_r^{(1)},$$
(5)
$$T_l^{(L3)}W_BT_r^{(1)}T_l^{(L2)}W_BT_r^{(2)}\mathrm{}T_l^{(L/21)}W_BT_r^{(L/21)}.$$
(6)
In the right sweep (5) the left blocks $`T_l`$ are calculated in the previous step of the sweep and the right blocks $`T_r`$ are taken from the previous left sweep (in the first right sweep from ISM); similarly for the left sweep.

The values of local thermodynamical quantities given by particular superblocks in the final sweep (after the steady state is reached) are spatially dependent. The values given by the superblock in the middle of the strip $`T_l^{(L/21)}W_BT_r^{(L/21)}`$ are the closest to the bulk ones. In this sense, the best transfer matrix eigenvalues as well as eigenvectors are those of the above-mentioned central superblock. The two largest eigenvalues are used for further finite-size scaling or extrapolation treatment.

## III DMRG with periodic boundary conditions

The translational invariance of the infinite lattice is preserved in finite strips with periodic boundary conditions when strip boundaries are connected with bulk intersite interactions. In this case the strip forms an infinitely long cylinder. If the radius of the cylinder is small enough, the model can be easily solved by exact numerical diagonalization methods.

In DMRG language, imposing periodic boundary conditions means that we have to introduce properly the connection of both ends of the superblock transfer matrix $`T`$. Thus, in distinction to open-boundary case the superblock is constructed from two Boltzmann weights connecting two blocks at both ends (see Fig. 1, the rightmost diagram).

$`T_{[2j+4]}(\sigma _1\xi _l\sigma _{j+4}\sigma _{j+3}\xi _r\sigma _2|\sigma _1^{}\xi _l^{}\sigma _{j+4}^{}\sigma _{j+3}^{}\xi _r^{}\sigma _2^{})=T_l^{(j)}(\sigma _1\xi _l\sigma _{j+4}|\sigma _1^{}\xi _l^{}\sigma _{j+4}^{})`$ (7)
$`\times W_B(\sigma _{j+4}\sigma _{j+3}|\sigma _{j+4}^{}\sigma _{j+3}^{})T_r^{(j)}(\sigma _{j+3}\xi _r\sigma _2|\sigma _{j+3}^{}\xi _r^{}\sigma _2^{})W_B(\sigma _2\sigma _1|\sigma _2^{}\sigma _1^{}),`$ (8)
where the block spin variable $`\xi =\{1,2,\mathrm{},m\}`$, and the primed variables are denoted by filled circles and ovals in Fig. 1.

In the first few steps the lattice is enlarged to the desired size; no degrees-of-freedom reduction is performed and the superblock transfer matrix remains equivalent to the exact one. As depicted in Fig. 1, the ISM starts with $`T^{(6)}=T_l^{(1)}W_BT_r^{(1)}W_B`$ defined on twelve sites where $`T_l^{(1)}=T_r^{(1)}=W_BW_B`$, and one Boltzmann weight, i.e. four new sites are added in each of further steps.

If $`2^j>m`$, the number of degrees of freedom should be reduced at each $`j`$th step to keep the order of the superblock matrix constant and equal to $`2^4m^2`$.

Summation in the equation for eigenvectors (1) of the transfer matrix (8) can be performed in two steps

$`\mathrm{\Phi }(\sigma _1\xi _l\sigma _{j+4}\sigma _{j+3}\sigma _{j+3}^{}\xi _r^{}\sigma _2^{}\sigma _1^{})`$ $`=`$ $`{\displaystyle \underset{\xi _r\sigma _2}{}}T_r(\sigma _{j+3}\xi _r\sigma _2|\sigma _{j+3}^{}\xi _r^{}\sigma _2^{})W_B(\sigma _2\sigma _1|\sigma _2^{}\sigma _1^{})`$ (10)
$`\times \mathrm{\Psi }(\sigma _1\xi _l\sigma _{j+4}\sigma _{j+3}\xi _r\sigma _2)`$
$`\mathrm{\Psi }(\sigma _1^{}\xi _l^{}\sigma _{j+4}^{}\sigma _{j+3}^{}\xi _r^{}\sigma _2^{})`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{\sigma _1\xi _l}{\sigma _{j+4}\sigma _{j+3}}}{}}T_l(\sigma _1\xi _l\sigma _{j+4}|\sigma _1^{}\xi _l^{}\sigma _{j+4}^{})W_B(\sigma _{j+4}\sigma _{j+3}|\sigma _{j+4}^{}\sigma _{j+3}^{})`$ (12)
$`\times \mathrm{\Phi }(\sigma _1\xi _l\sigma _{j+4}\sigma _{j+3}\sigma _{j+3}^{}\xi _r^{}\sigma _2^{}\sigma _1^{})`$
which is depicted graphically in Fig. 2.

This procedure uses the left and right transfer matrix blocks to calculate properly the left and right eigenvectors $`\mathrm{\Psi }_l`$ and $`\mathrm{\Psi }_r`$, respectively, of the whole superblock for the periodic boundary conditions. Once we have the $`\mathrm{\Psi }_l`$ and $`\mathrm{\Psi }_r`$, the left and right density matrices can be constructed.

$$\rho _l(\xi _l\sigma _{j+4}|\xi _l^{}\sigma _{j+4}^{})=\underset{\sigma _1\sigma _{j+3}\xi _r\sigma _2}{}\mathrm{\Psi }_l(\sigma _1\xi _l\sigma _{j+4}\sigma _{j+3}\xi _r\sigma _2)\mathrm{\Psi }_r(\sigma _1\xi _l^{}\sigma _{j+4}^{}\sigma _{j+3}\xi _r\sigma _2)$$
(13)
$$\rho _r(\sigma _{j+3}\xi _r|\sigma _{j+3}^{}\xi _r^{})=\underset{\sigma _1\xi _l\sigma _{j+4}\sigma _2}{}\mathrm{\Psi }_l(\sigma _1\xi _l\sigma _{j+4}\sigma _{j+3}\xi _r\sigma _2)\mathrm{\Psi }_r(\sigma _1\xi _l\sigma _{j+4}\sigma _{j+3}^{}\xi _r^{}\sigma _2),$$
(14)
and by its complete diagonalization

$$O_l(\xi _l^{new}|\xi _l\sigma _{j+4})\rho _l(\xi _l\sigma _{j+4}|\xi _l^{}\sigma _{j+4}^{})Q_l(\xi _l^{}\sigma _{j+4}^{}|\xi _l^{new})=\omega _i\delta _{ij}$$
(15)
sets of left and right eigenvectors stored in $`O_l`$ and $`Q_l`$ matrices, respectively, is obtained (analogously, for $`O_r`$ and $`Q_r`$). The indices $`i,j`$ ($`i,j=1,2,\mathrm{},2m`$) run over all states of $`m`$-state multi-spin variable $`\xi `$ and two-state variable $`\sigma `$. For the last steps of ISM and all FSM steps half of the eigenvectors (corresponding to their lowest eigenvalues) is discarded from the matrices $`O`$ and $`Q`$, and the information of the system carried by the density matrix is reduced. However, remaining eigenvectors (if $`m`$ is large enough) usually describe the system accurately because the truncation error $`\epsilon `$ defined as

$$\epsilon =\underset{\{\mathrm{discarded}\}}{}\omega $$
(16)
is very small ($`0\epsilon 1`$). $`_{\{all\}}\omega =1`$, as the eigenvectors are assumed to be normalized. The matrices $`O`$ and $`Q`$ enter the linear transformation as projectors mapping two blocks $`T_lW_B`$ onto one block $`T_l`$ through the following procedure

$`T_l^{(j+1)}(\sigma _1\xi _l^{new}\sigma _{j+5}|\sigma _1^{}\xi _l^{new}\sigma _{j+5}^{})={\displaystyle \underset{\genfrac{}{}{0pt}{}{\xi _l^{}\sigma _{j+4}^{}}{\xi _l\sigma _{j+4}}}{}}`$ $`O_l(\xi _l^{new}|\xi _l\sigma _{j+4})T_l^{(j)}(\sigma _1\xi _l\sigma _{j+4}|\sigma _1^{}\xi _l^{}\sigma _{j+4}^{})`$ (17)
$`\times `$ $`W_B(\sigma _{j+4}\sigma _{j+5}|\sigma _{j+4}^{}\sigma _{j+5}^{})Q_l(\xi _l^{}\sigma _{j+4}^{}|\xi _l^{new}).`$ (18)
Application to the right block $`T_r`$ is straightforward. As it is seen, we calculate the blocks $`T_l`$ and $`T_r`$ separately not using the standard mirror-reflection of $`T_l`$ to $`T_r`$. This procedure is necessary when dealing with anisotropic and/or inhomogeneous systems.

The calculated new blocks $`T_l^{(j+1)}`$ and $`T_r^{(j+1)}`$ are used in the next step of the ISM for construction of the new superblock

$$T_{[2j+6]}=T_l^{(j+1)}W_BT_r^{(j+1)}W_B.$$
(19)
Within the FSM, e.g., for a sweep to the right only the left blocks are calculated and $`T_r^{(L/2k)}`$ is taken from the previous left sweep

$$T_{[L]}=T_l^{(L/23+k)}W_BT_r^{(L/21k)}W_B.$$
(20)
The variable $`k`$ (indexing the steps within a sweep) runs over the values $`(k_0,k_0+1,\mathrm{},k_01,k_0)`$, where $`2^{k_0}2^{L/22}m`$. In the process of sweeping one of the Boltzmann weight is fixed (the upper one in Fig. 1) and the second one changes its position within the interval of $`2k_0`$ lattice sites. The local physical quantities are calculated at the lattice sites of the fixed Boltzmann weight and due to the rotational invariance of the problem are valid for all the rows of the periodic lattice.

## IV Results

It is well known that the DMRG describes better a strip with open boundary conditions than that with the periodic boundary conditions  because the precision of the largest eigenvalue of the superblock matrix is increasing proportionally to $`m`$ for open boundary conditions while for periodic boundary conditions only as $`\sqrt{m}`$.

However, if we are not interested in the largest eigenvalue of a finite-strip transfer matrix but in the estimation of the free energy of the whole 2D lattice (per spin), it is more effective to use a strip with periodic boundary conditions than that with open boundaries, as demonstrated in Table 3. The results with $`m=25`$ practically exactly reproduce the exact values for $`L=16`$. The estimation of the free energy for 2D models performed by DMRG can be improved by increasing the width of the strip. For a given $`m`$ the best results are obtained for $`L\mathrm{}`$, but in this case, for $`T=2.1`$ (i.e. below the critical temperature), the symmetry of the system is spontaneously broken. Exact free energy per site $`f_{\mathrm{Onsager}}^{(\mathrm{exact})}`$ was taken from .

The critical temperature and the properties of the infinite 2D system near the critical temperature should be derived from finite-size scaling ideas (as the finite-width strip is at criticality for $`T=0`$ only).

For calculation of the critical temperature, the phenomenological renormalization approach of Nightingale  have been used. Here, the scaling properties of the correlation length, found as logarithm of the ratio of two largest eigenvalues of the exact or superblock matrix, are exploited. The product of the inverse correlation length $`K_N`$ and the strip width $`L`$ should not depend on the $`L`$ at critical temperature $`T_\mathrm{C}^{}(L)`$

$$\frac{LK_L}{(L+2)K_{L+2}}=1.$$
(21)
The accuracy of the approximate critical temperature improves with size of the strip in the case of exact diagonalization. For DMRG calculations this statement is no longer valid, as for very large $`L`$ the symmetry of the system spontaneously breaks, and the phenomenological renormalization is not applicable any more. Thus, for given order of approximation $`m`$, there exists an optimum value of the strip width $`L^{\mathrm{opt}}`$. This can be estimated from the following considerations: For exact diagonalization or DMRG calculations with $`m`$ close to $`2^{L/22}`$, the difference of the approximate critical temperature from the exact critical temperature $`T_\mathrm{C}^{(\mathrm{exact})}=2\mathrm{ln}^1(1+\sqrt{2})`$  scales with the width of the strip as follows :

$$\frac{T_\mathrm{C}^{}(L)T_\mathrm{C}^{(\mathrm{exact})}}{T_\mathrm{C}^{(\mathrm{exact})}}L^{1/\nu },$$
(22)
i.e. the ratio

$$R\frac{\frac{\mathrm{d}}{\mathrm{d}L}T_\mathrm{C}^{}(L)}{\frac{\mathrm{d}^2}{\mathrm{d}L^2}T_\mathrm{C}^{}(L)}=\frac{\nu }{\nu +1}LL.$$
(23)
The optimum width $`L^{\mathrm{opt}}`$ should be less than $`L_\mathrm{C}`$ for which the ratio of the derivatives $`R(L_\mathrm{C})`$ (23) is substantially deviated from the originally linear behavior. In our calculations we have considered the DMRG results to be incorrect for $`R=0`$ or $`\mathrm{}`$. In the case of $`R=0`$, the precise value of $`L^{\mathrm{opt}}`$ is not too important as the first derivative or change of $`T_\mathrm{C}^{}(L)`$ is very small. Near $`R=\mathrm{}`$, a sharp drop of the second derivative of $`T_\mathrm{C}^{}(L)`$ to zero is required; indeed, the change of the distance from the line $`\frac{\nu }{\nu +1}L`$ by more than one order of magnitude takes place within one step of strip-width enlargement.

In Fig. 4 plots of strip-width-dependent critical temperatures $`T_\mathrm{C}^{}(L)`$ for two different boundary conditions and various block sizes $`m`$ are given. The estimations of the exact critical temperature for periodic and open boundary conditions were found as the values of $`T_\mathrm{C}^{}(L2)`$ if the first or second derivative of $`T_\mathrm{C}^{}(L)`$ changed their signs with respect to the value in the previous step. The curves for PBC cross the exact value of $`T_\mathrm{C}^{(\mathrm{exact})}`$. The curve maxima for OBC are quite far from it, and by increasing $`L`$, $`T_\mathrm{C}^{}`$ approaches the exact value very slowly.

The accuracy of the results for periodic boundary conditions (Fig. 5) is very high already at small values of $`m`$ and exceeds by an order the critical temperature estimation for maximum computer-accessible $`m`$ when using open boundary conditions. The critical temperature for not extremely large $`m=80`$ is given to seven digits. As the width of the strip can be increased only in discrete steps and the criterion of the optimum strip-width is somewhat vague, the accuracy of the critical temperature determination should be taken as large as a single step change of $`T_\mathrm{C}^{}(L)`$. These accuracy estimations together with deviations of our results from the exact critical temperature are given in Table 6. It should be noted that only ISM was performed in calculations of $`T_\mathrm{C}^{}(L)`$ in Figs. 4 and 5. The calculations with the FSM has been also done near the $`L^{\mathrm{opt}}`$ but only slight improvements of critical temperature were obtained. In calculation of the thermal critical exponent  a similar accuracy was reached; e.g., for $`m=44`$ within FSM the critical exponent $`\nu =1.0000016`$ is very close to the exact value $`\nu =1`$ .

## V Conclusion

The DMRG method for classical spin lattice strips with periodic boundaries was developed and applied to 2D Ising model. It was shown that this approach lead to more accurate results for 2D infinite lattice than DMRG with open boundary conditions. It was demonstrated that applying finite size scaling to strips treated by DMRG, an optimal width of the strip depending on the order of approximation existed, and a prescription how to find $`T_\mathrm{C}^{}(L^{\mathrm{opt}})`$ was given. For the Ising model it was shown by computations that for these $`L^{\mathrm{opt}}(m)`$ the value of the critical temperature, was for a given $`m`$ closest to the exact one. As our approach does not involve any information about the exact critical temperature $`T_\mathrm{C}^{(\mathrm{exact})}`$ or the universality class of the model, we believe that it is applicable to many different classes of spin lattice models. This belief is supported by analogous calculation for anisotropic triangular nearest-neighbor Ising model (ATNNI) with two different antiferromagnetic interactions $`J_1`$ and $`J_2`$ (model discussed in ). For this model the transfer matrix is non-symmetric and the phase diagram is quite different from that of the standard Ising model. For the periodic boundary conditions, the plot of critical temperatures is not monotonously decreasing as in the case of the Ising model (Fig. 5) but for large $`L`$ it turns up. Nevertheless, the accuracy of the critical temperature for the exactly solvable case (for external magnetic field $`H=0`$) is similar to the presented ones in this paper.

## Acknowledgments

This work has been supported by Slovak Grant Agency, Grant n. 2/4109/98. We would like to thank the organizers of the DMRG Seminar and Workshop in Dresden for the opportunity to participate in the meetings, especially for the useful discussion with T. Nishino.


text:
# SCALE DEPENDENT DIMENSIONALITY

## 1 Introduction

Is dimensionality dependent on the scale of resolution, or is it independent of this scale. This question becomes relevant in the light of some recent work (for example Cf.). It has ofcourse been pointed out that the spin half character of a collection of Fermions leads to the usual three dimensionality of our space, while the spin half itself is associated with the Compton wavelength as discussed in recent papers (Cf. for example). Further, it was argued that as we approach the Compton scale, we encounter lower dimensionality. In this paper we point out that indeed the dimensionality is scale dependent.

## 2 Scale Dependence

We first notice that at Planck scale $`l_P`$, we have

$$N^{3/4}l_PR$$
(1)
where $`N10^{80}`$ is the number of elementary particles and $`R10^{28}cm`$ the radius of the universe. This is not an empirical relation but rather can be deduced on the basis of a fluctuational creation of particles scheme recently discussed (Cf. for example). In this scheme, $`\sqrt{N}`$ particles are fluctuationally created and this happens in the Compton time $`\tau `$ of a typical elemental particle, a pion.
Further this corresponds to a fluctuational creation of $`N^{1/4}`$ Planck particles as recently argued, in the Planck time $`\tau _P`$. Indeed, we have

$$\dot{N}N^{1/4}/\tau _P\sqrt{N}/\tau $$
(2)
Equation (2) leads to (1).
(1) shows that at the Planck length, the fluctuational dimensionality is $`4/3`$. Interestingly this is the dimension of a Koch curve and a coastline . With this dimensionality we should have

$$MR^{4/3},$$

which indeed is true.
At the Compton scale of resolution, we have, as indeed can be deduced from (2) the well known Eddington formula,

$$R\sqrt{N}l$$
(3)
(3) shows the two dimensional character at the Compton length. Indeed as noted in the introduction three dimensionality is at scales much greater than the Compton wavelength - as we approach the Compton wavelength we encounter two dimensionality as can be seen from (3) - indeed this was the key to explain puzzling characteristics of quarks including their fractional charge and handedness.
Finally at scales $`L10cm`$, we have

$$N^{1/3}LR$$
(4)
(4) shows up the usual three dimensionality.
Interestingly, if we take the typical elementary particle the pion, and consider it successively as a $`4/3`$ dimensional object at the Planck scale, a two dimensional object at the Compton scale and three dimensional at our macro scale, and consider successive densities

$$\rho _Pm/(l_P)^{4/3},\rho _\pi m/l^2\text{and}\rho m/L^3,$$

we have,

$$M\rho _PR^{4/3}\rho _\pi R^2\rho R^3,$$

as required.