text:
# Irregular Magnetic Fields in Interstellar Clouds and Variations in the Observed Circular Polarization of Spectral Lines

## 1 Introduction

The breadths of spectral line radiation from interstellar gas clouds at radio wavelengths are recognized to be much greater than would be caused by Doppler shifts resulting from thermal motions alone at the temperatures of the atoms and molecules in the gas. These breadths are believed to be an indication of disturbances (most likely a form of turbulence or waves) which are fundamental to the structure and evolution of interstellar gas clouds and star formation (e.g., Vázquez-Semadeni et al. 2000, and reviews in Lada & Kylafis 1999 and in Franco & Carraminana 1999). Although the conclusions of earlier investigations that relate detailed information about the spectral lines specifically to the properties of turbulence have been questioned (see Dubinski, Narayan, & Phillips 1995), recent efforts along these lines are encouraging (Miesch & Scalo 1995; Lis et al. 1998). Of particular importance is the role of waves or turbulence in causing a magnetic pressure that acts parallel to the direction of the mean magnetic field lines in the clouds. This has long been viewed as a promising way to understand why the rate of collapse of interstellar clouds along the field lines is not more rapid than is indicated by the observational data. Magnetic disturbances have been thought to decay more slowly than ordinary disturbances, and hence more likely to persist and provide the component of the pressure that seems to be necessary (Arons & Max 1975). However, the premise that the magnetic disturbances are longer lived has been questioned in recent investigations (MacLow et al. 1998; Stone, Ostriker, & Gammie 1998; Padoan & Nordlund 1999). In any case, irregularities in the magnetic fields are expected to be associated with the bulk motions which are reflected by the non-thermal breadths of spectral lines. At some level, these should be evident in the observational data.

Information about the direction of the magnetic field in the interstellar gas has long been inferred from the linear polarization of starlight that results from selective absorption by dust grains. However, the polarizing power of the grains that absorb starlight apparently is greatly reduced in gas clouds of even modest densities (e.g., Goodman 1996). Though such polarization data are extensive, they do not thus seem to be sensitive indicators for the irregularities of the magnetic field in dense clouds. In contrast, there is no doubt that the magnetic field in interstellar clouds is evident in the circular polarization of spectral lines at radio frequencies through the Zeeman effect. Although such data are not so extensive as is desirable for extracting information on irregularities in the field, advances in the observational capabilities are continuing. It is thus useful to examine quantitatively the sensitivity of the circular polarization of spectral lines to plausible irregularities in the magnetic field in interstellar clouds. In particular, the polarization associated with a ray of radiation is a sum of contributions by the gas all along the path of the ray. The effects of deviations of the field from its mean value thus tend to cancel for the radiation that emerges from the cloud and to reduce its sensitivity as an indicator of irregularities. The goal of our investigation is to determine what information remains in the circular polarization of radiation that emerges from an irregular medium. Though not the emphasis of our investigation, an important challenge is to identify observational quantities for which the predictions distinguish between actual MHD phenomena and a simple generic description of the irregularities in the fields. Finally, we note that the emission of linearly polarized radiation has been observed from grains in dense clouds. This linearly polarized radiation holds promise to become a valuable tool (e.g., Hildebrand 1996; Rao et al. 1998) as more sensitive and extensive maps are obtained.

In Section 2, we describe how the irregularities are characterized and how the calculations of the circular polarization are performed. We focus on the quantity that commonly is obtained from the observational data—the single value $`B_\mathrm{f}`$ of the magnetic field that best fits the circular polarization (Stokes $`V`$ intensity) when the Zeeman splitting is much smaller than the breadth of the spectral line. By comparing the statistical properties of the $`B_\mathrm{f}`$ which are inferred from the observations with the properties of the $`B_\mathrm{f}`$ which we calculate, we propose to gain information about the actual strength of the irregular component of the field within interstellar clouds. As an example, data for the 21 cm spectral line toward Orion A are analyzed in this manner in Section 3.

## 2 Calculations

The calculations here are restricted to the limit that ordinarily is relevant for interstellar gas clouds—that in which the Zeeman splitting is much smaller than the breadth of the spectral line. The optical depth $`\tau (v)`$ for the intensity of the spectral line as a function of Doppler velocity $`v`$ can then be expressed as,

$$\tau (v)=f(s)\mathrm{exp}((vv_\mathrm{t})^2/v_{\mathrm{th}}^2)𝑑s,$$
(1)
where the integral is along a straight-line path of length $`s`$ through the medium in which the turbulent velocity $`v_\mathrm{t}`$ varies along the path. The integral can be understood as a sum of the contributions from components of the gas along the line of sight where the bulk velocity of the component at location $`s`$ is $`v_\mathrm{t}(s)`$. The exponential arises from the Maxwellian form for the distribution of atomic/molecular velocities at any location and depends upon the thermal velocity $`v_{\mathrm{th}}=\sqrt{2kT/m}`$. All of the necessary information about the excitation, the abundance of the relevant species, and the atomic/molecular data is incorporated into the opacity function $`f(s)`$. As measured by the Stokes $`V`$ intensity, the net circular polarization is then related to the difference $`\tau _V(v)1`$ between the optical depths associated with the two senses of circular polarization,

$$\tau _V(v)=2pv_{\mathrm{th}}^2f(s)B(s)(vv_\mathrm{t})\mathrm{exp}((vv_\mathrm{t})^2/v_{\mathrm{th}}^2)𝑑s,$$
(2)
where $`p`$ is a constant that depends upon the particular transition and $`B(s)`$ is the component of the magnetic field that is parallel to the path (and hence to the line of sight of the observer). Note that when $`B`$ is constant along the path of the radiation,

$$\tau _V=pB\tau (v)/v.$$
(3)
If the spectral line in the gas is being observed in absorption in front of a strong continuum source with intensity $`I_0`$ as occurs toward Orion A in the analysis in Section 3, the observed intensity $`I`$ is given by

$$I=I_0\mathrm{exp}(\tau )$$
(4)
and Stokes-$`V`$ is given by

$$V=\tau _VI,$$
(5)
so that if $`B`$ were constant

$$V=pBI/v$$
(6)
and $`B`$ could readily be determined from the observation of $`I`$ and $`V`$. Even though the magnetic field may not be constant along the path of the ray of radiation, a “best fit” (in the least squares sense) value $`B_\mathrm{f}`$ often is obtained from equation (3),

$$B_\mathrm{f}=(\tau _V\tau /v)𝑑v/p(\tau /v)^2𝑑v$$
(7)
or by the analogous fit to equation (6).

If there are irregularities in the magnetic field, the best fit values for $`B_\mathrm{f}`$ will not be the same for spectral lines that are emitted from different locations across the surface of the cloud. Even at a single location, irregularities will usually prevent equation (3) or (6) from being satisfied by a single value for the magnetic field at all Doppler velocities within the observed spectral line. These variations in the inferred $`B_\mathrm{f}`$ are related, though only indirectly through the emergent Stokes-$`V`$, to the variations in the magnetic field $`B(s)`$ along the path of the radiation. The purpose of the investigation here is to provide a statistical relationship between the variations in the values of $`B_\mathrm{f}`$ that are inferred from the observational data and the amplitudes of the actual irregularities in the magnetic fields that are the cause of such variations. While more sophisticated and involved analytical tools might be imagined, we will limit our attention here to two basic quantities. (i) The standard deviation of the values of $`B_\mathrm{f}`$ that are inferred from the radiation emerging at the various different locations which comprise the surface of the cloud

$$(B_\mathrm{f})_{\mathrm{sd}}=\left\{B_\mathrm{f}^2B_\mathrm{f}^2\right\}^{1/2}.$$
(8)
(ii) The average across the surface of the cloud of the residuals $`(B_{\mathrm{res}}^2)^{1/2}`$ that remain when best fits are obtained for $`B_\mathrm{f}`$ at the various locations where at a single location

$$B_{\mathrm{res}}^2=(\tau _VpB_\mathrm{f}\tau /v)^2𝑑v/(p\tau /v)^2𝑑v.$$
(9)
Irregularities in the velocity and magnetic fields in the gas are characterized here in a conventional manner by the exponent in the power law distribution for the amplitudes of the Fourier components, by the root mean square (rms) of the velocity and of the magnetic fields, and by the smallest wave number (corresponding to the largest scale length) to which the power law distribution extends. For turbulence, this scale length represents the largest scale length at which disturbances are injected into the gas. The scale length implies a “correlation length”, which we also use in describing the calculations. Representative configurations (or statistical “realizations”) for the fields are then created by statistical sampling with Gaussian distributions for the Fourier amplitudes. Such methods are standard and are described, for example, in our own previous investigations (Wallin, Watson, & Wyld 1998, 1999) as well as elsewhere (e.g., Dubinski, Narayan, & Phillips 1995) where turbulent velocity fields are created. The turbulent magnetic fields are created in the same way (e.g., Wiebe & Watson 1998). For comparison purposes, we present a few examples of the realizations created by statistical sampling where the slope of the power spectrum is exactly that given by the Kolmogorov exponent. However, we focus our main attention on calculations in which the spectrum is somewhat steeper. This steeper spectrum leads to velocities and magnetic fields in coordinate space for which the slope of the correlation function is in better agreement with that calculated utilizing the fields from the available MHD simulations (see below). We have computed the “structure function” for the magnetic and velocity fields, as well as for the mass density, from these MHD computations. We find that the variation of the structure function (as defined by, e.g., Frisch 1995) corresponds best to (separation)4/3, whereas the Kolmogorov variation is (separation)2/3. The variation of the structure function is related to the variation of the power spectrum (e.g., Frisch 1995). We thus focus on a distribution for the power spectrum in Fourier space in our statistical sampling $`(k^{7/3})`$ that is somewhat steeper than that of Kolmogorov $`(k^{5/3})`$. We have verified that the structure functions which are computed using the fields obtained from the steeper power spectrum do agree well with those that are computed directly from the fields of the MHD simulations. We note that the $`(k^{7/3})`$ dependence that we infer may not be a basic property of MHD turbulence, but may be a result of the numerical resolution in the MHD computations and the way in which the disturbances are injected. The representative configurations for the velocities and magnetic fields that are created according to the foregoing statistical procedure thus provide a “generic” description for the irregular fields in terms of a minimal set of basic parameters.

