text:
# Stochastic series expansion method for quantum Ising models with arbitrary interactions

## I introduction

Monte Carlo studies of classical and quantum many-body systems with long-range interactions are limited by time-consuming summations over the interacting particle pairs, the number of which grows quadratically with the system size. Many important problems in both basic and applied science can be mapped onto long-range interacting spin models, and hence it would be desirable to develop more efficient numerical techniques for tackling them. For classical Ising models, considerable progress has indeed been made on algorithms scaling almost linearly with the system size luijten . In the context of simulated annealing annealing , where the ground state of a classical system (typically with complicated interactions) is obtained through a simulation where the temperature is slowly lowered to zero, it has been suggested qann  that a more rapid convergence could be achieved by using a quantum model, e.g., the Ising model in a transverse (spin flipping) field. Even in an imaginary-time path-integral formulation, the quantum fluctuations can, at least in some cases santoro , relax the system towards its classical ground state more rapidly than thermal fluctuations. This is a strong motivation for developing more efficient simulation methods for quantum Ising models. Another important reason is the continued prominence of the transverse Ising model in the theory of magnetism, particularly in the context of quantum phase transitions sachdev ; guo ; rieger . Whereas transverse Ising models with short-range interactions have recently been actively studied using quantum Monte Carlo methods guo ; rieger , numerical work on long-range models has so far been limited to special cases rozenberg . In some of the best experimental realizations of the transverse Ising model the interactions are in fact long-ranged lihof .

Here a stochastic series expansion (SSE) sse  algorithm for transverse Ising models with long-range interactions is introduced in which the direct summation over the interacting spins is avoided. The computation time scales with the system size N as $`N\mathrm{ln}(N)`$ times the spatial integral of the absolute value of the interaction \[which normally converges as $`N\mathrm{}`$, or diverges only as $`\mathrm{ln}(N)`$\]. Both local and cluster-type updates are developed for the transverse Ising model with arbitrary interactions. The cluster update is a generalization of the classical Swendsen-Wang cluster method swendsen  to the transverse Ising model, and shares some features with a scheme previously used within the continuous-time world-line algorithm rieger . The way to treat the long-range interactions generalizes the scheme developed for the classical Ising model by Luijten and Blöte luijten ; rnote . The integration of these features in the SSE formalism should open new opportunities for detailed numerical studies of a wide range of important models. The algorithm is here tested on a ferromagnetic chain with interactions decaying as $`1/r^2`$, for which results in the classical limit are available for comparison anderson ; kosterlitz ; glumac ; cannas ; luijten2 .

In Sec. II the application of the SSE method to the transverse Ising model is described in detail. Local updates as well as classical and quantum-cluster updates are discussed. Results for the model with $`1/r^2`$ interactions are presented in Sec. III. Sec. IV concludes with a brief discussion.

## II STOCHASTIC SERIES EXPANSION

The SSE method sse  is an efficient alternative to worldline quantum Monte Carlo worldline . It is based on a generalization of the power-series scheme for the Heisenberg ferromagnet that was developed by Handscomb in the early 1960’s handscomb . Handscomb’s method was later extended to some other models hgen , but the requirement of analytically calculable traces of the terms of the expansion inhibited further progress. In the SSE method, a basis is instead chosen, and the traces are also evaluated stochastically, in combination with the sampling of the operator products in the series expansion of $`\mathrm{exp}(\beta \mathrm{H})`$. This starting point for quantum Monte Carlo is as generally applicable as the worldline (imaginary-time path-integral) approach. Recently, loop-type cluster updates evertz  have been developed and generalized for efficient SSE simulations of a wide range of models sseloop1 ; sseloop2 . However, since the loop updates rely heavily on the presence of off-diagonal pair (or multi-particle) interactions, they cannot be directly adapted to the transverse Ising model in the standard basis where the Ising term is diagonal. In the basis where the field is diagonal, loop updates can be easily implemented sseloop1 ; sseloop2  but then sign problems ssesign  appear when the interaction is frustrated. Here the SSE method is applied to an arbitrary transverse Ising model, i.e., with no limitations on the sign and range of the spin-spin interaction. Several types of local and cluster-type updates will be described.

### II.1 Configuration space

Consider the general Hamiltonian for the Ising model in a transverse field of strength $`h`$,

$$H=\underset{i,j}{}J_{ij}\sigma _i^z\sigma _j^zh\underset{i}{}\sigma _i^x,$$
(1)
where $`\sigma _i`$ is a Pauli spin matrix ($`\sigma _i^z=\pm 1`$) and $`J_{ij}`$ is the strength of the interaction between spins $`i`$ and $`j`$, which can be random or uniform and of any sign. The dimensionality is arbitrary. Define the operators

$`H_{0,0}`$ $`=`$ $`1,`$ (2)
$`H_{i,0}`$ $`=`$ $`h(\sigma _i^++\sigma _i^{}),i>0,`$ (3)
$`H_{i,i}`$ $`=`$ $`h,i>0,`$ (4)
$`H_{i,j}`$ $`=`$ $`|J_{ij}|J_{ij}\sigma _i^z\sigma _j^z,i,j>0,ij.`$ (5)
Up to a constant, the Hamiltonian can be written as

$$H=\underset{i=1}{\overset{N}{}}\underset{j=0}{\overset{N}{}}H_{i,j}.$$
(6)
The constants $`H_{i,i}`$ are introduced for purposes that will become clear below. Note that $`H_{0,0}`$ is not included as a term in the Hamiltonian (6) but will be important in the simulation scheme.

In the SSE approach sse  to finite-temperature quantum Monte Carlo, the partition function $`Z=\mathrm{Tr}\{\mathrm{exp}(\beta H)\}`$ is written as a power-series expansion, with the trace expressed as a sum over diagonal matrix elements in a suitably chosen basis. Using (6) then gives

$$Z=\underset{\alpha }{}\underset{n=0}{\overset{\mathrm{}}{}}\underset{S_n}{}\frac{\beta ^n}{n!}\alpha \left|\underset{l=1}{\overset{n}{}}H_{i(l),j(l)}\right|\alpha ,$$
(7)
where $`S_n`$ denotes a sequence of $`n`$ operator-index pairs (hereafter referred to as operators):

$$S_n=[i(1),j(1)],\mathrm{},[i(n),j(n)],$$
(8)
with $`i(l)\{1,\mathrm{},N\}`$ and $`j(l)\{0,\mathrm{},N\}`$. The standard basis $`\{|\alpha \}=\{|\sigma _1^z,\mathrm{},\sigma _N^z\}`$ is used.

Because of the constants added to $`H_{i,j}`$ in (5), the eigenvalues of these operators are $`2|J_{ij}|`$ and $`0`$. All non-zero terms in (7) are therefore positive and can be used as relative probabilities in an importance sampling scheme. A term is specified by a state $`|\alpha `$ and an operator sequence $`S_n`$. One can show that the total internal energy (including the constants added to $`H`$) is given by sse ; handscomb  $`E=n/\beta `$. Hence, the size of the operator sequence to be stored in computer memory scales as $`\beta NI_N(J)`$, where

$$I_N(J)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{N}{}}|J_{ij}|,$$
(9)
which converges or grows much slower than $`N`$ for most cases of interest.

In order to construct an efficient sampling scheme, it is useful to cut the expansion (7) at some power $`n=L`$, sufficiently high for the remaining truncation error to be exponentially small and completely negligible \[$`L`$ clearly has to be $`\beta NI_N(J)`$\]. One can then obtain an expansion for which the length of the operator sequence is constant, by considering random insertions of $`Ln`$ unit operators $`H_{0,0}`$ in the product in (7). Adjusting for the $`\left(\genfrac{}{}{0pt}{}{L}{n}\right)`$ possible insertions gives

$$Z=\frac{1}{L!}\underset{\alpha }{}\underset{S_L}{}\beta ^n(Ln)!\alpha \left|\underset{l=1}{\overset{L}{}}H_{i(l),j(l)}\right|\alpha ,$$
(10)
where $`[i(l),j(l)]=[0,0]`$ is now also an allowed operator in the sequence $`S_L`$, and $`n`$ denotes the number of non-$`[0,0]`$ operators. Note again that $`H_{0,0}`$ is not part of the Hamiltonian, but is introduced only for the purpose of constructing a computationally simpler updating scheme where the operator list has a fixed length.

It is useful to define states $`|\alpha (p)=|\sigma _1^z(p),\mathrm{},\sigma _N^z(p)`$ obtained by propagating $`|\alpha =|\alpha (0)`$ by the first $`p`$ operators in $`S_L`$;

$$|\alpha (p)=r\underset{l=1}{\overset{p}{}}H_{i(l),j(l)}|\alpha ,$$
(11)
where $`r`$ is a normalization factor. A non-vanishing matrix element in (10) then corresponds to the periodicity condition $`|\alpha (L)=|\alpha (0)`$, which requires that for each site $`i`$ there is an even number (or zero) of spin flipping operators $`[i,0]`$ in $`S_L`$. Definition (5) implies that the Ising operators $`[i,j]`$ may act only on states with $`\sigma _i^z=\sigma _j^z`$ if $`J_{ij}<0`$ (ferromagnetic), or $`\sigma _i^z=\sigma _j^z`$ if $`J_{ij}>0`$ (antiferromagnetic). There are no other constraints.

An SSE configuration is illustrated in Fig. 1. The vertical direction in this representation will be referred to as the SSE propagation direction. It can be related to the imaginary-time direction in standard path integral representations irsse . Note that this full configuration, including all the states $`|\alpha (p)`$ explicitly, does not have to be stored in the simulation. A single state and the operator sequence suffice for reproducing all the states, and such a representation is used in some stages of the simulation. For some updates it is convenient to generate other representations, as will be discussed below.

### II.2 Local updates

The sampling of Eq. (10) can be carried out using simple operator substitutions of the types

$`[0,0]_p[i,j]_p,i,j0,`$ (12)
$`[i,i]_{p_1}[i,i]_{p_2}[i,0]_{p_1}[i,0]_{p_2},i0,`$ (13)
where the subscript $`p`$ indicates the position ($`p=1,\mathrm{},L`$) of the operator in the sequence $`S_L`$. The power $`n`$ is changed by $`\pm 1`$ in the diagonal update (12) and is unchanged in the off-diagonal update (13). In the diagonal update the Ising terms $`[i,j]`$ and the constants $`[i,i]`$ are sampled. The constants are used in the off-diagonal update as a means of achieving easy insertions and removals of two spin-flipping operators $`[i,0]`$. With the value $`h`$ chosen for the constant in (4), the operator replacements do not change the weight of the SSE configuration. However, the off-diagonal update also leads to spin flips in the propagated states between $`p_1`$ and $`p_2`$; $`\sigma _i^z(p_1),\mathrm{},\sigma _i^z(p_21)`$ $``$ $`\sigma _i^z(p_1),\mathrm{},\sigma _i^z(p_21)`$. \[$`p_1>p_2`$ also has to be considered, leading to flipped $`\sigma _i^z(p_1),\mathrm{},\sigma _i^z(L1)`$ $`\sigma _i^z(0),\mathrm{},\sigma _i^z(p_21)`$\], which is allowed if (and only if) no Ising operators acting on site $`i`$ are present in $`S_L`$ between positions $`p_1`$ and $`p_2`$. Note that this constraint is completely local, regardless of the range of the interaction, and that the update requires no knowledge of the spin state. This is the reason for the advantage of this simulation scheme over worldline methods worldline ; rieger , where calculating the acceptance probability for every update requires a summation over all the spins interacting with those flipped. Here an allowed off-diagonal update (13) leaves the weight unchanged and can be carried out with probability one.

