text:
# The structure of small lead clusters

## I Introduction

The structure of a cluster is one of its primary properties and one which has been intensely studied, experimentally and theoretically Johnston (2002); Alonso (2000). However, there is still much to be learnt about the fundamentals of cluster structure and the possible structures that can be formed. For atomic clusters with pair interactions, it is relatively well-understood how the form of the potential determines the observed structure. For example, the effects of the width of the potential well Doye et al. (1995); Doye and Wales (1997a) and oscillations in the potential Doye et al. (2001); Doye and Wales (2001); Doye et al. (cond-mat/0205374) have been systematically studied. However, for the systems that are of most interest, the interatomic interactions are usually much more complex. In particular, metal clusters, which are of great technological relevance Argo et al. (2002), have a strong many-body character to their bonding.

This presents a number of challenges to our understanding of cluster structure in metals. First, there is the possibility that new types of structure could emerge as a result of many-body effects. Although, the structural types observed for pair potentials are also frequently observed for metals, e.g. the competition between icosahedral, decahedral and close-packed structures is also common for metals Cleveland and Landman (1991), there are an increasing number of intriguing exceptions. One seemingly common feature for small metal clusters is to exhibit structures with no discernible order Garzón et al. (1998); Michaelian et al. (1999); Garzón et al. (2000); Taneda et al. (2001); Michaelian et al. (2002). However, it might be that the disorder is a result of new structural principles that cannot be fully satisfied at the sizes considered (hence the disorder), but which could lead to novel high symmetry structures at certain magic sizes Soler et al. (2001). Indeed, there are general grounds to expect high symmetry structures to emerge irrespective of the potential Wales (1998). Studies that just reoptimize known cluster structures will of course miss such new features, and so it is important that efficient global optimization algorithms are used.

Secondly, the many-body character makes it increasingly difficult to relate the observed structure back to the interactions, even when the assumed form for the many-body potential is relatively simple. There has been some interesting progress recently in this area, namely into the causes of the disordered structures Soler et al. (2000), and the effect of the range of the attraction and repulsion on the competition between icosahedral, decahedral and close-packed clusters Baletto et al. (2002), but there is much still to be discovered. This task is particularly important because of the difficulty in producing good empirical metal potentials (it is not feasible to study the sizes in which we are interested in any other way). One needs this kind of physical insight to understand the strengths and deficiencies of a potential and how it could be improved. It would also help one to discriminate between different potentials that purport to model the same material but give rise to different structures.

Lead clusters illustrate some of these challenges. The first theoretical study on large clusters by Lim et al. using a glue potential seemed to indicate that the most stable clusters at relatively small sizes (from at least $`N55`$) are face-centred cubic (fcc) Lim et al. (1992), the preferred bulk structure Young (1991). This conclusion was based on a comparison of the energies of a series of Mackay icosahedra and fcc cuboctahedra. It is quite unusual to see bulk structures already being favoured at such small sizes, but this finding was rationalized on the basis of the particularly small value of $`\gamma `$, the ratio of the surface energies of the $`\{100\}`$ and $`\{111\}`$ faces Lim et al. (1992). Usually, the Mackay icosahedra have an energetic advantage at small sizes, because exclusively having $`\{111\}`$ facets gives them an appreciably lower surface energy. Furthermore, these results were not inconsistent with the experiments in the literature at that time, which were mass spectroscopic studies on very small lead clusters Mühlbach et al. (1982); LaiHing et al. (1987), and electron diffraction experiments on very large clusters, which were identified as fcc, but with possibly some vestiges of amorphous structure Yokozeki (1978).

This basic picture though has recently been challenged by both experimental Hyslop et al. (2001) and theoretical Hendy and Hall (2001) results. Electron diffraction of clusters from 3 to 7 nm indicates that the largest clusters are dominated by decahedra, but for the smaller clusters it was not possible to obtain an adequate fit to the diffraction pattern, suggesting that alternative structural models need to be considered Hyslop et al. (2001).

Simulations of the melting and freezing of large lead clusters (modelled by the glue potential used by Lim et al. Lim et al. (1992)) unexpectedly revealed that for a certain size range ($`600<N<4000`$, at least) fcc structures are not lowest in energy Hendy and Hall (2001). Instead, a new type of icosahedral structure, which is more stable than the fcc structures, spontaneously formed both on freezing and on heating at temperatures just below that for melting. Similar structures had been previously seen in some simulations of large lead clusters but it was not recognised that they could be lowest in energy Lim et al. (1994). They resemble anti-Mackay icosahedra Doye and Wales (1997b), which have a Mackay icosahedral core but with most of the outer layer in ‘hexagonal close-packed’ (hcp) surface sites rather than the ‘fcc’ sites that would continue the packing in the Mackay icosahedra.

These results naturally raise intriguing questions about the structures of small lead clusters. What alternative structural models might be needed to understand the experimental results? Do lead clusters, as modelled by the glue potential, really favour fcc structures at small sizes? Lim et al. clearly showed that other standard forms were not lower in energy, so the structures would have to be somewhat unusual. Interestingly, a simple analysis using macroscopic properties as inputs suggested that lead would be a particular likely candidate for disordered clusters to be low in energy Soler et al. (2000).

Here, we address some of these issues by performing global optimization for lead clusters with up to 160 atoms. We pay particular attention to characterizing the structures of these clusters and to understanding why the potential favours the lowest-energy structures.

## II Methods

To model the lead clusters we use a glue potential Ercolessi et al. (1988) of the form

$$E=\underset{i<j}{}\varphi \left(r_{ij}\right)+\underset{i}{}U\left(n_i\right),$$
(1)
where $`\varphi (r)`$ is a short-ranged pair potential, $`U(n)`$ is a many-body glue function and $`n_i`$ is a “generalized coordination number” for atom $`i`$. $`n_i`$ is defined as

$$n_i=\underset{j}{}\rho \left(r_{ij}\right),$$
(2)
where $`\rho (r)`$ is an “atomic density” function. These three functions have been fitted for lead using a variety of bulk and surface properties Lim et al. (1992). The use of surface energies is particularly important for the application of this potential to model clusters. As well as clusters, this potential has been also used to model the surface reconstructions and pre-melting of low-index lead surfaces Toh et al. (1994); Landa et al. (1995), and lead nanowires Gülseren et al. (1995, 1998).

The choice of this empirical potential is motivated by the need for computational efficiency in order that global optimization is feasible for the sizes we consider here, and by our intention to compare with previous results. The use of ab initio electronic structure methods for lead is prohibitively expensive, especially as relativistic effects would need to be included to obtain reasonable results Young (1991). For example, sophisticated density functional calculations are unable to reproduce the experimental surface energies and anisotropies Vitos et al. (1998); Fei .

The functions $`U(n)`$, $`\varphi (r)`$ and $`\rho (r)`$ are displayed in Figure 1. The pair potential has a very shallow well and so most of the binding energy comes from the glue term. The glue term has been chosen to have its minimum at $`n`$=12, consistent with the designation of $`n`$ as an effective coordination number. The form of $`\rho (r)`$ is particularly significant. As $`\rho (r)`$ decreases relatively slowly with increasing $`r`$ beyond the minimum in the pair potential, next-nearest neighbours make a significant contribution to $`n`$. Therefore, the difference in surface energies between the $`\{111\}`$ and $`\{100\}`$ faces is small because, although an atom on a $`\{111\}`$ face has fewer nearest neighbours, it has more next-nearest neigbours Lim et al. (1992). However, $`\rho (r)`$ then decreases relatively rapidly to zero at the cutoff at $`r=5.503`$Å, which typically occurs between the second and third neighbour shells.

For a pair potential the pair distances are the most important quantities. For example, the lowest-energy structure of a cluster involves a balance between maximization of the number of nearest neighbours, whilst minimizing the strain energy that results from nearest-neighbour pair distances deviating from the equilibrium pair value, $`r_{\mathrm{eq}}`$ Doye et al. (1995). However, for a glue potential, such as the current one, where the main contribution to the energy is from the glue function, the most important quantities are the $`n_i`$. Indeed, one of the key factors in generating a low-energy structure is to have the $`n_i`$ values as close as possible to $`n_{\mathrm{eq}}`$, the value of $`n`$ at the minimum of $`U`$. This can potentially lead to different ordering principles than for pair potentials. Only structures that have their nearest-neighbour pair distances close to $`r_{\mathrm{eq}}`$ are generally competitive for pair potentials. However, this constraint is relaxed for glue potentials, and particularly when, as in the current case, $`\rho (r)`$ initially falls off weakly with $`r`$.

For the atoms on the surface of a cluster $`n_i<n_{\mathrm{eq}}`$. Therefore, there will be a driving force for contraction of the surface to make the pair distances for the surface atoms smaller and hence their $`n_i`$ larger. At equilibrium the surface contraction will be balanced by the increase in energy due to the resulting compression of the cluster core.

These considerations represent a particular problem for some of the usual forms for atomic clusters, such as the Mackay icosahedra and to a lesser extent decahedra. The inherent strain in these clusters results in pair distances between surface atoms that are longer than $`r_{\mathrm{eq}}`$, and so the compression needed to increase $`n_i`$ for the surface atoms is particularly large. Therefore, these traditional structural forms are expected to become increasingly disfavoured by potentials for which the pair separation depends strongly on coordination number Soler et al. (2000). Instead, novel forms that are able to obtain large $`n_i`$ values for the surface atoms, whilst not having too large an energetic penalty for compression of the cluster interior, could potentially be lowest in energy.

The global optimization of the lead clusters was performed using the basin-hopping Wales and Doye (1997); Wales and Scheraga (1999) (or Monte Carlo minimization Li and Scheraga (1987)) approach. This method has proved particularly successful in locating putative global minima for a wide variety of cluster systems Web . The optimization task becomes rapidly more difficult with increasing $`N`$ (e.g. the number of minima on the potential energy surface is thought to scale exponentially with $`N`$ Tsai and Jordan (1993); Stillinger (1999); Doye and Wales (2002a)) and so, of course, the possibility that we have not been able to obtain the true global minimum increases. However, the structural principles and trends are clear from our results.

## III Global minima for $`N160`$

The energies and point groups for the putative global minima are given in Table 1. Point files will be made available online at the Cambridge Cluster Database Web . The energies of the global minima are represented in Figure 2 in such a way that makes particularly stable clusters stand out. All clusters in the range $`9N40`$ are depicted in Figure 3 and a selection of particularly stable larger clusters in Figure 4.