In detail, the foregoing generic description certainly is incomplete—most notably, it incorporates no correlation between the velocities and magnetic fields and it contains no information on the variations in the density of matter. To obtain a quantitative measure of the influence of these deficiencies, we also perform computations in which these quantities are provided by the results of numerical simulations by others for time-independent, compressible MHD turbulence. We refer the reader to the papers by these authors for a detailed description of the calculational methods (see Stone, Ostriker, & Gammie 1998). Certain basic aspects of these calculations should, however, be summarized. The turbulence is driven by continuously injecting disturbances at a range of scale lengths centered on one-sixth the length of an edge of the computational cube. This determines the minimum wavenumber in the spectrum of the turbulence, and hence the correlation length for the velocities. In the way that we measure the correlation length (see below), the resulting correlation length corresponds to about one-twelfth the length of the edge of the computational cube. Three cases are considered according to the relative strength of the mean magnetic field $`B_{\mathrm{avg}}`$ and the thermal gas pressure $`P_{\mathrm{th}}`$ as described by the usual “plasma parameter” $`\beta =(8\pi P_{\mathrm{th}}/B_{\mathrm{avg}}^2)=`$ 0.02, 0.2 and 2. For convenience, we label these as S, M, and W for strong, medium and weak, respectively. For all three cases, the rms turbulent velocity in each of the three orthogonal directions is close to the value of $`2.9c_\mathrm{s}`$ in the weak (W) field case where the velocities are essentially isotropic. With increasing strength of the mean magnetic field, the gas becomes somewhat anisotropic so that the rms turbulent velocities parallel and perpendicular to the mean field are $`3.4c_\mathrm{s}`$ and $`2.8c_\mathrm{s}`$, respectively, in the strong (S) case. Although we are limited to essentially a single ratio of the turbulent velocity to the sound velocity in the available MHD simulations, this ratio of approximately three is representative for relevant interstellar clouds. The chief difference between the three MHD cases is in the relative strengths of the mean magnetic fields which are in the ratio 1(W):3.2(M):10(S)—assuming that the same $`P_{\mathrm{th}}`$ is adopted in all three cases. The ratio of the standard deviation of the magnetic field in each of the three orthogonal directions to the mean magnetic field in these simulations is approximately 5/3 (W case), 2/3 (M case), and 0.15 \[parallel\]/0.25 \[perpendicular\] (S case). These ratios can be viewed as depending mainly on the Alfvenic Mach number—the ratio of the three dimensional rms turbulent velocity to the Alfven velocity. The Alfvenic Mach number is 5(W), 1.6(M) and 0.5(S) in the simulations (the sonic Mach number is approximately five in all three simulations). The turbulent magnetic fields are then similar in the three cases, as are the turbulent velocities. To relate the dimensionless results of the MHD computations to physical units, it is useful to recognize that there is freedom to choose the values for two additional scaling factors after $`\beta `$ has been specified. These can be the sound velocity $`c_\mathrm{s}`$ and the average matter density. Representative values for molecular clouds (e.g., Crutcher 1999) are $`\beta 0.01`$ to 0.1 and include $`3c_\mathrm{s}`$ within the range of rms turbulent velocities that are observed. Their gas temperatures typically are 10 to 30K, though the Orion A cloud that is analyzed in Section 3 is somewhat warmer ($`100`$ K).

Two idealizations will be considered for the opacity function $`f(s)`$—(i) $`f(s)`$ will be taken as a constant, and (ii) $`f(s)`$ will be taken as proportional to the matter density. The latter case is considered only in conjunction with the fields from the MHD simulations where we have information about the variations in the matter density.

In Figure 1, we show the results for the standard deviations $`(B_\mathrm{f})_{\mathrm{sd}}`$ that are obtained from the computed Stokes-$`V`$ as a function of the various parameters. For the statistically created fields, $`(B_\mathrm{f})_{\mathrm{sd}}`$ is independent of the average magnetic field in the gas. The results can thus be presented usefully in terms of the ratio $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ where $`B_{\mathrm{rms}}`$ is the rms of the irregular (or turbulent) component of the magnetic field in the gas. That is, $`B(s)=B_{\mathrm{avg}}+B_{\mathrm{irregular}}`$ for which $`B_{\mathrm{rms}}B_{\mathrm{irregular}}^2^{1/2}`$. Calculations are presented when $`B_\mathrm{f}`$ is determined separately for each of the (128)2 rays that emerge from the grid points on the surface and perpendicular to the surface of our computational cube consisting of (128)3 grid points. Calculations also are presented when the intensity and Stokes-$`V`$ for rays from separate areas of 8x8 surface grid points are first summed before the $`B_\mathrm{f}`$ are obtained from equation (7). The latter is intended to provide an indication of the influence of the finite angular resolution of a telescope. The calculations are performed as a function of the minimum wavenumber (our parameter $`k_{\mathrm{min}}`$ as defined in Wallin et al. 1998) at which the irregularities are injected into the medium. The correlation length is then computed from the structure function of the resulting fields in coordinate (that is, ordinary physical) space. We define the correlation length here as the separation at which the structure function “becomes flat”. Due (presumably) to computational imprecision, the structure function does not become completely flat but has a wavelike variation with an amplitude of a few percent. We thus take the correlation length to be the separation at which the structure function becomes flat after this wave has been subtracted. In Figure 1, the quantity $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ is then presented as a function of the number $`N_{\mathrm{corr}}`$ of such correlation lengths along an edge of the computational cube. Other than those at a single correlation length corresponding to $`N_{\mathrm{corr}}=12`$, all of the computations are performed with fields obtained from the power spectrum $`(k^{7/3})`$ that is somewhat steeper than Kolmogorov. Within the range of turbulent velocities in our computations (see Figure 2) which covers the range for interstellar clouds, any dependence of $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ upon the rms value of the turbulent velocities is negligible. For $`N_{\mathrm{corr}}12`$, the longest wavelengths in the turbulence are an appreciable fraction of the size of the cube and the purely statistical variations in $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ from one realization of the cube to another (with otherwise identical parameters) are appreciable. This is indicated in Figure 1 by the broader lines. These lines are obtained by averaging $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ for several statistical realizations. For $`N_{\mathrm{corr}}12`$, the variations in $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ from one statistical realization to another are negligible.

In the right-hand panel in Figure 1, we present $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ obtained by utilizing fields from the MHD computations of Stone et al. (1998). These computations are available only for a single correlation length which corresponds to $`N_{\mathrm{corr}}12`$ because the turbulence in these computations always is injected with the same distribution of scale lengths. In the strong field case, the correlation length is somewhat different in the directions perpendicular to the average magnetic field. For ease of comparison, the $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ based on the statistically created fields for $`N_{\mathrm{corr}}=12`$ in the left-hand panel are indicated by the two horizontal lines in the right-hand panel. In addition to the turbulent magnetic and velocity fields, there is a non-zero average magnetic field in the MHD simulations. This influences the medium, and potentially the results of our computations (a non-zero average magnetic field has no effect on the results of computations in the left-hand panel with the statistically created fields). As noted previously, observations suggest that molecular clouds fall between the medium and strong field cases in terms of the plasma parameter $`\beta `$. There is no reliable way to incorporate variations in the mass density into the computations that utilize the statistical fields. Direct comparisons can thus only be made with the computations for the MHD fields where variations in the mass density are ignored (open symbols) in computing the optical depths. The computations in which $`f(s)`$ in equation (1) is proportional to the mass density are indicated by filled symbols. For the computations with the finite (8x8) surface areas and when the optical depths do not depend upon mass density, $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ ranges from approximately 0.25 to 0.30 in the case of medium (M) field strength according to whether the line of sight is parallel or perpendicular to the direction of the average magnetic field. For the analogous correlation length $`(N_{\mathrm{corr}}=12)`$, the computation with the statistically created fields is midway between these values. When the optical depths are assumed to be proportional to mass density, $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ can be seen to range from 0.33 to 0.39 when the same fields from the MHD computations are utilized.

In Figure 2, we show the ratio $`(B_{\mathrm{res}}^2)^{1/2}/B_{\mathrm{rms}}`$ calculated with fields that are similar to those used for Figure 1. Unlike the results in Figure 1, $`(B_{\mathrm{res}}^2)^{1/2}/B_{\mathrm{rms}}`$ does depend upon the magnitude of the turbulent velocities. Thus, the curves in the left-hand panel in Figure 2 are labeled according to the ratio of the rms turbulent velocity to the rms thermal velocity. In addition to the curves which are computed with the steeper power spectrum, calculations are again presented for $`N_{\mathrm{corr}}=12`$ that utilize the Kolmogorov value for the exponent in the power spectrum (circles and crosses). As in Figure 1, the right-hand panel shows $`(B_{\mathrm{res}}^2)^{1/2}/B_{\mathrm{rms}}`$ computed with the MHD fields. The turbulent velocities of the MHD fields are slightly different from our choices for the curves in the left-hand panel. The results of the computations with the statistically created fields which are presented for purposes of comparison in the right-hand panel (horizontal lines) are thus computed with an rms turbulent velocity that is not exactly the same as that of any of the curves in the left-hand panel. The agreement between the results of the computations with the statistical and with the MHD fields is generally good. The calculations that are intended to reflect a finite (8x8) angular resolution coincide for the statistical and MHD fields when the optical depths do not involve mass density. When the optical depths depend upon mass density (specifically, when $`f`$ is proportional to density), $`(B_{\mathrm{res}}^2)^{1/2}/B_{\mathrm{rms}}`$ is seen to be increased for the 8x8 resolution from about 0.24 to about 0.30–0.33 for medium fields. The increase is greater for the strong fields. Because $`B_{\mathrm{rms}}`$ is much less than $`B_{\mathrm{avg}}`$ in strong field case, we suspect that these MHD fields are not so representative of interstellar clouds where these fields probably are more similar in strength.