If $`h0`$, the above updates of the operator sequence suffice for achieving ergodicity. If there are no Ising operators acting on a site $`i`$, $`\sigma _i^z(0),\mathrm{},\sigma _i^z(L1)`$ can also be flipped without changes in $`S_L`$. This update in principle makes simulations using the present scheme possible also for $`h=0`$, but in practice unconstrained spins occur frequently only at high temperatures, when $`n`$ is small. Other types of “classical” spin flips — flips of clusters — are also possible, and will be discussed in Sec. II C.

The simulation can be started with a random state $`|\alpha (0)`$ and a sequence $`S_L`$ containing only $`[0,0]`$ operators. The truncation $`L`$ can be chosen arbitrarily (small); it is adjusted during the equilibration part of the simulation, e.g., by requiring $`L>(4/3)n`$ after each update. This ensures than $`n`$ never reaches $`L`$ during the reminder of the simulation, and hence that there will be no detectable systematic errors arising from the truncation of the expansion sse . In the beginning of an updating cycle, the operator sequence $`S_L`$ and the state $`|\alpha (0)`$ is stored.

The diagonal update (12) is attempted successively for all $`p=1,\mathrm{},L`$. In the course of this process, the spin state is propagated by flipping spins $`\sigma _i^z`$ as off-diagonal operators $`[i,0]`$ are encountered in $`S_L`$, so that the states $`|\alpha (p)`$ are generated successively. For an $`[i,j][0,0]`$ update, i.e., removing a Hamiltonian operator, there are no constraints and the update should always be accepted with some non-zero probability. In the case of $`[0,0][i,j]`$, i.e., inserting an operator from the Hamiltonian, there are constraints, and the update may not be allowed for all $`i,j`$. However, initially the indices $`i,j`$ are left undetermined and it is assumed that any $`[i,j]`$ would be allowed. Under this assumption, the acceptance probabilities for the diagonal update are given by

$`P([0,0][i,j])`$ $`=`$ $`{\displaystyle \frac{\beta (Nh+2\underset{ij}{}|J_{ij}|)}{Ln}},`$ (14)
$`P([i,j][0,0])`$ $`=`$ $`{\displaystyle \frac{Ln+1}{\beta (Nh+2_{ij}|J_{ij}|)}},`$ (15)
where $`_{ij}`$ does not include $`i=j`$ and $`P>1`$ should be interpreted as probability one, as usual. These probabilities are simply obtained from the ratio of the new and old prefactor in (10) when $`nn\pm 1`$;

$$\beta ^{\pm 1}\frac{[L(n\pm 1)]!}{(Ln)!},$$
(16)
and the ratio between the matrix element $`1`$ of the $`[0,0]`$ operator and the sum $`Nh+2_{ij}|J_{ij}|`$ of the non-zero matrix elements of all $`[i,j]`$ operators. Staying with the assumption that any $`[i,j]`$ is allowed in the update $`[0,0][i,j]`$, the relative probability of an operator with the first index $`i`$ is $`P(i)=_jM_{ij}`$, where $`M_{ij}`$ is the non-zero matrix element corresponding to $`H_{ij}`$ (i.e., $`h`$ for $`i=j`$ and and $`2|J_{ij}|`$ else). The normalized cumulative probabilities $`P_c(k=1,\mathrm{},N)`$ are stored in a pre-generated table;

$$P_c(k)=\frac{\underset{i=1}{\overset{k}{}}P(i)}{_{i=1}^NP(i)}.$$
(17)
In order to select the first index $`i`$ of the operator $`[i,j]`$ to be inserted, a random number $`0R<1`$ is generated. The table $`P_c`$ is searched (using, e.g., a simple binary search) for the smallest $`k`$ for which $`P(k)R`$; the first index of the operator $`[i,j]`$ is then $`i=k`$. The second index can be chosen in a completely analogous way, with the relative probability for $`j`$, given $`i`$, being $`M_{ij}`$. For a random system with long-range interactions, a pregenerated table with $`N^2`$ elements is hence needed for storing all the cumulative probabilities for the second index. For non-random interactions in a translationally invariant system, the first index can be selected at random with equal probabilities without searching a table, and the size of the second table is reduced to $`N`$. For a short-range or truncated interaction the table size is smaller, corresponding to the number of spins within the range of the interaction; clearly, the whole selection process should then be reduced to a single step for obtaining both $`i`$ and $`j`$ (e.g., selecting one out of a total number $`N`$ of operators and reading the corresponding $`i,j`$ from a table). The two-step procedure is advantageous for non-random long-range interactions, where it allows for the reduction of the size of the probability table from $`N^2`$ to $`N`$. For random models, the storage requirement is always $`N^2`$, and it may then again be better to combine the first and second index searches, using a single size-$`N^2`$ table for all the cumulative probabilities of $`[i,j]`$. For short-range random interactions the size of the table is $`N`$ times the number of spins within the interaction range.

The operator $`[i,j]`$ generated as above may or may not be allowed in the current spin configuration $`|\alpha (p)`$. If $`\sigma _i^z(p)`$ and $`\sigma _j^z(p)`$ indeed are in an allowed state, $`[i,j]`$ is inserted at position $`p`$. Otherwise, the process for generating $`[i,j]`$ has to be repeated, until an allowed operator has been generated. The reject-and-repeat step leads to the correct probabilities for selecting among all the allowed diagonal operators $`[i,j]`$. Typically, an allowed operator is generated very quickly, since the interactions favor the allowed spin alignment. Note that the constants $`[i,i]`$ are always allowed (for $`h>0`$), so there is no risk of the search never terminating.

The off-diagonal update (13) can be efficiently carried out if $`S_L`$ is first partitioned into separate subsequences for each site $`i`$. Subsequence $`i`$ contains only spin-flipping operators $`[i,0]`$ and constants $`[i,i]`$. Their positions in $`S_L`$ are also stored, to be used for recombining the subsequences after the update. The constraints on modifications at site $`i`$ imposed by Ising operators $`[i,j]`$ or $`[j,i]`$ (for any $`j`$) can be stored as flags indicating the presence of one or several of these operators between neighboring subsequence operators. Updating a subsequence amounts to selecting two non-constrained neighboring operators at random from the subsequence, and carrying out the substitution (13) if the two operators are identical. If they are different, they can be permuted. A number proportional to the subsequence length of such pair updates are carried out for each subsequence, after which they are recombined into a new $`S_L`$.

The diagonal update (12) at all positions in $`S_L`$ require $`L\mathrm{ln}(N)\beta N\mathrm{ln}(N)I_N(J)`$ operations, where the factor $`\mathrm{ln}(N)`$ is the scaling of the average number of operations needed to search the cumulative probability table(s) in the case of long-range interactions. Partitioning $`S_L`$ into subsequences and updating all of them according to (13) requires on the order of $`L`$ operations. Hence, the number of operations for a full updating cycle of the degrees of freedom of the system (one Monte Carlo step) scales as $`\beta N\mathrm{ln}(N)I_N(J)`$. This should be compared to the $`\beta N^2`$ scaling in worldline methods worldline ; rieger , where one power of $`N`$ is due to the summation required to calculate the weight change when flipping a spin interacting with $`N`$ other spins. Here this summation has been circumvented by writing the interactions in the SSE formalism as fluctuating constraints that are purely local.

### II.3 Classical cluster update

In the Swendsen-Wang cluster algorithm swendsen  for the classical Ising model, i.e., with $`h=0`$ and a uniform nearest-neighbor interaction of strength $`J`$, auxiliary bond variables $`b_{ij}`$ are introduced in order to construct clusters of spins that can be flipped independently of each other. Given a spin configuration, and with initially all bond variables $`b_{ij}=0`$, for every interacting spin pair for which $`\sigma _i\sigma _j=J/|J|`$ (i.e., the orientation energetically favored) the bond variable is set, $`b_{ij}=1`$, with probability $`P=1\mathrm{e}^{2|J|\beta }`$. When all bonds have been visited, clusters of spins connected by $`b_{ij}=1`$ bonds are formed, and each of these clusters is flipped with probability $`1/2`$. Single spins not connected to any $`b_{ij}=1`$ bond are single-spin clusters. After the clusters have been flipped, all the bond variables are again set to zero and the process is repeated. This scheme can in fact be constructed using the SSE formalism, as an alternative to the Fortuin-Kastelin mapping ft , on which the Swendsen-Wang algorithm is based.

The relation to the Swendsen-Wang algorithm is shown as follows, by applying the SSE method to the classical Ising model, now again considering a general form of the interaction $`J_{ij}`$ and with the bond-operator $`H_{ij}=|J_{ij}|J_{ij}\sigma _i^z\sigma _j^z`$ as in Eq. (5). Since all operators $`H_{ij}`$ commute, the operator $`\mathrm{e}^{\beta H}`$ can be written as a product of operators $`\mathrm{e}^{\beta H_{ij}}=1+\beta H_{ij}+\mathrm{}`$. The uniqueness of the power-series expansion then implies that in the SSE, where $`\mathrm{e}^{\beta H}`$ is expanded directly, the probability of having one or more operators $`H_{ij}`$ on a bond $`i,j`$ when $`\sigma _i\sigma _j=J_{ij}/|J_{ij}|`$ is $`1\mathrm{e}^{2|J_{ij}|\beta }`$, i.e., exactly the probability of having the bond variable $`b_{ij}=1`$ in the Swendsen-Wang scheme. In a configuration $`\sigma _i\sigma _j=J_{ij}/|J_{ij}|`$ there can be no operators on the bond in the SSE, and the Swendsen-Wang $`b_{ij}=1`$ probability is also zero per construction. One can hence make the connection that one or more operators acting on a spin pair in the SSE scheme corresponds to a filled bond ($`b_{ij}=1`$) in the Swendsen-Wang algorithm. The definition of a cluster is then exactly the same in the two algorithms. Clearly, such a cluster in the SSE can also always be flipped, since the Ising operators only impose constraints on the relative orientations of connected spins, which is maintained when the cluster is flipped. Since the weight does not change, the flip should be done with probability $`1/2`$. The scheme is hence identical to the Swendsen-Wang algorithm, except that the filled bonds $`b_{ij}=1`$ in SSE are generated in a different way, using the diagonal update (12). Note that for a classical model, all the propagated SSE states (11) are identical, i.e., $`\sigma _i^z(p)=\sigma _i^z(0)`$ for all $`p=0,\mathrm{},L1`$, and hence no state propagations have to be considered as the diagonal update is carried out.