First, we will look at the global minima for $`N40`$ in detail, before surveying more briefly the results for larger clusters. For $`N8`$ the clusters exhibit the same structures as typically seen for pair potentials. However, Pb9 has a somewhat unusual form that can be described as two face-sharing octahedra, and so, as with Pb10, is the beginning of an hcp cluster. Then for $`N`$=11 and 12 more open structures with three-fold axes of symmetry are preferred.

Most of the global minima for $`13N33`$ are decahedral in origin. However, the growth sequence is not straightforward. The decahedra are generally asymmetric with the quasi-fivefold axis not passing through the centre of mass. So, although the growth sequence begins by adding atoms around the equator of the 13-atom Ino decahedron Ino (1969), before this shell is completed, asymmetric decahedra with a longer quasi-fivefold axis become lower in energy, starting at $`N`$=21. Furthermore, sometimes part of the structure is distorted away from the ideal decahedral positions, e.g. at $`N`$=18, 19, 26 and 32. There are also structures with two interpenetrating (Pb20) and face-sharing (Pb24) 13-atom decahedra, the latter with two additional shared capping atoms. It is noticeable that the decahedral global minima generally have a significant proportion of surface atoms in $`\{100\}`$ type environments. This feature reflects the small energy difference between fcc $`\{111\}`$ and $`\{100\}`$ faces noted earlier. For materials that more strongly favour $`\{111\}`$ faces, the most stable decahedral form is usually a Marks decahedron Marks (1984), because this structure maximizes the proportion of $`\{111\}`$ faces, whilst retaining a relatively spherical shape. However, for lead the most stable decahedral clusters occur at $`N`$=13 and 23 (Figure 2).

The other set of ordered global minima found for $`N40`$ are the hcp clusters at $`N`$=35-38. Again, these structures are somewhat unexpected, particularly as the fcc truncated octahedron is possible at $`N`$=38, but this is further evidence of a preference for structures with a significant proportion of $`\{100\}`$-like faces.

Of the other global minima for $`N40`$, Pb17 is related to the 11-atom global minimum, but it is hard to discern any overall order for those at $`N`$=28, 29 and 34. Pb39 and Pb40 are somewhat related to the preceding hcp structures, as is clear from the viewpoint chosen for Figure 3, but again there is little order apparent on the other side of the cluster.

Beyond $`N`$=38 none of the global minima that we have located can be assigned to any of the usual structural forms. However, it would be too simplistic just to characterize the clusters as disordered. From Table 1 we can see that high symmetry structures are still present. Furthermore, if the clusters were just disordered one would expect cluster properties to evolve fairly smoothly with size. However, it is clear from Figure 2 that there are “magic number” clusters that are particularly low in energy. Unsurprisingly, these magic numbers often correspond to the high symmetry clusters.

Although most of the clusters in this size range have no discernible overall order, there are common local surface motifs that are repeatedly visible. However, only at a relatively few sizes can these local preferences be assembled into a structure that has clear overall order.

As for the smaller clusters the surface structures of the global minima reflect the particularly low value of $`\gamma `$. However, this does not lead to structures with large $`\{100\}`$ faces, but rather to many surface atoms with $`\{100\}`$-like environments. The surfaces are typically covered with a patchwork of squares and triangles. So on the flat regions of the surface it is common to see atoms surrounded by three triangles and two squares (there are two ways this can be achieved), rather than the six triangles or four squares, that are typical of $`\{111\}`$ and $`\{100\}`$ surfaces, respectively.

Pb54 is somewhat related to the Mackay icosahedron. It has an uncentred 13-atom icosahedron at its centre and a clear five-fold axis of symmetry. Along this axis it looks similar to the the $`D_{5h}`$ structure that was found by Wolf and Landman as a low-energy isomer of the 55-atom Lennard-Jones clusters Wolf and Landman (1998), and which is related to the icosahedron by a single rearrangement in which the structure is twisted around a five-fold axis. However, there is a canted arrangement of squares and triangles around the equator of the cluster.

The axial configuration of Pb54 seems to be quite a common motif, and similar patterns can be seen in one of the chosen views for $`N`$=60, 78, 95 and 148, the last based on a six-fold rather than a five-fold symmetric version of the pattern, thus making the top surface flat, rather than pyramidal. As the size of these clusters increases the pattern is, of course, extended outwards. The resulting motif is clearest for the highly symmetric, 148-atom global minimum.

Pb148 is the most prominent magic number in this size range (Figure 2). In shape, the cluster is a hexagonal barrel. Although the outer surface has a clear sixfold symmetry, this is in fact broken by the octahedron at the centre of the cluster.

As flagged in the introduction, an important aim of this paper is not only to characterize the global minima for this lead potential, but to understand how the observed structures relate back to the form of the potential. We start by examining the decahedral 13-atom global minimum, for which the $`r_{ij}`$ and $`n_i`$ values have been included in Figure 1. It is noticeable that there is a significant dispersion of nearest-neighbour distances. In fact the longest distance is 11.7% longer than the shortest, which compares to a 2.3% difference for the same structure when optimized for the Lennard-Jones potential. As expected from the discussion in Section II, the structure distorts to move the $`n_i`$ values as close to $`n_{\mathrm{eq}}`$ as possible, rather than keeping the nearest-neighbour distances near to the minimum of the pair potential. This is achieved by an expansion along the fivefold axis and a contraction of the equator of the cluster. This reduces the $`n_i`$ values for the two vertex atoms on the fivefold axis, but increases the $`n_i`$ values for the other ten surface atoms, while maintaining the $`n_i`$ value for the central atom close to $`n_{\mathrm{eq}}`$. A similar anisotropy of the pair distances has previously been noted by Lim et al. in their analysis of the lead cuboctahedra Lim et al. (1992); there is a greater contraction for the $`\{100\}`$ faces of the cuboctahedra than the $`\{111\}`$ faces because of the enhanced contribution to $`n_i`$ from next neighbours across the diagonals of the squares on the $`\{100\}`$ faces.

To understand why novel structural forms are observed for this lead potential, we take Pb55 as an example and compare the contributions to the energy from a series of competing structures (Table 2 and Figure 5). The global minimum is based on the 54-atom structure illustrated in Figure 4 but with an additional surface atom. Also possible at this size are a fcc cuboctahedron, an Ino decahedron and a Mackay icosahedron.

In Table 2 we have decomposed the pair energy into two components:

$$E_{\mathrm{pair}}=n_{\mathrm{nn}}ϵ+E_{\mathrm{strain}},$$
(3)
where $`n_{\mathrm{nn}}`$ is the number of nearest neighbours, $`ϵ`$ is the depth of the pair potential, and $`E_{\mathrm{strain}}`$ is the energetic penalty for distances that deviate from $`r_{\mathrm{eq}}`$, the distance corresponding to the minimum of the pair potential. More formally,

$$E_{\mathrm{strain}}=\underset{i<j,r_{ij}<r_0}{}ϵ\varphi (r_{ij}),$$
(4)
where $`r_0`$ is a distance criterion that distinguishes nearest from next-nearest neighbours. For all the structures we consider, there is a clear separation between these coordination shells.

As expected the pair energy only contributes a small fraction of the total energy. It is also noticeable that $`E_{\mathrm{strain}}`$ is of similar magnitude to $`E_{\mathrm{pair}}`$. This is in marked contrast to what occurs for pair potentials, where minimization of the strain energy is a key element of a structure’s stability. $`E_{\mathrm{strain}}`$ for the global minimum is particularly large.

Although the pair energy is small in magnitude, it is structure sensitive and so it can still determine the relative stabilities of structures when the energies from the glue term, $`E_{\mathrm{glue}}`$, are similar. For example, the major component of the difference in energy between the 55-atom cuboctahedron and decahedron is the greater pair energy of the cuboctahedron.

It is clear from Table 2 that the global minimum’s stability is a result of its significantly lower glue energy, which is a result of the atoms being able to achieve $`n_i`$ values that are closer to the ideal value, $`n_{\mathrm{eq}}`$. However, this lower glue energy is partially offset by the higher pair energy resulting from the distortion of the pair distances that is necessary to achieve an increase in $`n_i`$.

If we look at the atomic contributions to the energy (Table 2 and Figure 5(a)) it is clear that the lower energy results from a lower average energy for the surface atoms (particularly atoms 13–33 in Figure 5(a)), which outweighs the somewhat less favourable energies for the atoms in the interior of the cluster.

It is also interesting to understand what structural features of the global minimum lead to the larger value of $`n_i`$. We analyse this in Figure 5(b) and (c), first by looking at the cumulative contribution to $`n_i`$ from pairs with distances less than $`r`$:

$$n_i(r_{ij}<r)=\frac{1}{N}\underset{ij,r_{ij}<r}{}\rho (r_{ij}).$$
(5)
It is particularly interesting to note that $`n_i^<(r)`$ for the global minimum only becomes largest beyond 4.926Å. Therefore, although the contribution to $`n_i`$ from distances beyond this distance is small in magnitude, it is key in stabilizing the global minimum relative to the more conventional forms.

We can analyse this further by considering $`n_{\mathrm{pairs}}(r)`$, the number of pairs of atoms that are separated by less than $`r`$. It can be seen from Figure 5(c) that the number of pair distances within the radius of the cutoff distance for $`\rho `$ is significantly larger for the global minimum than for the competing structures, which in turn correlates with the larger value of $`n_i`$. However, this is only true because the cutoff is located between the second and third coordination shell. The cuboctahedron, decahedron and icosahedron all have a relatively narrow distribution of nearest-neighbour distances, which leads to a clear distinction between the second and third coordination shell. By contrast, the global minimum has a much more disperse nearest-neighbour shell and hence there is no clear distinction between a second and third neighbour shell. Instead, there is a steady increase in $`n_{\mathrm{pairs}}(r)`$ beyond the start of the second neighbour shell.

Although the above analysis has been presented for a single cluster, repeating this procedure for other sizes has confirmed the generality of the conclusions.

## IV Conclusions

We have shown by the locating the global minima for small lead clusters interacting with a many-body potential of the glue form that, contrary to the original conclusion of Lim et al., these clusters do not adopt fcc geometries for $`N160`$. Instead, they form a series of novel structures that are a consequence of the many-body character of the potential. These results naturally lead one to wonder at what size bulk-like fcc structures will develop. To help us answer this question we have plotted in Figure 6 the energies of the global minima, alongside those for a number of sequences of high-symmetry structures and the novel icosahedral forms that Hendy obtained by simulations of freezing Hendy and Hall (2001) and by construction Hendy and Doye (to be submitted).