## 3 An Application to Observational Data

A statistical analysis such as the foregoing requires measurements of Stokes-$`V`$ in the radiation from a large number of locations on the surface of an interstellar gas cloud. There is only one existing data set that meets this requirement—an H I absorption-line Stokes $`I`$ and $`V`$ map across the surface of a region toward Orion A (T. H. Troland, R. M. Crutcher, D. A. Roberts, W. M. Goss, & C. Brogan, in preparation). The observations were performed with the VLA in the “C” configuration resulting in a spatial resolution of approximately 0.04 pc (15 arcseconds). The VLA synthesis mapping technique provides uniform sensitivity to Stokes $`I`$ and $`V`$ over the mapped region with no missing spatial points, which is a desirable characteristic for our analysis. Although observational data are available for a larger region, we limit the analysis to the region in which the background source is strong enough that the detections of the magnetic field are generally at a “three sigma” confidence level or better. With the angular resolution of the instrument, the region that meets our criteria for analysis consists of approximately 100 non-overlapping beams. Hence, we can obtain that number of independent observational values for $`B_\mathrm{f}`$ and $`B_{\mathrm{res}}^2`$ as a statistical sample for analysis. Within this area, there are no positions where the magnetic field was not detected at our sensitivity cutoff. Thus, the data set has no bias for positions with stronger magnetic fields or against positions with weaker magnetic fields. However, this data set does have significant limitations for an astrophysically meaningful analysis. Toward Orion A the 21 cm spectrum consists of a pair of partially overlapping absorption lines as indicated in the spectral data displayed in Figure 3 from a representative location on the surface. These two H I velocity components arise in a feature which has been called the Orion lid. An extensive analysis of H I absorption due to the Orion lid has been carried out by van der Werf & Goss (1989) based on earlier Stokes $`I`$ mapping. It consists primarily of neutral, atomic gas that probably is located close to and on our side of the H II region Orion A. The lid extends over at least 1.6 pc in the plane of the sky. Van der Werf & Goss argue that the more negative velocity H I is primarily photodissociated H2, some of which has been shocked, and that the more positive velocity component arises in the near envelope of the molecular cloud that is behind the H II region. We only utilize data from the positive velocity side of the intensity minimum that occurs in Figure 3 at approximately 5 km s-1. This is primarily because the line shape of the other component leads to typically a factor of two lower sensitivity to $`B_\mathrm{f}`$, so fewer positions are available for statistical analysis. In addition, the more positive velocity component is likely to be more relevant to the general interstellar medium, and we wish to avoid confusion in the analysis due to the overlap of these two components. Specifically, only spectral channels on the positive velocity side and where the absorption is less than ninety percent of the maximum absorption are included. The data give only the optical depth of the H I line. In order to infer column densities, one must know the spin temperature of the line. Toward the main exciting star of the H II region, $`\theta ^1`$C Orioni, the H I column density has been measured with the ultraviolet Lyman-$`\alpha `$ absorption line; when compared with the 21-cm line optical depth, this yields a spin temperature of about 100 K. However, there is no information about possible variations in spin temperature over the mapped area. There is only the upper limit H2/H I $`<10^4`$ available; at least toward $`\theta ^1`$C Orioni, the Orion lid is not a molecular cloud. Also, there is no direct information about the volume density of H I. An assumption of a spherical cloud and the fact that absorption is seen over at least 1.6 pc in the plane of the sky leads to $`n_\mathrm{H}10^3`$ cm-3. However, it seems likely that the Orion lid has a more sheet-like geometry, which would raise the volume density, but by an unknown amount. Because of these limitations, our analysis of this Orion A data set should be considered to be primarily an example of the application of our analysis technique rather than one from which firm astrophysical conclusions should be drawn.

The formal mean square error associated with the least squares fit to obtain $`B_\mathrm{f}`$ in equation (7) from a single, observed spectral line profile is

$$\sigma _\mathrm{f}^2=B_{\mathrm{res}}^2/(N1),$$
(10)
where $`N`$ is the number of intensities within the spectral line that are used to find the value of $`B_\mathrm{f}`$ from that spectral line. For the spectra included in Figure 3, the confidence level thus corresponds to $`B_\mathrm{f}3\sigma _\mathrm{f}`$.

Differences between Stokes-$`V`$ and $`I/v`$ that lead to $`\sigma _\mathrm{f}^2`$ can result from noise in the background source and in the instruments, as well as from the irregularities in the magnetic field that we are seeking to understand. In the relevant regime ($`\tau _V1`$ and $`V/I1`$), the noise makes a contribution to $`\sigma _f^2`$ that can be expressed as,

$$\sigma _{\mathrm{th}}^2\underset{i}{}\sigma _{V,\mathrm{th}}^2(B_\mathrm{f}/V_i)^2,$$
(11)
where $`\sigma _{V,\mathrm{th}}^2`$ is the mean square variation in the measured $`V`$ that would occur even if the magnetic field were constant along the path of a ray. The sum is over the index $`i`$ that designates the Stokes parameters $`V_i`$ which are used in the computation of a particular value of $`B_\mathrm{f}`$. The quantity $`\sigma _{V,\mathrm{th}}^2`$ is computed for an adjacent interval of velocities that is the same size as that used in computing $`B_\mathrm{f}`$, but which is located outside of the spectral line. At these velocities, measurements of $`V`$ are not influenced by variations in the magnetic field whereas the variations due to the noise should be essentially the same as within the spectral line. We thus utilize

$$\sigma _{V,\mathrm{th}}^2=\underset{j}{}V_j^2/(N1),$$
(12)
where the index $`j`$ now designates locations in velocity that are spaced equivalently to those that are used in computing $`B_\mathrm{f}`$, but which are located in the interval outside of the spectral line. The quantities $`\sigma _\mathrm{f}^2`$ and $`\sigma _{\mathrm{th}}^2`$ are computed for each profile and are then averaged over the face of the cloud. We find for these data toward Orion A

$$\sigma _{\mathrm{th}}^2/\sigma _\mathrm{f}^20.1.$$
(13)
Since $`\sigma _{\mathrm{th}}^2/\sigma _\mathrm{f}^2`$ is relatively small, for simplicity we will ignore the contribution of the noise in our further discussions and make the simplifying approximation that $`\sigma _\mathrm{f}`$ and $`(B_\mathrm{f})_{\mathrm{sd}}`$ are due entirely to variations in the magnetic field.

From these observational data toward Orion A, we find (in agreement with Troland et al.) for the mean value of $`B_\mathrm{f}`$ obtained by averaging over the approximately 100 independent telescope beams that are included in our analysis

$$B_\mathrm{f}^{\mathrm{obs}}=170\mu \mathrm{G},$$
(14)
for the standard deviation of these values

$$(B_\mathrm{f})_{\mathrm{sd}}^{\mathrm{obs}}=130\mu \mathrm{G},$$
(15)
and for the average of the residuals

$$(B_{\mathrm{res}}^2)^{1/2}^{\mathrm{obs}}=83\mu \mathrm{G}.$$
(16)
A histogram of the approximately 100 values of $`B_{\mathrm{f}}^{}{}_{}{}^{\mathrm{obs}}`$ used to find $`B_\mathrm{f}^{\mathrm{obs}}`$ is given in Figure 4.

To apply our calculations in Figures 1 and 2 for the purpose of inferring the strength of the irregular component of the actual magnetic field in the Orion A cloud as parameterized by $`B_{\mathrm{rms}}`$, the size of the actual gas cloud measured in correlation lengths $`N_{\mathrm{corr}}`$ enters. Of course, real interstellar gas clouds are unlikely to be statistically uniform entities as are our idealized cubes. We can nevertheless estimate a correlation length from the variation of the measured quantities across the face of the Orion A cloud. That is, the structure function is computed for the 21 cm line and a correlation length is inferred in the same way as was done for our theoretical velocities in obtaining Figures 1 and 2. In this way, we infer that there are approximately ten correlation lengths across the face of the region that is included in our analysis. The size of each telescope beam then represents about one correlation length. Hence, the relationship between the beam size and the correlation length is known, but the relevant size of the cloud (which is needed to obtain $`N_{\mathrm{corr}}`$) is uncertain. Observations indicate only that the face of the entire cloud is greater than about 1.6 pc or about forty correlation lengths. The dimension of the cloud along the line of sight is unknown, but most likely is considerably smaller. Computations are thus performed for the quantities in Figures 1 and 2 by first summing the intensities within non-overlapping surface areas that are one correlation length on a side. These are then averaged over the surface of the cube. Previously, we described first summing over fixed areas of the surface of the cube to obtain the curves for “8x8 rays” in these Figures. Now the sum is over a number of rays which varies with $`N_{\mathrm{corr}}`$. This seems to be the most useful way to relate our calculations to the available observational data for Orion A. The results are represented by the lines labeled “Obs” in Figures 1 and 2. The “Obs” curve in Figure 1 is somewhat less sensitive to the uncertainty about $`N_{\mathrm{corr}}`$ than are the other curves. In Figure 2, the “Obs” curve is quite insensitive to $`N_{\mathrm{corr}}`$ and has a value that is close to 0.4 over the entire range in Figure 2. The actual $`B_{\mathrm{rms}}`$ that is then implied by the value of 83 $`\mu `$G for $`(B_{\mathrm{res}}^2)^{1/2}^{\mathrm{obs}}`$ is