It is interesting to note that the SSE scheme for the classical Ising model should in fact be more efficient than the standard Swendsen-Wang algorithm at high temperatures. This is because the number of operators in the SSE operator list scales as $`E(T)/T`$, where $`E(T)`$ is the total energy at temperature $`T`$ ($`EN`$) and for large $`T`$ the construction of the clusters based on the operator list should then be faster than visiting all the bonds, as is done in the Swendesn-Wang algorithm. However, in practice the interesting physics occurs when the number of SSE operators per interacting spin pair is of the order of one or larger, and then there are no advantages of the SSE classical cluster algorithm relative to Swendesen-Wang.

The classical SSE cluster update can also be used in the presence of the transverse field ($`h>0`$). The clusters are defined in terms of bonds signifying the presence of one or more Ising operator, as above, without regard for the single-spin flipping operators $`H_{i,0}`$ and constants $`H_{i,i}`$. These operators can be neglected because when a cluster is flipped, all spins $`\sigma _i^z`$ belonging to the cluster are implicitly flipped in all propagated states (11), i.e., $`\sigma _i^z(p)`$ $``$ $`\sigma _i^z(p)`$ for all $`p=0,\mathrm{},L1`$ (this is the reason for the term “classical cluster” even when $`h>0`$) and hence all operations with the single-spin operators remain valid and produce the same factors in the weight before and after the cluster flips. Note again that only the first state, i.e., $`\sigma _i^z(0)`$, $`i=1,\mathrm{},N`$, has to be stored when constructing the classical clusters.

In the case of long-range interactions, a cluster can consist of several intertwined pieces on the lattice, as illustrated for a two-dimensional case in Fig. 2. Regardless of the range of the interaction, the construction of the clusters, given an SSE operator list, can be easily carried out using a number of operations scaling as the number of operators in the list.

Since the classical SSE cluster update is equivalent to the Swendsen-Wang algorithm in the classical limit and only takes the Ising terms into account also in the quantum case, it cannot be expected to be efficient much beyond the classical limit $`h=0`$. For a non-random system that undergoes a phase transition at $`T_c(0)`$ when $`h=0`$, the critical temperature is reduced by the transverse field; $`T_c(h)<T_c(0)`$. Hence, the classical clusters will percolate for $`T>T_c`$ and this update will not be efficient close to $`T_c`$. The primary reason to introduce the classical cluster update here was to demonstrate the relationship between SSE and the Swendsen-Wang algorithm. In the case of long-range interactions, the scheme becomes very similar to the Luijten-Blöte algorithm luijten , again just differing in the way the bonds are generated.

### II.4 Quantum-cluster update

The purpose of the quantum-cluster update is to effect flips of spins $`\sigma _i^z(p)`$ only in a limited number of propagated states $`p`$, in different states for different sites $`i`$. In other words, these clusters will be finite and irregularly-shaped both in the space and SSE propagation (imaginary time) direction. In the process, operator substitutions $`H_{i,i}H_{i,0}`$ (constant to spin-flip, and vice versa) will also be accomplished. This update hence replaces the local off-diagonal update (13).

To discuss the quantum-cluster update, it is useful to introduce the notion of vertices sseloop1 ; sseloop2 . Looking at the graphical representation of a configuration in Fig. 1, it can be noted that the vertical “lines” of same spins between two operators acting on a given site constitute redundant information. The full configuration can be represented by a list of positions (on the lattice) of the operators, and the spin states (on one or two sites for the model considered here) before and after the operators act. These relevant spins are called legs of the 2-spin vertices (corresponding to constant and spin-flip operators) or 4-spin vertices (corresponding to Ising bond-operators). All possible vertices for the transverse Ising model are shown in Fig. 3. Note that only those Ising vertices that are compatible with the sign of the interaction between a given pair of spins are allowed for those spins; again, this is due to the choice of constant in the bond-operator (5). In the computer, the vertices are linked to each other by pointers, so that from a given vertex-leg one can reach the next or previous vertex that has a leg on the same site (i.e., there are links that replace the segments of vertical lines of same spins in Fig. 1). A detailed discussion of the practical implementation of a linked vertex list has been given in Ref. sseloop2, .

To construct and flip a quantum-cluster, one of the legs of one of the $`n`$ vertices is picked at random, and the corresponding spin is flipped. Depending on the type of the vertex, different actions are taken, examples of which are given in Fig. 4. The arrow pointing into the vertex indicates the entrance leg. In the case of an Ising vertex, all the four spins are flipped and the cluster building process branches out from all the legs, as indicated by the arrows pointing out from the vertex. Using the pointers of the linked vertex list, the arrows point to legs of other vertices; these become new entrance legs which are put on a stack and subsequently processed one-by one. If the entrance leg is on a constant or spin-flip vertex, only the entrance spin is flipped. The vertex type then also changes, in terms of operators from $`H_{i,0}`$ to $`H_{i,i}`$, and vice versa. In these cases there is no branching-out and no new legs are put on the stack, i.e., this particular branch of the cluster terminates. If a link points to a spin that has already been flipped (i.e., two arrows point toward each other), that leg should not be used again as an entrance and is hence not put on the stack. Therefore, each vertex-leg can be visited at most once (each spin can be flipped at most once) and the cluster is completed when there are no more entrance-legs on the stack. The reason that the cluster can always be flipped is again that the SSE weight is not affected; the matrix element of the Ising bond-operator is not affected when both spins are flipped (in the absence of an external field in the $`z`$-direction, which would necessitate a modified approach), and the matrix elements for the constant and spin-flip operators are both equal to $`h`$.

The construction of a single cluster, which is flipped with probability one, is a quantum-mechanical analogue of the classical Wolff algorithm wolff ; in the absence of the transverse field the clusters are identical to those of the Wolff algorithm. Note, however, that there is a difference when constructing more than one cluster: The number of operators in the SSE operator list, and their positions on the lattice, do not change in the quantum-cluster update. The clusters are therefore completely deterministic once the operator list is given. Hence, when constructing several clusters using the same SSE operator list, it is quite likely that the same cluster is constructed and flipped multiple times. This is clearly not efficient. However, one can also construct all clusters, as in the Swendsen-Wang scheme, and only flip them with probability $`1/2`$. This is done by always starting a new cluster from a vertex-leg which has not yet been visited. Every vertex-leg belongs uniquely to one cluster, and clearly the number of operations required to complete this updates then scales as $`L`$, i.e., typically as $`\beta N`$.

A natural definition of a Monte Carlo step including the quantum-cluster update is a full sweep of diagonal updates, followed by the construction of the linked list of vertices, in which all clusters are constructed and flipped with probability $`1/2`$. After that, the updated vertex list is mapped back into a state $`|\alpha (0)`$ and an operator sequence $`S_L`$. Free spins, i.e., those that are not acted on by any operators, can again be considered as single-spin clusters and should also be flipped with probability $`1/2`$. No local off-diagonal updates (13) are needed.

Since the quantum-cluster update explicitly includes the quantum mechanical features of the configurations (i.e., the presence of spin-flip operators), it can be expected to work well also close to a quantum phase transition ($`T_c=0`$) driven by varying $`h`$.

## III 1D INVERSE-SQUARE FERROMAGNET

As a non-trivial demonstration of the method, a ferromagnetic chain with interactions decaying as $`1/r^2`$ is considered next. The interaction is summed over all $`i,j`$ in (1), i.e., each pair is counted twice. Periodic boundary conditions are used. $`J_{ij}`$ includes both distances in the periodic system, i.e.,

$$J_{ij}=J_{ji}=\frac{J}{2}\left(\frac{1}{|ij|^2}+\frac{1}{(N|ij|)^2}\right),$$
(18)
where $`J`$ sets the over-all energy scale.

The classical $`1/r^2`$ Ising chain has been the subject of numerous studies anderson ; kosterlitz ; glumac ; cannas ; luijten2 . The long-range interaction allows for a finite-$`T`$ phase transition even in one dimension. The transition is of an unusual kind, with the correlation length exponent $`\nu =\mathrm{}`$, and a discontinuous jump in the magnetization at $`T_c`$. It can be thought of as a one-dimensional analogue of the Kosterlitz-Thouless transition, with the topological excitations being kink solitons kosterlitz . The model is also important because it can be mapped onto the Kondo problem anderson .

For small $`h/J`$, one can expect a behavior similar to the classical case, i.e., a finite-$`T`$ phase transition to a ferromagnetic state. For $`h\mathrm{}`$ the system becomes disordered, and there should therefore be a finite $`h_\mathrm{c}`$ for which the system undergoes a quantum phase transition (i.e., $`T_\mathrm{c}=0`$). For $`h<h_c`$, $`T_c>0`$ and one can expect the same universality class as in the classical case, since the quantum fluctuations become irrelevant at $`T_\mathrm{c}`$. Here only a single field-strength $`h/J=0.5`$ is considered. The simulations show that $`T_\mathrm{c}>0`$ in this case.

The model is invariant with respect to flipping all spins, which means that for any finite system the average magnetization vanishes. The squared magnetization,

$$M^2=\left(\frac{1}{N}\underset{i}{}\sigma _i^z\right)^2,$$
(19)
is therefore calculated. Results for $`M^2`$ with statistical errors in the fifth decimal place can easily be obtained for systems with several hundred spins (and there are no problems in going to considerably larger systems). For small systems the results are in perfect agreement with exact diagonalization data.

A “tempering” scheme, where $`\beta `$ is considered as an additional discretized dimension of the configuration space tempering , was also implemented in the simulations. Transitions satisfying detailed balance are carried out between neighboring $`\beta `$ values. This way, results can be obtained on a dense temperature grid with much less effort than by several fixed-$`\beta `$ simulations. A temperature spacing $`\mathrm{\Delta }T/J=0.010.02`$ was used.

Figure 5(a) shows results for systems with $`N`$ up to $`512`$. At high temperatures, $`M^2`$ decreases with increasing $`N`$, as expected, and there is a slight increase with $`N`$ at low $`T`$. The curves intersect at $`T/J1.4`$. A discontinuous magnetization jump at $`T_c`$ in the thermodynamic limit implies that $`M^2`$ should become size independent at $`T_\mathrm{c}`$ for sufficiently large $`N`$. A notable difference between the finite-size behavior of $`M(T)`$ seen in Fig. 5 and the magnetization curves for the classical system is that in the latter case the curves do not intersect, but the infinite-size value $`M(T_c)`$ is approached with a logarithmic correction luijten2 . The reason for the different form of the finite-size scaling for $`h>0`$ should be clarified.

Figure 5(b) shows in more detail the behavior in the region where the curves intersect. The point of intersection moves slowly towards higher $`T`$ as $`N`$ increases, and larger $`N`$ would be needed to extract $`T_c`$ accurately. Based on the data presented here $`T_\mathrm{c}/J=1.42\pm 0.01`$. This can be compared with $`T_\mathrm{c}(h=0)1.53J`$ for the classical model luijten2 . A reduction of $`T_\mathrm{c}`$ is expected on account of quantum fluctuations for $`h>0`$. The quite small reduction for $`h/J=0.5`$ is consistent with the $`T0`$ magnetization being only slightly reduced from the classical value $`M(0)=1`$. It would clearly also be interesting to study the quantum phase transition, but that problem is beyond the scope of this paper.