The figure confirms that the icosahedra and decahedra are always higher in energy than the best fcc structures (except at the smaller sizes considered in the last section) and that the fcc truncated octahedra become slightly lower in energy than the cuboctahedra Lim et al. (1992); Hendy and Hall (2001). More interestingly, we can clearly see that both the global minima we have found and the new icosahedra produced by Hendy Hendy and Hall (2001); Hendy and Doye (to be submitted) are significantly lower in energy than the best fcc clusters. Thus, extrapolations between these two sizes ranges Hen  suggest that fcc clusters might well not be the lowest in energy for intermediate sizes, and hence that fcc clusters are not lowest in energy until at least $`N\mathrm{15\hspace{0.17em}000}`$ Hendy and Doye (to be submitted). This is a surprising result, for although it is not uncommon to find small ($`N<100`$) metal clusters that do not exhibit any of the usual cluster structures, it is unprecedented for this behaviour to persist up to such large sizes.

The results are also relevant to the ongoing issue of disordered metal clusters. Like recent theoretical results for gold Garzón et al. (1998); Michaelian et al. (1999); Garzón et al. (2000), cadmium, zinc Michaelian et al. (2002) and vanadium Taneda et al. (2001) many of our global minima do not fit with the fcc, hcp, decahedral and icosahedral structures that are often found for close-packed materials. However, to call these lead clusters disordered would be too strong because, although most of the clusters for $`N>40`$ have no overall stuctural order, there are common local structural preferences which at a few sizes result in highly symmetric ordered structures that are particularly stable.

For pair potentials, clusters tend to retain a lattice structure away from the magic numbers, because it is unfavourable for the pair distances to deviate significantly from the equilibrium value. For example, most small Lennard-Jones clusters can be considered to be based upon Mackay icosahedra, either with an incomplete outer layer or covered by an ordered overlayer Northby (1987). By contrast, for metal clusters the many-body character of the bonding can make it favourable to break the lattice structure away from the magic numbers so that the atoms can (in the language of the current potential) increase their effective coordination numbers ($`n_i`$’s). A structure with no overall order results. Therefore, if one only examines a few cluster sizes the presence of particularly stable ordered structures may be missed.

From our analysis of the energetics of the competing structural forms we have seen that the shape of $`\rho (r)`$, in particular the shoulder and the position of the cutoff, is key to the stability of the novel structures that we find to be lowest in energy. This dependence on the cutoff is somewhat worrying both because it is a rather long-range feature of the potential and because its position is not physically motivated, but chosen more for computational convenience. These results illustrate how sensitively cluster structure depends on the potential; the correct determination of the relative energies of competing clusters is a stringent test of any potential Doye and Wales (2002b).

Our results also help us to understand the structures exhibited by lead nanowires Gülseren et al. (1998) modelled using the same potential. Gülseren et al. were surprised at the apparent contradiction between the non-fcc character of their nanowires and the cluster results of Lim et al. Lim et al. (1992), and so suggested a number of reasons for the differing structural tendencies. However, our results show that the non-fcc character is common to both systems.

###### Acknowledgements.

J.P.K.D is grateful to the Royal Society for the award of a University Research Fellowship. S.C.H. would like to acknowledge the support of the ISAT linkages fund administered by the Royal Society of New Zealand.


format: text

text:
# Correlation effects on the Fermi surface of the two-dimensional Hubbard model

## Abstract

Effects of electron correlation on the Fermi surface is investigated for the two-dimensional Hubbard model by the quantum Monte Carlo method. At first, an infinitesimal doping from the half filling is focused on and the momentum dependent charge susceptibility $`\kappa (k)=\frac{dn(k)}{d\mu }`$ is calculated at a finite temperature. At the temperature $`T\frac{t^2}{U}`$, it shows peak structure at $`(\pm \pi /2,\pm \pi /2)`$ on the Fermi surface (line). It is consistent with the mean-field prediction of the $`d`$-wave pairing state or the staggerd flux state. This momentum dependent structure disappears at the high temperature $`TU`$. After summarizing the results of the half filling case, we also discuss the effects of the doping on the momentum dependent charge susceptibility. The anisotropic structure at half filling fades out with sufficient doping.

Effects of electron-electron interaction have been a subject of intense investigation in condensed matter physics. One of the most remarkable phenomena is the Mott transition, which is a quantum phase transition driven by the interaction. At a rational filling, the ground state of the interacting electron systems sometimes belongs to the Mott insulator. Furthermore, various kinds of instabilities toward, e.g. superconductivity or antiferromagnetism exist near the Mott transition. The electronic properties near the Mott transition can deviate from a simple Fermi liquid and it leads to interacting phenomena. In particular, the behavior of the quasiparticles near the Fermi surface is an important topic to be investigated.

Recently singular momentum dependence of the low energy excitations has been focused on in the strongly correlated electron system. In this context, the spectral weight of the two-dimensional Hubbard model has been investigated . Moreover, the deformation of the Fermi surface due to the interaction is observed in the t-J model . Also the angle resolved photoemission spectroscopy (ARPES) experiments suggest anisotropic properties in the low energy excitations .

We investigate the momentum dependent charge susceptibility $`\kappa (k)=\frac{dn(k)}{d\mu }`$ in the two-dimensional Hubbard model by the Quantum Monte Carlo (QMC) simulations. The $`\kappa (k)=\frac{dn(k)}{d\mu }`$ measures a charge fluctuation with momentum resolution. Since the infinitesimal shift of the chemical potential $`\mu `$ corresponds to an infinitesimal doping at a finite temperature, the $`\kappa (k)`$ also shows the distribution of an infinitesimally doped carriers in momentum space, $`\delta n(𝒌)\kappa (𝒌)\delta \mu .`$ Without interaction, the $`\kappa (𝒌)`$ is peaked uniformly on the Fermi surface (line). At half filling, it a square in the Brillouin zone, $`|k_x|+|k_y|=\pi `$. Even in the presence of the interaction, we call this square the Fermi surface (FS) at half filling in this paper.

The ground state of the two-dimensional Hubbard model at half filling belongs to an insulating fix point (antiferromagnetic Mott insulator). Since the Mott insulator has a finite charge gap, the charge fluctuation does not exist at zero temperature. However, it is interesting to discuss the effects of the interaction on the charge fluctuation at a finite temperature. We study correlation effects on the momentum dependent charge susceptibility the $`\kappa (𝒌)`$, focusing on the difference between ($`\pi `$,0) and ($`\pi /2`$, $`\pi /2`$) points. Many kinds of fixed points can play a role in the charge properties at a finite temperature even if they can not be the long-range ordered at the zero temperature. Therefore, comparing the QMC results with possible mean-field solutions is interesting . After summarizing the discussion at half filling, we also study the effects of finite doping on the $`\kappa (𝒌)`$ in this paper.

The two-dimensional Hubbard model on a square lattice is

$``$ $`=`$ $`t{\displaystyle \underset{j,k,\sigma }{}}(c_{j\sigma }^{}c_{k\sigma }+c_{k\sigma }^{}c_{j\sigma })`$
$`+U{\displaystyle \underset{i}{}}(n_i1/2)(n_i1/2)\mu {\displaystyle \underset{i,\sigma }{}}n_{i\sigma },`$
where $`j,k`$ denotes nearest-neighbor links, $`U`$ the Coulomb interaction, and $`t`$ hopping amplitude. Periodic boundary condition is imposed. We obtain approximation-free results by using the finite temperature auxiliary field quantum Monte Carlo (QMC) method . To calculate the momentum dependent charge susceptibility $`\kappa (𝒌)`$, we directly evaluated the $`\kappa (𝒌)`$ in the QMC simulations as follows,

$`\kappa (𝒌)={\displaystyle \frac{\beta }{L^2}}\left(n(𝒌){\displaystyle \underset{𝒌}{}}n(𝒌)n(𝒌){\displaystyle \underset{𝒌}{}}n(𝒌)\right),`$
where $`\beta `$ denotes an inverse temperature and $`L`$ the linear system size.

At half filling, the sign problem does not occur due to the particle-hole symmetry. Since the simulation is carried out in the grand canonical ensemble, information on the infinitesimally doped systems are statistically taken into account. even though the averaged system is half-filled. Away from half filling, we focus on the relatively large doping case to avoid the difficulty of the negative sign problem.

At first, let us discuss the half filling case . Figure 2 shows the $`\kappa (𝒌)`$ for $`U/t=0`$ and $`U/t=4`$ at $`T/t=0.2`$. Without interaction, the $`\kappa (𝒌)`$ is peaked on the FS and its value is constant there. On the other hand, with the interaction $`U/t=4`$, the $`\kappa (𝒌)`$ shows peak structure around ($`\pm \pi /2`$,$`\pm \pi /2`$) points. It gives an anisotropic structure on the FS and indicates that the ($`\pm \pi /2`$,$`\pm \pi /2`$) points are more sensitive to the shift of the chemical potential than ($`\pm \pi `$,0) or ($`0,\pm \pi `$) points. One may say that the holes are doped near the ($`\pm \pi /2`$,$`\pm \pi /2`$) points upon a small doping. Figure 2 shows the temperature dependence of the $`\kappa (𝒌)`$ for $`U/t=4`$ at several temperatures. The peak structure at the $`(\pm \pi /2,\pm \pi /2)`$ points is clear at $`Tt^2/U`$ and becomes ambiguous as the temperature increases and finally vanishes above $`TU`$. Here we note that, at these temperatures, the antiferromagnetic correlation length is smaller than the linear system size and the size dependence is negligible. Our results suggest that the interaction brings about the anisotropy in the charge fluctuation. The anisotropy is evident at the temperatures of order of the effective exchange interaction $`Jt^2/U`$. For the high temperatures $`T>U`$ where the Coulomb interaction $`U`$ becomes irrelevant. the anisotropic feature vanishes .

Although the ground state belongs to the antiferromagnetic Mott insulator, a different picture can appear as a crossover at an intermediate temperature $`TJ`$. Comparing the QMC results with possible three different mean-field solutions ( $`d`$-wave pairing state, staggered flux state, and Néel state), the peak structure is qualitatively reproduced by the $`d`$-wave pairing state or the staggered flux state .