$$B_{\mathrm{rms}}210\mu \mathrm{G},$$
(17)
which indicates that $`B_{\mathrm{rms}}\left|B_\mathrm{f}^{\mathrm{obs}}\right|`$. The observed value $`(B_\mathrm{f})_{\mathrm{sd}}^{\mathrm{obs}}=130\mu `$G, together with the “Obs” curve in Figure 1, seem to imply a similar, though somewhat larger, value for $`B_{\mathrm{rms}}`$. The ordinate for the “Obs” curve ranges from approximately 0.35 at $`N_{\mathrm{corr}}=10`$ to slightly less than 0.2 at $`N_{\mathrm{corr}}=50`$. If $`N_{\mathrm{corr}}10`$ as seems most appropriate, the inferred random magnetic field is

$$B_{\mathrm{rms}}370\mu \mathrm{G}.$$
(18)
In view of the involved procedures that are used in these two assessments of $`B_{\mathrm{rms}}`$, the apparent discrepancy is not surprising. If the optical depths increase with mass density, the MHD calculations in Figures 1 and 2 indicate that the sense will be to raise the calculated values in both Figures. The increase will be a larger fraction in Figure 1 than in Figure 2—an adjustment that is in the correct direction to reduce the discrepancy between the values for $`B_{\mathrm{rms}}`$ that are inferred in the two approaches. An adjustment of this nature will also decrease the value for $`B_{\mathrm{rms}}`$ that is inferred from $`(B_{\mathrm{res}}^2)^{1/2}^{\mathrm{obs}}`$ to about 130 $`\mu `$G if one simply scales linearly from the results in the right hand panel in Figure 2. Our best estimates thus yield the conclusion that $`B_{\mathrm{rms}}`$ is similar in magnitude to the mean value of the magnetic field $`B_\mathrm{f}^{\mathrm{obs}}`$ for this gas toward Orion A.

We cannot exclude quantitatively the possibility that large scale changes in the average magnetic field make a significant contribution to the variations that are inferred for the magnetic field in the Orion A cloud. An indication that such changes are small in directions along the line of sight is evident in Figure 3. The mean values of the magnetic fields are essentially the same in the feature on the left-hand side (the feature that we are analyzing) and in the feature on right-hand side of the profile. Since these two features probably are separated along the line of sight, we can conclude that the magnetic field probably does not change significantly in this direction. Deviations between nearby channels—which presumably are due to changes over much shorter distances—tend to be greater than the difference between the mean values in the two features. The Troland et al. (in preparation) map of $`B_\mathrm{f}`$ does show evidence of some systematic gradients across the map, though these variations are not simple and do not appear to be large enough to dominate. Note that $`(B_{\mathrm{res}}^2)^{1/2}/B_{\mathrm{rms}}`$ is sensitive to changes in the average magnetic field along the line of sight, but not to changes across the face of the cloud. In contrast, $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ is sensitive to changes that occur in directions across the face of the cloud, but not in the direction along the line of sight.

## 4 Discussion

A number of deficiencies can be imagined in the simplifying idealizations of our calculations. Nevertheless, the conclusion appears to be firm that plausible turbulent magnetic fields are likely to have a significant effect on the magnetic field strengths that are inferred from the observation of the circular polarization of spectral lines at radio wavelengths. That is, their effect is unlikely to “average out”. Clearly the exact ratio depends upon specific considerations outlined in the foregoing. One-third is a representative estimate for the ratio $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ of the rms deviations of the observationally inferred field to the rms of the irregular (or turbulent) magnetic field that actually exists in the cloud. As expected, there is a tendency for deviations from the average to cancel when the contributions are summed along the paths of the rays of radiation. When the opacity function $`f(s)`$ is a function of location $`s`$ because it depends upon matter density or other considerations, $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ will naturally tend to increase. This increase can be viewed as a result of effectively reducing the number of elements along the path that contribute which, in turn increases the scatter simply due to statistics. In MHD simulations, there should be some correlation between the matter density and the strength of the magnetic field. This also tends to increase $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$, though we suspect that the effect is small based on our examination of the correlation between density and field strength. The results for the strong field (S) case in Figure 1 might seem to be in conflict with this conclusion. We suspect, however, that the relatively larger values of $`(B_\mathrm{f})_{\mathrm{sd}}/B_{\mathrm{rms}}`$ in this case are due mainly to the “statistical” effect of greater contrasts due to the anisotropic variations in the density.

In a larger sense, the idealizations here are an incomplete description for star forming clouds because gravity is absent in these idealizations. Gravitational attractions will tend to cause clustering of matter, and hence more extreme variations in $`f(s)`$ than considered in our calculations. In contrast, the results when $`f(s)`$ is constant in Figures 1 and 2 suggest that purely statistical sampling provides velocity and magnetic fields that may be adequate to understand the circular polarization of clouds where gravity is not a major effect and when the finite resolution of the telescope is considered. This is not an altogether happy conclusion since it potentially limits our prospects for identifying the physical phenomena for which the non-thermal line breadths and variations in the magnetic fields are indicators.

Our interpretation that the irregular component of the magnetic field in the Orion A cloud is similar in magnitude to the average magnetic field ($`B_{\mathrm{rms}}=\left|B_\mathrm{f}^{\mathrm{obs}}\right|`$) implies that the magnetic field is quite disordered in this cloud. That the bulk motions of the gas are able to deform the field to this degree tends to indicate that the turbulent kinetic energy is at least comparable with the energy in the magnetic field, and thus that the turbulent velocity is at least comparable with the Alfven velocity. In numerical computations of MHD turbulence (including the MHD fields that we are using), the energy of the irregular component of the magnetic field tends to be similar to the kinetic energy of the turbulence (typically 30 to 60%; Vázquez-Semadeni et al. 2000). For Alfven waves, these two components of the energy are equal. From this, it also follows that the turbulent and Alfven velocities should be comparable when $`B_{\mathrm{rms}}=\left|B_\mathrm{f}^{\mathrm{obs}}\right|`$. Padoan & Nordlund(1999) have emphasized that MHD computations in which the turbulent velocities are super-Alfvenic lead to predicted features for molecular clouds that are different in important ways from computations with lower turbulent velocities. In view of the uncertainties in our present analysis for $`B_{\mathrm{rms}}`$, we are unable to say whether the turbulent velocities are somewhat larger or somewhat smaller than the Alfven velocity. Alternatively, the ratio of the turbulent and Alfven velocities could be obtained from the spectral linebreadth in Figure 3 and the average field strength in equation (14) if the mass density were known reliably. However, as discussed in Section 3, the mass density is poorly known. A mass density of about $`10^4`$ H-atoms cm-3 is required to yield an Alfven velocity that is equal to the turbulent velocity. As discussed in Section 3, the Orion A cloud probably is “sheet-like”. At least along one line of sight, its molecular fraction is small. It may not, thus, be representative of star-forming, molecular clouds. We note that in other molecular clouds where the matter density is better known, the velocities of the internal motions have been interpreted as equal to the Alfven velocities (Crutcher 1999).

We are especially grateful to J. M. Stone for kindly making available the results of the MHD simulations. We also acknowledge valuable contributions by A. B. Men’shchikov in developing the initial computational methods, and by H. W. Wyld in a number of helpful conversations. This work has been supported by NSF Grants AST98-20641 and AST99-88104.


format: text

text:
# Putting Proteins back into Water.

\[

## Abstract

We introduce a simplified protein model where the solvent (water) degrees of freedom appear explicitly (although in an extremely simplified fashion). Using this model we are able to recover the thermodynamic phenomenology of proteins over a wide range of temperatures. In particular we describe both the warm and the cold protein denaturation within a single framework, while addressing important issues about the structure of model proteins.

\]

One of the main goals of statistical physics in the last decade has been to understand the ”folding code”: how the amino-acid sequence of a protein (coded in DNA, the ”genetic code”), uniquely determines its functional (”native”) structure, or fold. Understanding the principles that drive a protein to fold to its native structure is of great conceptual and practical relevance, since it could lead, for example, to high specificity drugs.

Proteins are extremely complex structures: they are long heteropolymers made of up to $`20`$ different amino-acids species, each of them with its own chemical, electrostatic and steric properties; the physiological solvent, an aqueous solution, and its characteristics play a fundamental role both in the dynamics and in the thermodynamics of folding. It is therefore not surprising that only in recent times statistical physicists have begun working on this problem, mainly after the introduction of the so-called HP model, where the above mentioned richness has been reduced to a manageable level. In the HP model, proteins are modeled as self-avoiding polymers on a lattice (two or three dimensional), greatly reducing the number of accessible conformations. The chemical and electrostatic properties of amino-acids have also been simplified: indeed, it has been recognized since long that the main force stabilizing the native conformations of globular proteins is the hydrophobicity of non-polar amino-acids. Consequently, the important properties of amino-acids are reduced to two: they are either polar (ions or dipoles, labeled with P) or non-polar (H).

Hydrophobicity can be described as the tendency of hydrophobic molecules to reduce as much as possible their surface of contact with water: two hydrophobic molecules try to stick together in order to hide from water their mutual surface of contact. Consequently, hydrophobicity has been introduced in the HP model as an effective attractive interaction between H amino-acids. Then, the solvent degrees of freedom can be neglected. Here we show that such a simplification can be removed, and water can be taken into account, keeping the complexity of the model at a still manageable level: the benefits are a better description of the protein phenomenology (namely, cold destabilization and eventually denaturation) and some insights on the structure of the protein core.

In the last fifteen years there has been a growing body of evidence for the so called cold destabilisation of proteins: the free energy difference $`\mathrm{\Delta }F_D^N`$ between denaturate and native conformations of proteins has parabolic shape, with a maximum at temperatures of the order of $`1525^0`$C, or lower, implying that at lower temperatures the native conformation is less and less stable. In some cases, even the cold denaturation of proteins has been obtained .