The high accuracy of these simulations demonstrate that the algorithm indeed is very efficient. The computer resources used for this work were quite modest; on the order of $`200`$ CPU hours on an SGI Origin2000. The scaling of the CPU time is close to linear in $`N`$ for the $`1/r^2`$ interaction, for which the interaction sum (9) converges rapidly. Only the local updates discussed in Sec. II B were used in these simulations. The cluster updates have been tested as well and improve the performance . The quantum-cluster update should be particularly useful for studying the quantum phase transition, where there will be a broad distribution of the sizes of the clusters constructed in this update.

## IV DISCUSSION

A new approach to long-range interacting quantum models has here been developed within the framework of transverse Ising models. It is important to note that the technique can also be generalized to other types of systems, with the usual caveat of sign problems ssesign . What is particular about the Ising interaction is that it can be written so that a spin-spin term either gives zero or a constant when acting on an arbitrary basis state. This is what is needed in order to reduce the interactions to local constraints in the SSE formalism. However, the algorithm can easily be modified to cases where the diagonal interaction can take several non-zero values. The first modification is in the diagonal update. For the Ising model, the probability of selecting a given bond (5) is given by a matrix element corresponding to the spin pair being in a configuration energetically favored by the interaction. If the spins are in a non-favored configuration (correspinding here to a vanishing matrix element) the update is simply rejected. In the general case, the probability to use in this update should correspond to the largest diagonal matrix element on a given bond, and if the actual configuration corresponds to a smaller matrix element the update should be accepted only with a probability reflecting this smaller value (i.e., the ratio between the actual value and the largest matrix element). The quantum-cluster update can be modified by using ideas developed within the “directed-loop” algorithm sseloop2 . For example, there could be $`4`$-particle vertex processes where the whole vertex is not necessarily reversed as in Fig. 4(a). The process could instead either go straight through the vertex (modifying the vertex only at the entrance and exit legs) or “bounce” back without modifying the vertex at all. The details of how this is done in practice will of course depend on the types of diagonal and off-diagonal terms in the Hamiltonian. The main point to note is that in the SSE approach all the information needed to update the vertices is contained in the vertices themselves, which are always local and can be generated in the diagonal update based purely on local decisions.

The transverse Ising simulation algorithm has here been tested on a one-dimensional model with long-range interactions decaying as $`1/r^2`$. The program requires almost no modifications for higher-dimensional systems, and random interactions are also very easy to implement. Future studies will have to address how well the method works in practice for a variety of systems that are more challenging because of frustrated interactions, long-range frustrated interactions, or even randomly frustrated long-range interactions. For short-range interactions, it would also be interesting to see how the SSE quantum-cluster method constructed here compares to the transverse Ising cluster method previously developed for continuous-time worldline simulations rieger .

###### Acknowledgements.

I would like to thank Patrik Henelius for numerous discussions. This work was supported by the Academy of Finland (project 26175).


format: text

text:
# Fractionalization in Spin Liquid Mott Insulators: Vison Wavefunctions and Gaps

## Abstract

We determine conditions under which a spin-liquid Mott insulator $`|0`$ defined by a Gutzwiller projected BCS state at half-filling is $`Z_2`$ fractionalized. We construct a trial vison ($`Z_2`$ vortex) state $`|V`$ by projecting an $`hc/2e`$ vortex threading the hole of a cylinder/torus and examine its overlap with $`|0`$ using analytical and numerical calculations. We find that generically the overlap vanishes in the thermodynamic limit so the spin-liquid is $`Z_2`$ fractionalized. However, for microscopic parameters appropriate for high Tc cuprates, we estimate the vison gap $`E_vJ`$, consistent with recent experimental bounds, due to the proximity to the bipartite symmetric point where $`E_v=0`$.

Spin liquid insulators doped with holes have long been regarded as natural candidates for high temperature superconductivity since the resonating valence bond (RVB) proposal of Anderson anderson . There has been a revival of interest in this idea since such wavefunctions have recently been shown to provide a remarkably good description paramekanti  of the superconducting (SC) state of the high Tc cuprates over a wide doping range. In addition, spin liquid states for frustrated magnets are of great interest both theoretically spin-liquids  and experimentally RVBexpts .

It has also long been suspected that RVB insulators may exhibit quantum number fractionalization or spin-charge separation kivelson ; read ; sachdev ; wen , i.e., support neutral, spin-$`1/2`$ and charged, spin-$`0`$ excitations. Senthil and Fisher senthil\_1  have recently emphasized that such fractionalized states in $`d=2,3`$ may be described as deconfined phases of $`Z_2`$ gauge theories, and necessarily imply a novel gapped excitation read ; wen ; senthil\_1 : the $`Z_2`$ vortex or “vison”. The presence/absence of visons threading noncontractible loops leads to a topological ground state degeneracy on cylinders or torii read ; wen ; senthil\_1 . They further proposed senthil\_2  an ingenious flux-trapping experiment to detect the vison. However, a series of experiments to detect visons in highly underdoped cuprates have led to negative results, setting bounds on the vison gap $`E_v150180`$K moler , an order-of-magnitude smaller than the natural scale $`J10^3`$K.

Motivated by the above developments we address the following issues: (a) How can one think about the vison in terms of electronic coordinates? (b) Under what conditions is the vison well-defined, leading to $`Z_2`$-fractionalization and topological degeneracy in spin liquid insulators? (c) How large is the vison gap, the energy scale below which fractionalization is apparent? (d) What are the implications for cuprates and for other frustrated magnetic materials?

Vison Wavefunction: For conceptual clarity, consider first a square lattice of $`L^2`$ sites “wrapped” into a cylinder along the $`\widehat{x}`$ direction. The spin-liquid ground state $`|0`$ of interest to us is given by the Gutzwiller projection $`𝒫=_𝐫(1n_𝐫n_𝐫)`$ of the $`N`$-particle d-wave BCS wavefunction at half-filling ($`N=L^2`$):