Next, we discuss doping effects on $`\kappa (𝒌)`$. Due to a severe negative sign problem, it is difficult to treat systems with very small doping (near the half-filled). The result of the $`\kappa (𝒌)`$ for the doped system is shown in Figure 3 with the filling $`_𝒌n(𝒌)/L^20.64`$. The correlation effects suppress the charge fluctuation and the sum of the $`\kappa (𝒌)`$ over the Brillouin zone is substantially reduced. Further the $`\kappa (𝒌)`$ is spread over the total Brillouin zone for the interacting case ($`U/t=4`$) in comparison with the non-interacting case ($`U/t=0`$). However, as for the anisotropy on the Fermi surface (line), a clear peak structure does not observed. This is clearly different from the case of the half filling. It implies that, with sufficient doping, instabilities which is driven by the interaction may become irrelevant and the deviation from a simple Fermi liquid is negligible.

We are grateful to M. Imada, Y. Kato and S. Ryu for discussions. Y.H is supported in part by Grant-in-Aid from the Ministry of Education, Science and Culture of Japan. The computation in this work has been done in part using the facilities of the Supercomputer Center, ISSP, University of Tokyo.


text:
# Mixed configuration-interaction and many-body perturbation theory calculations of energies and oscillator strengths of J=1 odd states of neon

## I Introduction

A combined configuration-interaction (CI) many-body-perturbation-theory (MBPT) method, applied previously to divalent atoms Savukov and Johnson (2002), is extended to particle-hole states of closed-shell atoms. After derivation of CI+MBPT expressions for particle-hole states, we will apply the theory to calculations of energies and electric-dipole transition probabilities for neon.

For neon, many accurate measurements of transition rates are available, providing important tests of theory. Reciprocally, the theory might help resolve existing discrepancies among oscillator strengths ($`f`$-values) for transitions from the ground state to several excited states, for which experiments disagree. There is also a certain deficiency in existing ab-initio theories in neon, for which discrepancies among many measurements and theoretical calculations are unsettled. For example, the only other elaborate ab-initio calculations (Avgoustoglou and Beck (1998)) give an oscillator strength for the $`\left[2p_{3/2}^13s_{1/2}\right]_1`$ neon state larger than most experimental values by more than two standard deviations. Extensive calculations performed by Hibbert et al. (1993) for many transition rates along the neon isoelectronic sequence use a general configuration-interaction code (CIV3) Hibbert (1975). The calculations utilize parametric adjustments with measured fine structures, but do not completely agree with experiments in neon and have an accuracy similar to other semiempirical calculations of Seaton (1998). However, the two calculations disagree with each other for several transitions. We hope that our calculations may help to understand better the theoretical problems in neon and provide guidance for the analysis of experimental data.

Some possible applications of the present CI+MBPT method include the study of neonlike ions, Ne I – Si V, S VII, Ar IX, Ca XI, and Fe XVII that have astrophysical interest and have been included in the Opacity Project (Seaton (1987)). The transition data in neon and other noble gases are also used in plasma physics, and in studying discharges that find many industrial applications in lamps and gas lasers. The methods presented here might be also used for improving the accuracy of MBPT or for extending CI+MBPT to more complicated open-shell atoms.

The principal theoretical difficulty arises from the sensitivity of transition amplitudes to the interaction between closely spaced fine-structure components. Although it is possible to obtain energies which are reasonably precise on an absolute scale using coupled-cluster methods (Ilyabaev and Kaldor (1992)), accurate fine-structure splittings seem very difficult to obtain without semiempirical adjustments. This is why semiempirical approaches, which have fine-structure intervals carefully adjusted, are more successful in neon than are ab-initio calculations. However, as we will demonstrate in this paper, CI calculations corrected with MBPT are also capable of accurately predicting fine-structure splittings and, consequently, transition amplitudes. In this paper, we will demonstrate the excellent precision of CI plus second-order MBPT. Third-order corrections, for which numerical codes already exist Avgoustoglou et al. (1992), can also be included, providing even further improvement in accuracy.

In the following section, we use the effective Hamiltonian formalism and particle-hole single-double coupled equations to derive expressions for the second-order Hamiltonian matrix of the CI+MBPT method. In the final expressions, we present a quite accurate new MBPT that can predict energies of hole states and can describe appropriately the interactions in particle-hole atoms. The accuracy of hole energies obtained with the new MBPT will be illustrated for neon and low-$`Z`$ neon-like ions. Our CI+MBPT energies and $`f`$-values for many states of neon are tabulated. Their agreement with experiment and other theories are shown.

## II CI+MBPT method

The accuracy of the Rayleigh-Schrödinger variant of second-order MBPT given in Safronova et al. (2001) is insufficient for purpose, so that more accurate single-double equations must be used. The formulas for the correlation operator and a system of coupled equations for the correlation coefficients are given in Avgoustoglou et al. (1995); we follow the notation of Avgoustoglou et al. (1995) in the the paragraphs below. Under certain conditions, those equations can be further simplified and rewritten in the following form:

$`\left(\epsilon _b\epsilon _\alpha \right)\chi _b^\alpha `$ $`=`$ $`R_b^\alpha `$
$`\left(\epsilon _b+\epsilon _c\epsilon _\alpha \epsilon _\beta \stackrel{~}{g}_{bcbc}\right)\chi _{bc}^{\alpha \beta }`$ $`=`$ $`R_{bc}^{\alpha \beta }\stackrel{~}{g}_{bcbc}\chi _{bc}^{\alpha \beta }`$
$`\left(\epsilon _v\epsilon _r\right)\chi _v^r`$ $`=`$ $`R_v^r`$
$`\left(\epsilon _a\epsilon _v\right)\chi _v^a`$ $`=`$ $`R_v^a`$
$`\left(\epsilon _v+\epsilon _b\epsilon _r\epsilon _s\right)\chi _{vb}^{rs}`$ $`=`$ $`R_{vb}^{rs}`$
$`\left(\epsilon _v+\epsilon _b\epsilon _a\epsilon _s\right)\chi _{vb}^{as}`$ $`=`$ $`R_{vb}^{as}\text{.}`$ (1)
In the second equation of this set, the term $`\stackrel{~}{g}_{bcbc}\chi _{bc}^{\alpha \beta }`$ is subtracted from both sides of this equation to make the right-hand side small. Since large random-phase approximation (RPA) corrections in the particle-hole CI+MBPT are treated by CI, the quantities $`W_{va^{}}^{v^{}a}`$ entering this set of equations on the right-hand side in Ref. Avgoustoglou et al. (1995) are small and have been neglected here. The concern might be raised for the correlation coefficients $`\chi _{vb}^{as}`$ and $`\chi _v^r`$, which generally would have small factors $`\left(\epsilon _v+\epsilon _b\epsilon _a\epsilon _s\right)`$ or $`\left(\epsilon _v\epsilon _r\right)`$ in front. However, for the large CI model space, energies of the core-virtual orbitals $`bs`$ are well separated from the energies of the valence-hole orbitals $`av`$. The quantities $`R`$ in zero approximation can be set to:

$`R_j^i`$ $`=`$ $`\mathrm{\Delta }_{ij}`$
$`R_{bc}^{ij}\stackrel{~}{g}_{bcbc}\chi _{bc}^{ij}`$ $`=`$ $`g_{ijbc}`$
$`R_{vb}^{is}`$ $`=`$ $`g_{isvb}`$ (2)
to obtain the first-order effective Hamiltonian,

$$H_{v^{}a^{},va}^{\text{eff}}=(\epsilon _v\epsilon _a)\delta _{v^{}v}\delta _{a^{}a}+H_{v^{}aa^{}v}^{(1)}$$
(3)
and the correlation coefficients $`\chi `$. Here we define the first-order correction $`H_{v^{}aa^{}v}^{(1)}=\mathrm{\Delta }_{v^{}v}\delta _{a^{}a}+\stackrel{~}{g}_{v^{}aa^{}v}`$ to the effective Hamiltonian. For faster convergence of CI and for subtraction of the dominant monopole contributions in RPA diagrams, a $`V^{(N1)}`$ Hartree-Fock (HF) model potential for which $`\mathrm{\Delta }_{nm}=\stackrel{~}{g}_{nama}`$, $`\mathrm{\Delta }_{na}=\mathrm{\Delta }_{an}=\mathrm{\Delta }_{ab}=0`$ is introduced.

Further improvement of accuracy can be achieved through iterations. After one iteration we obtain the second-order contribution to the effective Hamiltonian,

$$H_{v^{}aa^{}v}^{(2)}=\delta R_v^v^{}\delta _{a^{}a}+\delta R_a^a^{}\delta _{v^{}v}+\delta \stackrel{~}{R}_{va}^{av^{}}\text{,}$$
(4)
where

