text:
# Bayesian curve fitting for lattice gauge theorists

## 1 INTRODUCTION

The physical consequences of a quantum field theory can be deduced from the correlation functions of the fields. Usually this is done by matching the correlation functions (or suitable combinations of them) to known behavior which contains the observables of interest. The path integrals which yield the correlation functions are often computed using Monte Carlo estimates.

To estimate the low-lying energy spectrum of the stationary states of a theory, one studies a matrix of correlation functions

$$h^{ij}(t)=\underset{t_0}{}\varphi _i(t+t_0)\varphi _j^{}(t_0),$$
(1)
where $`\varphi _i(t)=0`$ and $`\varphi _i^{}(t)`$ are operators capable of creating the states of interest. As the temporal separation $`t`$ becomes large, the elements of this correlation matrix tend (neglecting boundary conditions) to a handful of damped exponentials whose decay rates yield the lowest-lying energies of the system. To reliably determine these decay rates, the correlators $`h^{ij}(t)`$ must be calculated for large enough $`t`$ such that they are well approximated by only small numbers of exponentials. The Monte Carlo estimates of the correlation matrix $`h_D^{ij}(t)`$, where now $`D`$ denotes Monte Carlo estimates, must then be fit to known asymptotic forms $`h_M^{ij}(t)`$, where $`M`$ denotes the model fitting functions. The fit is carried out by adjusting the fit parameters, collectively denoted by $`𝒖`$, in $`h_M^{ij}(t)`$ in order to maximize the likelihood of finding the results $`h_D^{ij}(t)`$ given $`h_M^{ij}(t)`$. This amounts to minimizing

$$\chi ^2=\underset{\alpha \beta }{}\left(M_\alpha (𝒖)D_\alpha \right)C_{\alpha \beta }^1\left(M_\beta (𝒖)D_\beta \right),$$
(2)
where $`M_\alpha (𝒖)`$ refers to $`h_M^{ij}(t)`$ and $`D_\alpha `$ refers to $`h_D^{ij}(t)`$. The indices $`\alpha `$ and $`\beta `$ each include the indices of the correlation matrix as well as the instants in time at which the correlation matrix is sampled. $`C_{\alpha \beta }`$ is the covariance matrix associated with $`D_\alpha `$. One commonly-used method is to use the fitting function (which neglects boundary conditions)

$$h_M^{ij}(t)=\underset{p=0}{\overset{n1}{}}Z^{ip}Z^{jp}\mathrm{exp}(tE_p),$$
(3)
with $`E_{p+1}>E_p`$, and fit the $`n\times n`$ matrix $`h_D^{ij}(t)`$ for $`t_{\mathrm{min}}tt_{\mathrm{max}}`$. If the $`n`$ operators have been constructed using a variational approach applied to a much larger set of operators, the off-diagonal matrix elements $`h_D^{ij}(t)`$ are often statistically consistent with zero. In such cases, little is gained by including these correlation functions in the fitting, and hence, one may fit only to the $`n`$ diagonal correlators using the fitting function

$$h_M^{(i)}(t)=\underset{p=0}{\overset{n1}{}}Z^{ip}\mathrm{exp}(tE_p),$$
(4)
where now the $`Z^{ip}`$ are usually real and positive.

When using either Eq. (3) or (4) as the fitting function, it is crucial that $`t_{\mathrm{min}}`$ be chosen such that contributions from all higher energy eigenstates are negligible. Generally one is guided by the quality of fit to determine this, but the guidelines for doing so are not very clear. For example, one could start with $`t_{\mathrm{min}}=0`$, for which the fit quality is usually very bad, and increase $`t_{\mathrm{min}}`$ by one temporal lattice spacing at a time, keeping $`t_{\mathrm{max}}=t_{\mathrm{min}}+t_{\mathrm{window}}`$ for fixed $`t_{\mathrm{window}}`$, until some acceptable fit quality, such as $`Q=0.2`$, is achieved. But the acceptable value of $`Q`$ and the value of $`t_{\mathrm{window}}`$ are somewhat subjective, and it is not clear how to incorporate the uncertainties in choosing these values into the final best-fit values for the energy levels. More stringent values of $`Q_{\mathrm{acc}}`$ and $`t_{\mathrm{window}}`$ generally lead to larger values of $`t_{\mathrm{min}}`$, and hence, larger errors in the best-fit energies. The subjectivity of these errors and the loss of information associated with the discarded $`t<t_{\mathrm{min}}`$ measurements of the correlators are undesirable features of this approach.

For illustrative purposes, consider two examples, one with an excellent signal (Example 1), the other with a questionable signal (Example 2). Both examples are energies of a gauge field in the presence of a static positive and a static negative charge in 2+1-dimensional compact U(1) lattice gauge theory using an improved action on an anisotropic lattice. Typical fit results for the lowest few energies using the standard analysis with the fitting function from Eq. (4) are shown in Figures 1-3. $`S_g`$ denotes the chosen symmetry channel and $`R`$ is the number of lattice sites between the charges. Tables 1 and 2 detail the fit results for the lowest-lying energies $`E_0`$ and $`E_1`$ in each example. In Example 1, a suitable choice for $`t_{\mathrm{min}}`$ and $`t_{\mathrm{max}}`$ is easily guided by the fit quality $`Q`$ and the final results are reasonably insensitive to minor changes in the fit range. In Example 2, one might agonize over whether $`t_{\mathrm{min}}=4`$ or $`5`$ should be used. The standard analysis provides no means of incorporating such an uncertainty into the error estimate.

In order to incorporate information from small time separations, many exponentials must be retained in the model fitting functions. Standard fits using many exponentials quickly encounter severe problems: often, fitting instabilities occur and uncertainties in the estimates of the parameters of interest become very large. The source of these problems is the fact that unconstrained fits allow physically insensible or impossible values for those parameters which are not adequately constrained by the data. The solution to this problem is to introduce constraints. Bayesian statistics allows the introduction of such constraints in a natural way. The Bayesian approach is now widely used throughout many disciplines, such as economics, medicine, astrophysics, and condensed matter physics.

This talk is a report on progress being made in exploring methods of extracting physical observables from Monte Carlo estimates of correlation functions which incorporate Bayesian inference into the statistical analysis. Such procedures can fully exploit the information present in the correlation functions at small temporal separations and use this information in a way consistent with fundamental probabilistic hypotheses. The computed errors on the best-fit parameters include both statistical uncertainties and systematic errors associated with the treatment of the excited-state contaminations. Example 1 will be used to verify that such new methods work in cases where the signal is excellent. Example 2 will then be used to see if the new methods provide an improved means of extracting the information in cases when the signal is more obscured. An earlier report on such explorations is given in Ref. . Bayesian methods have been previously discussed in the context of lattice QCD in Ref. , and more recently, using maximum entropy methods in Refs. .

Since these methods employ Bayesian inference, a brief review of the Bayesian approach is given in Sec. 2. The application of the Bayesian approach to the particular problem of interest here, namely, extraction of energies from Monte Carlo determinations of correlation functions, is then detailed in Sec. 3.

## 2 BAYESIAN REGRESSION

Bayesian regression is different from the classical (frequentist) approach. In the classical approach, the data are the only source of information explicitly taken into account in constructing an estimate or test. In the Bayesian approach, an estimate or test is produced by combining the current data with information from past experience and/or theoretical constraints (prior information). In this section, Bayesian regression is described in general terms. Some introductory textbooks on Bayesian statistics are given in Refs. . For particular emphasis on Bayesian regression, see, for example, Refs. , and for a few other discussions of Bayesian methods in physics, see Refs. .

The Reverend Thomas Bayes was a Presbyterian minister born in 1702 in London, England, and who died April 17, 1761. His theory of probability was published posthumously by Richard Price in 1763 in an article entitled an Essay towards solving a problem in the doctrine of chances in the Philosophical Transactions of the Royal Society. The theorem presented by Bayes was restricted to the binomial distribution, but its generality was most likely recognized by Bayes. Bayes’ theorem was generalized beyond the binomial distribution by Laplace in 1774, most likely without knowledge of Bayes’ work. Standard (frequentist) statistical methods were developed later than Bayesian methods. Although linear regression and goodness of fit first appeared in the late 1800’s, the field blossomed during the 1920’s with the works of Fisher, Neyman, and Pearson, with a flurry of research and applications during World War II. Bayesian methods are much older, but were largely ignored or actively opposed until the 1950’s when they were championed by prominent non-statisticians, such as physicist H. Jeffreys and economist A. Bowley. Their popularity grew rapidly during the 1970’s with the advent of affordable computers, and heated debates over the superiority of either method have raged ever since.

The cornerstone of Bayesian statistics is Bayes’ theorem:

$$P(M|DI)=\frac{P(D|MI)P(M|I)}{𝑑MP(D|MI)P(M|I)},$$
(5)
where $`D`$ represents the data, $`M`$ denotes the model which one believes should describe the data, and $`I`$ represents our prior knowledge or background information about the system. The conditional probability distribution $`P(D|MI)`$ is known as the likelihood of the data, $`P(M|I)`$ is the prior probability distribution, and $`P(M|DI)`$ is called the posterior probability distribution. The term in the denominator of Eq. (5) is independent of $`M`$ and so can often simply be absorbed into an overall normalization. Bayes’ theorem essentially states that

$$\mathrm{posterior}\mathrm{likelihood}\times \mathrm{prior}.$$

Bayesian regression uses the posterior distribution for all statistical inference. Model parameters are estimated using one’s favorite statistic but evaluated using the posterior distribution. Common measures of central tendency include the mode, mean, and median, and common measures of dispersion include the variance, skewness, and kurtosis. For example, if the model is specified by a set of parameters $`𝒖`$, then the mean value and variance of a model parameter $`u_j`$ are given by