The HP model is unable to deal with cold destabilization since its low temperature state is compact and more and more stable down to $`T=0`$: is a good description of cold destabilisation and eventually denaturation relevant for protein folding? We think that the answer is affirmative for at least two reasons.

In order to describe protein folding with a simple model, it is important to capture the essential physics of the process, at the temperatures at which it takes place. If the stability of native conformations of proteins begins to decrease below $`1525^o`$C, it is unlikely, at least a priori, that the physics responsible for such a behavior is not important around the maximal stability temperature, in a range relevant for in vivo protein folding. A further reason to believe that a good model for protein folding should also agree with the cold destabilisation phenomenology is that, actually, there is no clear-cut distinction between the physics that stabilises proteins, and the one that destabilises them. In both cases a re-analysis of the concept of hydrophobicity and of hydrophobic hydration is necessary.

Already Frank and Evans  identified the origin of hydrophobicity in the partial ordering of water around non-polar molecules (such as, for example, pentane, benzene and some amino-acids). Water molecules tend to build ice-like cages around non-polar molecules. Although a detailed analysis of these structures is, to our knowledge, still lacking (actually recently some better understanding and consensus are emerging), we can guess their energetic and entropic properties. Indeed, water molecules forming these cages are highly hydrogen bonded, much as in ice; consequently, their formation is energetically favorable with respect to bulk liquid water. Yet, the possible molecular arrangements in the cages are a small number compared to all the disordered molecular conformations typical of liquid water. The latter are energetically unfavorable with respect to bulk water because water molecules fail to form hydrogen bonds with hydrophobic amino-acids. Therefore the free energy of formation of a cage ($`F_{cage}F_{nocage}=\mathrm{\Delta }F`$) is a balance between an enthalpy gain/loss and an entropy loss/gain: ordered cages give an enthalpy gain ($`\mathrm{\Delta }H<0`$) and an entropy loss ($`\mathrm{\Delta }S<0`$); the scenario is the opposite for disordered states. All of the above arguments call for a model able to reproduce (at least qualitatively) such a rich phenomenology.

The model we propose here borrows two of the simplifications from the HP model: proteins are still modeled as heteropolimers on a lattice, made of just two different amino-acid species: polar (P) and non-polar (H). Then, we put proteins back into water: every site of the lattice that is not occupied by the polymer is occupied by water (in general, by a group of water molecules that can be arranged in $`q`$ states). Water is described using the Muller-Lee-Graziano (MLG) two-states model (Fig.1a) , that Silverstein et al. have recently shown to be consistent with a molecular model of the water-amino acid system . The energy of each $`H`$ amino-acid depends on the states of the water sites it is in contact with: as a simplifying assumption (see Fig. 1b), we say that out of the $`q`$ possible states of a water site, one can be singled out to be a cage conformations (labeled $`s=0`$), energetically favorable with energy $`J`$ ($`J>0`$), and the remaining $`q1`$ ($`s=1,\mathrm{},q1`$) states are energetically unfavorable with energy $`K>0`$ (they represent the disordered states of reduced hydrogen-bond coordination). We stress that the term (un)favorable is always with respect to bulk liquid water. Water sites that are not in contact with $`H`$ amino-acids (that is, bulk water sites) do not contribute to the energy (whereas they would have an energetic description according to the MLG model, that yet has five free parameters, too many for a simple theoretical model). $`P`$ amino-acids do not interact with water so that their energy is always $`0`$: such a crude approximation is made with the idea that hydrophobicity is the leading effect stabilizing the native conformation of proteins. Some better description of the water-P interaction would be welcome, but such ingredient is unnecessary for our present purposes.

Given a protein of $`N`$ amino-acids, with the sequence $`a_1,a_2,\mathrm{},a_N`$ ($`a_i=P`$ or $`H`$), the energy of the protein is then

$$E=\underset{<i,H>}{}(J\delta _{s_j,0}+K(1\delta _{s_j,0}))$$
(1)
where the sum is over the water sites that are nearest neighbors of some $`H`$ amino-acid. Starting from (1) we can write the partition function of the system as

$$Z_N=\underset{C}{}Z_N(C)$$
(2)
where $`Z_N(C)`$ is the partition function associated to a single conformation $`C`$:

$$Z_N(C)=q^{n_0(C)}\left((q1)e^{\beta K}+e^{\beta J}\right)^{n_1(C)}$$
(3)
where the dependence on the water degrees of freedom has been explicitly calculated. $`n_1(C)`$ is the number of water sites nearest neighbors of some $`H`$ amino-acid, $`n_0`$ is the number of bulk water sites.

We deal with model proteins of length up to $`N=17`$ on the square lattice, and compute the partition function, and all the thermodynamic quantities and averages by exact enumeration of the $`2155667`$ different conformations. We show the results for the particular sequence $`PHPPHPPHPHPPHPPHH`$. We choose $`J=1`$ (actually, both $`K`$ and the temperature $`T`$ can be normalized with respect to $`J`$), $`K=2`$ and $`q=10^5`$ (a better determination of these values could come from molecular dynamics and structural studies). We take the Boltzmann constant $`k_B=1`$.

In Fig.2 the specific heat $`C_v`$, and the average number of monomer-monomer contacts, $`n_c`$, are shown. The low-temperature peak in the specific heat coincides with a jump of $`n_c`$: at lower temperatures the protein is swollen, and maximizes the number of water-H contacts, in agreement with cold denaturation. The number of contacts, $`n_c`$, begins decreasing coinciding with the high-temperature peak of the specific heat, that therefore coincides with the usual warm denaturation phenomenon. Between $`T_c`$ and $`T_w`$ there is a region where the most probable conformation is the one represented in the inset of Fig.2: as it can be seen, it is compact with a hydrophobic core, out of reach for water (we also checked that this native state is unique, in that its Boltzmann weight is the largest above $`T_c`$). We have analyzed the behavior of different protein lengths and of different sequences, and we have always found the same qualitative behavior of $`C_v`$ and $`n_c`$. Our model is therefore able to describe, within a single framework, both cold and warm denaturation. Moreover, it shows a native state with a mostly hydrophobic core.

Although the ratio between $`T_c`$ and $`T_w`$ in Fig.2 is unphysical, using the full MLG model it is possible to come closer to real values: the price to be paid is the larger number of parameters to adjust. In this Letter we address the physical principles responsible for the thermodynamic behavior of proteins on a broad range of temperatures: we believe that the differences between the bimodal model and the MLG model (and other possible more refined models) govern the details of the behavior more than the essential features.

We next compare the free energy, enthalpy and entropy variations of folding of our model with those from the literature. Indeed, such a comparison is a difficult one, since it is hard to define what a denatured state is in our theoretical calculations. Therefore, as a simple approximation, we consider as denaturate those conformations with at most $`4`$ monomer-monomer contacts (a polymer of $`17`$ monomers over a square lattice has at most $`9`$ monomer-monomer contacts). The native state has $`8`$ monomer-monomer contacts.

In Fig.3 we show $`F_{Denaturate}F_{Native}=\mathrm{\Delta }F_N^D`$, $`\mathrm{\Delta }H_N^D`$ and $`T\mathrm{\Delta }S_N^D`$. They coincide qualitatively with the ones from experiments. We point out the presence of two temperatures below and above which $`\mathrm{\Delta }F_N^D<0`$: the denatured state of our model protein is more stable than the native state. Between these two temperatures, instead, $`\mathrm{\Delta }F_N^D>0`$, and the native state is the most stable. In the same temperature range where $`\mathrm{\Delta }F_N^D>0`$, $`\mathrm{\Delta }H_N^D`$ and $`T\mathrm{\Delta }S_N^D`$ have a strong temperature dependence: they even change sign, a signature of the rich physics behind the water-protein system. At high temperatures we find that both $`\mathrm{\Delta }H_N^D`$ and $`\mathrm{\Delta }S_N^D`$ saturate ($`T\mathrm{\Delta }S`$ grows linearly, therefore $`\mathrm{\Delta }S`$ saturates), as experimentally observed. Some particular care should be paid to the low temperature behavior of $`\mathrm{\Delta }H_N^D`$ and $`T\mathrm{\Delta }S_N^D`$. Indeed, $`\mathrm{\Delta }H_N^D`$ goes to a constant value, which is consistent with a lower bound for the energies, and $`T\mathrm{\Delta }S_N^D`$ tends to $`0`$ with $`T`$. Experiments should be made below $`T_c`$ to assess such a behavior (although a recent model suggests such scenario ). We find therefore that our model reproduces qualitatively the known calorimetric data of protein denaturation over a broad range of temperatures.

The hydrophobic effect is often modeled through attractive effective $`HH`$ interactions. Within our framework, we consider a system of two $`H`$ amino-acids in solution and we compare the partition function of the system when the two amino-acids are in contact, $`Z_c`$, with the one when the two amino-acids have no mutual contacts $`Z_0`$. The effective attractive interaction is defined as $`ϵ=T\mathrm{ln}(Z_c/Z_0)`$ ($`ϵ`$ is positive if attractive, with this definition). The $`T\mathrm{}`$ limit is

$$ϵ(T\mathrm{})=2K\frac{2}{q}(K+J)$$
(4)
and is attractive for large values of $`q`$: it is the usual hydrophobic effective interaction. Yet, the $`T=0`$ limit is $`ϵ=2J`$, repulsive. A meaningful effective interaction should at least include such a temperature dependence. Actually, the strong temperature dependence of $`ϵ`$ is not the only limitation to a definition of an attractive $`HH`$ interaction. Indeed, such an interaction can be meaningful only for amino-acids surrounded by water molecules, but it cannot be defined in the core of proteins, where water is absent. As a consequence, in the absence of some true interactions between amino-acids, the hydrophobic interaction alone is not able to favor thermodynamically the native state against different compact states obtained by reordering only the core of the protein. As an example, the two conformations in Fig.4, corresponding to the sequence $`PPHPPPPPPPPPHHHP`$ have the same probability to occur in our model, since they hide and expose to water the same number of H amino-acids.