$`\delta R_v^v^{}`$ $`=`$ $`{\displaystyle \underset{sCI}{}}{\displaystyle \frac{\mathrm{\Delta }_{v^{}s}\mathrm{\Delta }_{sv}}{\epsilon _v\epsilon _s}}`$ (5)
$`{\displaystyle \underset{scd}{}}{\displaystyle \frac{g_{cdvs}\stackrel{~}{g}_{v^{}scd}}{\epsilon _c+\epsilon _d\epsilon _v^{}\epsilon _s\stackrel{~}{g}_{cdcd}}}`$
$`+{\displaystyle \underset{stc}{}}{\displaystyle \frac{g_{v^{}cst}\stackrel{~}{g}_{stvc}}{\epsilon _v+\epsilon _c\epsilon _s\epsilon _t}}\text{,}`$
$`\delta R_a^a^{}`$ $`=`$ $`{\displaystyle \underset{scd}{}}{\displaystyle \frac{g_{cda^{}s}\stackrel{~}{g}_{ascd}}{\epsilon _c+\epsilon _d\epsilon _a\epsilon _s\stackrel{~}{g}_{cdcd}}}`$ (6)
$`+{\displaystyle \underset{scd}{}}{\displaystyle \frac{g_{cda^{}s}\stackrel{~}{g}_{ascd}}{\epsilon _a^{}+\epsilon _c\epsilon _s\epsilon _t\stackrel{~}{g}_{a^{}ca^{}c}}}\text{,}`$
$`\delta \stackrel{~}{R}_{va}^{av^{}}`$ $`=`$ $`{\displaystyle \underset{tu}{}}{\displaystyle \frac{g_{av^{}tu}\stackrel{~}{g}_{tuva^{}}}{\epsilon _v+\epsilon _a^{}\epsilon _t\epsilon _u}}`$ (7)
$`+{\displaystyle \underset{cd}{}}{\displaystyle \frac{\stackrel{~}{g}_{cdva^{}}g_{av^{}cd}}{\epsilon _c+\epsilon _d\epsilon _a\epsilon _v^{}\stackrel{~}{g}_{cdcd}}}`$
$`+{\displaystyle \underset{tCI}{}}{\displaystyle \frac{\stackrel{~}{g}_{av^{}ta^{}}\mathrm{\Delta }_{tv}}{\epsilon _v\epsilon _t}}+{\displaystyle \underset{dtCI}{}}{\displaystyle \frac{\stackrel{~}{g}_{dv^{}ta^{}}\stackrel{~}{g}_{atvd}}{\epsilon _v+\epsilon _d\epsilon _a\epsilon _t}}`$
$`{\displaystyle \underset{dt}{}}{\displaystyle \frac{\stackrel{~}{g}_{dv^{}tv}\stackrel{~}{g}_{ata^{}d}}{\epsilon _a^{}+\epsilon _d\epsilon _a\epsilon _t\stackrel{~}{g}_{a^{}da^{}d}}}`$
$`+{\displaystyle \underset{tCI}{}}{\displaystyle \frac{\mathrm{\Delta }_{v^{}t}\stackrel{~}{g}_{taa^{}v}}{\epsilon _v+\epsilon _a^{}\epsilon _t\epsilon _a}}`$
$`+{\displaystyle \underset{dt}{}}{\displaystyle \frac{\stackrel{~}{g}_{datv}\stackrel{~}{g}_{v^{}ta^{}d}}{\epsilon _a^{}+\epsilon _d\epsilon _v^{}\epsilon _t\stackrel{~}{g}_{a^{}da^{}d}}}`$
$`{\displaystyle \underset{dt}{}}{\displaystyle \frac{\stackrel{~}{g}_{data^{}}\stackrel{~}{g}_{v^{}tvd}}{\epsilon _v+\epsilon _d\epsilon _v^{}\epsilon _t\stackrel{~}{g}_{adad}}}.`$
Note that in the last equation we have extended the single-double method. The last term entering $`\delta \stackrel{~}{R}_{va}^{av^{}}`$ in the single-double formalism would normally not contain $`\stackrel{~}{g}_{adad}`$ in the denominator. However, if we do not modify this denominator, we find that in the third-order MBPT, large terms proportional to $`\stackrel{~}{g}_{adad}`$ will appear leading to a decrease in accuracy. A physical reason for modifying the denominator of this term is that the process described by this term contains two holes in the intermediate states with large interaction energy. This interaction should be treated nonperturbatively, for example, by inclusion of $`\stackrel{~}{g}_{adad}`$ into the denominator as we have done on the basis of the single-double equations in other terms. Finally, this term is almost equal to the seventh term (they are complex conjugates and their Goldstone diagrams are related by a reflection through a horizontal axis), and for convenience they are set equal in numerical calculations. The angular reduction for $`\delta \stackrel{~}{R}_{va}^{av^{}}`$ can be easily obtained using the second-order particle-hole formulas given in Ref. Safronova et al. (2001).

## III A solution of the hole-energy problem

### III.1 Breit corrections

Apart from Coulomb correlation corrections, the Breit magnetic interaction is also important in neon and the isoelectronic ions. The breakdown of various Coulomb and relativistic contributions to the energy of $`3s`$ states of neon are given in Ref. Avgoustoglou et al. (1995). Breit corrections cancel, but for higher excited states they may not. Hence, to improve the accuracy of fine-structure splittings, we include the Hartree-Fock hole Breit correction $`B_{aa}^{\text{(HF)}}`$ in our calculations,

$$B_{aa}^{\text{(HF)}}=\underset{c}{}\stackrel{~}{b}_{acca}\text{.}$$
(8)
We have checked that the first-order corrections $`B^{(1)}`$ to the energies of $`J=2`$ and $`J=1`$ states given in Table I of Ref. Avgoustoglou et al. (1995) agree with our $`B_{aa}^{\text{(HF)}}`$ contributions, 0.00062 and 0.00090 a.u., for $`2p_{3/2}`$ and $`2p_{1/2}`$ states, respectively. We omit the small frequency-dependent Breit, quantum-electrodynamic, reduced-mass, and mass-polarization corrections. Small as they are, those corrections are further reduced after subtraction for the fine-structure intervals. More careful treatment of relativistic corrections is needed in calculations of high-$`Z`$ neon-like ions.

### III.2 Calculations of hole energies for neonlike ions

Since we propose a new variant of the MBPT expansion, we would like first to demonstrate that this expansion is convergent for hole states. The theoretical hole energies shown in Table 1 have been obtained in the $`V^{(N)}`$ HF potential using Eq. (6) for $`\delta R_a^a`$ to calculate second-order corrections. The extra term in the denominator is important and is necessary for convergence of the perturbation expansion. Experimental hole energies in the National Institute of Standards and Technology (NIST) database Ref. nis  are found as the limit energies for the neon isoelectronic sequence. For neutral neon only one limit, the p3/2 energy is given in NIST nis . The 2p1/2-2p3/2 splitting 780.4269(36) cm-1 has been measured in Ref. Harth et al. (1985), and using this value we find the experimental p1/2 energy. Table 1 demonstrates the good agreement of our theoretical p3/2, p1/2 energies as well as the same fine structure interval for neon-like ions. Our fine structure interval, whose correctness is crucial for transition amplitude calculations, differs from experiment just by about 10 cm-1. Note that the HF value 187175 cm-1 for the 2p3/2 state is 8.5% higher than the experimental value 173930 cm-1, and, after adding correlation corrections, we obtain improvement by a factor of ten. For the fine structure, the HF value 1001cm-1 disagrees even more, by 28%. If we use Rayleigh-Schrödinger perturbation theory, the corrections are twice as large as our results, and the agreement with experiment does not improve.

## IV Neon energies and oscillator strengths of J=1 odd states

To test the accuracy of the CI+MBPT method, we first calculated energies of several lowest odd J=1 neon states, Table 2. The number of configurations in CI was chosen to be 52. The order of eigenstates obtained in CI+MBPT is the same as the order of the experimental levels. We abbreviate long NIST designations since the levels are uniquely specified by energy or by order.

The pure ab-initio energies differ from experimental energies by 0.0069 a.u., but after subtraction of the systematic shift (which does not make much difference in transition calculations), the agreement is at the level of 0.0001 a.u. for almost all states. Therefore, we consider the accuracy of CI+MBPT adequate for correct prediction of level mixing and oscillator strengths. For the 3s states, agreement with experiment for the fine structure interval is much better than that obtained by Avgoustoglou et al. (1995), 0.0002 versus 0.0012 a.u.; a possible explanation for this could be that single-double equations miss important corrections which we included by modifying the denominators. In Ref. Avgoustoglou et al. (1995), however, the systematic shift is small.

Finally, we present our CI+MBPT oscillator strengths in neon. After diagonalization of the second-order effective Hamiltonian, we obtain wave functions in the form of expansion coefficients in the CI space and use them to calculate oscillator strengths. Size-consistent formulas for dipole matrix elements for transitions decaying into the ground state are provided in Ref. Avgoustoglou and Beck (1998), where the absorption oscillator strength $`f`$ is also defined. We give in this table ab-initio values of the oscillator strengths $`f`$. The dominant part of the RPA corrections is included at the level of CI. Small normalization corrections are omitted.

Many experiments have disagreements in oscillator strengths far exceeding the cited errors (see Fig. 1 and Table 4): hence, for comparison, we give in Table 3 two statistical averages: the first is a weighted according to cited standard deviations and the second is an unweighted average. For the 3s levels, the experimental data compiled in Ref. Avgoustoglou and Beck (1998) and for the higher excited levels in Ref. Zhong et al. (1997) have been included in the averaging. Average values obtained here are not necessarily the most accurate, but they serve well for comparison and for a test of our probably less accurate calculated values.

A more careful analysis of experimental techniques to exclude systematic errors, which are definitely present, is necessary; our values can provide some guidance. For $`p_{3/2}^13d`$ states, since the energy separation of the two states is small, experiments give the sum of the two oscillator strengths, and the value 0.0196 rather than 0.0130 should be compared with the experimental values 0.0187 (0.0199). In this table, we also compare our theory with other semiempirical theories. Surprisingly, early calculations by Aleksandrov et al. (1983a) agree well with our calculations. A fair agreement, considering the high sensitivity of these transitions to correlation correction, is also obtained with the other theories in the table.

## V Conclusions

In this paper, we have introduced CI+MBPT theory for particle-hole states of closed-shell atoms. A difficulty that the hole energy has poor convergence is overcome with modifications of denominators in MBPT. Good precision for hole states and for particle-hole states is illustrated for many energy levels of neon. Apart from energies, our theory is tested in calculations of oscillator strengths. Agreement with averaged experimental values is achieved.

###### Acknowledgements.

The work of W. R. J. and I. M. S. was supported in part by National Science Foundation Grant No. PHY-01-39928.


text:
# Generation of Pulsed Polarization Entangled Two-Photon State via Temporal and Spectral Engineering

## I Introduction

Quantum entanglement eprb , once discussed only in the context of the foundations of quantum mechanics, is now at the heart of the rapidly developing field of quantum information science nielson . Many researchers are hoping to exploit the unique features of the entangled states in order to surpass the “classical limit” in applications such as quantum lithography litho , the quantum optical gyroscope gyro , quantum clock synchronization and positioning clock , etc.

A particularly convenient and reliable source of entangled particles is the process of type-II spontaneous parametric down-conversion (SPDC) klyshkobook ; kiess ; shih1 ; rubin ; kwiat1 , in which an incident pump photon is split into two orthogonally polarized lower energy daughter photons inside a crystal with a $`\chi ^{(2)}`$ nonlinearity. Initially, the photon pairs are entangled in energy, time, and momentum due to the energy and momentum conservation conditions that govern the process. In addition, polarization entanglement may be obtained as a result of specific local operations on the photon pair kiess ; shih1 ; rubin ; kwiat1 . An optical process such as this has the advantage that the photons, once generated, interact rather weakly with the environment, thus making it possible to maintain entanglement for relatively long periods of time.

In this paper, we propose and analyze efficient generation schemes for pulsed two-photon polarization entangled states. Through temporal or spectral engineering of the two-photon state produced in the type-II SPDC process pumped by a ultrafast laser pulse (ultrafast type-II SPDC), it is possible to generate pulsed polarization entangled states that are free of any post-selection assumptions. Pulsed polarization entangled states are an essential ingredient in many experiments in quantum optics. They are useful, for example, as building blocks for entangled states of three or more photons keller . (In general, the SPDC process only results in two-photon entanglement.) Pulsed two-photon entangled states are also useful in practical quantum cryptography systems, since the well-known arrival times permit gated detection.