$`u_j`$ $`=`$ $`{\displaystyle 𝒅𝒖u_jP(M(𝒖)|DI)},`$ (6)
$`\mathrm{var}(u_j)`$ $`=`$ $`{\displaystyle 𝒅𝒖(u_ju_j)^2P(M(𝒖)|DI)}.`$ (7)
The mode of $`u_j`$ is found by maximizing the posterior probability, that is, maximizing the probability of the model $`M`$ being correct given data $`D`$ and subject to the background information $`I`$.

The likelihood $`P(D|MI)`$ is the probability distribution of the finding the data $`D`$ given the particular set of model parameters $`M`$ and background information $`I`$, and from the Central Limit theorem has the form

$$P(D|MI)\mathrm{exp}(\chi ^2/2),$$
(8)
where $`\chi ^2`$ is the familiar form given in Eq. (2).

The new and key feature of Bayesian statistics is the prior $`P(M|I)`$. In standard maximum likelihood fits, the Bayesian prior $`P(M|I)`$ is taken to be a constant and thus, maximizing the likelihood becomes equivalent to minimizing $`\chi ^2`$. When the number of fit parameters is small compared to the number of data points, this is a reasonable thing to do. However, when the number of fit parameters becomes comparable or larger than the number of data points, the above procedure can become unstable and/or yield parameter estimates with overly large uncertainties, as mentioned previously. The role of the Bayesian prior $`P(M|I)`$ is to filter out parameter values that are impossible and improbable, given some prior knowledge about the solutions, and hence, introduce constraints into the fits.

Prior distributions are specified based on information accumulated from past studies (earlier data), the opinions of subject-area experts, and/or from theoretical constraints. One often considers logical connections between the parameters and symmetries which the prior distribution must obey. To simplify the computational burden, practitioners often restrict the choice of prior to some familiar distributional form. In some cases, the prior can be endowed with very little informative content (such as maximum entropy which amounts to letting a monkey throw balls into bins). Strictly speaking, using the observed data to influence the choice of prior is against the Bayesian philosophy. However, it is legitimate to use a small subset of the data to guide the construction of the prior and to use the rest of the data in the regression itself. It is very important that one should avoid putting in more information than is truly known, since results do and should depend on the prior. The prior is often viewed both as a nuisance and as an opportunity.

## 3 CHOICE OF PRIOR AND TESTING

By introducing a Bayesian prior, we can now use the model function of Eq. (4) but summing over many more exponentials:

$$h_M^{(i)}(t)=\underset{p=0}{\overset{N_{\mathrm{exp}}1}{}}A_p^{(i)}\mathrm{exp}(tE_p),$$
(9)
where $`N_{\mathrm{exp}}`$ is the number of exponentials, $`N_{\mathrm{cor}}`$ is the number of correlators to be simultaneously fit, and $`i=0,\mathrm{},N_{\mathrm{cor}}1`$. To ensure positivity of the coefficients and to order the energies, define

$$A_n^{(i)}=\left(b_n^{(i)}\right)^2,E_n=E_{n1}+ϵ_n^2,$$
(10)
so that the actual parameters used in the fitting are $`u_\alpha =\{E_0,ϵ_n,b_n^{(i)}\}`$. In total, there are $`(N_{\mathrm{cor}}+1)N_{\mathrm{exp}}`$ fitting parameters. We then choose a Gaussian prior

$$P(M|I)\mathrm{exp}\left(\underset{\alpha }{}\frac{(u_\alpha \eta _\alpha )^2}{2\sigma _\alpha ^2}\right).$$
(11)
In this work, we shall not include a prior for the first $`N_{\mathrm{cor}}`$ energies, so that the sum over $`\alpha `$ in Eq. (11) includes all parameters except $`E_0,ϵ_1,ϵ_2,\mathrm{},ϵ_{N_{\mathrm{cor}}1}`$. Hence, there are $`N_{\mathrm{cor}}N_{\mathrm{exp}}`$ terms in the summation over $`\alpha `$. In using such a prior, we incorporate our expectation that the parameter corresponding to the index $`\alpha `$ most likely has a value between $`\eta _\alpha \sigma _\alpha `$ and $`\eta _\alpha +\sigma _\alpha `$.

To simplify matters, we assume that the energies of the excited-state contamination above the $`N_{\mathrm{cor}}`$-th level are most likely to be equally spaced. Also, due to the variational construction of our operators, we expect that correlator $`j`$ is dominated by the $`E_j`$ exponential. We further simplify matters by taking the same $`\eta _\alpha `$ and $`\alpha _\alpha `$ for all the $`b_j^{(j)}`$. All other coefficients are expected to be small and are given the same expected range of values. To summarize, our prior is specified using