Therefore this model suggests that it is improper to define interactions of hydrophobic origin inside proteins, and that the detailed structure of the cores of proteins should be stabilized by other mechanisms. Indeed, in the biochemistry literature the debate is still strong whether the hydrophobic interaction alone is able to enforce the full native state of proteins or other interactions should also be taken into account. Effective interactions can be safely defined whenever they substitute some non-changing environment. When a protein is folding, its amino-acids find instead an ever changing environment that depends on water and on the other amino-acids. Even the reliability of two-body effective interactions vs. many-body ones is an open issue still to be settled. It is therefore intrinsically difficult to define effective potentials of some general validity between amino-acids: our model points out such a problem for hydrophobic interactions.

In conclusion, we have introduced a model of proteins in water that is able to reproduce the known features of proteins (namely, cold destabilisation and warm denaturation, a native state with a mostly hydrophobic core, and the correct free energy, enthalpy and entropy of folding). We also checked our results for different protein lengths, sequences, parameter values and even implementing the full MLG model for the description of water. Although some details may change, the overall behavior is consistent and robust. Moreover, lattice models are intended to be only qualitatively instructive, whereas a quantitative description can be given only by more detailed off-lattice models.


text:
# Decoherence in a superconducting ring

## I Introduction and results

A superconducting ring is a good candidate for macroscopic quantum memory . Its ability to store quantum information stems from the complicated “vacuum structure”, viz. the presence of different sectors of states corresponding to different values of magnetic flux through the ring.

It is easy to suppress spontaneous transitions between these different sectors by making the ring sufficiently thick . In the absence of such transitions, the main mechanism by which quantum information stored in the ring deteriorates is decoherence, i.e. a decrease, in time, of matrix elements between states with different values of the flux. Decoherence is due to interaction of the macroscopic variable, in our case the supercurrent in the ring, with the “environment” comprised by local degrees of freedom—fluctuations of the electromagnetic field.

The supercurrent density required to support a single flux quantum through a macroscopic ring is very small. This makes that state locally almost indistinguishable from the state with no flux at all. Based on this observation, we have suggested  that in this case decoherence may not be as strong as it is generally believed to be for macroscopic systems. The purpose of the present paper is to discuss decoherence in the ring quantitatively.

In general, the environment starts out in an arbitrary initial state with some density matrix $`\rho _0`$. Two different values of the macroscopic variable give rise to two different histories of the environment, which are represented by evolution operators $`U_1(t)`$ and $`U_2(t)`$. We define the amount of coherence (at time $`t`$) between these two histories as the “overlap”

$$𝒞(t)=\text{Tr}[U_1(t)\rho _0U_2^{}(t)].$$
(1)
Next, we define a decoherence exponent $`D(t)`$ by

$$\mathrm{exp}[D(t)]=\left|𝒞(t)\right|^2.$$
(2)
For an environment comprised by linear oscillators, which in addition interact linearly with the macroscopic variable, decoherence had been widely discussed in connection with macroscopic tunneling . That problem is similar but not identical to ours: in quantum memory tunneling is suppressed. Indeed, our problem is more analogous to the “error-free” case of ref. ; however, we will be able to supply definite values for all of our parameters.

In principle, the linear approximation can break down at large times even for weakly nonlinear environments. In addition, whatever we designate as the environment will interact with a still larger system and will loose energy to that system at some characteristic dissipation rate. For example, in our case, the electromagnetic field can be absorbed on the surface of the superconductor, as well as on the surfaces of any surrounding objects. In such cases, the evolution of the environment can no longer be represented by unitary operators.

It is certainly of interest to understand how nonlinearity and dissipation in the environment affect quantum coherence in the original system. However, for our present purposes the linear approximation will be sufficient. A mechanism for decoherence exists already at the linear level . Dissipation will only be discussed phenomenologically, to make sure that it does not significantly alter the evolution of the oscillators at relevant times.

Specifically, we consider the following problem. Suppose at the initial time $`t=0`$ a supercurrent has been set up in the ring but the electromagnetic field is still in the thermal state corresponding to zero current, with a Gibbs distribution $`w_m`$ of probabilities for the energy eigenstates $`|m`$. Because the field is displaced from its true equilibrium (corresponding to a nonzero average magnetic field), it will oscillate. Consider two oscillation histories, corresponding to two different values of the current (one of which may be zero), and call their evolution operators $`U_1(t)`$ and $`U_2(t)`$. The decoherence exponent in this case will be determined by

$$\mathrm{exp}[D(t)]=\underset{m}{}\left|m|U_2^{}U_1|m\right|^2w_m.$$
(3)
We describe a calculation of $`D(t)`$ below.

The problem just posed is admittedly somewhat artificial. Imagine that we switch from the state with no flux (and no current) to the state with unit flux by dragging a single flux line through the bulk of the superconductor. There is no reason to expect that the current will equilibrate first, and the field next, as in our model problem. Nevertheless, we think that it is good to know the answer to this model problem before turning to more realistic setups.

The answer is this. Consider a circular ring of radius $`R_0`$ made out of a London-type superconducting round wire of radius $`R_1`$, in the limit

$$\delta R_1R_0,$$
(4)
where $`\delta `$ is the London penetration depth. We find that at intermediate times $`R_1/ctR_0/c`$ ($`c`$ is the speed of light), the exponent $`D(t)`$ grows as $`t`$ for $`t\mathrm{}/k_BT`$ and as $`t^2`$ for $`t\mathrm{}/k_BT`$, up to some logarithms. At large times, $`tR_0/c`$, the exponent saturates at a certain limiting value, for which we obtain the following estimate:

$$D_{\mathrm{lim}}\frac{f(u)}{16\alpha _{\mathrm{EM}}\mathrm{ln}^4(R_0/R_1)}\left(\frac{\pi c}{R_0\omega _{\mathrm{min}}}\right)\left(\frac{\delta }{R_1}\right)^2,$$
(5)
where $`\alpha _{\mathrm{EM}}`$ is the fine-structure constant, and $`u=k_BT/\mathrm{}\omega _{\mathrm{min}}`$ with $`\omega _{\mathrm{min}}c/R_0`$; $`f(u)=1`$ for small $`u`$, and $`f(u)=u`$ for large $`u`$.

The main result is the macroscopic suppression factor $`(\delta /R_1)^2`$ in (5). Taking for estimates $`\delta 10^5`$ cm, $`R_0=1`$ cm, $`R_1=1`$ mm, $`\omega _{\mathrm{min}}=\pi c/R_010^{11}`$ s<sup>-1</sup>, and replacing the logarithm (which controls logarithmic accuracy) with unity, we obtain $`xT/(1\mathrm{K})`$ and $`D_{\mathrm{lim}}10^7f(x)`$.

These results apply not only to a ring made out of a solid superconductor, but also to a ring coated with a superconducting film, as long as the thickness of the film is much larger than the London depth $`\delta `$. Similar calculations can be done for any type of persistent circular current, of which a supercurrent in a macroscopic ring is but one limiting case. The opposite limit is a circular current in a single atom—e.g. a circular Rydberg state , which has already been proposed as a basis for quantum computation . It would be interesting to do the calculation for that case and see how the decoherence rate changes as the current loop becomes microscopic.

We now outline the main steps leading to estimate (5).

## II The oscillator Hamiltonian

We will only consider decoherence at times

$$tR_1/c,$$
(6)
when, as we will see, the main contribution to the decoherence exponent comes from electromagnetic modes with low frequencies:

$$\omega c/R_1.$$
(7)
For such modes, the response of a superconductor of the London type (the only type we consider) is determined, to a good accuracy, by the London penetration depth $`\delta `$.

To count the oscillator modes, we enclose the ring in a large normalization sphere of radius $`R_0`$, centered at the center of the ring. We define a system of spherical coordinates $`(r,\theta ,\varphi )`$, with angle $`\theta `$ measured form the axis of the ring’s rotational symmetry.

Only the azimuthal component of the vector potential, indeed only the $`\varphi `$-independent modes of it, interact with the current. For these modes, the vector potential takes the form

$$\text{A}(r,\theta ,\varphi ;t)=\text{e}_\varphi (\varphi )A(r,\theta ;t),$$
(8)
where $`\text{e}_\varphi `$ is the azimuthal unit vector. Then,

$$^2\text{A}=\text{e}_\varphi \left(^2\frac{1}{r^2\mathrm{sin}^2\theta }\right)A,$$
(9)
and the wave equation for $`𝐀`$ takes the form

$$\ddot{A}^2A+\left[\frac{1}{r^2\mathrm{sin}^2\theta }+V(r,\theta )\right]A=\frac{4\pi }{c}j.$$
(10)
Here $`j`$ is the constant (in time) supercurrent density, and $`V(r,\theta )`$ is the “potential” that represents the Meissner effect: inside the superconductor

$$V(r,\theta )=1/\delta ^2,$$
(11)
while outside the superconductor both $`j`$ and $`V`$ are zero.

The equation of motion (10) can be obtained, after substitution (8), from the effective Hamiltonian

$$H=d^3x\left(\frac{1}{8\pi c^2}\dot{\text{A}}^2+\frac{1}{8\pi }\text{B}^2+\frac{V}{8\pi }\text{A}^2\frac{1}{c}\text{A}\text{j}\right),$$
(12)
where $`\text{B}=\times \text{A}`$ is the magnetic field, and $`\text{j}=\text{e}_\varphi j`$. Such Hamiltonian formulation neglects any dissipation (absorption) of the electromagnetic field. This approximation is justified in Appendix A. We stress again that decoherence we calculate in this paper is not due to any such dissipation, but is a result of a “measurement” of the supercurrent, performed by fluctuations of the electromagnetic field.