We begin in section II with a discussion of the limitations of the state-of-the-art techniques for the generation of polarization-entangled photon pairs. In section III, we present a universal Bell-state synthesizer, which makes use of a novel interferometric method to temporally engineer the two-photon wavefunction. We follow up in section IV with an analysis of the spectral properties of the two-photon wavefunction produced in ultrafast type-II SPDC. We then propose a method for the efficient generation of pulsed polarization entanglement via the spectral engineering of the two-photon wavefunction. We also discuss two two-photon polarization entangled states with novel spectral characteristics: the frequency-correlated state; and the frequency-uncorrelated state. Entangled states with such spectral properties might be useful for quantum-enhanced positioning and in multi-source interference experiments.

## II Two-Photon Entanglement in ultrafast type-II SPDC

Let us briefly review one of the standard techniques for generating polarization entangled two-photon state via type-II SPDC shih1 ; rubin ; kwiat1 . A typical experimental setup is shown in Fig. 1. A type-II nonlinear crystal is pumped with a UV laser beam and the orthogonally polarized signal and idler photons, which are in the near infrared, travel collinearly with the pump. After passing through a prism sequence to remove the pump, the signal-idler photon pair is passed through a quartz delay circuit (e-o delay $`\tau `$) before the beamsplitter (BS) splits the SPDC beam into two spatial modes. A polarization analyzer and a single-photon detector are placed at each output port of the beamsplitter for polarization correlation measurements. In this case, the two quantum mechanical amplitudes in which both photons end up at the same detector are not registered since only coincidence events are considered. That is, a state post-selection has been made shih1 ; rubin ; decaro . A noncollinear type-II SPDC method developed later resolved this state post-selection problem and is usually regarded as the “standard” method for generating polarization entangled photon pairs kwiat1 .

Both the collinear and noncollinear type-II SPDC methods work very well for the generation of polarization entangled photon pairs when the UV pump laser is continuous wave (cw). In this case, however, there is no information available regarding the photons’ arrival times at the detectors. Such timing information can be quite useful in certain applications, as discussed in section I. If the UV pump has the form of an ultrafast optical pulse ($`100`$ fsec), then the photon pair arrival times can be known within a time interval on the order of the pump pulse duration. However, it has been theoretically and experimentally shown that, in general, type-II SPDC suffers the loss of quantum interference if the pump is delivered in the form of an ultrafast pulse pulsedspdctheory ; pulsedspdcexp .

We will first briefly review the theoretical treatment of the ultrafast type-II SPDC process and discuss the physical mechanism of the loss of quantum interference (or decoherence). With this understanding in hand, we then discuss in the subsequent sections two methods for eliminating this decoherence through spectral and temporal engineering of the two-photon state.

From first-order perturbation theory, the quantum state of type-II SPDC may be expressed as rubin ; pulsedspdctheory