$$\eta _\alpha =\{\begin{array}{ccc}\epsilon ,\hfill & \alpha ϵ_j,\hfill & j=N_{\mathrm{cor}}\mathrm{}N_{\mathrm{exp}}1,\hfill \\ \mathrm{\Gamma },\hfill & \alpha b_j^{(j)},\hfill & j=0\mathrm{}N_{\mathrm{cor}}1,\hfill \\ \gamma ,\hfill & \alpha b_j^{(i)},\hfill & ij,\hfill \end{array}$$
(12)
and for the Gaussian widths, we similarly use $`\sigma _\epsilon `$, $`\sigma _\mathrm{\Gamma }`$, and $`\sigma _\gamma `$. Choosing values for $`N_{\mathrm{exp}}`$, $`\epsilon `$, $`\mathrm{\Gamma }`$, $`\gamma `$, $`\sigma _\epsilon `$, $`\sigma _\mathrm{\Gamma }`$, and $`\sigma _\gamma `$ completes the specification of our prior. $`N_{\mathrm{exp}}`$ is easily chosen by increasing its value until the energies of interest stabilize. The widths must be chosen sufficiently large so that the fits are not overly constrained beyond our knowledge, but overly large widths could lead to overly large uncertainties in the best-fit parameters. Note that describing the prior is very straightforward, in contrast to describing how one subjectively chooses fitting ranges $`t_{\mathrm{min}}`$ and $`t_{\mathrm{max}}`$ in the standard analysis.

The results of applying this Bayesian curve fitting method in the two examples previously described are shown in Figs. 4 and 5. The fit values shown correspond to the mode of each quantity, that is, the value of the quantity at the maximum of the posterior distribution. The errors have been computed assuming a normal distribution about this maximum for reasons discussed below. These errors automatically combine the statistical errors in the data with the systematic uncertainties in our prior knowledge. The method works remarkably well in both examples, eliminating the need to subjectively determine $`t_{\mathrm{min}}`$. The insensitivity of the results to the parameters in the prior probability is demonstrated in Fig. 6. Notice that the error estimates increase significantly only in case (e) which is a drastic change in the prior. This insensitivity indicates that $`E_0`$ and $`E_1`$ are largely being determined by the data. Of course, the data contain little information about many of the other parameters, and hence, such parameters depend strongly on the prior. We have found that this new method works well as long as $`\epsilon `$ is not set much larger than the expected level spacings in the system and the allowed ranges of the coefficients are not set too large, as in case (e).

Currently, there remains one fly in the ointment, namely, the computation of the error estimates. We have evaluated the multi-dimensional integrals needed to calculate the variances in the fit parameters using both the Metropolis method and HMC and have encountered severely long auto-correlations. The presence of many poorly constrained parameters leads to long, narrow ridges in the probability distribution which hamper the Markov chain Monte Carlo integration. For the Metropolis method, we have commonly observed auto-correlation lengths in the tens of thousands. The problem is somewhat ameliorated with HMC, reducing the auto-correlation lengths by factors of twenty or so, but clearly, this is still not sufficient. Work continues on this topic. An alternative method of computing the errors (which deviates from a strict Bayesian philosophy) is suggested in Ref. .

## 4 CONCLUSION

Bayesian regression techniques are a viable alternative method for extracting physical observables from stochastically-determined correlation functions. By allowing the regression to take into account prior knowledge of the system from theoretical considerations and/or previous experience, Bayesian methods can make use of information obtained at short temporal separations and can incorporate systematic effects into error estimates of the observables. However, Bayesian methods are not a cure for bad data.

Much of this work arose from a collaborative effort to explore Bayesian methods. Others involved in this effort are Peter Lepage, Bryan Clark, Christine Davies, Kent Hornbostel, Paul Mackenzie, and Howard Trottier. The author wishes to acknowledge especially fruitful discussions with Peter Lepage, and also with Robert Swendsen, Julius Kuti, Jimmy Juge, David Richards, and Mike Peardon. This work was supported by the U.S. National Science Foundation under Award PHY-0099450.


format: text

text:
# Compact phases of polymers with hydrogen bonding

## Abstract

We propose an off-lattice model for a self-avoiding homopolymer chain with two different competing attractive interactions, mimicking the hydrophobic effect and the hydrogen bond formation respectively. By means of Monte Carlo simulations, we are able to trace out the complete phase diagram for different values of the relative strength of the two competing interactions. For strong enough hydrogen bonding, the ground state is a helical conformation, whereas with decreasing hydrogen bonding strength, helices get eventually destabilized at low temperature in favor of more compact conformations resembling $`\beta `$-sheets appearing in native structures of proteins. For weaker hydrogen bonding helices are not thermodynamically relevant anymore.

The collapse of a self-avoiding flexible polymer chain in a “bad” solvent has been studied for many years Gro . Following de Gennes seminal work on showing the intimate connection of polymer collapse with tricritical systems DesG , most of the theoretical effort has been concerned with the universal features of the $`\theta `$-point, the thermodynamic second order transition between the swollen and the compact phase vander , this last phase being usually regarded as a structureless liquid globule phase Gro . The possibility of more complex behavior in the compact phase has been investigated only recently, revealing the existence, at lower temperatures than the collapse gas-to-liquid transition, of a liquid-to-solid and a solid-to-solid transition Zhou .

On the other hand, protein molecules undergo similar transitions between denatured, molten globule, and native states, which are solid-like structures with a well defined three-dimensional conformation Fersht . The main driving force of protein collapse is believed to be the hydrophobic effect, which shields most of the non-polar side chains in the core of the native protein structure from water Dill . This could indeed be grossly described as a “bad” solvent effect. Yet, native structures of proteins are very peculiar, when compared to typical compact conformations of self-avoiding polymer chains. The benchmark of protein nativeness is perhaps the ubiquitous presence of highly ordered local motifs, called secondary structures, known to be stabilized by hydrogen bonding Pauling .

In this Letter, we propose a minimal off-lattice homopolymer model, where a usual two-body isotropic attractive interaction -mimicking the hydrophobic effect- is competing with a directed attractive interaction mimicking the angular dependence of hydrogen bonding Cre . We consider a chain of $`N`$ beads at positions $`\stackrel{}{r}_i`$ in the three-dimensional continuum space $`𝐑^\mathrm{𝟑}`$, with $`1iN`$. The chain constraint is enforced strictly, by keeping the distance between consecutive beads along the chain constant and unitary, $`\left|\stackrel{}{r}_i\stackrel{}{r}_{i1}\right|=1`$, for $`2iN`$, while no other constraint is considered. We model the hydrophobic effect water  by considering a pair-wise attractive square well potential with a hard wall, $`E_{hp}=_{0i1<jN}V_{hp}\left(\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|\right)`$, with:

$$V_{hp}\left(r\right)=\{\begin{array}{ccc}\mathrm{}\hfill & \mathrm{for}r\sigma \hfill & \\ 1\hfill & \mathrm{for}\sigma r\lambda \hfill & \\ 0\hfill & \mathrm{for}r\lambda \hfill & \end{array},$$
(1)
where $`\sigma `$ is the hard-core radius of each bead, and $`\lambda `$ is the range of the attractive interaction. In the following we will always consider the case $`\sigma =1`$, $`\lambda =1.5`$, as in Zhou .

In order to model hydrogen bonding, we need to break isotropy and favor a preferred direction between the two ‘hydrogen-bonded’ beads. We use the same type of directed interaction proposed by Chen and Kemp Jeff , so that the two planes, each containing one of the two hydrogen-bonded beads and its nearest neighbors along the chain, will both be preferibly orthogonal to the contact vector between them Multibody : $`E_{hb}=_{2i<j+1N}V_{hb}(\stackrel{}{r}_i\stackrel{}{r}_j,\stackrel{}{u}_i,\stackrel{}{u}_j)`$, where $`\stackrel{}{u}_i=\left(\stackrel{}{r}_{i+1}\stackrel{}{r}_i\right)\times \left(\stackrel{}{r}_i\stackrel{}{r}_{i1}\right)`$, and:

$$V_{hb}(\stackrel{}{r},\stackrel{}{u}_i,\stackrel{}{u}_j)=0.5\left(\left|\widehat{r}\widehat{u}_i\right|^m+\left|\widehat{r}\widehat{u}_j\right|^m\right)V_{hp}\left(\left|\stackrel{}{r}\right|\right),$$
(2)
where $`\widehat{}`$ denotes normalized vectors. The directionality degree of hydrogen bonding is controlled by the exponent $`m`$; a large value corresponds to a strong ‘directionality’. We have mainly studied the $`m=12`$ case, since lower values of $`m`$ do not favor protein-like secondary structures in this parametrization.

The interplay between hydrophobic collapse and hydrogen bonding is controlled by the relative strength $`\alpha `$ between the two interactions when the following total Hamiltonian is considered:

$$_\alpha =E_{hp}+\alpha E_{hb}.$$
(3)
Whenever two beads come into contact, i.e. their mutual distance falls within the well, they always gain a negative unitary energy contribution from $`E_{hp}`$. A further negative contribution could come from $`E_{hb}`$, depending on how well the hydrogen bond is formed between them, ranging from $`0`$, in the worst case, to $`\alpha `$, in the best one. Thus, from the microscopic point of view the two energy terms are cooperative. If hydrogen bonding is switched off, $`\alpha =0`$, we are back to the usual case of isotropic pairwise attraction considered in Zhou , which yields a compact groundstate with no secondary structures. In the other limit of no hydrophobic interaction, $`\alpha =\mathrm{}`$, the ground-state has already been shown to be a long straight helix, when $`m6`$ Jeff . Since the ground state differs significantly in the two limiting cases, one should actually expect a non-trivial competition between the two energy terms for intermediate values of $`\alpha `$, despite the microscopic cooperativity. This competition is induced at a global macroscopic level as a consequence of chain connectivity and excluded volume constraints Glgl . In this Letter, we will focus on its thermodynamic implications.

Our results qualitatively agree with previous work on an analogous lattice model, where hydrogen bonding was mimicked via the introduction of rotating spins BMST . We remark, nonetheless, that the extension of such results to our off-lattice model is higly non trivial, since the geometrical order implicit in the lattice structure could mask or enhance artificially secondary structure formation. As an example, while isotropic compaction of a homopolymer chain on a cubic lattice is sufficient to produce some amount of secondary structure Chan , this is not true for an off-lattice homopolymer Onu ; Greg .

Our aim is to determine, by means of Monte Carlo simulations, the density of states $`\rho \left(E\right)`$ of a polymer chain with the Hamiltonian (3), so that the partition function $`Z_N\left(T\right)`$ of a $`N`$-bead chain at reduced temperature $`T`$ can be easily reconstructed: $`Z_N\left(T\right)=_E\rho \left(E\right)\mathrm{exp}\left[E/T\right]`$. We have employed a set of standard moves currently used in simulations of polymer chain; pivot, crankshaft, and reptation moves Sokal . In order to avoid trapping in local energy minima, we have employed a novel simulation method, based on generalized ensemble techniques Jesper . The key notion, using generalized ensembles, is that a proper reweighting of temperature as a function of energy should allow the chain to escape from such energy minima Mult . The method lends itself in a natural way to be formulated within an iterative convergence scheme, and the possibility of properly employing the statistical information from more different steps of such scheme greatly increases its effectiveness Jesper .

Nevertheless, the presence of frustration provides an inherent limitation to such method, since it is based on the knowledge of local properties of the phase space, and the competition between different energy terms results in different regions of the phase space sharing the same total energy but having different local densities of states. This turns out to be the case within our model, causing a very slow convergence to equilibrium. Therefore, we have introduced a finer, two-dimensional representation of the full multi-dimensional phase-space, by identifying a conformation through both its hydrophobic energy $`E_{hp}`$ and hydrogen bond energy $`E_{hb}`$. Our simulation method can be easily adapted in order to compute the density of states $`\rho (E_{hp},E_{hb})`$ as a function of both energy terms. Details on the employed reweighting scheme will be published elsewhere long . In this way, not only convergence to equilibrium is more easily obtained, but the partition function $`Z_N(T,\alpha )=_{E_{hp},E_{hb}}\rho (E_{hp},E_{hb})\mathrm{exp}\left[\left(E_{hp}+\alpha E_{hb}\right)/T\right]`$, and hence any other relevant thermodynamic quantity, can now be reconstructed for any given value of the relative strength $`\alpha `$ between the two competing energy terms. The effectiveness of this sampling strategy shows that ther two energy terms serve as relevant order parameters.

The main result of this work, obtained for a chain with $`N=17`$ beads, is shown in Fig. $`1`$. The logarithmic density of states $`S(E_{hp},E_{hb})=\mathrm{ln}\left[\rho (E_{hp},E_{hb})\right]`$, the microcanonical entropy, is represented as a surface plot in the employed two-dimensional representation of the conformational space. Effective free energy landscapes can easily be reconstructed within the same representation:

$$F_\alpha (T,E_{hp},E_{hb})=\left(E_{hp}+\alpha E_{hb}\right)/TS(E_{hp},E_{hb}),$$
(4)
where the free energy $`F_\alpha (T,E_{hp},E_{hb})`$ is given in dimensionless units. The surface plot in Fig. $`1`$ can thus be interpreted as the opposite of the free energy landscape at infinite temperature, so that free energy valleys are seen as entropic ridges. In the phase space region with the lowest values of both energy terms, the entropy surface exhibits a rich yet regular structure, which is going to play a crucial role in determining the thermodynamic properties at low temperature. Three different ridges, separated by non-convex regions of the entropy surface corresponding to free energy barriers, can be identified.

The properties of the conformation ensembles populating such entropic ridges, or free energy valleys, can be readily identified by computing several order parameters, which measure the compactness degree and the amount of secondary structure content. Compactness is usually measured by means of the square gyration radius:

$$R_g^2=\underset{i=1}{\overset{N}{}}\left(\stackrel{}{r}_i\stackrel{}{r}_{cm}\right)^2/N,$$
(5)
where $`\stackrel{}{r}_{cm}=_{i=1}^N\stackrel{}{r}_i/N`$ is the center of mass vector. As for secondary structures, we define the helical content of a conformation as:

$$\mathrm{\Sigma }_h=\underset{ji=3}{\overset{5}{}}\left[\left(V_{i1,j1}+V_{i,j}+V_{i+1,j+1}\right)/3\right]^m,$$
(6)
and the parallel and antiparallel sheet content similarly:

$`\mathrm{\Sigma }_{ps}`$ $`=`$ $`{\displaystyle \underset{ji6}{}}\left[\left(V_{i1,j1}+V_{i,j}+V_{i+1,j+1}\right)/3\right]^m,`$ (7)
$`\mathrm{\Sigma }_{as}`$ $`=`$ $`{\displaystyle \underset{ji5}{}}\left[\left(V_{i1,j+1}+V_{i,j}+V_{i+1,j1}\right)/3\right]^m,`$ (8)
where $`V_{i,j}=V_{hb}(\stackrel{}{r}_i\stackrel{}{r}_j,\stackrel{}{u}_i,\stackrel{}{u}_j)`$, ($`0V_{i,j}1`$) measures to what extent a hydrogen bond is formed between beads $`i`$ and $`j`$. Each term in the above sums (6),(7),(8), again between $`0`$ and $`1`$, measures to what extent the $`(i,j)`$ pair can be considered the center of a local portion of a given secondary structure. Within our definitions a simple hydrophobic contact, which is not a good hydrogen bond, does not contribute to secondary structure counting.

As is shown in Fig. $`1`$, the ridges in the entropy surface are associated, with increasing number of hydrophobic contacts and decreasing number of hydrogen bonds, to helices with $`4`$ beads per turn, to helices with $`5`$ beads per turn, and to sheet-like conformations, respectively. The gyration radius decreases accordingly long , since long straight helices are extended objects. Whereas helices (a) and (b) are indeed representative of the two helical ridges, the sheet-like conformation (c) is just one among many possible different representatives. In this region we expect the occurrance of many different free energy minima, possibly giving rise to glassy behavior. Such frustration is of course not resolved within our bivariate parametrization.

We now discuss in detail the case $`\alpha =3`$. The mean square gyration radius $`<R_g^2>`$, and the specific heat per monomer $`C=T^2\left(<_\alpha ^2><_\alpha >^2\right)/N`$ are shown in Fig. $`2`$. The behavior of the secondary structure order parameters, $`<\mathrm{\Sigma }_h>`$ and $`<\mathrm{\Sigma }_{ps}>+<\mathrm{\Sigma }_{as}>`$, is shown in the inset. The curves resulting from two different independent simulations are shown; the accuracy is quite good down to temperatures as low as 0.1. The specific heat curve exhibits, with decreasing temperature, one shoulder and two higher and sharper peaks. Specific heat peaks are usually related to a phase transition, but care should of course be taken in generalizing results from such a small system (see e.g. Grass  for a detailed discussion of problems arising in finite-size scaling of $`\theta `$-collapse). For a finite-size analysis of such transitions, we show the free energy contour plots at the corresponding temperatures in Fig. $`3`$.

As signalled by the decrease of the gyration radius, the specific heat shoulder is related to the collapse of the chain from the swollen high temperature phase. The first peak close to the shoulder corresponds to a sharp increase of the gyration radius, and is related to the formation of many hydrogen bonds and to the appearance of helical structure, whereas the hydrophobic energy is almost not changed. The free energy contour plot clearly shows the existence of two competing minima, so this globule-to-helix transition is first order. The second peak is marked by a sharp decrease in the gyration radius, and is related to the breaking of helices and to the formation of sheet-like structures, as (c) in Fig. $`1`$. Sheet-like conformation are compact objects, having less hydrogen bonds but more hydrophobic contacts than helices. This last transition is again first order, as shown in the free energy contour plot. Three different minima are present, originating from the entropic ridges identified in Fig. $`1`$, but helices with $`5`$ beads per turn (b) never get efficiently populated, since they suffer competition from either side.

The very structure of the specific heat, one shoulder and then two peaks with decreasing temperature, is similar to what is found in the usual $`\alpha =0`$ case, even if the thermodynamically stable phases are completely different. It is tempting to interpret our results within the same overall framework proposed in Zhou , that is to say with decreasing temperature the chain first undergoes a gas-to-liquid collapse, then a first order liquid-to-solid transition, and finally a solid-to-solid transition which is again first order (in the absence of hydrogen bonding, the last transition is a continuous polymorphic transition Zhou ). In our model, the possibility of hydrogen bonding simply acts in selecting helices and sheets among all possible solid crystalline conformations.

In Fig. $`1`$ we have also summarized the different thermodynamic static properties of the polymer chain when $`\alpha `$ is varied. The yellow lines can be thought of as dynamical trajectories only in the infinitely slow cooling case. Actual dynamics does not take place within the effective free energy landscape (4), since kinetic barriers in the full phase space are smoothed over by the coarse-graining of our representation. This is most likely the case for the helix-to-sheet transition, where we expect the underlying energy landscape to be much more roughed.

All trajectories in Fig. $`1`$ start from a common point at infinite temperature, but then explore different regions of the phase space, according to different strengths of hydrogen bonding. Nevertheless a collapse transition, related to the shoulder in the specific heat, is common to all $`\alpha `$ values, and moreover takes place in all cases for similar values of the competing energies. At lower temperatures, hydrogen bond strength greatly affects the thermodynamic behavior. When $`\alpha =0`$, the conformational ensemble populated at low temperature does not include, as expected, the regions of the phase space with a high content of secondary structure. As in Zhou , two further transitions are present, liquid-to-solid and solid-to-solid. If $`\alpha =1`$, after the last solid-to-solid transition the sheet-like region of the phase space becomes efficiently sampled at low temperatures. If $`\alpha =3`$, as we have already seen, the chain first undergoes a transition from the liquid globule phase to the helical region and then to the sheet-like region. Both transitions are first order. Finally, if $`\alpha =4`$, hydrogen bonds are strong enough to produce a helical ground state, as in the $`\alpha =\mathrm{}`$ case Jeff .

Note that only the marginal border of the ‘green’ sheet-like region is thermodynamically relevant at low temperatures (see also the small increase of the order parameter in Fig. $`2`$). It is believed that $`\alpha `$-helices are more likely to be formed by residues with small side chain groups, whereas the loss in conformational entropy suffered by bigger side chain groups, when arrenged in helical conformation, favors the formation of $`\beta `$-sheets Rose . This general picture is consistent with our results. In fact, no side groups are present and helices are indeed entropically favored, since they sit on the top of a ridge in the entropy surface, whereas sheet-like conformations do not.

To summarize, we have introduced a simple model for an off-lattice self-avoiding polymer chain with two competing attractive inteactions, isotropic and directionalized. By means of Monte Carlo simulations we have determined the density of states of the chain within a two-dimensional representation of the phase space, and hence the phase diagram for different values of the relative strength of the two competing energies. If the directionalized interaction is strong enough, different conformational ensembles compete closely with each other at low temperature, which have peculiar proteinlike features, such as helices and sheets.

It is a pleasure to acknowledge enlightening discussions with A. Maritan, C. Rischel, K. Sneppen and G. Tiana.


text:
# The boundary between long-range and short-range critical behavior

## Abstract

We investigate phase transitions of two-dimensional Ising models with power-law interactions, using an efficient Monte Carlo algorithm. For slow decay, the transition is of the mean-field type; for fast decay, it belongs to the short-range Ising universality class. We focus on the intermediate range, where the critical exponents depend continuously on the power law. We find that the boundary with short-range critical behavior occurs for interactions depending on distance $`r`$ as $`r^{15/4}`$. This answers a long-standing controversy between mutually conflicting renormalization-group analyses.

Many of the intermolecular forces that play a central role in large areas of chemistry, physics, and biology have a long-ranged nature. Well-known examples are electrostatic interactions, polarization forces, and van der Waals forces. Remarkably, there are still considerable deficiencies in our knowledge of the critical behavior induced by these interactions, which include the prominent case of Coulombic criticality (see Ref. mef-story  for a review).

But also for the simpler case of purely attractive, long-range interactions, there exists a smothering controversy, which we aim to resolve in this work. The present understanding of critical behavior in systems with such algebraically decaying interactions is largely based on the renormalization-group (RG) calculations for the O($`n`$) model by Fisher, Ma, and Nickel FMN . Their analysis revealed that different regimes can be identified for the universal critical properties, characterized by the decay power of the interactions. In view of the small number of global parameters determining the universality class, the location of the boundaries between these regimes is of considerable interest. It is, therefore, disturbing to note that there appears to be still no consensus on the theoretical side regarding the precise location of the boundary between short-range and long-range critical behavior. Here, we address this issue by means of numerical calculations, which allow us, with minimal prior assumptions, to decide between contradictory RG scenarios.

Our approach specializes to the Ising model, $`n=1`$, in $`d`$ dimensions, described by the reduced Hamiltonian

$$/k_\mathrm{B}T=K\underset{ij}{}\frac{s_is_j}{r_{ij}^{d+\sigma }},$$
(1)
where the spins $`s_k=\pm 1`$ are labeled by the lattice site $`k`$, the sum runs over all spin pairs, and the pair coupling depends on the distance $`r_{ij}=|\stackrel{}{r}_i\stackrel{}{r}_j|`$ between the spins. According to the analysis of Fisher *et al.* FMN , universality classes are parametrized by $`\sigma `$, and the following three distinct regimes were identified: a) The classical regime. The upper critical dimension is given by $`d_u=2\sigma `$, so that mean-field-type critical behavior occurs for $`\sigma d/2`$. (b) The intermediate regime $`d/2<\sigma <2`$; here the critical exponents are *continuous functions* of $`\sigma `$. (c) The short-range regime. For $`\sigma 2`$ the universal properties are those of the model with short-range interactions, e.g., between nearest-neighbors only; thus one observes that for $`d=3`$ van der Waals interactions (decaying as $`1/r^6`$) actually lie quite close to the boundary between regimes (b) and (c).

Although the general outline of these results has been widely accepted, one part of this picture has become the scene of a debate that appears, up to now, to be unresolved. This concerns the situation close to $`\sigma =2`$. In Ref. FMN  is has been conjectured that, throughout the intermediate regime (b), the correlation-function exponent $`\eta `$ is *exactly* given by $`2\sigma `$. On the other hand, in the short-range regime (c), $`\eta `$ takes a constant (but $`d`$-dependent) value $`\eta _{\mathrm{sr}}>0`$ for all $`d<4`$, resulting in a *jump discontinuity* in $`\eta `$ at the decay power $`\sigma =2`$. While this remarkable phenomenon is not forbidden by thermodynamic arguments (which only require $`\eta 2+\sigma `$), it has attracted considerable attention over the past decades and various efforts have been undertaken to reinvestigate the corresponding RG scenario Sak ; Yama ; GT ; HN . In addition, it may be remarked that this scenario does not capture the one-dimensional case, where a phase transition is rigorously known to be absent Rue  for $`\sigma >1`$, rather than for $`\sigma >2`$ (see also Ref. invsqr  and references therein). The issue was first addressed by Sak Sak , who pointed out that higher-order terms in the RG equations considered in Ref. FMN  generate additional short-range interactions in the renormalization process, which affects the competition between the long-range and the short-range fixed points of the RG transformation. As a consequence, for $`d<4`$ the boundary between the intermediate and the short-range regime was found to shift from $`\sigma =2`$ to $`\sigma =2\stackrel{~}{\eta }`$, where the $`\epsilon `$ expansion for $`\stackrel{~}{\eta }`$ agrees to lowest orders with that of the short-range exponent $`\eta _{\mathrm{sr}}`$. In a field-theoretic approach, Honkonen and Nalimov HN  proved, to all orders in perturbation theory, the stability of the short-range fixed point for $`\sigma >2\eta _{\mathrm{sr}}`$ and of its long-range counterpart for $`\sigma <2\eta _{\mathrm{lr}}`$, where $`\eta _{\mathrm{lr}}`$ is the anomalous dimension of the field, evaluated at the *long-range* fixed point. Note that these authors also pointed out that the former result can be obtained from simple scaling arguments, but the latter *cannot*. Appealing aspects of these findings are firstly the continuous and monotonic $`\sigma `$ dependence of the correlation-function exponent (provided that $`\eta _{\mathrm{lr}}`$ and $`\eta _{\mathrm{sr}}`$ coincide at $`\sigma =2\eta _{\mathrm{sr}}`$), and secondly that the theory has now attained consistency with the exact results for the one-dimensional case.

However, the analysis of Sak  has also been the subject of criticism. Van Enter vE  pointed out an apparent inconsistency in Sak’s results. Three-dimensional *short-range* O($`n`$) models exhibit a nonzero spontaneous magnetization below the critical temperature. Van Enter vE  derived that, for $`n2`$, this broken symmetry is inconsistent with the presence of a long-range perturbation with $`\sigma <2`$. This contradicts Sak’s finding that such perturbations are irrelevant in a non-vanishing range $`2\eta _{\mathrm{sr}}<\sigma <2`$. While this inconsistency emerges for $`n2`$ and antiferromagnetic perturbations only, it may be taken as a sign that the renormalization scheme Sak ; HN  is incomplete, which could thus also affect the purely ferromagnetic $`n=1`$ case. Another conflict with Sak’s results arose from the analyses of Yamazaki Yama . A qualitative difference between his and Sak’s result for the correlation-function exponent as a function of $`\sigma `$ is the absence of a kink at $`\sigma =2\eta _{\mathrm{sr}}`$. However, we note a problem with the internal consistency of his results: The smoothly varying $`\eta (\sigma )`$ was obtained by using a long-range (i.e., $`k^\sigma `$-type, where $`k`$ is the wave vector used in the Fourier representation) instead of a short-range ($`k^2`$-type) propagator in the RG calculations, despite the fact that, in the same calculation, it was found that the long-range term in the Landau-Ginzburg-Wilson Hamiltonian is irrelevant, and only the short-range term survives. A more serious objection was raised by Gusmão and Theumann GT  (but ignored in Ref. HN ), who argued that the parameter $`2\sigma `$ (i.e., $`2\alpha `$ in the notation of Ref. HN ) is *not* a valid expansion parameter. They reconsidered the problem using renormalized perturbation theory in $`\epsilon ^{}=2\sigma d`$, and found that the long-range expansion in $`\epsilon ^{}`$ is stable with respect to short-range perturbations up to $`\sigma =2`$. This implies a restoration of the early result of Fisher *et al.* FMN  that the boundary between the intermediate and the short-range regime lies at $`\sigma =2`$. While, at least for $`d=3`$, the difference between this value and $`\sigma =2\eta _{\mathrm{sr}}`$ may be small, it indicates a fundamental dichotomy, which is also relevant for other systems with competing fixed points.

Until recently, it appeared not to be feasible to shed some light on this unsatisfactory situation by means of Monte Carlo simulations, for the following reasons. Since the problems reside in the transition region between the intermediate and the short-range regime, one expects corrections to scaling to converge only slowly with increasing system size. This necessitates the simulation of large systems. Since each spin interacts with all other spins in the system, the amount of work per Metropolis-type spin visit then becomes prohibitive. Use of a cluster Monte Carlo algorithm LB95  that employs a number of operations per spin flip that is *independent of the system size*, and moreover suppresses critical slowing down, has now enabled us to obtain highly accurate data for sufficiently large system sizes.

We present numerical results for the critical exponent $`\eta `$ and the Binder cumulant Binder  as a function of $`\sigma `$, that allow us to distinguish between the different theories. The simulated systems are defined on $`L\times L`$ lattices with periodic boundaries and sizes from $`L=4`$ to $`L=1000`$. We chose to study two-dimensional systems not only to maximize the attainable linear system size, but in particular because there the exponent $`\eta _{\mathrm{sr}}=\frac{1}{4}`$ has a much larger value than for $`d=3`$ ($`\eta _{\mathrm{sr}}0.037`$); this maximizes both the size of the disputed region $`2\eta _{\mathrm{sr}},2`$ and the magnitude of the supposed jump in $`\eta (\sigma )`$. The length of the simulations was such, that (for the largest system sizes) typically a relative uncertainty of one part in a thousand was reached for the Binder cumulant. The precise form of the pair interaction was taken as

$$\stackrel{~}{K}(|\stackrel{}{r}|)=K_{r_x\frac{1}{2}}^{r_x+\frac{1}{2}}dx_{r_y\frac{1}{2}}^{r_y+\frac{1}{2}}dy\frac{1}{(x^2+y^2)^{(d+\sigma )/2}},$$
(2)
where $`\stackrel{}{r}=(r_x,r_y)`$ denotes the difference between the integer coordinates of the two interacting spins. It is important to note that this differs from the interaction in Eq. (1) only in powers of $`r`$ that decay faster than $`r^{d\sigma }`$; these terms are irrelevant and will only affect nonuniversal terms, such as the location of the critical point. The critical exponents and boundaries of the regimes (a)–(c) will not be altered. As mentioned, we have adopted periodic boundary conditions, which are expected to speed up the approach to the thermodynamic limit. No artificial cutoff was imposed on the interaction range: Interactions take place between all periodic images.

We have carried out simulations for decay powers $`\sigma =1.2`$, $`1.4`$, $`1.6`$, $`1.75`$, $`1.85`$, $`1.95`$, $`2.0`$, $`2.05`$, $`2.25`$, $`2.5`$, $`2.75`$, and $`3.0`$. For each value of $`\sigma `$, we have determined the critical coupling $`K_c`$ on the basis of a finite-size scaling analysis of the susceptibility-like quantity $`m^2`$ and of the dimensionless amplitude ratio $`Qm^2^2/m^4`$, where $`m`$ is the magnetization density. $`Q`$ is a variant of the fourth-order cumulant introduced by Binder Binder . For Ising-like systems above the upper critical dimension, i.e., in the classical regime, it has been predicted brez-zj  to take the universal critical value $`8\pi ^2/\mathrm{\Gamma }^4(\frac{1}{4})0.456947`$, the same value as in the mean-field model LB95 ; this has been confirmed numerically in Ref. LB9697 . In the short-range regime it takes a universal value close to $`Q=0.856216\mathrm{}`$ KB ; SS . In the intermediate regime, finally, $`Q`$ has been calculated by means of a singular expansion in $`\epsilon ^{}`$, up to second order in $`\sqrt{\epsilon ^{}}`$ EL99 ; CT01 . From the magnetic susceptibility, which diverges as $`L^{\gamma /\nu }`$ at criticality, $`\eta `$ can be extracted using the scaling law $`\gamma =(2\eta )\nu `$. An important issue in these finite-size analyses is the fact that the correction-to-scaling exponents are essentially unknown, and actually depend on the borderline value of $`\sigma `$. In particular, there will be logarithmically decaying finite-size corrections at the borderline itself. Considerable effort has been exercised to encompass all such uncertainties in the quoted margins for our estimates for $`K_c`$, $`Q`$, and $`\eta `$ (see Table 1).

While the detailed least-squares analysis depended on the value of $`\sigma `$, the following general procedure was applied. An (effective) leading correction-to-scaling exponent $`y_i`$ was obtained from a fit of the amplitude ratio $`Q`$, in which also temperature-dependent terms were taken into account. Subsequently, $`\gamma /\nu `$ was extracted from an analysis of $`m^2`$ in which $`y_i`$ was kept fixed. This allowed us to include higher-order corrections. The analysis was then repeated with $`y_i`$ varied within its uncertainty margins, in order to determine the uncertainty in $`\gamma /\nu `$. Finally, the entire analysis was repeated for different minimum system sizes. For $`\sigma =1.75`$, a consistent analysis could only be obtained by allowing for logarithmic corrections to scaling, as would be expected if the crossover from long-range to short-range criticality occurs at $`\sigma =2\eta _{\mathrm{sr}}`$.

The critical coupling exhibits a very smooth but also remarkable, almost parabola-shaped dependence on $`\sigma `$, with a maximum around $`\sigma =2.17`$. This nonmonotonic variation may be due to our use of the integrated coupling (2) which tends to enhance the close-neighbor couplings, in particular for large $`\sigma `$. For $`\sigma >2`$, both $`Q`$ and $`\eta `$ are fully consistent with the known values for the two-dimensional Ising universality class. On the other hand, for $`\sigma 1.60`$, $`\eta `$ is in excellent agreement with the RG prediction $`2\sigma `$ FMN ; $`Q`$ exhibits a remarkably linear deviation from its mean-field value EL99 . However, in the disputed region $`1.75\sigma <2.0`$, where we invested most of our computing efforts, the uncertainties are still considerable; this is by itself already indicative of the occurrence of a crossover between two fixed points.

In the absence of accurate predictions for the amplitude ratio $`Q`$ (depicted in Fig. 1), this quantity does not provide much further evidence to distinguish between the different theoretical scenarios. However, there is no evidence for a jump in $`Q`$ at $`\sigma =2`$; this fails to support the Gusmão and Theumann GT  analysis. For $`\sigma =1.85`$ and $`\sigma =1.95`$, we find that $`Q`$ is compatible with the 2D Ising value, and the deviation of about 1.5 standard errors at $`\sigma =1.75`$ is still within acceptable limits.

On the other hand, the exponent $`\eta `$ offers clear evidence in order to distinguish between the different scenarios, cf. Fig. 2. While also this quantity suffers from relatively large uncertainties in the crossover region, these uncertainties are small enough to exclude the scenario of Ref. GT  which implies that a transition to short-range critical behavior occurs at $`\sigma =2`$ and that $`\eta =2\sigma `$ for $`\sigma <2`$. Indeed, instead of the 2D Ising value $`\eta =\frac{1}{4}`$, this exponent should then take the values $`0.10`$ at $`\sigma =1.85`$ and $`0.05`$ at $`\sigma =1.95`$. However, at these two values of $`\sigma `$, our results for $`\eta `$ lie within one standard deviation from the 2D Ising value, and three respectively almost five standard deviations above $`2\sigma `$.

In summary, we have demonstrated that both the correlation-function exponent $`\eta `$ and the fourth-order amplitude ratio $`Q`$ take their (universal) short-range Ising values for $`\sigma >2\eta _{\mathrm{sr}}`$. For $`\sigma >2`$ this could be shown with high numerical accuracy, but also for $`2\eta _{\mathrm{sr}}<\sigma <2`$ the results are precise enough to convincingly exclude a crossover from long-range to short-range critical behavior at $`\sigma =2`$. Instead, they support a crossover at $`\sigma =2\eta _{\mathrm{sr}}`$, in agreement with Refs. Sak ; HN . This provides a clear answer to the dilemma posed by conflicting renormalization-group scenarios.

Furthermore, our findings for $`\eta `$ are in excellent agreement with $`\eta =2\sigma `$ in the intermediate range $`d/2<\sigma <2\eta _{\mathrm{sr}}`$, supporting the conjecture that all higher-order contributions vanish in the $`\epsilon ^{}`$ expansion for $`\eta `$ FMN  and in contrast with the scenario proposed in Ref. Yama .

The amplitude ratio $`Q`$ depends approximately linearly on $`\sigma `$ for $`d/2<\sigma <2\eta _{sr}`$. Remarkably, the most prominent deviations from linearity occur near $`\sigma =2\eta _{sr}`$, while the $`\epsilon ^{}`$ expansion predicts a square-root-like singularity at the opposite end of the intermediate range EL99 .

Finally, we quote a plausible explanation for the breakdown of the description of Ref. GT  in the range $`2\eta _{\mathrm{sr}}<\sigma <2`$ Sakpc . If the validity of the $`\epsilon ^{}`$-expansion is not uniform in $`\sigma `$, the region of convergence may shrink to zero when $`\sigma `$ approaches 2. Then one cannot fix $`\epsilon ^{}`$ and take the limit $`\sigma 2`$, as is done in Ref. GT . This would also explain why the long-range ($`\epsilon ^{}`$) expansion does not provide any indication of a change to short-range behavior at $`\sigma =2\eta _{\mathrm{sr}}`$.

###### Acknowledgements.

We are indebted to Professor J. Sak for providing his valuable comments on the status of the theoretical analyses of the border between the regions of the intermediate- and the short-range interactions. This research is supported in part by the Dutch FOM foundation (“Stichting voor Fundamenteel Onderzoek der Materie”) which is financially supported by the NWO (“Nederlandse Organisatie voor Wetenschappelijk Onderzoek”).


text:
# X-RAY SPECTROSCOPY OF CATACLYSMIC VARIABLES

## MOTIVATION

Cataclysmic variables (CVs) are interacting binaries in which the accreting object (the primary) is a white dwarf (see Warner 1995 for a comprehensive review), unlike in X-ray binaries in which the accreting star is a neutron star or a black hole. CVs are numerous, as the following line of reasoning illustrates. In the Galaxy, there are about 35 classical novae per year (Shafter 1997). Since classical nova eruptions require the primary to have accreted $`1\times 10^4`$ M of hydrogen-rich fuel, we infer the total accretion rate in CVs throughout our Galaxy to be about 3.5$`\times 10^3`$ M yr-1. If the mean accretion rate per CV is 3.5$`\times 10^{10}`$ M yr-1, then the total number of CVs in our Galaxy is about 10 million. Although such calculations are fraught with uncertainties, particularly in the mean accretion rate, and hence the total Galactic population may differ from this estimate by a factor of 10, there is no doubt that CVs are commonplace in our Galaxy. If we further assume a mean X-ray luminosity of CVs to be about 4$`\times 10^{31}`$ ergs s-1 (Mukai & Shiokawa 1993), the Galactic CVs might contribute $`4\times 10^{38}`$ ergs s-1 of unresolved X-ray emission.

Moreover, the relatively high space density also means that the nearest members of the class are less than 100 pc from the Earth, close enough for detailed studies of the component stars and the accretion processes in many wavelength ranges: thus the CVs have proved to be an excellent laboratory for the study of accretion processes. Moreover, the relative shallowness of the potential well of a white dwarf means that the accretion disks of CVs are too cool to emit X-rays. Rather, CVs generate X-ray photons when material accretes onto the white dwarf. Because of the modest luminosities, these X-rays do not strongly influence the appearances of the accretion disks or alter the secular evolution of the systems. In many instances, therefore, we are able to study the accretion disk in the optical and in the UV, and the accretion process from the disk to the white dwarf in the X-rays, independent of each other. Such decoupling in CVs can be used to our advantage in understanding these systems, which can then be applied to a wider range of accretion powered X-ray sources.

## CV SUBCLASSES

Historically, optical astronomers have used the eruptive properties of CVs to classify them into dwarf novae, old novae, and nova-like systems. These are based on whether they have been seen to exhibit either of the two distinct types of outbursts known in CVs: dwarf nova outbursts, which are sudden brightenings of the accretion disk, and classical nova eruptions, which are thermonuclear runaways on the surface of the white dwarf. It is believed that dwarf nova outbursts are caused by an instability of the disk at a low accretion rate; systems that accrete at a high rate, and systems without disks, do not show dwarf nova outbursts. It is likely that all CVs undergo classical nova eruptions once a sufficient amount of fresh fuel is accreted. This results in a typical recurrence time of $`>`$10,000 years. Because high accretion rate systems with massive white dwarfs have the most frequent runaways, they therefore comprise the majority of old novae (=systems seen to undergo a classical nova eruption).

The X-ray properties of a CV is strongly influenced by the magnetic field of the primary, leading to an orthogonal scheme of CV subclassification. In Polars, or AM Her type systems, the primary is so magnetic (10 to $`>`$200 MG) that no accretion disk can form, optical/infrared cyclotron emission is strongly detected, and the white dwarf spins synchronously with the binary orbit. The intermediate polars (IPs, or DQ Her type systems) are generally less magnetic, partial accretion disks are probably present, cyclotron emission is weak or undetected, and the white dwarf spin period is shorter than the orbital period. In both these subclasses, accretion proceeds quasi-radially along field lines onto the magnetic polar regions, and forms accretion “columns” or “curtains.” There a strong shock forms, heating the plasma to X-ray temperatures (kT$`>`$ 10 keV). The white dwarf surface below the shock is heated to kT$``$20 eV by these X-rays and by direct injection of kinetic energy by dense blobs. The IPs are the strongest source of 2–10 keV X-rays among CVs, with estimated luminosities in excess of 1033 ergs s-1, while the Polars are generally strong ($`10^{32}`$ ergs s-1) soft X-ray sources.

Although less conspicuous to X-ray observers, the majority of CVs are the non-magnetic, low accretion rate dwarf novae. Dwarf novae in quiescence are relatively weak (a few times $`10^{31}`$ ergs s-1) sources of medium energy X-rays. These X-rays are generated in a compact equatorial boundary layer between the disk and the white dwarf surface (see, e.g., Patterson and Raymond 1985a; Mukai et al 1997). There is no doubt that the material in the accretion disk just beyond the boundary layer has sufficient kinetic energy stored in its Keplerian motion to produce these X-rays. However, it is not clear exactly how the gas is heated to X-ray emitting temperatures. The very fact that medium energy X-rays are observed require that the emitting region is optically thin and the shocks are sufficiently strong (Pringle & Savonije 1979). Numerical solutions obtained by Narayan & Popham (1993) are consistent with such a picture; however, such studies to date suffer from a limited dimensionality and a dependence on the unknown viscosity parameter, $`\alpha `$. At higher accretion rates (dwarf novae in outburst, nova-like systems, and old novae), the boundary layer should become optically thick and radiate predominantly in the soft X-rays/EUV range (Patterson and Raymond 1985b). Although some systems do show the expected soft component, observations show that the high accretion rate systems also emit medium energy X-rays, perhaps from the optically thin surface of the boundary layer.

The magnetic CVs are better studied in X-rays than the non-magnetic systems, because they are brighter, and because the theory is more straightforward. I, too, will concentrate mostly on magnetic CVs in the rest of this paper, with occasional references to non-magnetic systems.

## X-RAY EMISSION FROM MAGNETIC CVS

### The Aizu Model

The Aizu model was originally developed for non-magnetic white dwarfs with a spherically symmetric accretion flow (Aizu 1973). However, its current applications are for magnetic CVs, which locally resemble the spherical case. In Aizu’s original paper, numerical solutions are given in the limit of small shock height, $`h`$. There are more modern treatments, either with different approximations, or with refinements: of the latter, the addition of a second cooling function (Wu et al 1994) representing cyclotron cooling is the most important. I will use the term “Aizu model” inclusively, including all these variations.

In the Aizu model for magnetic CVs, material follows the field lines to the surface, forming a strong shock. If the shock height is small, the shock temperature $`kT_s`$ is determined by the white dwarf mass $`M_{wd}`$. For accreting material of solar abundances, $`kT_s`$21 keV for $`M_{wd}=0.6M_{}`$ and $`kT_s`$57 keV for $`M_{wd}=1.0M_{}`$. The shock height in turn is determined by equating the settling time of the post-shock flow with its cooling time. Aizu assumed optically thin thermal bremsstrahlung to be the sole mechanism for cooling. Although Compton cooling probably is not relevant to most magnetic CVs (Frank et al 1983), we will return to the subject of an additional cooling mechanism in the next section.

In the Aizu model, the emergent X-ray spectrum is the sum of emissions from plasmas at different temperatures, from $`kT_s`$ to the white dwarf photospheric temperature. As shown in Fig.1, the bulk of the post-shock region (note the logarithmic scale of the X axis) is hot and low-density. At an accretion rate per unit area of 1 g cm-2s-1 onto a 0.6$`M_{}`$ primary, the shock height is about 0.053 white dwarf radii (Rwd, or 460 km), and the density just below the shock is about $`1.0\times 10^{16}`$ cm-3. Thus, it is reasonable to ignore optical depth effects in modelling the X-ray continuum emission from just below the shock, as long as the accretion column radius is no more than, say, 500 km. Although the local accretion rate and the accretion column radius are both poorly known, these numbers can be considered typical: They result in a reasonable total luminosity, which is somewhat better constrained to be typically a few times 1032 ergs s-1 for Polars and a few times 1033 ergs s-1 for IPs.

### Simple X-ray spectral models

For low signal-to-noise, low spectral resolution spectra obtained with early X-ray satellite, however, a single temperature Bremsstrahlung model with a simple absorber was sufficient to fit the data. This was mostly the case in the 1970s, except that the need to include a reflection component was already seen in the HEAO-1 data on the prototype Polar, AM Her (Rothschild et al 1981). In magnetic CVs, a solid angle of $`2\pi `$ steradians is subtended by the white dwarf surface, as seen from the post-shock plasma. Therefore, roughly half the original X-ray photons strike the surface, resulting in electron scattering, photoelectric absorption, and fluorescent line emission whose relative efficiencies are a strong function of the photon energy. This causes a characteristic bump in the continuum above 10 keV, as seen with HEAO-1, and the Fe K$`\alpha `$ line at 6.4 keV. Indeed, emission line features in the 6–7 keV range are ubiquitous among magnetic CVs. However, in the proportional counter data, it was difficult to separate the contributions of the 6.4 keV reflection feature from those of thermal plasma at 6.7 keV (He-like Fe) and at 6.97 keV (H-like Fe).

Norton and Watson (1989) have presented their study of the spin modulation of IPs using EXOSAT data. They find that the spin modulation is stronger at lower energies, generally validating the “accretion curtain” model, in which the spin modulation is caused by photoelectric absorption in the accretion flow towards the magnetic poles. However, they also find that the energy dependence was not as strong as would be expected from pure photoelectric absorption. This suggests partial covering absorber: unfortunately, the statistics of the spin phase resolved spectra of magnetic CVs are usually not good enough to quantify this in detail, and one cannot distinguish spatial and temporal partial covering when analyzing phase-averaged spectra.

The Ginga LAC provided good data on magnetic CVs, often extending up to 30 keV. Using these, Ishida (1991) found typical Bremsstrahlung temperatures of $``$15 keV for Polars and $``$30 keV for IPs.

### Thermal lines from the post-shock plasma

Ishida et al (1994) observed the unusual IP, EX Hya, with the ASCA SIS. In addition to the Fe K$`\alpha `$ lines, they discovered K$`\alpha `$ lines of Ne, Mg, Si, S, and Ar for the first time, directly establishing the presence of multi-temperature plasma in this system (see Fig.2). Moreover, something about the temperature distribution can be inferred from the line ratios: in fact, Fujimoto and Ishida (1997) have converted these into the shock temperature, and estimated the primary of EX Hya to be 0.48 M. They also infer the temperature at “the bottom of the shock” to be $``$0.65 keV, yet do not discuss the fact that the true photospheric temperature must be much lower in this system (the EUVE spectrum extends to $``$180Å, without showing any signs of the soft, blackbody-like component; see Hurwitz et al 1997). The Fe K$`\alpha `$ lines in EX Hya are resolved into 3 components, although the 6.4 keV line is rather weak in this system. Finally, the spectrum of EX Hya around 1 keV is dominated by the Fe L features, which is a characteristic of thermal plasma with kT$``$1 keV.

The strength of the Fe L complex in EX Hya (also seen in the X-ray spectrum of stellar coronae) points out an inadequacy of the Aizu model: the assumption of the Bremsstrahlung cooling breaks down for optically thin, kT$``$1 keV plasma. Such a plasma cools, instead, predominantly by emitting Fe L photons (see, e.g., Gehrels and Williams 1993). Thus, the Aizu model overestimates the cooling time appropriate for the kT$``$1 keV region, and overestimates its height. More importantly, integrating plasma emission models over the emission measure distribution calculated using the Aizu model would result in a significant overprediction of the Fe L feature in magnetic CVs.

However, this turns out to be the least of problems, when applied to most magnetic CVs other than EX Hya. Among ASCA archival spectra of magnetic CVs, strong Fe L features are the exception, not the rule. The case for this is particularly strong for the typical IP, V1223 Sgr (Fig.3). Even though the continuum is strongly detected down to $``$0.4 keV, there is no sign of the Fe L features. This, therefore, cannot be explained as due to photoelectric absorption (which would remove continuum and line photons equally; any photoelectric absorption which removes the Fe L photons will not let $``$0.4 keV continuum photons through), and points to a real problem in applying the standard spectral model to real magnetic CVs. Of the 20 or so magnetic CVs observed with ASCA, only EX Hya, YY Dra, V884 Her, and AE Aqr show strong Fe L features. Similarly, many magnetic CVs show only weak evidence for K$`\alpha `$ lines of elements other than Fe. Several lines of arguments suggest that EX Hya has an exceptionally low accretion rate for an IP, and hence a relatively tall, low density shock; this may be somehow connected to the ease with which low energy lines have been detected in this system.

A similar situation exists for non-magnetic CVs. They, too, should exhibit multi-temperature plasma emission. There are several cases where strong Fe K and L features are seen (e.g., V603 Aql), which is a direct evidence for at least two temperature components; there are other cases where the Fe L features are absent (e.g., HT Cas). For both non-magnetic and magnetic CVs, it is highly unlikely that an X-ray emitting plasma can cool to white dwarf photospheric temperatures without emitting Fe L photons. Rather, we may try to modify the standard model by considering a geometry in which the lower temperature component is selectively hidden, while allowing the higher temperature component (which also emits low energy continuum photons) to escape.

The study of Fe K$`\alpha `$ lines have turned out to be more fruitful: the fluorescent component can often be separated from the thermal components, and evidence has been found for a significant broadening of the thermal lines in some systems, while they appear narrow in others (Hellier et al 1998). For example, the 6.4 keV line in the IP AO Psc appears narrow, while the thermal lines are blended into a broad bump; in contrast, V1223 Sgr (which is otherwise very similar to AO Psc) appears to have three narrow Fe K$`\alpha `$ lines in the ASCA data. Having considered, and discarded, several possibilities such as bulk Doppler motion, Hellier et al consider Compton scattering to be the most likely broadening mechanism. If the opacity for the thermal line core is high enough, line photons will be trapped, and they can only escape after Compton scattering, which degrades the photon energy. The possibility that resonant scattering in the post-shock region can also enhance the equivalent width of the Fe K$`\alpha `$ feature, when viewed along the magnetic axis, has been discussed by Terada et al (1999). The true test of this attractive idea would be a phase-resolved observation of a system whose geometry is established from other means. Similarly, better data as well as more detailed calculations are required to test if the resonant trapping/Compton scattering mechanism can quantitatively explain the broadened line profiles.

## PRE-SHOCK FLOW

If the shock height is small, then a substantial fraction of the X-rays emitted in the post-shock region must traverse the pre-shock flow.

### Absorption Effects

As already mentioned, Norton and Watson (1989) discovered that partial-covering absorption was necessary to model the spectrum of IPs. The true complexity of the intrinsic absorption is revealed by the following exercise, divsed by Osborne (1994, private communication). In this, one fits a simple (Bremsstrahlung plus a single, cold, absorber) model to the ASCA spectrum of a magnetic CV, first in the entire ASCA range ($``$0.4–10 keV). Next, raise the lower limit of energy range to 1.0, 2.0, 3.0, etc. keV and repeat the fit. One finds that the best fit kT and NH values change as a function of the lower limit, revealing that the observed spectrum cannot be adequately modeled by a simple absorber.

This is not surprising: the lines-of-sight to different parts of the post-shock gas pass through different amounts of pre-shock flow. One can generalize the partial covering absorber model to, say, a power-law distribution of NH (Done and Magdziarz 1998). Unfortunately, even this is a phenomenological model; unless we can determine the complex geometry of the accretion flow, one cannot create a realistic model of the complex absorption.

Moreover, the pre-shock flow should be ionized. The degree of ionization is, again, geometry-dependent: the pre-shock flow is less ionized if more radiation can escape from the sides of the post-shock region. If little radiation can escape from the sides (if accretion column radius is much larger than the shock height, this would be the case just above the shock), the ionization parameter $`\xi `$ can be several hundred. Indeed, ASCA and ROSAT spectra of several magnetic CVs do suggest the presence of warm absorbers, but do not uniquely prove it: an alternative fit involving multiple partial covering absorber can usually be found (see, e.g., Fig. 4). To reach a final conclusion on this issue, one would need to attempt a direct detection of warm absorber edges in higher spectral resolution, high signal-to-noise spectra.

### Scattering and re-emission

However, absorption is not the only physical process that takes place in the pre-shock flow. Some electron scattering will take place, particularly if the plasma is significantly ionized; also, photoelectric absorption will be followed by continuum and/or line emission at lower energies. Because the geometry of the pre-shock flow is far from spherically symmetric around the emission region, the relative importance of absorption vs. emission depends strongly on the viewing geometry, which in turn depends on the binary inclination angle and the spin phase.

In most magnetic CVs, the case for photoionized plasma emission has been ambiguous. For example, Kallman et al (1993) included such a component in their interpretation of the BBXRT spectrum of BY Cam. However, Done and Magdziarz (1998) consider that low energy features in this system, seen with ASCA, are due to thermal emission from multi-temperature plasma. The best case so far of photoionized plasma emission may be found in the spin minimum of AM Her: although the predominant magnetic pole is hidden behind the white dwarf limb at this phase, a weak X-ray source remains in view (Ishida et al 1997). The pre-shock flow is a plausible origin of this component. However, emission from the second pole cannot be excluded.

## CHANDRA HETG OBSERVATION OF V1223 SGR

In previous sections, I have reviewed the state of the art in X-ray spectroscopy of (mostly magnetic) CVs in 1999. Although many recent data (notably those obtained with ASCA and BeppoSAX) remain to be exploited fully, there are several important questions that require higher resolution data. Perhaps the instrument that would have best matched the scientific objectives was the XRS on-board ASTRO-E; nevertheless, grating instruments on XMM-Newton and Chandra both have the potential to advance our knowledge of CV X-ray spectra significantly. In the rest of this article, I will describe preliminary results from a Chandra HETG observation of X-ray bright IP, V1223 Sgr, obtained on 2000 April 30/May 1 for approximately 52 ksec. Analysis and interpretation of the data are still ongoing.

First, the Fe K$`\alpha `$ feature of this system clearly consists of 3 separate components (see Fig. 5). The detected line widths for the 6.7 and 6.97 keV thermal lines appear consistent with instrumental resolution, confirming the ASCA result. Thus, V1223 Sgr should serve as a baseline; a Chandra HETG observation of AO Psc, the system for which the line broadening is most securely detected in ASCA data, is scheduled in Cycle 2.

The 6.4 keV (fluorescent) component also largely appears to be narrow, however a shoulder, extending down to $``$6.2 keV, is tentatively detected. This might indicate that some 6.4 keV line photons (but not the thermal line photons) suffer Compton degradation. Alternatively, some of the fluorescence line may originate in the pre-shock flow, which may have sufficient bulk velocity to Doppler shift the line energy by up to $``$200 eV.

The MEG data (Fig.6) reveal a rich emission-line spectrum, comparable to that of EX Hya observed with the ASCA SIS. Lines are detected from H- and He-like ions of O, Ne, Mg, Si, S, and probably Ar. Moreover, the MEG resolves the resonant, intercombination, and forbidden triplets from He-like ions of O, Ne, and Mg. The forbidden lines are very weak, suggesting relatively high density, and the intercombination lines are comparable to the resonant components, suggesting photoionization origin. Moreover, we do not detect strong Fe L features; their weakness is another evidence for photoionization origin. Finally, when plotted in 10 eV bins, the presence of the OVII warm-absorber edge at 0.739 keV is strongly suggested.

It is therefore likely that these features originate in the immediate pre-shock flow. The only other body near the emission region, the white dwarf itself, has a much higher density and therefore a lower ionization parameter value than is required to produce these features. As we have already speculated, the pre-shock flow on the other hand is a natural origin both for the warm absorber and of the photoionized plasma emission. Detailed analyses are in progress to quantify the conditions of the pre-shock flow using these new diagnostics opened up by the Chandra grating instrument.

We are left with a major question, however: Where are the expected thermal lines from the post-shock, cooling plasma? The presence of a strong Fe L bump in the ASCA spectrum of EX Hya, among other clues, suggests that it is the lines from coronal plasma that we have detected with the SIS. Analyses of Chandra and XMM-Newton grating observations of this system are expected to confirm this interpretation. Even though similar K shell features are seen, V1223 Sgr appears completely different in this respect. Comparison of ASCA spectra suggest that V1223 Sgr is more typical of magnetic CVs, whereas EX Hya belongs to a smaller group perhaps with 3 other systems, characterized perhaps by low accretion rates.

One possible explanation for the difference between the two magnetic CVs (perhaps two groups of magnetic CVs) is the following. In the Aizu model, the plasma density increases sharply towards the bottom of the shock. Eventually, the post-shock region becomes optically thick in the continuum. In EX Hya, this happens at about 0.65 keV, as suggested by the analysis of Fujimoto and Ishida (1997), while in V1223 Sgr, this happens at a much higher temperature, say kT$``$2 keV as has been suggested by Beardmore et al (2000). Earlier works calculated the optical depth across the accretion column just below the shock; their conclusion that that region is optically thin remains valid. However, most of the observed medium energy (0.5–10 keV) X-rays actually originate much lower down in the post-shock region (Fig.7). For example, this plot shows that, for a $`kT_s`$=20 keV shock, we need to use E$``$40 keV continuum photons to probe the upper, low-density part of the post-shock region. Most of the photons we observe with ASCA, for example, must come from the bottom $``$1% of the shock! The situation is more extreme for low energy line photons: the emissivity of K$`\alpha `$ lines of elements up to Si and S, as well as those of the Fe L lines, peaks at temperatures below 2 keV in a collisionally ionized plasma. In an Aizu shock with $`kT_s`$=20 keV, this is the bottom $``$0.1% of the shock, where the density is $``$20 times greater than immediately below the shock.

These line photons cannot escape freely by travelling vertically, either, partly because of the photoelectric absorption in the cool, pre-shock, flow. Moreover, the vertical electron scattering optical depth within the post-shock region is non-negligible: the characteristic cooling time of the post-shock plasma is inversely proportional to the local accretion rate, and is $``$1 s at 1 g cm-2s-1. Thus we expect $``$1 g cm-2 of vertical column in the post-shock region; in fact, a numerical integration which properly accounts for the deceleration shows that the actual number is $``$2 g cm-2, with a weak dependence on the primary mass. It therefore appears that most of the low energy line photons from the bottom can be hidden by a significant optical depth; exactly how this works depends on the shock geometry. In addition, the white dwarf atmosphere has a scale height of several kilometers, when heated to $`kT20`$ eV (this applies only to the area surrounding the accretion column; underneath the column, the ram pressure will reduce the scale height). Thus, there may be a puffed up ring of white dwarf atmosphere, further reducing our ability to see the low energy line emissions. Thermal Fe K$`\alpha `$ lines, which come from higher up in the post-shock region, is probably responsible for the 6.7 and 6.97 keV features seen in V1223 Sgr. The shock models would suggest that the Fe K$`\alpha `$ emitting region is optically thin in the continuum. This is consistent with the Hellier et al. (1998) interpretation of the broad Fe K$`\alpha `$ lines, which requires a continuum optical depth of 0.05–1.0 in the line emitting region, where the temperature is kT$``$5 keV.

The accretion rate also controls the efficiency of pre-shock flow heating. Low accretion rate shocks are so tall that most of the radiation escapes from the side, lowering the degree of photoionization just above the shock. Thus, EX Hya could have very little emission lines from this region while V1223 Sgr is dominated by it.

## CONCLUSIONS SO FAR AND FUTURE PROSPECTS

The Aizu shock model attempts to connect the optically thick region immediately below the shock to the extremely optically thick white dwarf atmosphere. With hindsight, it is therefor natural that optical depth plays a significant role in shaping the observed X-ray spectra of magnetic CVs. It appears that the low energy X-ray line features of V1223 Sgr and EX Hya have different origins, despite some superficial similarities. There is evidence that the low energy lines seen in V1223 Sgr are from a dense, photoionized plasma, probably in the pre-shock flow. The lines in EX Hya is from the cooling post-shock plasma; these are weaker or hidden in V1223 Sgr. It is still likely that the Fe K$`\alpha `$ lines share the same origins (fluorescence off the white dwarf surface for the 6.4 keV line, thermal plasma emission for the other two components), since photoionization is unlikely to be strong enough to produce a significant Fe K$`\alpha `$ feature in any CVs. It has been suggested that resonant trapping of these line photons is the cause for the broadening of these thermal Fe K$`\alpha `$ lines in some systems; We await the recently approved Chandra grating observation of AO Psc to explore the details of this mechanism. Along with more observations, a more detailed modelling is also necessary. However, if optical depth does play an important roll in shaping the X-ray spectrum, as I have argued in this paper, then a realistic geometry must be used in numerical simulations. Therefore, it appears essential to improve our knowledge of the detailed shape of the post-shock region through, e.g., the observations of eclipses and self-occultations of the X-ray emitting region.

What can we expect for non-magnetic CVs? Despite the differences in geometry, X-ray emitting regions in magnetic and non-magnetic systems share some essential features: boundary layer emission in the latter must originate in a shock-heated plasma which cools down, and becomes denser, before it can settle down onto the primary surface. The depth of the gravitational potential, the size of the emitting region, and the local accretion rate are all similar between magnetic and non-magnetic CVs. Therefore, it is possible that the optical depth effects become important below some critical temperature, hiding any intrinsic line features, in non-magnetic as well as in magnetic systems. With future high quality X-ray spectra, we may begin to detect broadened lines, warm absorbers, and photoionized plasma emissions in non-magnetic CVs, a very exciting prospect indeed.

## ACKNOWLEDGEMENTS

I thank Chris Done and Coel Hellier for their critical reading of this manuscript. I also thank Manabu Ishida, and Tim Kallman for the insights that they have shared with me over the years.