Now, expand the field $`A`$ in a complete set of modes $`f_n(r,\theta )`$:

$$A(r,\theta ;t)=\underset{n}{}f_n(r,\theta )X_n(t).$$
(13)
This expansion defines a set of oscillator coordinates $`X_n`$. The modes $`f_n`$ satisfy the following equation:

$$\left[^2+\frac{1}{r^2\mathrm{sin}^2\theta }+V(r,\theta )\right]f_n=\frac{\omega _n^2}{c^2}f_n,$$
(14)
where $`\omega _n`$ are real eigenfrequencies, together with the continuity conditions at the surface of the superconductor and some boundary condition at the boundary of the normalization volume. They are normalized in that volume by

$$d^3xf_n(r,\theta )f_n^{}(r,\theta )=\delta _{nn^{}}.$$
(15)
A typical low-frequency mode is shown schematically in Fig. 1. Notice that its main support is outside the superconductor.

Substituting the expansion (13) in the Hamiltonian (12) we obtain the Hamiltonian for the oscillators:

$$H=\frac{1}{8\pi c^2}\underset{n}{}\left(\dot{X}_n^2+\omega _n^2X_n^2\right)\frac{1}{c}\underset{n}{}X_nd^3xj(r,\theta )f_n(r,\theta ),$$
(16)
which is of the Caldeira-Leggett  type. The role of a macroscopic variable is played by the total supercurrent $`I`$. Because in our case the total supercurrent has no dynamics of its own, it is not necessary to split each integral in (16) into a product of $`qI`$ and a coupling $`C_n`$. Nevertheless, to conform to the established notation, we will denote these integrals as

$$qC_n=\frac{1}{c}d^3xj(r,\theta )f_n(r,\theta ).$$
(17)
The quantum-mechanical problem specified by the Hamiltonian (16) is exactly solvable, and the overlap (1) is easily calculable. For the decoherence exponent of two states corresponding to two different values $`q=q_1`$ and $`q=q_2`$, with the initial state being a thermal ensemble at $`q=0`$, we obtain

$$D(t)=\frac{2(\mathrm{\Delta }q)^2}{\pi \mathrm{}}_0^{\mathrm{}}𝑑\omega \frac{J(\omega )}{\omega ^2}[1\mathrm{cos}\omega t]\mathrm{coth}(\mathrm{}\beta \omega /2),$$
(18)
where $`\mathrm{\Delta }q=q_1q_2`$, $`\beta =1/k_BT`$, and $`J(\omega )`$ is the spectral density defined as in refs. :

$$J(\omega )=\frac{\pi }{2}\underset{n}{}\frac{C_n^2}{m_n\omega _n}\delta (\omega _n\omega );$$
(19)
according to (16) the masses of all oscillators are

$$m_n=\frac{1}{4\pi c^2}.$$
(20)
Eq. (18) is similar to expressions found in the literature , but we should stress that ours is not a macroscopic tunneling problem: as we pointed out before, in quantum memory, we strive to prevent spontaneous transitions (tunneling or otherwise) between states with different values of $`q`$. These different values now simply label various macroscopically distinct histories.

Expression (18) applies at both zero and finite temperatures. Two limiting cases can be considered: $`k_BT\mathrm{}\omega `$, and $`k_BT\mathrm{}\omega `$, where $`\omega `$ is a typical frequency at which the integral in (18) saturates. In the first case (“low” temperatures),

$$D(t)=\frac{2(\mathrm{\Delta }q)^2}{\pi \mathrm{}}_0^{\mathrm{}}𝑑\omega \frac{J(\omega )}{\omega ^2}[1\mathrm{cos}\omega t];$$
(21)
in the second case (“high” temperatures),

$$D(t)=\frac{4(\mathrm{\Delta }q)^2}{\pi \mathrm{}^2\beta }_0^{\mathrm{}}𝑑\omega \frac{J(\omega )}{\omega ^3}[1\mathrm{cos}\omega t].$$
(22)
To complete the calculation of $`D`$ we now need to find the couplings $`C_n`$, for the low-frequency modes of interest.

## III Decoherence at intermediate times

Our strategy will be as follows. We first compute the spectral density $`J(\omega )`$ for $`\omega `$ in the range

$$c/R_0\omega c/R_1.$$
(23)
The form of the spectral density at these $`\omega `$ determines the decoherence exponent at intermediate times

$$R_1/ctR_0/c.$$
(24)
Then, in the next section, we establish that at smaller $`\omega `$, $`\omega c/R_0`$, $`J(\omega )`$ is proportional to $`\omega ^3`$. According to (18), that means that at $`tR_0/c`$ the exponent $`D(t)`$ does not grow beyond a certain limiting value $`D_{\mathrm{lim}}`$. Finally, we estimate $`D_{\mathrm{lim}}`$.

For each mode of oscillation, we define a wavenumber $`k_n`$ by

$$k_n=\omega _n/c.$$
(25)
At any time in the range (24), the main contribution to the decoherence exponent comes from modes with

$$k_n1/ct.$$
(26)
For counting these modes, the normalization volume can be replaced by a tube, coaxial with the wire, of radius $`R_{\mathrm{norm}}`$ that satisfies

$$ctR_{\mathrm{norm}}R_0.$$
(27)
This new normalization volume is more convenient than the original sphere because the problem now acquires an approximate cylindrical symmetry: we can view the wire as a straight cylinder, with fields subject to periodic boundary conditions at the ends.

This cylindrical problem has a natural set of cylindrical coordinates, $`(\rho ,\phi ,z)`$. The expression (8) for the vector potential can now be approximated by

$$\text{A}(\rho ,\phi ;t)=\text{e}_zA(\rho ,\phi ;t),$$
(28)
and the mode equation (14) by

$$\left[^2+V(\rho )\right]f_n(\rho ,\phi )=\frac{\omega _n^2}{c^2}f_n(\rho ,\phi ).$$
(29)
As seen from (29), the mode functions $`f_n`$ are eigenfunctions of $`i/\phi `$, the angular momentum about the $`z`$ axis. We concentrate on functions with zero angular momentum (i.e. those independent of the angle $`\phi `$), because for other values the coefficients $`C_n`$ in (19) are suppressed by additional powers of $`k_nR_1`$. Thus, we are left with $`f_n`$ that depend only on the radial coordinate $`\rho `$, the distance to the wire’s axis.

By assumption, the wire’s radius $`R_1`$ is much larger than the penetration depth $`\delta `$. Then, the solution to equation (14) inside the superconductor is approximately

$$f_n(\rho )A_n\mathrm{exp}[(\rho R_1)/\delta ].$$
(30)
On the outside, the solution is a combination of Bessel functions; near the boundary, for $`k_n`$ satisfying (26), it reduces to

$$f_n(\rho )A_n+B_n\mathrm{ln}(\rho /R_1).$$
(31)
Eq. (31) already takes into account the continuity of the solution on the surface; the continuity of the first derivative requires $`B_n=A_nR_1/\delta `$, so (31) becomes

$$f_n(\rho )A_n\left(1+\frac{R_1}{\delta }\mathrm{ln}(\rho /R_1)\right).$$
(32)
Notice the enhancement of the coefficient of the logarithm; it comes from matching the first derivatives.

The coefficient $`A_n`$ in (30), (32) is determined from the normalization condition (15). To compute the normalization integral, we need the solution further away from the surface of the superconductor, where the form (32) no longer applies. But (32) does allow us to pick up the coefficients of the Bessel functions. Writing $`\mathrm{ln}(\rho /R_1)=\mathrm{ln}(k_n\rho )\mathrm{ln}(k_nR_1)`$, we observe that the coefficient of $`Y_0(k_n\rho )`$ (the function of the second kind) is proportional to $`A_nR_1/\delta `$, while the coefficient of $`J_0(k_n\rho )`$ (the function of the first kind) is additionally enhanced by $`\mathrm{ln}(k_nR_1)`$.

We will be content with logarithmic accuracy; to such accuracy, away from the surface we have

$$f_n(\rho )=A_n\frac{R_1}{\delta }\mathrm{ln}(k_nR_1)J_0(k_n\rho ).$$
(33)
To count the modes, we use the Dirichlet boundary condition at $`\rho =R_{\mathrm{norm}}`$:

$$f_n(R_{\mathrm{norm}})=0.$$
(34)
Normalizing (33) in a cylinder of radius $`R_{\mathrm{norm}}`$ and length $`L`$, we find

$$A_n=\left(\frac{k_n}{2LR_{\mathrm{norm}}}\right)^{1/2}\frac{\delta }{R_1}\frac{1}{|\mathrm{ln}k_nR_1|}.$$
(35)
Finally, using for $`j`$ the equilibrium current density

$$j(\rho )=\frac{I}{2\pi R_1\delta }\mathrm{exp}[(\rho R_1)/\delta ],$$
(36)
where $`I`$ is the total supercurrent, we obtain for the spectral density

$$q^2J(\omega )=\frac{\pi }{4c^2}\left(\frac{\delta }{R_1}\right)^2\frac{I^2L}{\mathrm{ln}^2(\omega R_1/c)}.$$
(37)
We now specialize to the case when $`q_1=0`$, and $`q_2=I_s`$ is the current corresponding to one quantum of flux through the ring: to logarithmic accuracy

$$I_s=\frac{c\mathrm{\Phi }_0}{2L\mathrm{ln}(L/R_1)},$$
(38)
where $`\mathrm{\Phi }_0=\pi \mathrm{}c/e`$ is the flux quantum. We obtain