$$|\psi =\frac{i}{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑑t|0,$$
(1)
where

$$=ϵ_0d^3\stackrel{}{r}\chi ^{(2)}E_p(z,t)E_o^{()}E_e^{()}$$

is the Hamiltonian governing the SPDC process. The pump electric field, $`E_p(z,t)`$, is considered classical and is assumed to have a Gaussian shape in the direction of propagation. The operator $`E_o^{()}`$ ($`E_e^{()}`$) is the negative frequency part of the quantized electric field of o-polarized (e-polarized) photon inside the crystal. Integrating over the length of the crystal, $`L`$, Eq. (1) can be written as

$$|\psi =C𝑑\omega _e𝑑\omega _o\mathrm{sinc}\left(\frac{\mathrm{\Delta }L}{2}\right)_p(\omega _e+\omega _o)a_e^{}(\omega _e)a_o^{}(\omega _o)|0,$$
(2)
where $`C`$ is a constant and $`\mathrm{\Delta }k_p(\omega _p)k_o(\omega _o)k_e(\omega _e)`$. The pump pulse is described by $`_p(\omega _e+\omega _o)=\mathrm{exp}\{(\omega _e+\omega _o\mathrm{\Omega }_p)^2/\sigma _p^2\}`$, where $`\sigma _p`$ and $`\mathrm{\Omega }_p`$ are the bandwidth and the central frequency of the pump pulse, respectively.

For the experimental set-up shown in Fig. 1, the electric field operators at the detectors are

$`E_1^{(+)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle 𝑑\omega ^{}\{\mathrm{cos}\theta _1e^{i\omega ^{}(t_1+\tau )}a_e(\omega ^{})\mathrm{sin}\theta _1e^{i\omega ^{}t_1}a_o(\omega ^{})\}},`$
$`E_2^{(+)}`$ $`=`$ $`{\displaystyle \frac{i}{\sqrt{2}}}{\displaystyle 𝑑\omega ^{}\{\mathrm{cos}\theta _2e^{i\omega ^{}(t_2+\tau )}a_e(\omega ^{})+\mathrm{sin}\theta _2e^{i\omega ^{}t_2}a_o(\omega ^{})\}},`$
where $`\theta _1`$ and $`\theta _2`$ are the angles of the polarization analyzers $`A_1`$ and $`A_2`$ and $`\tau `$ is the e-o delay. Here we have assumed no spectral filtering before detection. The coincidence count rate at the detectors is proportional to

$`R_c`$ $``$ $`{\displaystyle 𝑑t_1𝑑t_2|0|E_2^{(+)}E_1^{(+)}|\psi |^2}`$ (3)
$`=`$ $`{\displaystyle 𝑑t_+𝑑t_{}|𝒜(t_+,t_{})|^2},`$
where $`t_+=(t_1+t_2)/2`$ and $`t_{}=t_1t_2`$.

The two-photon amplitude $`𝒜(t_+,t_{})`$ has the form

$$𝒜(t_+,t_{})=\mathrm{cos}\theta _1\mathrm{sin}\theta _2\mathrm{\Pi }(t_+,t_{}+\tau )\mathrm{sin}\theta _1\mathrm{cos}\theta _2\mathrm{\Pi }(t_+,t_{}+\tau ),$$
(4)
where

$`\mathrm{\Pi }(t_+,t_{})=\{\begin{array}{cc}e^{i\mathrm{\Omega }_pt_+}e^{\sigma _p^2\{t_+[D_+/D]t_{}\}^2}\hfill & \text{for }0<t_{}<DL\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}`$ (7)
The parameters $`D_+`$ and $`D`$ are defined to be

$`D_+`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{u_o(\mathrm{\Omega }_o)}}+{\displaystyle \frac{1}{u_e(\mathrm{\Omega }_e)}}\right){\displaystyle \frac{1}{u_p(\mathrm{\Omega }_p)}},`$
$`D`$ $`=`$ $`{\displaystyle \frac{1}{u_o(\mathrm{\Omega }_o)}}{\displaystyle \frac{1}{u_e(\mathrm{\Omega }_e)}},`$
where, for example, $`u_o(\mathrm{\Omega }_o)`$ is the group velocity of o-polarized photon of frequency $`\mathrm{\Omega }_o`$ in the crystal.

It is clear from Eqs. (3) and (4) that the degree of quantum interference is directly related to the amount of overlap between the two two-photon wavefunctions $`\mathrm{\Pi }(t_+,t_{}+\tau )`$ and $`\mathrm{\Pi }(t_+,t_{}+\tau )`$. Figure 2 shows calculated two-photon wavefunctions for type-II SPDC in a 2 mm Beta-Barium Borate (BBO) crystal with a central pump wavelength of 400 nm. The two-photon wavefunction for the cw-pumped case is symmetric in $`t_+`$ and $`t_{}`$. As shown in Fig. 2(a), it has a rectangular shape in the $`t_{}`$ direction and extends to infinity in the $`t_+`$ direction. In the case of ultrafast type-II SPDC, the two-photon wavefunction is strongly asymmetric, see Fig. 2(b). As we shall show shortly, this is the origin of loss of quantum interference in ultrafast type-II SPDC.

Having learned the exact shape of the two-photon wavefunction $`\mathrm{\Pi }(t_+,t_{})`$, we are now in a position to study the overlap of the wavefunctions in Eq. (4). The situation is well illustrated in Fig. 3. When $`\tau =0`$, there is no overlap at all and hence no quantum interference. Increasing $`\tau `$ brings the two two-photon wavefunctions together and they begin to overlap. This increases the indistinguishability of the two wavefunctions, which represent the two-photon polarization states $`|V_1,H_2`$ and $`|H_1,V_2`$. Because the two wavefunctions are symmetric in cw-pumped type-II SPDC, perfect overlap may be achieved with the proper value of $`\tau `$, see Fig. 2(a). With an ultrafast pump, however, the wavefunctions are asymmetric and perfect overlap is not possible. The amount of relative overlap may be increased in one of two ways: (i) by decreasing the thickness of the crystal; or (ii) through the use of narrowband spectral filters to expand the two-photon wavefunction, see Fig. 2(c). Both methods increase the relative amount of overlap, which in turn increase the overall indistinguishability of the system. The drawback of these methods, obviously, is the reduced number of available entangled photon pairs.

## III Temporal Engineering of the Two-Photon State

In this section, we present a method in which complete overlap of two two-photon wavefunctions may be achieved with no change in the spectral properties of the photons. Consider the experimental setup shown in Fig. 4. A type-II BBO crystal is pumped either by a cw- or by a pulsed- pump laser. As in Ref. kwiat1 , we restrict our attention to the photons found in the intersections of the cones made by the e- and o-rays exiting the crystal note . A $`\lambda /2`$ plate rotates the polarization in one arm before the two photons are brought together at a polarizing beam splitter (PBS). Note that since the two photons always have the same polarizations when they reach the PBS, they always exit different ports. Thus, the state post-selection assumption is not necessary.

The two Feynman alternatives (reflected-reflected and transmitted-transmitted) leading to coincidence detection are shown in Fig. 5. Note that there are only two two-photon amplitudes and in both cases the e-ray (o-ray) of the crystal is always detected by $`D_1`$ ($`D_2`$). This means that, unlike the experiment discussed in the previous section, any spectral or temporal differences between the o- and e-photons provides no information at the detectors that might make it possible to distinguish between the two wavefunctions. The two biphoton wavefunctions, which represent $`|H_1|H_2`$ and $`|V_1|V_2`$, are therefore quantum mechanically indistinguishable so that the state exiting the PBS is

$$|\mathrm{\Phi }=\frac{1}{\sqrt{2}}\left(|H_1|H_2+e^{i\phi }|V_1|V_2\right),$$

where $`\phi `$ is the phase between the two terms. This phase may be adjusted by tilting the quartz plates QP1 and QP2. By setting $`\phi =0`$, the Bell state

$$|\mathrm{\Phi }^{(+)}=\frac{1}{\sqrt{2}}\left(|H_1|H_2+|V_1|V_2\right)$$

is attained. From this state, only simple linear operations are needed to transform to any of the other Bell states.

Let us now analyze this interferometer more formally. If we assume that the phase difference between the two amplitudes $`\phi `$ is set to zero, then the electric field operators that reach the detectors may be written as

$`E_1^{(+)}`$ $`=`$ $`{\displaystyle 𝑑\omega ^{}\{\mathrm{cos}\theta _1e^{i\omega ^{}(t_1+\tau )}a_{Ve}(\omega ^{})\mathrm{sin}\theta _1e^{i\omega ^{}t_1}a_{He}(\omega ^{})\}},`$
$`E_2^{(+)}`$ $`=`$ $`{\displaystyle 𝑑\omega ^{}\{\mathrm{cos}\theta _2e^{i\omega ^{}t_2}a_{Vo}(\omega ^{})\mathrm{sin}\theta _2e^{i\omega ^{}(t_2+\tau )}a_{Ho}(\omega ^{})\}},`$
where, for example, $`a_{Vo}(\omega ^{})`$ is the annihilation operator for a photon of frequency $`\omega ^{}`$ with vertical polarization which was originally created as the o-ray of the crystal.

Then the two-photon amplitude $`𝒜(t_+,t_{})`$ may be expressed as

$$𝒜(t_+,t_{})=\mathrm{cos}\theta _1\mathrm{cos}\theta _2\mathrm{\Pi }(t_++\tau /2,t_{}+\tau )+\mathrm{sin}\theta _1\mathrm{sin}\theta _2\mathrm{\Pi }(t_++\tau /2,t_{}\tau ).$$
(8)
In contrast to the amplitude represented by Eq. (4), the two-photon wavefunctions on the right-hand side of Eq. (8) overlap completely when $`\tau =0`$, regardless of their shape, see Fig. 6. Since complete overlap is possible for any shape of the two-photon wavefunction, this method should work for both cw-pumped and pulse-pumped SPDC. As such, this interferometer may be considered as a universal Bell-state synthesizer $``$ perfect quantum interference may be observed regardless of pump bandwidth, crystal thickness, or SPDC wavelength, with no need for spectral filters.

The typical peak-dip effect may be observed by varying the delay $`\tau `$, i.e.,

$$R_c=\{\begin{array}{cc}\frac{1}{2}+\frac{1}{2}\mathrm{exp}(D_+^2\sigma _p^2\tau ^2/2D^2)\hfill & \text{for }\theta _1=\theta _2=45^{}\hfill \\ \frac{1}{2}\frac{1}{2}\mathrm{exp}(D_+^2\sigma _p^2\tau ^2/2D^2)\hfill & \text{for }\theta _1=\theta _2=45^{}.\hfill \end{array}$$

Note that the thickness of the crystal $`L`$ does not appear in the above equation and all other experimental parameters, $`D_+`$, $`\sigma _p`$, and $`D`$, only affect the width of the quantum interference, not the maximum visibility. Figure 7 shows the calculated coincidence rate as a function of the delay $`\tau `$. It is easy to see that the coincidence rate $`R_c`$ is a monotonically varying function of $`\tau `$, with complete quantum interference expected at $`\tau =0`$. This insensitivity to small changes in path length difference provides for a stable source of high-quality Bell-states. This robustness is not found in two-crystal or double-pulse pump type-II SPDC schemes, where a small phase change in the interferometer results in sinusoidal modulations at the pump central wavelength. Such schemes require active phase stabilization and, therefore, are not practical kim3 . Our scheme does not suffer this disadvantage.

In addition, if $`\tau 0`$, the amplitudes do not overlap completely and a more general state is generated:

$$\rho =\epsilon \rho _{ent}+(1\epsilon )\rho _{mix},$$

where $`\rho _{ent}=|\mathrm{\Phi }^{(+)}\mathrm{\Phi }^{(+)}|`$ and $`0\epsilon 1`$. Such partially entangled states are called Werner states and have recently been the subject of experimental studies in the cw domain werner . Our scheme provides an efficient way to access a broad range of two-qubit states in both the pulsed and cw domains.

## IV Spectral Engineering of the Two-Photon State

Using the temporal engineering technique discussed in the preceding section, it is possible to generate polarization entangled states with either cw or pulsed pumping schemes. It should be pointed out that neither the asymmetry of the two-photon wavefunction nor the spectral properties of the photon pairs are affected by this technique. In this section, we are interested in removing the asymmetry shown in Fig. 2(b) through a careful choice of wavelengths of photons involved in the interaction. In this way, it is possible to obtain pulsed polarization-entangled photon pairs directly.

Recall that the source of poor wavefunction overlap in ultrafast type-II SPDC is the asymmetry introduced by a coupling between $`t_+`$ and $`t_{}`$ in the two-photon wavefunction given in Eq. (7). The approach here is to remove the coupling term so that the two-photon wavefunction becomes symmetric, as originally proposed by Keller and Rubin pulsedspdctheory . An example of such a wavefunction is shown in Fig. 8. Note that it has a rectangular shape in $`t_{}`$ (just as in cw-pumped type-II SPDC) and has a Gaussian shape in $`t_+`$ due to the Gaussian pulse envelope. It is not difficult to see that two-photon wavefunctions of this type may be completely overlapped, thus yielding full quantum interference. Therefore, pulsed polarization-entangled photon pairs may be generated using the well-known techniques described in Ref. shih1 ; kwiat1 . It is also expected that a triangular shaped correlation function, which has been observed in cw-pumping SPDC, should be observed in ultrafast type-II SPDC, as well.

Referring to Eq. (7), it is clear that if $`D_+/D`$ is made to vanish, then $`\mathrm{\Pi }(t_+,t_{})`$ becomes symmetric. Figure 9 shows the value of $`D_+`$ and $`D`$ for a type-II BBO crystal as a function of pump pulse central wavelength. Note that $`D_+=0`$ when the central wavelength of the pump pulse is 757 nm. The two-photon wavefunction for this case is shown in Fig. 8 for a crystal thickness of 2 mm and a pump bandwidth of 8 nm.

In this example, this approach has an additional advantage: the entangled photon pairs are emitted with central wavelengths of 1514 nm, which is within the standard fiber communication band. Such pulsed entangled photons at communication wavelengths may be useful for building practical quantum key distribution systems using commercially installed optical fibers. At this wavelength, of course, single photon detection is more problematic. However, the recent development of single photon counting techniques using InGaAs avalanche photodiodes may soon provide good single photon counters at this wavelength counter .

We have shown that by making the $`D_+`$ term vanish through a careful choice of the pump wavelength, an initially asymmetric two-photon wavefunction can be made symmetric. Since the temporal and spectral properties of the photons are related by simple Fourier Transforms, it should come as no surprise that $`D_+=0`$ leads to a spectrally symmetric state, as well.

The spectral properties of the two-photon state from ultrafast type-II SPDC are best illustrated in plots of the two-photon joint spectrum, which can be regarded as a probability distribution for the photon frequencies. Recall that the two-photon state $`|\psi `$ given in Eq. (2) contains the phase mismatch term $`\mathrm{sinc}(\mathrm{\Delta }L/2)`$ and the pump envelope term $`_p(\omega _e+\omega _o)`$. The joint spectrum function is simply the square modulus of the product of these two terms:

$$S(\omega _e,\omega _o)=|\text{sinc}(\mathrm{\Delta }L/2)_p(\omega _e+\omega _o)|^2.$$

The joint spectrum function for a typical ultrafast type-II SPDC configuration (type-II BBO pumped with an ultrafast UV pulse) is shown in Fig. 10(a). Here we have assumed $`L=2`$ mm, the pump wavelength is 400 nm, the pump bandwidth is 2 nm, and the SPDC wavelength is 800 nm. Note that the joint spectrum function is asymmetric: the frequency range of the idler photon is much larger than that of the signal. The joint spectrum becomes symmetric, however, when $`D_+=0`$, as shown in Fig. 10(b). Note that the signal and the idler photons have identical spectra. Thus, the photon pair is distinguishable only in the polarizations of the photons, as required for a polarization entangled state.

The spectral properties of the photon pairs represented in Fig. 10(b) display the tendency of frequency-anticorrelation in the sense that a positive detuning for one photon is accompanied by a negative detuning for the other. This effect follows from the energy conservation condition that constrains the SPDC process. In ultrafast type-II SPDC, the anticorrelation is not as strong as in the cw-pumped case due to the broad bandwidth of the pump pulse. However, the general tendency of anti-correlation, $`\omega _s=\omega _p/2\pm \omega `$ and $`\omega _i=\omega _p/2\omega `$ where $`\omega `$ is the detuning frequency, is clearly visible in Fig. 10. This frequency-anticorrelation, however, is not a required feature of the two-photon state. In ultrafast type-II SPDC, two-photon states with novel spectral characteristics $``$ the frequency-correlated state and the frequency-uncorrelated state $``$ may be generated through appropriate choices of the parameters that affect the joint spectrum.

The idea behind two-photon states with novel spectral characteristics is not new. The output characteristics of a beamsplitter and a Mach-Zehnder interferometer for frequency-correlated and frequency-uncorrelated photon pairs are theoretically studied in Ref. campos  with no discussion of how such states might be generated. Frequency-correlated states are also studied in Ref. giov , but the discussions are limited to the generation of frequency-correlated states using quasi-phase matching in a periodically poled crystal. Here, we discuss how frequency-correlated states and frequency-uncorrelated states may be generated via appropriate spectral engineering of the two-photon state in ordinary bulk crystals. Since our main focus is polarization entangled photon pairs with novel spectral characteristics, we restrict our attention to the case in which the two-photon joint spectrum is symmetric, i.e., $`D_+=0`$.

Frequency-correlated state generation via the SPDC process is somewhat counter-intuitive since conservation of energy requires the frequencies of the photon pair to sum to the frequency of the pump photon. This requirement imposes a strong constraint in cw-pumped SPDC, where the pump the pump field is monochromatic. In ultrafast SPDC, however, the pump field has a bandwidth of several terahertz frequency (several nanometers in wavelength) and so the frequencies of the photon pair need only sum to some value within the range of pump frequencies. This extra freedom makes it possible to generate the frequency-correlated state.

Examination of the joint spectrum function $`S(\omega _s,\omega _i)`$ suggests that if the bandwidth of the pump envelope $`_p(\omega _s+\omega _i)`$ is large enough and if the parameters are chosen such that the $`\text{sin}(\mathrm{\Delta }L/2)`$ may be approximated by the delta function $`\delta (\omega _s\omega _i)`$, then the two-photon state becomes

$$|\psi =C𝑑\omega _p(2\omega )a_s^{}(\omega )a_i^{}(\omega )|0,$$
(9)
where $`C`$ is a constant. Eq. (9) shows the signature of the frequency-correlated state: $`\omega _s=\omega _p/2\pm \omega `$ and $`\omega _i=\omega _p/2\pm \omega `$. Figure 11(a) shows the calculated joint spectrum function for the frequency-correlated two-photon state. In this example, the pump central wavelength was set to 757 nm in order to satisfy the condition $`D_+=0`$, i.e., to assure direct generation of polarization entangled photon pairs. This also has the effect of properly aligning the $`\text{sinc}(\mathrm{\Delta }L/2)`$ function. The only other requirement is that this function should be narrow enough that it may be approximated as $`\delta (\omega _s\omega _i)`$. This is achieved by increasing the value of $`L`$, the crystal thickness, since the width of the sinc function is inversely proportional to the crystal length. In Fig. 11(a), it is assumed that a 12 mm thick type-II BBO crystal is pumped by a 20 nm bandwidth ultrafast pulse centered at 757 nm and that, at zero detuning, the wavelengths of the SPDC photons are 1514 nm. The structure of frequency-correlation is clearly illustrated: as the signal detuning increases, the idler detuning is also increased. Such frequency-correlated states have been shown to be useful for certain quantum metrology applications clock .

By slightly modifying the condition for the generation of the frequency-correlated state, it is possible to generate the frequency-uncorrelated state. By frequency-uncorrelated state, we mean that the frequencies of the two photons are uncorrelated, in the sense that the range of the available frequencies for a particular photon is completely independent of the frequency of its conjugate. As far as the photons are concerned, this means that the spectral properties of one photon are in no way correlated with the spectral properties of the other. The joint spectrum of such a state is shown in Fig. 11(b) where the joint spectrum function is calculated for type-II SPDC in a 5-mm thick BBO crystal pumped by a 10 nm bandwidth ultrafast pump centered at 757 nm. Again, the pump pulse central wavelength is set to 757 nm to ensure symmetry in the joint spectrum. The more general case is discussed in Ref. grice , where it is also shown that such states are essential in experiments involving interference between photons from different SPDC sources.

## V Summary

The generation of polarization entangled states requires the coherent superposition of two two-photon wavefunctions. An additional requirement is that the two wavefunctions must be identical in all respects except polarization. This is not the case, in general, when the photon pairs originate in a type-II SPDC process, although only a simple delay is required when the process is pumped by a cw laser. When an ultrafast pumping scheme is employed, however, these techniques typically do not result in the complete restoration of quantum interference. Some improvement is seen with (i) thin nonlinear crystals and/or (ii) narrow spectral filters before the detectors. Unfortunately, both of these techniques result in greatly reduced count rates. By engineering the two-photon state of ultrafast type-II SPDC in time or in frequency, it is possible to recover the quantum interference without discarding any photon pairs. In addition, spectral engineering holds the promise of two-photon states with novel spectral properties which have not been available previously.

The temporal engineering of the two-photon state is accomplished through a novel interferometric technique which changes the temporal distribution of the two-photon amplitude $`𝒜(t_+,t_{})`$, which is comprised of two two-photon wavefunctions $`\mathrm{\Pi }(t_+,t_{})`$. The shapes of the wavefunctions are not altered in this process $``$ only the way in which they overlap in $`𝒜(t_+,t_{})`$. The temporally engineered source does not suffer the loss of quantum interference common in ultrafast type-II SPDC and the interferometer can be viewed as a universal Bell-state synthesizer since the quantum interference is independent of the crystal properties, the bandwidth of spectral filters, the bandwidth of the pump laser, and the wavelengths.

A different approach is taken in the spectral engineering technique. Here, the asymmetric two-photon wavefunction $`\mathrm{\Pi }(t_+,t_{})`$ from ultrafast type-II SPDC is symmetrized through careful control of the crystal and pump properties. Such spectrally engineered two-photon states exhibit complete quantum interference with no need for auxiliary interferometric techniques. The spectral engineering techniques may also be employed in the generation of two-photon states with novel spectral characteristics: the frequency-correlated state and the frequency-uncorrelated state. Unlike typical two-photon states from the SPDC process, which exhibit strong frequency-anticorrelation, the frequency-correlated state is characterized by a strong positive frequency correlation, while the frequencies of the two photons are completely uncorrelated in the frequency-uncorrelated state.

The temporal and spectral engineering techniques studied in this paper are expected play an important role in generating the pulsed entangled photon pairs essential in applications such as practical quantum cryptography, multi-photon entangled state generation, multi-photon interference experiments, etc. Since temporal engineering allows one to control the decoherence in a stable way, it also facilitates the study of the effects of decoherence in entangled multi-qubit systems. In addition, frequency-correlated discussed in this paper may be useful for certain quantum metrology applications.

###### Acknowledgements.

We would like to thank M.V. Chekhova, S.P. Kulik, M.H. Rubin, and Y. Shih for helpful discussions. This research was supported in part by the U.S. Department of Energy, Office of Basic Energy Sciences. The Oak Ridge National Laboratory is managed for the U.S. DOE by UT-Battelle, LLC, under contract No. DE-AC05-00OR22725.


text:
# Structure and diffusion time scales of disordered clusters

## Abstract

The eigenvalue spectra of the transition probability matrix for random walks traversing critically disordered clusters in three different types of percolation problems show that the random walker sees a developing Euclidean signature for short time scales as the local, full-coordination constraint is iteratively applied.

,

Physical processes such as diffusion in porous media, electrical and thermal conduction in composite materials, and information flow in random networks may be modeled using random walks on disordered clusters. For random walks on uniform Euclidean systems, their mean-square displacement goes as $`<R^2(t)>t`$, where $`t`$ is the number of steps. In contrast, constraining the random walker to traverse a critically disordered cluster such as an incipient infinite percolation cluster makes its mean-square displacement go as $`<R^2(t)>t^{2/d_w}`$, where $`d_w>2`$ is known as the walk dimension. This slowing down of the random walker, called anomalous diffusion (see, for reviews, ), is due to the delays it experiences by being trapped in or reflected from the irregular features of the cluster.

The disordered cluster where the random walker is constrained to traverse can be of many different types. For this article, we have chosen to investigate three closely-related models. These models are ordinary , fully coordinated (FC) , and iterated fully coordinated (IFC) percolation , all on the square lattice. Although the disordered clusters produced by each of these models look qualitatively different, it has been found previously that they exhibit the same static and dynamic critical properties . This work extends the analysis by investigating the entire eigenvalue spectra of their corresponding transition probability matrices . Furthermore, we have investigated the effects of culling either fully-coordinated (FC) sites or partially-coordinated (PC) sites on the structure of the disordered clusters and the characteristics of the random walk on them.

The eigenvalues $`\lambda `$ of the transition probability matrix $`W`$ directly provide the time-dependent behavior of the random walker . Eigenvalues near $`\lambda =1`$ characterize the long time behavior of the random walker while those at the region near $`\lambda =0`$ only enter the characterization of the random walker’s short time behavior.

In Fig. 1 the eigenvalue spectra for a blind ant random walker traversing a critically disordered cluster constructed from ordinary, FC, and IFC percolation are shown together with the sketch for the eigenvalue spectrum for a random walk on a uniform Euclidean lattice. At the region near $`\lambda =1`$, all three percolation models exhibit the same power-law behavior. Quantitatively, the eigenvalue density $`n(\lambda )`$ near $`\lambda =1`$ scales according to $`n(\lambda )|\mathrm{ln}\lambda |^{d_s/21}`$ . For the above three models, therefore, the long time behavior of the random walker is the same. This agrees with the previous result that the three models exhibit the same walk dimension $`d_w`$ .

The region near $`\lambda =0`$ contributes to the behavior of the random walker only for short time scales. In this region the data for IFC percolation peaks at $`\lambda =0`$, similar to the way the solid curve for the uniform Euclidean case peaks in that region . Also, as we iterate from ordinary to FC and then to IFC percolation, this closeness to the uniform Euclidean case improves. Since the random walker perceives mostly the local structure of the cluster for short time scales, the local structure of the disordered cluster formed becomes increasingly more Euclidean as we iterate from ordinary to FC and then to IFC percolation. In particular, although critically disordered clusters formed in IFC percolation exhibit a fractal nature for long time scales they closely mimic the uniform Euclidean case for short time scales.

We further investigate the nature of the critically disordered clusters by culling either the fully-coordinated (FC) sites or the partially-coordinated (PC) sites. FC sites are sites with four neighbors also occupied while PC sites have at least one neighbor that is empty. Culling PC sites would produce clusters consisting only of the FC sites and thus leaving only the part of the cluster that is highly compact.

Shown in Fig. 2 are the eigenvalue spectrum when either FC or PC sites are culled. For ordinary percolation, when PC sites are culled the eigenvalue spectrum is mostly zero. However, when FC sites are culled the eigenvalue spectrum still mimics the spectrum with no sites culled. This shows that PC sites are the dominant sites in clusters created from ordinary percolation in determining diffusion time scales, as they are also dominant structurally .

For FC percolation, culling either FC or PC sites will produce eigenvalue spectra that behave in a similar fashion. This also agrees with the previous result that neither FC nor PC sites dominate in the cluster’s structure produced by FC percolation .

For IFC percolation, the eigenvalue spectrum when FC sites are culled is mostly zero while the spectrum when PC sites are culled produces most of the features for the original spectrum. In particular, consider the features near $`\lambda =0`$ and near $`\lambda =1`$. The spectrum when PC sites are culled behaves in a similar fashion to those regions of the original spectrum where no culling is done (see Fig. 1). Thus, in IFC percolation FC sites are dominant not only structurally but also in dictating the dynamics of the Brownian motion.

In summary, although FC and IFC percolation exhibit the same fractal structure, short time diffusion characteristics become increasingly Euclidean as we iterate the full coordination requirement. Critically disordered clusters produced from IFC percolation, in particular, exhibit both a fractal structure for long time scales and a Euclidean structure for short time scales.