$$|0𝒫|BCS=𝒫\left[\underset{𝐫,𝐫^{}}{}\phi ^{}(𝐫𝐫^{})c_𝐫^{}c_𝐫^{}^{}\right]^{N/2}|\mathrm{vac}.$$
(1)
Here $`\phi (\rho )=L^2_𝐤\mathrm{exp}(i𝐤\rho )[\mathrm{\Delta }_𝐤/(\xi _𝐤+\sqrt{\xi _𝐤^2+\mathrm{\Delta }_𝐤^2})]`$ is the internal pair wavefunction with the variational parameters notation  $`\mu `$ and $`\mathrm{\Delta }`$ which determine $`\xi _𝐤=ϵ(𝐤)\mu `$, with $`ϵ(𝐤)=2t\left(\mathrm{cos}k_x+\mathrm{cos}k_y\right)4t^{}\mathrm{cos}k_x\mathrm{cos}k_y`$ and $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }\left(\mathrm{cos}k_x\mathrm{cos}k_y\right)/2`$. The $`𝐤`$’s are chosen consistent with periodic boundary condition (PBC) along $`\widehat{x}`$.

With optimized variational parameters, $`𝒫|BCS`$ is an energetically good variational ground state for the large-$`U`$ Hubbard model gros ; paramekanti . However, we focus more generally on understanding topological degeneracy and fractionalization in the phases described by this class of states rather than just the energetically optimal wavefunction. Toward this end, we define a trial vison state $`|V_y`$ by Gutzwiller projecting an $`hc/2e`$ vortex threading the cylinder along $`\widehat{y}`$ by $`|V_y𝒫|\left(hc/2e\right)_y`$, following Senthil and Ivanov ivanov  and building on earlier ideas footnote1 .

To construct the $`hc/2e`$ vortex we simply modify the pair $`\phi ^{}(𝐫𝐫^{})`$ in the BCS ground state to $`\mathrm{exp}[i𝐐(𝐫+𝐫^{})/2]\phi _A(𝐫𝐫^{})`$, where $`𝐐=2\pi \widehat{x}/L`$ twists the phase of the condensate by $`2\pi `$ going around the cylinder. However, the requirement of a single-valued wavefunction, when $`𝐫`$ (or $`𝐫^{}`$) goes around the cylinder, constrains $`𝐤`$’s to satisfy antiperiodic boundary conditions (APBC) along $`\widehat{x}`$ in the Fourier transform defining $`\phi _A(𝐫𝐫^{})`$; hence the subscript $`A`$. For the fully projected state at half-filling one can show that the phase factors drop out; in first quantized form, these are independent of the spin configuration. Then the trial vison state $`|V_y`$ simplifies to a form identical to eq. (1) with $`\phi ^{}\phi _A`$.

Note that a projected $`hc/e`$ vortex with $`𝐐=4\pi \widehat{x}/L`$ does not require any change in boundary conditions on $`\phi `$ and is therefore identical to the ground state $`|0`$. The manifest $`Z_2`$ character (two visons is the same as none) of projected vortices makes them attractive candidates for visons. On a torus $`|0`$ corresponds to PBC along both $`\widehat{x}`$ and $`\widehat{y}`$, the vison states $`|V_x`$ and $`|V_y`$ correspond to PBC/APBC and $`|V_{xy}`$ to APBC along both directions.

Phase diagram and Symmetries: A superconducting state with a vortex threading the hole of the cylinder/torus carries current and is orthogonal to the ground state. However, Gutzwiller projecting the BCS wavefunction at half-filling destroys SC order paramekanti  and results in a spin-liquid insulator footnote2 . The question we would like to address is: does any remnance of vorticity survive projection? If it does, then $`V_\alpha |0=V_\alpha |V_\beta =0`$, for $`\alpha ,\beta =x,y,xy`$, the vison is well-defined leading to topological degeneracy footnote3 , and one is in a $`Z_2`$ fractionalized phase. If, on the other hand, the proposed vison state has nonzero overlap with the state $`|0`$ in the thermodynamic limit, one is in a conventional non-fractionalized phase. Our main task then is to compute these overlaps and establish the existence of disconnected topological sectors as evidence for fractionalization.

To determine the conditions under which the vison survives, it is useful to consider the full parameter space for projected BCS states at half-filling, which is the $`(t^{},\mu )`$ plane with $`t=1`$ and $`\mathrm{\Delta }`$ held fixed at some non-zero value in order to describe a RVB liquid of singlet pairs.

We can use symmetry arguments to show that not all states in this space are distinct upon projection and $`|0(t^{},\mu )=|0(t^{},\mu )`$. In brief, we can change $`t^{}t^{}`$, $`\mu \mu `$ by a global particle-hole transformation $`c_{𝐫\sigma }^{}(1)^{x+y}c_{𝐫\sigma }^{}`$ in the wavefunction, also redefining the $`|\mathrm{vac}`$ since empty sites transform into doubly occupied ones. However, with exactly one particle per site such a transformation interchanges $``$-spins and $``$-spins at every site for any configuration. Since the BCS wavefunctions are spin-singlets this leaves the state invariant. We thus restrict attention to $`t^{}0`$ in Fig. 1(A).

Short-range RVB limit: Note that large negative $`\mu `$ with $`\mu <4(|t|+|t^{}|)`$ is the bosonic limit randeria  of the BCS wavefunction, with $`\phi (𝐫𝐫^{})`$ decaying exponentially in real space. In this limit $`|0`$ may be viewed as a short-range RVB state. Earlier studies of closely related dimer models kivelson ; read  indicate that, for short-range singlet bonds there exist four topological sectors on a torus, which may be straightforwardly related to sectors with/without visons in the $`\widehat{x}`$, $`\widehat{y}`$ directions. Further, two dimer configurations from different sectors have exponentially small overlap $`\mathrm{exp}(\alpha L)`$. A state is a superposition of many configurations, and assuming that individually small overlaps do not add up coherently, the overlap of two states from different sectors would vanish as $`L\mathrm{}`$.

We have numerically checked (see below) that the overlap vanishes for large $`L`$ validating the above assumption, and implying fractionalization in the short range RVB limit. Using the symmetry discussed above, we conclude that states in the region with (positive) $`\mu >4(|t|+|t^{}|)`$ are also fractionalized, even though $`\phi (𝐫𝐫^{})`$ is not obviously short-ranged here.

Analytical insights for $`\mu =\mu _{\mathrm{BCS}}`$: We next focus on states for which $`\mu =\mu _{\mathrm{BCS}}(t^{})`$ (see the curve in Fig. 1(A)), where $`\mu _{\mathrm{BCS}}`$ is the value at which the unprojected BCS state is at half-filling. The unprojected state may then be viewed as the ground state of a d-wave BCS mean field (MF) Hamiltonian and it is a coherent superposition of number eigenstates, sharply peaked at the correct mean density. Gutzwiller projection picks out the $`N=L^2`$ contribution with no double occupancy.

Define local $`SU(2)`$ gauge transformations affleck  $`𝒰`$, generated by $`T_i^+=c_i^{}c_i^{}`$, $`T_i^{}=c_ic_i`$, $`T_i^z=_\sigma c_{i\sigma }^{}c_{i\sigma }^{}1`$, which mix empty and doubly occupied sites as an $`SU(2)`$ doublet, but act trivially on the projected subspace with exactly one particle per site. Gutzwiller projection is then equivalent to projection onto the $`SU(2)`$ singlet subspace, and we may write any state $`𝒫|\mathrm{\Phi }=_𝒰𝒰|\mathrm{\Phi }`$, where the integral is over all group elements $`𝒰_\theta =\mathrm{exp}(i_i\stackrel{}{T}_i\stackrel{}{\theta }_i)`$. The overlap of interest can then be written as

$$V|0=_𝒰(hc/2e)|𝒰|BCS,$$
(2)
which now reduces the problem to computing overlaps of unprojected states.

For a nonbipartite lattice (e.g., with $`t^{}0`$), one cannot gauge away the off-diagonal term in the d-wave BCS MF Hamiltonian and $`𝒰|BCS`$ is a SC state (or vortex vacuum) for arbitrary $`𝒰`$ ivanov . Then $`(hc/2e)|𝒰|BCS=0`$ for arbitrary $`𝒰`$ and $`V|0`$ vanishes. Thus non-bipartiteness is a sufficient condition for fractionalization in d-wave states when $`\mu =\mu _{\mathrm{BCS}}`$ previous .

On the other hand, for a bipartite system (in our case $`t^{}=0`$), it is well known affleck  that one can gauge away the off-diagonal part of the MF Hamiltonian (provided it doesn’t break time reversal) and transform the BCS state into a staggered flux state. Now $`(hc/2e)|𝒰|BCS`$ is non-zero for some choice $`𝒰`$ and one cannot use eq. (2) to argue that $`V|0`$ vanishes.

Overlaps within $`Z_2`$ gauge theories: Before proceeding to the numerical results, it is useful to understand how overlaps like $`V|0`$ are expected to scale with system size $`L^2`$. In a non-fractionalized phase the overlap will remain non-zero as $`L\mathrm{}`$. However, in a fractionalized phase we argue that it should vanish exponentially for large $`L`$. This is easy to see when the matter fields are gapped (e.g., the short-range RVB limit) and may be integrated out to obtain a $`2+1`$ dimensional $`Z_2`$ gauge theory. Deep in the deconfining phase of a pure gauge theory, we have evaluated the overlap perturbatively and it scales as $`V|0\mathrm{exp}(L/\xi )`$. To understand this result note that the overlap in the deconfined phase behaves similar to the spin-spin correlation of the dual Ising model at distance $`L`$, where the dual spins represent the gapped visons. Since the dual Ising model is disordered when the gauge theory is deconfining, the exponential decay follows. Another way to understand this result is that the overlap is nonzero due to tunneling of the vison (a “string” of length $`L`$) out of the hole of the cylinder. Such tunneling events involve intermediate states in which the vison pierces the surface of the cylinder. This costs an energy proportional to the vison gap, and leads to an overlap exponentially small in $`L`$. The situation with gapless matter fields is less well studied, even though numerically (see below) we find clear evidence for exponential decay.

Numerical results: We finally discuss the results of numerical evaluation of vison overlaps using the variational Monte Carlo method. For computational simplicity numerics  we focus on the overlap $`V_x|V_y`$. We choose $`t=1`$ and fix $`\mathrm{\Delta }=1.25`$ for all the calculations reported below. The following numerical results are obtained in regimes where we have given analytical arguments above. (1) We have explicitly checked that the overlaps are invariant under $`(t^{},\mu )(t^{},\mu )`$. (2) In the regime of short-range RVB we find that the system is $`Z_2`$ fractionalized; see Fig. 2(A) for large negative $`\mu `$. (3) On the curve $`\mu =\mu _{\mathrm{BCS}}`$ we find that the non-bipartite systems are fractionalized. (4) For the special point $`(t^{},\mu )=(0,0)`$ with bipartite symmetry, we indeed find that the overlaps are non-zero, independent of $`L`$ (see Fig. 2(A)), the vison is not well defined and the system is not fractionalized. Although consistent with the analytic arguments, the numerical results are important because the arguments were not rigorous for (2) and (3), and silent on point (4).

We next evaluate overlaps at $`t^{}=0.0,0.25,0.5`$ for $`1.5\mu 1.5`$ which, as we discuss below, covers the region of interest for possible spin-liquid insulators in the vicinity of the high Tc cuprates. We show in Fig. 2 the overlap $`V_x|V_y`$ as a function of system size $`L`$ for a range of $`t^{}`$ and $`\mu `$. We find that at small $`L`$ the overlaps are finite and then cross over on a scale $`\xi ^{}`$ to an asymptotic decay $`\mathrm{exp}(2L/\xi )`$. This behavior can be simply described by the functional form $`V_x|V_y(L)=a(1\mathrm{tanh}((L\xi ^{})/\xi ))`$ which appears to fit the data well. (The factor of two is included in the asymptotic decay definition since we are looking at overlap of states with two visons in different directions). In the regime of parameters where we can reliably extract both $`\xi ^{}`$ and $`\xi `$, we find $`\xi ^{}2\xi `$. More generally, it is hard to access the asymptotic behavior in the region of interest, and we can extract crossover scale $`\xi ^{}`$ more accurately than $`\xi `$.

We plot the inverse length scale $`1/\xi ^{}`$ for $`t^{}=0,0.25,0.5`$ at different $`\mu `$ in Fig. 3. For $`t^{}=0`$, we see that $`\xi ^{}\mathrm{}`$ as $`\mu 0`$, fully consistent with the finite overlap independent of $`L`$ on accessible system sizes for $`(t^{},\mu )=(0,0)`$ seen in Fig. 2(A). It also clear that $`\xi ^{}`$ remains finite everywhere away from the origin though it may become quite large in its vicinity.

Since fractionalization manifests itself at length scales larger than $`\xi ^{}`$, our numerical results strongly suggest that all points in the phase diagram of Fig. 1(A) are fractionalized, except for the origin where $`\xi ^{}`$ diverges. $`(t^{},\mu )=(0,0)`$, with its special bipartite symmetry, is a non-fractionalized “singular point” in the space of spin-liquid insulators that we study.

Vison gaps: Now the question arises: what regime in parameter space is relevant for the cuprates? One possibility is that there are no fractionalized states in the vicinity of the observed phases in real materials. An alternative worth exploring in view of the success of RVB wavefunctions in understanding the SC state paramekanti , is that we take the insulating limit (hole doping $`x0`$) of $`𝒫|BCS`$. Optimizing variational parameters $`\mu ,\mathrm{\Delta }`$ for the large-$`U`$ Hubbard model ($`U=12,t^{}=1/4`$), we find $`\mathrm{\Delta }1.25`$ while $`0.3\mu 0.3`$ as shown in Fig. 1(B). For this region (the shaded ellipse in Fig. 1(A)) we conclude from Fig. 3(B) that $`\xi ^{}25`$ lattice spacings.

We now convert this to an estimate of the vison gap $`E_v`$. We expect $`E_v=\alpha J/(\xi ^{})^z\alpha J/\xi ^{}`$ which vanishes at the “singular point” with bipartite symmetry. Here $`J`$ is the nearest-neighbor superexchange, $`\alpha 𝒪(1)`$ is a dimensionless constant and the dynamical exponent $`z1`$. For $`\alpha =1`$ and $`J=1200K`$, the estimated $`E_v50`$K (but this energy could be higher if $`\alpha `$ is larger). More concretely, we emphasize that proximity to the ‘bipartite point’ $`(t^{}=0,\mu =0)`$ can lead to a very small vison gap for fractionalized RVB states. This is consistent with existing experimental bounds moler  and suggests that fractionalization is apparent only at very low $`T`$ and does not play a role in pseudogap anomalies.

The ideas and methods developed here should be explored further in highly frustrated magnets, with large $`t^{}`$ or ring-exchange interactions ring-exchange , which may be more promising candidates for observing fractionalization.

Acknowledgments: We thank T. Senthil for extensive discussions and for generously sharing his ideas and unpublished notes with us. We thank also L. Balents, D. M. Ceperley, M. P. A. Fisher, E. Fradkin and A. Melikidze for helpful discussions. MR and NT gratefully acknowledge the hospitality of the Physics Department at University of Illiois and support through DOE grant DEFG02-91ER45439 and DARPA grant N0014-01-1-1062. AP was supported by NSF DMR-9985255 and PHY99-07949 and the Sloan and Packard foundations. We acknowledge the use of computational facilities at TIFR including those provided by the DST Swarnajayanti Fellowship.


text:
# Matter-wave gap solitons in atomic bandgap structures

## Abstract

We demonstrate that a Bose-Einstein condensate in an optical lattice forms a reconfigurable matter-wave structure with a band-gap spectrum, which resembles a nonlinear photonic crystal for light waves. We study in detail the case of a two-dimensional square optical lattice and show that this atomic bandgap structure allows nonlinear localization of atomic Bloch waves in the form of two-dimensional matter-wave gap solitons.

Recent experiments Christiani et al. (2002); Greiner et al. (2001, 2002); Denschlag et al. (2002) have demonstrated that a Bose-Einstein condensate (BEC) loaded into an optical lattice is a perfect test ground for a range of fascinating physical effects. In particular, this is because the effective potential of an optical lattice can be easily manipulated by changing the geometry, polarization, phase, or intensity of the laser beams. Due to inherent nonlinearity of the coherent matter waves introduced by inter-atomic interactions, BEC in a lattice potential can form a periodic nonlinear system, which is expected to display rich and complex dynamics.

On the other hand, propagation of light in dielectric structures with a periodic variation of the refractive index is receiving growing attention. Photonic bandgap (PBG) materials Joannopoulos et al. (1995) – artificial periodic structures with a high index contrast – can be used to effectively control the flow of light. The study of photonic crystals made of a Kerr nonlinear material, the so-called nonlinear photonic crystals, has revealed that such structures exhibit a wealth of nonlinear optical phenomena and, in particular, they can support self-trapped nonlinear localized modes of the electromagnetic field in the form of optical gap solitons Mingaleev and Kivshar (2001); Slusher and Eggleton (2003). Dynamically reconfigurable PBG structures - optically-induced refractive index gratings in nonlinear materials - are now offering new ways to control light propagation and localization Evfremidis et al. (2002).

As demonstrated in this letter, BEC in an optical lattice can be regarded as a fully reconfigurable analog of a nonlinear photonic crystal for matter waves - an ”atomic band gap” (ABG) structure. A deep analogy between coherent light and matter waves suggests that the concepts employed in the study of nonlinear PBG structures can be borrowed for analyzing BEC in optical lattices. Many properties of ABG structures can potentially be exploited for high-precision control and manipulation of coherent matter waves in a similar way to how PBG structures are used to manipulate light. The effect of nonlinear localization in bandgaps is one of these properties.

Localized nonlinear excitations of coherent matter waves - bright atomic solitons \- could be very useful for applications such as atomic interferometry due to their robust nature. However, their generation has so far been experimentally achieved only in BECs with attractive atomic interactions (analogous to self-focusing optical media) which are unstable against collapse above a small ($`10^3`$) critical number of particles.

In theory, a shallow 1D optical lattice can support bright matter-wave solitons even in repulsive BEC Zobay et al. (1999) with large atom numbers. These solitons, described in a framework a coupled-mode theory Zobay et al. (1999), are localized on a large number of lattice wells, and are predicted to exist only in atomic bandgaps. In the opposite case of a deep lattice, the condensate can be described by the superposition of ground state modes in the individual wells. The mean-field treatment of the condensate in this regime leads to a discrete equation which admits solutions in the form of stationary modes localized on a few lattice sites - discrete solitons Trombettoni and Smerzi (2001); Abdullaev et al. (2001), in complete parallel with spatial solitons in periodic optical structures Kivshar and Agrawal (2003) and localized modes of atomic lattices Flach and Willis (1998).

The theory of the nonlinear localized matter waves in optical lattices is mostly limited to 1D case 111After submission of this manuscript, the generation of 2D arrays of localized states via modulational instability of Bloch waves was theoretically demonstrated in B.B. Baizakov, V.V. Konotop, and M. Salerno, J. Phys. B: At. Mol. Opt. Phys. 35, 5105 (2002). However, as can be deducted from the analogous studies of 2D PBG structures Mingaleev and Kivshar (2001), the nonlinear localization of in BEC in higher-dimensional lattices is qualitatively different because both the symmetry and dimensionality of the lattice start to play an important role in the formation and properties of the bandgap structure and corresponding nonlinear modes. In this Letter, we analyze the as yet unexplored problem of the existence and stability of 2D matter-wave gap solitons of BEC with repulsive interatomic interactions loaded into optical lattices, which are analogous to localized states of light waves in nonlinear photonic crystals. We show that their accurate description is only possible within a full mean-field model of BEC in a periodic lattice potential, beyond the tight-binding approximation or coupled-mode theory often employed for the study of BEC in 1D optical lattices.

The dynamics of a three-dimensional BEC cloud loaded into a two-dimensional optical lattice can be described by the Gross-Pitaevskii (GP) equation for the macroscopic wavefunction, $`\mathrm{\Psi }(𝒓,t)`$ of the condensate,

$$i\frac{\mathrm{\Psi }}{t}=\left[\frac{1}{2}^2+V(𝒓)+\gamma _{3\mathrm{D}}|\mathrm{\Psi }|^2\right]\mathrm{\Psi }.$$
(1)
This equation has been made dimensionless by using the natural length, time, and energy scales of the optical lattice. The characteristic length scale in this model, $`a_L=k_L^1=d/\pi `$, defines the size of a single well of the lattice potential through the lattice constant that depends on the laser wavelength as $`d=\lambda /2`$. The depth of the lattice potential, measured in the lattice recoil energy $`E_L=\mathrm{}^2/(2ma_L^2)`$, is proportional to the laser intensity and can be varied from $`0`$ to $`20E_L`$ in experiment. The characteristic frequency scale is $`\omega _L=E_L/\mathrm{}`$. The interatomic interactions are characterized by the coefficient $`\gamma _{3\mathrm{D}}=8\pi (a_s/a_L)`$, where $`a_s`$ is the $`s`$-wave scattering length, positive for repulsive interactions.

In the simplest case of a square optical lattice, the combined potential of the lattice and magnetic trap, $`V(𝒓)=V_L(𝒓)+V_T(𝒓)`$, can be written as $`V(𝒓)=V_0\left(\mathrm{sin}^2x+\mathrm{sin}^2y\right)+\omega _{i}^{}{}_{}{}^{2}𝒓^2/2,`$ where $`V_0`$ is the amplitude of the optical lattice, and $`\omega _i`$ are the (normalized) trapping frequencies in the corresponding directions. The current experiments with 2D optical lattices are performed with a strongly anisotropic BEC cloud, with the ratio of the trapping strengths $`\mathrm{\Omega }=\omega _x/\omega _{y,z}10^1`$ Greiner et al. (2001). To simplify our analysis, we assume that the relatively weak magnetic confinement characterized by trap frequencies $`\omega _i`$ has small effect on the stationary states of the condensate in the lattice. Under this assumption, the trap component of the confining potential in the lattice plane can be neglected, and the model can be reduced to a two-dimensional GP equation by assuming separable solutions of the form, $`\psi (𝒓;t)=\mathrm{\Phi }(z)\psi (x,y;t)`$ where the function $`\mathrm{\Phi }(z)`$, describing the condensate in the direction perpendicular to the lattice, is a solution of the 1D quantum harmonic oscillator problem, with the normalization condition $`\mathrm{\Phi }^2(z)𝑑z=1`$. Then the standard dimensionality reduction procedure Louis et al. (2003), leads to the following equation for the 2D condensate wavefunction in the 2D lattice potential:

$$i\frac{\psi }{t}=\left[\frac{1}{2}_{}^2+V_\mathrm{L}(𝒓)+|\psi |^2\right]\psi ,$$
(2)
where $`_{}=^2/x^2+^2/y^2`$, $`V_\mathrm{L}(𝒓)V_\mathrm{L}(x,y)`$ is the periodic potential of the optical lattice, and the wavefunction is rescaled as $`\psi \psi \sqrt{\gamma _{2\mathrm{D}}}`$, with $`\gamma _{2\mathrm{D}}=\gamma _{3\mathrm{D}}/\sqrt{2}`$.

Spectrum of atomic Bloch waves in the optical lattice can be found by exploiting the analogies with the theory of single-electron states in crystalline solids. Stationary states of the condensate in an infinite periodic potential of a 2D optical lattice are described by solutions of Eq. (2) of the form: $`\psi (𝒓,t)=\varphi (𝒓)\mathrm{exp}(i\mu t)`$, where $`\mu `$ is the chemical potential, corresponding to the energy level of the stationary state in the lattice potential. The case of noninteracting condensate formally corresponds to Eq. (2) being linear in $`\psi `$. According to the Bloch theorem, the stationary wavefunction can then be sought in the form $`\varphi (𝒓)=u_𝒌(𝒓)\mathrm{exp}(i𝒌𝒓)`$, where the wavevector $`𝒌`$ belongs to a Brillouin zone of the square lattice, and $`u_𝒌(𝒓)=u_𝒌(𝒓+𝒅)`$ is a periodic (Bloch) function with the periodicity of the lattice. For the values of $`𝒌`$ within an $`n`$-th Brillouin zone, the dispersion relation for the 2D Bloch waves, $`\mu _n(𝒌)`$, is found by solving the following linear eigenvalue problem

$$\left[\frac{1}{2}(i\mathbf{}_{}+𝒌)^2+V_L(𝒓)\right]u_{n,𝒌}=\mu _n(𝒌)u_{n,𝒌}.$$
(3)
The eigenvalue problem (3) is simpler than a general 2D problem due to the separability of the lattice potential. The typical dispersion relation for the energy of the Bloch states in the lowest Brillouin zones is shown in Fig. 1 for a moderate value of the lattice depth $`V_0=5`$. The dispersion diagram is presented in the reciprocal lattice space, and the dispersion relations are calculated along the characteristic high-symmetry directions of the irreducible Brillouin zone (see Fig. 1, bottom right). It can be seen that the absolute bands and gaps are determined by the spectra of the Bloch waves in the middle (point $`\mathrm{\Gamma }`$) and on the edge (point $`\mathrm{M}`$) of the Brillouin zone.

To understand how the variation in the lattice parameters affects the band structure, we note that, owing to the scaling properties of the model, any increase in the well spacing, $`a_L`$, translates into decreasing the well depth $`V_0`$. Thus the global behaviour of the band structure can be understood by examining the Bloch wave spectrum as a function of the lattice depth. Figure 2 presents the chemical potentials corresponding to the atomic Bloch waves at the edges of the absolute bands plotted as functions of $`V_0`$. The dashed line $`\mu =V_0`$ (dashed line in Fig. 1) separates the quasi-unbound ($`\mu >V_0`$) and strongly bound ($`\mu <V_0`$) condensate states in the lattice (see also Ref. Drese and Holthaus (1997)). By varying the interwell separation or the amplitude of the lattice, different regimes of the atomic bandgap structure can be accessed. The tight-binding regime of the condensate dynamics corresponds to the domain $`V_0V_0^{}=\mu `$, where the bands ”collapse” to discrete levels of bound states in a single isolated well (see Fig. 2). The narrow-gap regime of a shallow lattice, where the coupled-mode theory applies, is found for the opposite case, $`V_0V_0^{}`$.

The study of the 1D atomic bandgap structures Alfimov et al. (2002a); Louis et al. (2003) indicates that the correct structure and dynamical properties of the nonlinear localized modes can only be revealed through the analysis of the full continuous mean-field GP model that bridges the gap between the coupled-mode theory and discrete tight-binding approach. One promising development in approximate analysis would be generalization of the Wannier function formalism Alfimov et al. (2002b) to the 2D case. However, the complexity involved in the generation of optimally localized higher-dimensional Wannier functions Marzari and Vanderbilt (1997) may not provide significant advantage compared with the direct solution of the model equation (2).

Here, we find spatially localized stationary solutions of Eq. (2) numerically. Our numerical procedure involves minimization of the functional $`=f^{}fd𝒓`$, where $`f(\varphi )=[_{}^2+\mu V_L(𝒓)|\varphi |^2]\varphi `$, by following a descent technique with Sobolev preconditioning García-Ripoll and Pérez-García (2001). The minimization procedure yields a stationary state when $`(\varphi )=0`$. Previous applications of this method to the analysis of optical solitons has shown that it is suitable for tackling both the fundamental and higher-order nonlinear localized modes. Our numerical method is independent of the lattice symmetry; it offers a versatile technique for finding and analyzing nonlinear localized modes in a range of different bandgap structures, and for making accurate predictions about their density, degree of localization, and spatial structure in different areas of the parameter space.

We have identified different families of matter-wave gap solitons of the repulsive BEC in the 2D atomic bandgap structure. These solitons can exist in all gaps, excluding the semi-infinite gap of the spectrum (marked by 1 in Figs 1, 2) below the first band. We note in passing that only conventional nonlinear localization of BEC with attractive interactions is possible in the semi-infinite gap. The change in the peak amplitude of the the lowest-order ”fundamental” gap solitons across the first gap, along with their spatial localization properties, are demonstrated in Fig. 3, for $`V_0=5`$. Gap solitons of repulsive BEC branch off the lower-energy (i.e. smaller values of $`\mu `$) gap edge (see Figs. 3, top). Near this edge, i.e. at low particle numbers and amplitudes, their central peak is localized around a single well \[Figs. 3 (a)\]. Further within the gap, the localized mode of the BEC contains larger numbers of atoms, and has a strongly localized central peak. Near the linear tunnelling threshold, $`\mu V_0`$, the nonlinear localized states develop extended tails in the orthogonal directions $`(x,y)`$, along which the tunnelling is assisted by the lower inter-well barrier heights (compared to $`2V_0`$ in diagonal directions). Although the stability analysis of the states localized on the extended background (for $`\mu >V_0`$) is difficult due to the large computational domain, we have confirmed that strongly localized 2D modes are dynamically stable. In the lattice geometry under consideration, the line $`\mu =V_0`$ lies within the first gap (marked 2 in Figs. 1, 2), and only localized modes within this gap have been analyzed.

The higher-amplitude gap solitons have a nontrivial structure of the tails \[see Fig. 3 (b)\], i.e. the zeros of the matter wave are centered within individual lattice wells (see Figs. 3, insets), rather than on the potential maxima. This structure is determined by the spatial structure of the Bloch wave at the lower-energy edge of the second band (point $`M`$ in Fig. 1), and cannot be described within the framework of the discrete model.

We have also found families of higher-order gap modes that can be identified as bound states of the fundamental gap solitons. These states can be centered on the lattice potential minima or maxima, similarly to higher-order ”odd” and ”even” states in 1D lattices Louis et al. (2003), and can exhibit symmetry-breaking instabilities. Two examples of the lowest-order in- and out-of-phase odd modes, centered on a lattice maximum, are shown in Fig. 4.

The crucial issue of the potential observation of 2D gap solitons is the stability of the 2D localized state in a 3D BEC cloud, i.e. dynamical stability of the 3D state $`\mathrm{\Psi }(x,y,z;t)`$ with the initial condition given by $`\mathrm{\Phi }(z)\varphi (x,y)`$, where $`\varphi (x,y)`$ is the two-dimensional stationary gap soliton. The analysis of stability of 2D solitons in a 3D BEC cloud, wich will provide clues for possible experimental observation of multi-dimensional gap solitons, is the subject of our separate study.

In conclusion, we have demonstrated that the interaction of BEC with the lattice potential is analogous to the light scattering by a nonlinear photonic bandgap structure. We have studied the properties of two-dimensional atomic bandgap structures and demonstrated the existence of gap solitons, the spatially localized states of BEC existing in the gaps of the matter-wave spectrum. We believe the analogy between the physics of BEC in optical lattices and photonic crystals can be useful for revealing many novel features of the matter-wave dynamics in reconfigurable atomic bandgap structures.

The authors thank C.M. Savage and J.J. García-Ripoll for invaluable advice. This work is supported by the Australian Partnership for Advanced Computing and Australian Research Council. The authors are part of the ARC Centre for Quantum-Atom Optics.


text:
# Star Formation at the Twilight of the Dark Ages: Which Stars Reionized the Universe?

## 1 Introduction

When and how the Universe became reionized is one of the fundamental questions in cosmology. There is currently a great deal of interest in this question due to recent observational and theoretical developments. The discovery of several $`z6`$ quasars in the SDSS whose spectra are consistent with showing zero flux below Lyman-$`\alpha `$ (a ‘Gunn-Peterson’ trough) may indicate that the IGM was predominantly neutral at $`z6`$ (Fan 2001; Becker et al. 2001). The ionization history of the Universe may also be constrained via observations of the cosmic microwave background (CMB). The recently released first-year results from the WMAP satellite constrain the reionization redshift to the range $`z_{\mathrm{reion}}=17\pm 5`$ (Kogut et al. 2003).

What is the nature of the sources that produced the photons responsible for this transition from neutral to ionized? In Cold Dark Matter (CDM) models, halos large enough to form significant amounts of molecular hydrogen start appearing around redshift $`z30`$ (Tegmark et al. 1997). The stars that formed in these halos, out of pristine, metal-free and magnetic-field-free gas, with molecular hydrogen as the only coolant, were almost certainly quite different from the stars we see around us today. Numerical simulations suggest that these first stars were extremely massive, on the order of a few hundred solar masses (Abel et al. 2000, 2002; Bromm et al. 2002). This first generation of very massive, metal free stars is often referred to as ‘Pop III’. Because of their high temperatures and low metallicities, Pop III stars may produce up to $`20`$ times as many hydrogen-ionizing photons per solar mass as Pop II stars (Bromm et al. 2001). It is therefore natural to think that these stars may play an important role in early reionization. Ironically, less is known about the *second* generation of stars, precisely because the metals, magnetic fields, and photons produced by the first stars complicate the situation. The key processes that determine how efficiently Pop III stars can form and when the transition to Pop II occurs — destruction of H2 by Lyman-Werner photons, catalysis of H2 by X-rays, and the production, dispersal, and mixing of heavy elements — remain poorly understood (Machacek et al. 2001, 2003; Ricotti et al. 2001, 2002; Cen 2002; Yoshida et al. 2003).

Around $`z30`$–20, halos that are large enough to cool by atomic processes start to collapse. Two processes are likely to be responsible for regulating star formation in these halos. If the Universe has already been reionized, the UV background will prevent gas from collapsing into halos with temperatures smaller than the effective Jeans mass, $`2\times 10^5`$ K, corresponding to circular velocities $`V_c30`$–50 km/s (e.g. Thoul & Weinberg 1996; Gnedin 2000). We hereafter refer to the suppression of star formation resulting from this effect as photoionization ‘squelching’. Supernovae and massive stars also regulate star formation by heating the ISM and driving winds that remove the gas from the galaxy (Dekel & Silk 1986; Mac Low & Ferrara 1999).

In this Letter, we present predictions for the cosmic star formation rate density (SFRD) due to both Pop II and Pop III stars from $`z30`$–3. We compare our results with observations at the low redshift end ($`z3`$–6). We use a semi-analytic model to explore the effects of photoionization squelching and supernova feedback on these predictions, and consider several possible reionization histories (including a model with multiple reionizations). We then assess the relative contribution from Pop II and Pop III stars to the ionizing photon budget during the redshift range relevant to the epoch of reionization as constrained by WMAP and the SDSS quasar observations ($`z6`$–30). Throughout, we assume cosmological parameters consistent with the WMAP data (Spergel et al. 2003): matter density $`\mathrm{\Omega }_m=0.3`$, baryon density $`\mathrm{\Omega }_b=0.044`$, dark energy $`\mathrm{\Omega }_\mathrm{\Lambda }=0.70`$, Hubble parameter $`H_0=70`$ km/s/Mpc, fluctuation amplitude $`\sigma _8=0.9`$, and a scale-free primordial power spectrum $`n_s=1`$.

## 2 The halo-collapse model

In order for gas to be converted into stars, the first-order condition is that it must be inside a collapsed halo of sufficient mass to allow the gas to cool and become dense. Following this line of argument, we may model the global star formation rate density by assuming that it is proportional to the rate at which gas collapses into halos in a given mass range:

$$\dot{\rho }_{}=e_{}\rho _b\frac{\mathrm{d}F_h}{\mathrm{dt}}(M>M_{\mathrm{crit}}),$$
(1)
where $`\frac{\mathrm{d}F_h}{\mathrm{dt}}(M>M_{\mathrm{crit}})`$ is the time derivative of the fraction of the total mass in collapsed halos with masses greater than $`M_{\mathrm{crit}}`$, obtained from the halo mass function $`dn_h/dM(M,z)`$ given by the Press-Schechter model (Press & Schechter 1974) or one of its variants, and $`\rho _b`$ is the mean density of baryons. One may then assign the efficiency of conversion of gas into stars $`e_{}`$ and the critical mass $`M_{\mathrm{crit}}`$ for separate populations depending on the main coolant and the mode of star formation. For example, $`M_{\mathrm{crit}}`$ corresponds to a halo virial temperature of about $`10^4`$ K for halos that cool via atomic processes, while $`T_{\mathrm{crit}}100`$ K for molecular cooling. Pop III stars are generally assumed to form in the lower temperature, H2-cooled halos, with a much lower efficiency than Pop II stars, which are associated with larger, HI-cooled halos (for a more detailed discussion, see e.g. Loeb & Barkana 2001). We hereafter refer to this as the ‘halo-collapse’ model.

## 3 Merger Tree models

Variants on simple models like the one presented above have been used in many recent analytic studies of early star formation and reionization (e.g. Cen 2002; Wyithe & Loeb 2002, 2003; Venkatesan et al. 2003; Haiman & Holder 2003). However, there is a well-developed semi-empirical approach to modeling the physics of atomic cooling, Pop II star formation and chemical enrichment, and supernova feedback, within the framework of hierarchical merging predicted by CDM models. The effect of photoionization squelching on Pop II star formation has also been included in some semi-analytic models (Kauffmann et al. 1993; Somerville 2002; Benson et al. 2002). This approach has been used in a large number of studies of galaxy formation at lower redshift $`z0`$–5 (e.g. Kauffmann et al. 1993; Cole et al. 1994; Kauffmann et al. 1998; Cole et al. 2000; Somerville & Primack 1999; Somerville et al. 2001). It is interesting to see how the results based on these more realistic recipes compare with the simple ‘halo-collapse’ model described above, and to study the relative importance of the various processes that are expected to regulate star formation at very high redshift.

Here we use the models developed in Somerville & Primack (1999, SP) and Somerville et al. (2001, SPF), with photoionization squelching added as described in Somerville (2002), using a recipe based on the numerical results of (Gnedin 2000). We follow halo merger histories down to halos with temperature $`T_{\mathrm{vir}}=10^4`$ K, where atomic cooling becomes possible. We shall refer to these models as the ‘merger tree’ models.

## 4 Results

### 4.1 The Cosmic Star Formation History

We show the predicted star formation rate density (SFRD) for Pop II and Pop III stars in Fig. 1. For Pop III stars, we have used the halo collapse model (Eqn. 1) with $`e_{}^{\mathrm{III}}=0.001`$ and $`M_{\mathrm{crit}}^{\mathrm{III}}=1.0\times 10^6h^1M_{}`$. Also shown in Fig. 1a are the results from detailed numerical hydrodynamic simulations of Pop III star formation in a cosmological volume by Yoshida et al. (2003), for an assumed characteristic Pop III mass of $`100M_{}`$. We see that with these parameter choices, our simple recipe reproduces the simulation results fairly well. Of course, we do not really expect Pop III star formation to continue to $`z3`$. However, in the absence of more detailed modeling of metal production and mixing, we do not know when the Pop III mode will ‘shut off’.

Fig. 1a also shows the prediction of the halo collapse model for Pop II stars, with $`M_{\mathrm{crit}}^{\mathrm{II}}`$ corresponding to a temperature111See SP for the conversions between virial mass, velocity, and temperature as a function of redshift of $`10^4`$ K and efficiencies of $`e_{}^{\mathrm{II}}=0.1`$ and $`e_{}^{\mathrm{II}}=1.0`$. The $`e_{}^{\mathrm{II}}=1`$ halo collapse model represents a sort of upper limit for star formation in halos that cool via HI, as it would require all new gas flowing into these halos to cool and form stars instantaneously. The predicted SFRD may be compared with observational estimates and limits on this quantity at ‘low’ redshift $`3z6`$ (see figure caption). The assumed efficiencies of $`e_{}^{\mathrm{II}}=0.1`$ and $`e_{}^{\mathrm{II}}=1`$ approximately bracket the range in observational estimates. We note in passing that the Lanzetta et al. high estimate exceeds even the extreme case of $`e_{}^{\mathrm{II}}=1`$ at $`z=8`$–10.

We also show the star formation history produced in the more realistic semi-analytic merger tree models described in Section 3. To study the sensitivity of our results to our recipes for star formation, photoionization squelching, and supernova feedback, we investigate different choices of ingredients and parameters. The merger tree models are summarized in Table 1. The reader only interested in the main result may skip directly to Section 4.2.

Models MT-1a and MT-1b do not include supernova feedback or photoionization squelching. Star formation is then regulated only by the rate at which gas can cool via atomic processes and collapse, and by the star formation timescale $`\tau _{}=m_{\mathrm{cold}}/\dot{m}_{}`$. We consider two choices for this timescale, which we refer to as ‘constant efficiency’ and ‘accelerating’. In the constant efficiency (CE) recipe, $`\tau _{}`$ is a constant (we take $`\tau _{}=0.1`$). In the ‘accelerating’ recipe, $`\tau _{}t_{\mathrm{dyn}}`$, where $`t_{\mathrm{dyn}}r_{\mathrm{vir}}/v_{\mathrm{vir}}`$ is the dynamical time of the halo. This is similar to the scaling observed in nearby galaxies (Kennicutt 1983, 1998) and is commonly used as an empirical basis for star formation recipes in semi-analytic models and hydrodynamic simulations. Referring to Fig. 1a, we see that the ‘accelerating’ scaling leads to dramatically more star formation at high redshift compared with the constant efficiency scaling. We use this recipe in the rest of the models.

In models MT-2a–c, we include the effect of photoionization squelching (but no SN feedback), with different assumed reionization histories, motivated by the joint WMAP and SDSS constraints discussed above. We show in Fig. 1a a model in which hydrogen is reionized at $`z_{\mathrm{reion}}=20`$ or $`z_{\mathrm{reion}}=15`$ and remains reionized, and a model with a ‘double reionization’ based on the model of Cen (2002), in which the Universe is reionized at $`z_{\mathrm{reion}}^{(1)}=15.5`$, recombines by $`z_{\mathrm{rec}}=13.5`$, and experiences a second reionization at $`z_{\mathrm{reion}}^{(2)}=6`$. The ‘dips’ in the star formation history are due to the suppression of gas infall by the photoionizing background, and demonstrate that a substantial fraction of the total star formation in the Universe would be taking place in small mass halos in the absence any kind of feedback.

Models MT-3a–c, shown in Fig. 1b, all assume the Cen (2002) double-reionization history described above, and also include supernova feedback. Note that the ‘dips’ are now much less dramatic, because much of the gas has already been removed from the small halos that are affected by photoionization squelching. An uncertain aspect of implementing supernova feedback in semi-analytic models is the fate of the reheated gas. In model MT-3a, the gas reheated by SN is removed from the disk but retained in the halo (‘retention’ feedback), in MT-3b the reheated gas is removed from the disk *and* dark matter halo of all galaxies (‘ejection’ feedback), and in model MT-3c the gas is ejected from the halo only if the halo virial velocity is less than 100 km/s. This threshold for ejection of gas by super-winds is motivated by theoretical arguments (Dekel & Silk 1986) and observations of nearby galaxies (Martin 1999). The global star formation rate changes by as much as a factor of six at $`z3`$ depending on these choices, but by less than a factor of two at $`z15`$. Note that the ‘dip’ following the first reionization is considerably more pronounced in the model with ejection feedback, as gas which has been accreted before reionization and ejected is not allowed to re-collapse in small halos while the photoionizing background is ‘switched on’. We consider models MT-3a and MT-3c to be the most realistic of the models considered here. Models with similar ingredients and parameter values have been shown to reproduce the luminosity function of galaxies at $`z0`$ (SP) as well as of Lyman break galaxies at $`z3`$ (SPF).

### 4.2 Production of Ionizing Photons

In Fig. 2 we show the cumulative number of hydrogen-ionizing photons per hydrogen atom in the Universe, produced by Pop II stars in our fiducial models (MT-3a and MT-3c), and by our model for Pop III star formation. For Pop II, we use the results of Leitherer et al. (1999) for the number of $`\lambda <912`$ Å photons produced by low metallicity stars with a Salpeter IMF. For Pop III, we assume that each star produces $`1.6\times 10^{48}`$ photons s$`{}_{}{}^{1}M_{}^{1}`$ for a lifetime of 3 million years (Bromm et al. 2001). We emphasize that we have shown the number of ionizing photons *produced*, without attempting to fold in the fraction of these photons that manage to escape from the galaxy or to propagate through the IGM. The ionization fraction $`x_e`$ is expected to scale as this quantity times $`f_{\mathrm{esc}}f_{\mathrm{ion}}/C_{\mathrm{clump}}`$ (e.g., Spergel et al. 2003), where $`f_{\mathrm{ion}}`$ is the number of ionizations per UV photon, $`f_{\mathrm{esc}}`$ is the fraction of ionizing photons that escape from the galaxy, and $`C_{\mathrm{clump}}`$ is the clumping factor, reflecting the clumpyness of the IGM. The recombination time at $`z1025`$ is $`0.2`$–0.7 times the Hubble time (Cen 2002). Observational constraints on $`f_{\mathrm{esc}}`$ for both nearby and high redshift galaxies range from a few percent to $``$ 50 percent (Leitherer et al. 1995; Hurwitz et al. 1997; Steidel et al. 2001). Clumping factors at high redshift $`z10`$ are expected to be on the order of $`C_{\mathrm{clump}}=2`$–10 (Cen 2002). Therefore the combination $`f_{\mathrm{esc}}f_{\mathrm{ion}}/C_{\mathrm{clump}}`$ is plausibly of order $`0.1`$–0.2 in the relevant redshift range.

In Fig. 2, we see that the total number of ionizing photons produced by Pop II stars in our fiducial model overtakes the contribution from Pop III stars at around $`z16`$–12. If reionization occurred near the lower redshift end of the WMAP range (which seems easier to reconcile with our model), then our results suggest that Pop II stars contributed at least half of the ionizing photons. If reionization occurred as late as $`z_{\mathrm{reion}}12`$, then Pop II stars (which are known observationally to exist) may even have been able to reionize the Universe on their own. One should also keep in mind that given the non-negligible rate of star formation (both Pop III and Pop II) even at $`z20`$, we might expect there to have been sufficient metal pollution to have halted the formation of very massive Pop III stars at quite a high redshift. However, the redshift of this transition is highly uncertain, and is almost certainly not sharp, but will occur at different epochs in different environments.

## 5 Discussion

We have shown predictions of the global star formation history from the time when the first stars began to form $`z30`$ until the epoch of the most distant observed galaxies $`z3`$–6. We summarize our main conclusions as follows:

* There are two main competing effects determining the Pop II star formation efficiency (SFE) at high redshift. (i) The star formation activity shifts toward smaller halos at earlier times, leading to decreased effective SFE as supernova feedback and photoionization squelching reduce the available cold gas supply in these halos. (ii) If the SFE is higher when the gas density is larger, as suggested by observations in nearby galaxies, then this is partially counteracted by the increased efficiency expected due to the higher characteristic gas densities.
* Our results suggest that the background of hydrogen-ionizing photons during reionization at $`z15`$–20 is roughly equally divided between Pop II and Pop III stars. While Pop III stars produce a larger number of ionizing photons per baryon, this is compensated by the lower SFRD of Pop III stars in the relevant redshift range $`z15`$–20.
* If $`z_{\mathrm{reion}}`$ is pushed to the lower limit of the range allowed by the WMAP results, the background produced by Pop II stars alone may be sufficient to reionize the universe, thus removing the only compelling *observational* argument for the existence of Pop III stars.

Using similar semi-analytic recipes for star formation applied within N-body simulations with radiative transfer, Ciardi et al. (2003) also concluded that Pop II stars alone can produce reionization histories that are within the WMAP constraints.

There are several uncertain factors that could shift our results by factors of a few in different directions. The characteristic masses of the Pop III stars could be a factor of 2–3 higher than we have assumed here, resulting in a corresponding scaling in the star formation and ionizing photon production rate. However, we may have underestimated the number of ionizing photons produced by Pop II stars at early times, as they may have lower metallicities than the Leitherer et al. (1999) models used here, and may have a somewhat top-heavy IMF (Larson 1998). This could lead to a factor of $`2`$–3 increase in the number of ionizing photons predicted by the Pop II stars (Tumlinson et al. 2003; Schaerer 2003). It is clear that more detailed modeling of chemical evolution is needed to better constrain these uncertainties, as well as the transition redshift between massive Pop III formation and Pop II formation with a normal IMF.

If correct, the implications of an ionizing background that is composed nearly equally of photons originating from Pop II and Pop III stars are interesting in several respects. The halos hosting the Pop II galaxies are rare peaks and will be highly clustered, while the lower-mass halos hosting the Pop III stars will be much less clustered, leading potentially to a rather complex topology for reionization. Because the Pop III stars will produce an even larger excess of helium ionizing photons, there are important implications for helium reionization as well (e.g. Venkatesan et al. 2003).

## Acknowledgments

We thank L. Moustakas for a careful reading of the manuscript and useful suggestions.