$$(\mathrm{\Delta }q)^2J(\omega )=\frac{\pi ^3c^2\mathrm{}^2}{16e^2L\mathrm{ln}^2(L/R_1)}\left(\frac{\delta }{R_1}\right)^2\frac{1}{\mathrm{ln}^2(\omega R_1/c)}$$
(39)
($`e`$ is the absolute value of the electron charge).

With this form of $`J(\omega )`$, the integral in eq. (18) gets the main contribution from $`\omega 1/t`$. So, at $`t\mathrm{}/k_BT`$, we apply eq. (21), to obtain (to logarithmic accuracy)

$$D(t)=\frac{t\pi ^3c^2\mathrm{}}{16e^2L\mathrm{ln}^2(L/R_1)}\left(\frac{\delta }{R_1}\right)^2\frac{1}{\mathrm{ln}^2(ct/R_1)}.$$
(40)
At $`t\mathrm{}/k_BT`$, we apply the “high”-temperature limit (22). Although the integral in (22) in our case is convergent, we impose a low-frequency cutoff at $`\omega =c/R_0`$, a representative value below which eq. (37) is not applicable. This will give us an estimate of the accuracy of the result. We obtain

$$D(t)=\frac{t^2\pi ^2c^2k_BT}{8e^2L\mathrm{ln}^2(L/R_1)}\left(\frac{\delta }{R_1}\right)^2\left[\frac{1}{\mathrm{ln}(ct/R_1)}\frac{1}{\mathrm{ln}(R_0/R_1)}+O\left(\frac{1}{\mathrm{ln}^2(ct/R_1)}\right)\right].$$
(41)
The second term in the bracket is a correction due to the low-frequency cutoff. Under the condition (24), this term is small compared to the first term (the one from $`\omega 1/t`$), but as $`t`$ approaches values of order $`R_0/c`$ it becomes comparable to the first term. Thus, at $`tR_0/c`$ the expansion indicated in (41) breaks down, and we need a different method to establish the form of $`D(t)`$ at large $`t`$.

## IV Decoherence at large times

We now turn to modes with

$$k_n1/R_0.$$
(42)
For counting these, we use the original normalization sphere and the Dirichlet boundary condition

$$f_n()=0.$$
(43)
The modes are normalized in the sphere by the condition (15).

Let us expand each mode function $`f_n`$ in associated Legendre functions $`P_l^1(\mathrm{cos}\theta )`$:

$$f_n(r,\theta )=(2\pi N_l)^{1/2}\underset{l1}{}S_{nl}(r)P_l^1(\mathrm{cos}\theta ),$$
(44)
where $`S_{nl}`$ are some radial functions, and $`N_l`$ is the normalization integral for $`P_l^1`$:

$$N_l=\frac{2l(l+1)}{2l+1}.$$
(45)
Note the absence of $`l=0`$: there is no (nonzero) $`P_0^1`$.

Eq. (15) now translates into the following normalization condition for the radial functions $`S_{nl}`$:

$$\underset{l1}{}_0^{}S_{nl}(r)S_{n^{}l}(r)r^2𝑑r=\delta _{nn^{}},$$
(46)
while the mode equation (14) leads to the following equation for $`S_{nl}`$:

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}S_{nl}\right)+\frac{l(l+1)}{r^2}S_{nl}k_n^2S_{nl}+\underset{l^{}}{}V_{ll^{}}S_{nl^{}}=0.$$
(47)
The “potential” $`V_{ll^{}}`$ here is given by

$$V_{ll^{}}(r)=(N_lN_l^{})^{1/2}_0^\pi V(r,\theta )P_l^1(\mathrm{cos}\theta )P_l^{}^1(\mathrm{cos}\theta )\mathrm{sin}\theta d\theta ,$$
(48)
in terms of the “potential” $`V(r,\theta )`$.

The potential (48) is nonzero only for $`R_{\mathrm{in}}<r<R_{\mathrm{out}}`$, where $`R_{\mathrm{in}}`$ and $`R_{\mathrm{out}}`$ are the inner and outer radii of the ring. In the region $`r<R_{\mathrm{in}}`$, eq. (47) is a free radial equation, in which we can now neglect $`k_n^2`$. A solution regular at $`r=0`$ is then

$$S_{nl}(r)=D_{nl}r^l,$$
(49)
where $`D_{nl}`$ are some constants. In the outside region, $`r>R_{\mathrm{out}}`$, where eq. (47) is also free, the solution is a linear combination of two spherical waves. As long as $`r1/k_n`$, so that $`k_n^2`$ is still negligible, we can write

$$S_{nl}(r)=F_{nl}\left(\frac{r}{R_{\mathrm{out}}}\right)^l+G_{nl}\left(\frac{R_{\mathrm{out}}}{r}\right)^{l+1},$$
(50)
with some constants $`F_{nl}`$ and $`G_{nl}`$. On the other hand, at large distances $`r1/k_n`$, we have (cf. ref. )

$$S_{nl}(r)\frac{F_{nl}(2l+1)!!}{rk_n^{l+1}R_{\mathrm{out}}^l\mathrm{cos}\delta _n}\mathrm{sin}(k_nr\pi l/2+\delta _n),$$
(51)
where $`\delta _n`$ is a phase shift.

Because the smallest possible $`l`$ is $`l=1`$, from the normalization condition (46) we see that

$$F_{nl}=O(k_n^2)$$
(52)
(the phase shift $`\delta _n`$ is $`O(k_n^3)`$). For each $`n`$, the sets of coefficients $`F_{nl}`$, $`G_{nl}`$, and $`D_{nl}`$ are related through scattering on the potential $`V_{ll^{}}`$. Because this scattering occurs at $`rR_01/k_n`$, it is insensitive to $`k_n`$, so $`G_{nl}`$ and $`D_{nl}`$, and therefore also $`f_n`$ on the surface of the wire, are all $`O(k_n^2)`$. Then, according to eq. (19), the spectral density $`J(\omega )`$ is $`O(\omega ^3)`$.

This suppression of the spectral density at small $`\omega `$ can be explained by noting that the large wavelength fluctuations of the field “see” simultaneously two diametrically opposite segments of the wire, with currents that add up to zero. Such fluctuations therefore have only derivative interactions with the current; hence the extra powers of $`k_n`$, or $`\omega `$.

From eq. (18), we now see that the low frequency region $`\omega c/R_0`$ does not significantly contribute to decoherence, and the main contribution at $`tR_0/c`$ comes from modes with $`\omega c/R_0`$. We can estimate decoherence at these large times by using (18) with the expression (39), which correctly describes higher frequencies, and an infrared cutoff at $`\omega _{\mathrm{min}}c/R_0`$. In this way, we arrive at the estimate (5).

## A The role of dissipation

Here we list some estimates for absorption rates of low-frequency electromagnetic field on the surface of a superconductor. The purpose is to show that the corresponding dissipation time is much larger than the timescale $`R_0/c`$, at which the decoherence exponent saturates.

The response kernel $`Q`$, which determines the current induced in the superconductor by the electromagnetic field, is defined by

$$j_{\mathrm{ind}}(\text{x},t)=d^3\text{x}^{}𝑑t^{}Q(\text{x},\text{x}^{};tt^{})A(\text{x}^{},t^{}).$$
(A1)
The low-frequency expansion of the Fourier transform of $`Q`$ (inside the superconductor) is

$$Q(\omega )=\frac{1}{4\pi \delta ^2}i\omega \sigma /c^2+\mathrm{},$$
(A2)
where the first term is the one taken into account in eq. (12), while the second term describes dissipation due to a finite conductivity $`\sigma `$.

We only consider low frequencies,

$$\omega c/R_1,$$
(A3)
and conventional ($`s`$-wave) superconductors with critical temperatures $`T_c`$ of order of a few Kelvin. For $`R_11`$ mm, we have $`\mathrm{}\omega 2\mathrm{\Delta }_0`$, where $`\mathrm{\Delta }_0`$ is the zero-temperature gap. So, at sufficiently low temperatures, a single photon cannot break a Cooper pair, and the dissipation is entirely due to thermally excited quasiparticles. Therefore, we estimate the conductivity as

$$\sigma \sigma _N\mathrm{exp}[\mathrm{\Delta }(T)/k_BT],$$
(A4)
where $`\sigma _N`$ is the normal-state conductivity at $`T`$ around $`T_c`$, and $`\mathrm{\Delta }(T)`$ is the temperature-dependent gap.

Using (A4) with $`\sigma _N10^{18}\mathrm{s}^1`$ (corresponding to $`10^8\mathrm{\Omega }^1\mathrm{m}^1`$ in SI) and $`\delta 10^5`$ cm, we find that for our range of frequencies the second term in (A2) is much smaller than the first. This allows us to expand the surface impedance as

$$\zeta (\omega )=\omega [4\pi c^2Q(\omega )]^{1/2}=i\frac{\omega \delta }{c}\left[1+2\pi i\omega \sigma \delta ^2/c^2+O(\omega ^2)\right].$$
(A5)
The dissipative effect (absorption) is represented by the real part of $`\zeta `$,

$$\zeta _R(\omega )=2\pi \omega ^2\sigma \delta ^3/c^310^6(\omega /10^{11}\mathrm{s}^1)^2\mathrm{exp}[\mathrm{\Delta }(T)/k_BT],$$
(A6)
where the estimate is for the same values of $`\sigma _N`$ and $`\delta `$ as before.

For example, suppose that the device is inside a cavity, and the cavity and the ring are made of the same superconducting material. Then, assuming that radiation does not leak out, dissipation is mainly due to absorption on the walls of the cavity. The rate of dissipation is estimated as

$$\tau ^1\zeta _Rc/R_{\mathrm{cav}},$$
(A7)
where $`R_{\mathrm{cav}}>R_0`$ is the radius of the cavity. The timescale $`\tau `$ is thus $`1/\zeta _R`$ times larger than $`R_{\mathrm{cav}}/c`$, which in turn is larger than $`R_0/c`$, the timescale at which decoherence saturates.