text:
# Generation of ultrabright tunable polarization entanglement without spatial, spectral, or temporal constraints

## Abstract

The need for spatial and spectral filtering in the generation of polarization entanglement is eliminated by combining two coherently-driven type-II spontaneous parametric downconverters. The resulting ultrabright source emits photon pairs that are polarization entangled over the entire spatial cone and spectrum of emission. We detect a flux of $``$12 000 polarization-entangled pairs/s per mW of pump power at 90% quantum-interference visibility, and the source can be temperature tuned for 5 nm around frequency degeneracy. The output state is actively controlled by precisely adjusting the relative phase of the two coherent pumps.

Polarization entanglement has been used to demonstrate a variety of quantum effects from quantum teleportation telep  to quantum cryptographic protocols crypto . The quality of polarization-entangled photons sources can be characterized by their pair flux and the purity of the entangled state they generate kwiat94 ; kwiat95 ; kwiat99 ; shih ; Brilliant . For the existing sources the requirements of high flux and high purity are often in conflict. Consider, for example, type-II spontaneous parametric downconversion (SPDC) in a noncollinearly phase-matched beta-barium borate (BBO) crystal. Here kwiat95  spatial and spectral filtering are necessary to eliminate non-entangled photons that would reduce the purity of the output state. A source of polarization-entangled photons has been proposed kwiat94  and demonstrated shih  in which the outputs of two different SPDC crystals are combined interferometrically. It was recognized that such a setup would generate entangled photons independent of their wavelengths and angles of emission kwiat94 . The two-crystal interferometer, however, did not show the promised high flux and high visibilities shih ; this was attributed to technical difficulties in the alignment.

Our group has investigated the use of a collinearly-propagating geometry and long periodically-poled crystals to simplify alignment and to increase the output flux in both type-I elliott  and type-II SPDC chris . In the case of type-II SPDC in periodically poled potassium titanyl phosphate (PPKTP) chris , we have obtained *post-selected* polarization-entangled photons. However, spatial and spectral filtering are still required to obtain a high-purity entangled state and the post-selection process involves a 3-dB loss. In this Letter we report on a robust implementation of the coherent addition of two SPDC sources based on a single PPKTP crystal. Our scheme fully exploits, for the first time to our knowledge, the properties of interferometric combining of two coherent SPDC sources kwiat94  to yield an ultrabright source of polarization entanglement that has no spatial or spectral constraints. Moreover, collinear operation allows us to control the output state by locking the pump phase of the same interferometer. This setup produces $``$10 times more polarization-entangled pairs/s per mW of pump than any other continuous-wave (cw) source in the literature kwiat99 ; Brilliant .

Figure 1 illustrates a source that coherently combines the outputs of two SPDC crystals. A laser is split by a 50–50 beam splitter (BS) and pumps the two crystals that are phase matched for collinear type-II SPDC. In the low-gain regime, the bi-photon state just after the crystals is given by

$$|\mathrm{\Psi }=\frac{1}{\sqrt{2}}\left(|H_A(\omega _s)V_A(\omega _i)+e^{i\varphi _p}|H_B(\omega _s)V_B(\omega _i)\right),$$
(1)
where $`A`$ and $`B`$ refer to the two arms of the interferometer, $`\omega _s`$ and $`\omega _i`$ are the signal and idler frequencies, respectively, and $`\varphi _p=k_p(L_BL_A)`$ is the difference of the delays accumulated by the pump (with wave vector $`k_p`$) in the paths $`L_A`$ and $`L_B`$ between the 50–50 BS and the crystals. A half-wave plate (HWP) is used to rotate the polarizations by 90 in arm $`B`$, so that the output state after the polarizing beam splitter (PBS) is

$$|\mathrm{\Psi }=\frac{1}{\sqrt{2}}\left(|H_1(\omega _s)V_2(\omega _i)+e^{i\varphi }|V_1(\omega _s)H_2(\omega _i)\right),$$
(2)
where 1 and 2 refer to the two PBS output ports. The overall phase $`\varphi =\varphi _p+\varphi _s+\varphi _i`$ is determined by the pump phase $`\varphi _p`$ and the phase delays accumulated by the signal and idler, respectively, with $`\varphi _{s,i}=k_{s,i}(L_B^{}L_A^{})\mathrm{\Delta }\varphi _{\lambda /2}(\omega _{s,i})`$. The first term of $`\varphi _{s,i}`$ is the delay due to the arm lengths $`L_A^{}`$ and $`L_B^{}`$ between the crystals and the PBS, and the second term is the phase difference introduced by the HWP. Note that the phase delays introduced by the identical crystals in the two arms cancel. Under collinear phase matching $`k_p=k_s+k_i`$, and $`\varphi `$ is equal to the phase difference accumulated by the pump in the Mach-Zehnder interferometer formed between the 50–50 BS and the PBS except for a fixed offset due to the HWP. The phase of the output bi-photon state in Eq. 2 can therefore be precisely controlled by locking the Mach-Zehnder interferometer as seen by the pump alone: one can generate the triplet (for $`\varphi =0`$) or the singlet state ($`\varphi =\pi `$), as well as intermediate states that are linear combinations thereof.

The HWP constrains the signal (idler) fields from the two crystals to exit at output 1 (2) in Fig. 1, ensuring kwiat94 ; shapiro00  that the two sources are indistinguishable so that all the photons are polarization entangled regardless of their wavelengths and directions of emission. Spatial and spectral filtering is unnecessary in this two-crystal configuration, thus promising a source that has a much higher photon-pair flux, plus a larger bandwidth and spatial extension than BBO sources. Due to energy and momentum conservation, one expects the emitted photon pairs from this broadband spatially-extended source to show spectral and spatial entanglement. Additional advantages of this scheme include automatic erasure of timing distinguishability, non-degenerate operation, and source tunability.

To implement the interferometric source described above it is crucial that the two SPDC sources be identical. Source differences introduce an element of distinguishability between the two paths that would lead to a mixed state output. We therefore implemented the scheme based on a single crystal with counter-propagating pump beams derived from a single laser. The single-crystal approach is particularly useful with periodically-poled crystals, as it mitigates imperfections in the profiles of the periodic gratings.

We used a 10-mm-long ($`X`$ crystallographic axis), 1-mm-thick ($`Z`$ axis), and 4-mm-wide ($`Y`$ axis) hydrothermally-grown PPKTP crystal with a grating period of 9.0 $`\mu `$m. At a temperature of $`32^{}`$C this grating period phase matches type-II collinear degenerate downconversion of a 398.5-nm pump polarized along the $`Y`$ axis and propagating along the crystal’s $`X`$ axis. The crystal was housed in an oven and was maintained at its operating temperature with $`\pm `$0.1C precision. This crystal was previously characterized and used in type-II collinear SPDC to yield single-beam quantum interference with a 99% visibility chris . We used second harmonic generation to measure the temperature and wavelength tuning behavior in PPKTP using a cw tunable laser centered around 797 nm. The second-harmonic measurements are well described by the Sellmeier phase-matching equations for PPKTP sellmeier , which allow us to calculate the spatial and spectral properties of the downconverted photons, as well as the phase-matching angles’ dependence on the crystal temperature. The latter predictions have been verified by imaging the emitted photons with a CCD camera and narrow spectral filters.

The experimental setup is shown in Fig. 2. The frequency-doubled cw Ti:sapphire pump laser at 398.5 nm was split by a BS that had a splitting ratio of $``$50–50. To balance the powers of the two pump beams we inserted a half-wave plate (HWP1) to vary the horizontally polarized pump power in the (counter-clockwise propagating) brighter path. The crystal was not phase matched for a vertically polarized pump. Each pump beam focussed to a waist of $``$150 $`\mu `$m at the center of the PPKTP crystal. The generated beams were collimated with 300-mm radius-of-curvature dichroic mirrors (DM1,2) and combined at a PBS after the polarization of one of the beams was rotated by $`90^{}`$ with a HWP. The dichroic mirrors were coated for high reflectivity (HR) at 797 nm and for high transmission (HT) at 398.5 nm, with a residual reflectivity of 0.2% at the pump wavelength. The pump beams, which propagated collinearly with the downconverted beams, were weakly reflected by the four mirrors (DM1-4) and recombined on the PBS, which had a $``$20–80 splitting ratio at the pump wavelength. The resultant pump beam from port 2 of the PBS was directed by a dichroic mirror (DM5, HR at 398.5 nm and HT at 797 nm) for detection with an ultraviolet photodiode (UV-PD). The BS and PBS in Fig. 2 formed a Mach-Zehnder interferometer for the pump and the detected fringes were used to stabilize the interferometer with a side-locking technique. This provided a convenient and robust way to control the phase $`\varphi `$ of the output state in Eq. 2. By inserting a dispersive medium (such as a thin glass plate) in one of the arms of the Mach-Zehnder interferometer we introduced a fixed but variable offset between the phase of the pump fringes and the phase of the output state (the overall offset phase includes other dispersive elements in the interferometer). By varying this phase offset, we were able to lock the phase of the output state at an arbitrary value while optimizing the side-locking feedback signal.

We placed two irises in the output beam paths to control the acceptance angle of the detection system. We estimate that an iris diameter of 1 mm corresponded to an internal emission solid angle of $``$3.5 $`\times 10^5`$ sr at the crystal. Flat dichroic mirrors (not shown in Fig. 2) similar to DM1 were used to eliminate residual pump light. The output photons were detected with single-photon Si detectors (Perkin-Elmer SPCM-AQR-14) through polarization analyzers (composed of a half-wave plate and a polarizer). The outputs of the single-photon detectors were counted and also sent to an AND gate (TTL logic family 74F) for coincidence counting. The coincidence window for this configuration was measured to be $`39.4`$ ns. This parameter allowed us to correct for the rate of accidental coincidences in all of the data reported. For example, when 12 000 coincidences/s were measured an average of 67 000 singles/s were detected at each single-photon detector and $`250`$ coincidences/s were due to accidental Poisson processes.

A summary of our experimental results, with the accidentals removed, is shown in Figs. 35 for $`\varphi =\pi `$ (singlet). The temperature of the crystal was set to $``$32C to ensure frequency degenerate operation. Figure 3 shows the quantum interference in the coincidence counts when the analyzer angle in arm 2 was varied for a fixed angle in arm 1 with no narrowband interference filter. We observed a visibility of 100 $`\pm `$ 3% (85 $`\pm `$ 3%) when analyzer 1 is set to 0 (45). In what follows we will use the 45-visibility as an indication of the quality of the state generated.

In Fig. 4 we report the 45-visibility for the singlet state as a function of the iris diameter. Two sets of data are shown, one using a 3-nm interference filters centered at 797 nm placed in front of the detectors (diamonds) and one in which the interference filter was removed (squares). In both cases the visibility is almost constant as a function of the iris diameter. This allows us to increase the pair flux (Fig. 4 inset) while preserving the purity of the output state. With the 3-nm filter we observed a visibility of 90% and a flux $``$ 12 000 pairs s-1 mW-1 with a 4-mm iris. Under this condition, following Ref. Aspect , we tested Bell’s inequality and obtained $`S=2.599\pm 0.004`$, violating the classical limit by more than 100 $`\sigma `$.

Figure 4 can be compared with data obtained in a single-pass configuration with similar collection geometry properties reported in Ref. chris . The visibility of quantum interference in the single-pass experiment drops much faster as the iris diameter increases than in this interferometric configuration. The nearly constant visibility in Fig. 4 arises from effective spatial and spectral indistinguishability in this dual-pumped interferometric configuration.

Two main factors limited the visibility: wavefront distortion and diffraction caused by the components of the interferometer, and defects in the electric-field poling of the crystal. Wavefront distortion and diffraction lead to spatial distinguishability between the two downconverted beams. Inhomogeneity in the crystal grating introduces a temporal mismatch between the two paths. Both these effects were mitigated somewhat by closing the iris and by adding spectral filters. To investigate the effects of wavefront distortion and diffraction we measured the interferometer visibility directly by injecting a laser beam at 797 nm through arm 2 of the PBS in Fig. 2 and observing the fringe signal in arm 1. The input and output beam diameters could be varied with irises. When we changed the diameter of the output beam for a fixed input beam diameter, the visibility showed the same plateau for small iris diameters as in Fig. 4. When we decreased the input beam diameter for a fixed output beam diameter the visibility increased linearly, approaching 100%. This suggests that the diffraction inside the interferometer was responsible for the flat plateau in the visibility of Fig. 4. We note that a slight mismatch in the length of the two interferometric arms can also degrade the visibility.

No interference filter was used in obtaining the data shown in Fig. 5 and the iris diameter was fixed at 2.2 mm. The temperature of the crystal was then scanned between 20 and 50C and the 45-visibility was measured. We used our knowledge of the Sellmeier equations to calculate the phase-matched signal and idler wavelengths for each temperature setting, and hence obtain the abscissas shown in this figure. Figure 5 shows that the 45-visibility is essentially independent of the signal and idler emission wavelengths for a range of $``$5 nm around degeneracy.

In conclusion we have demonstrated a source of polarization-entangled photon pairs with high flux and state purity. The cw source is based on the interferometric combination of two coherently-driven type-II sources of spontaneous parametric downconversion from a single PPKTP crystal. This dual-pumped source is uniquely characterized by the fact that all the emitted photon pairs are polarization entangled, regardless of their wavelengths and directions of emission. Therefore it can be tuned, has a wide bandwidth and an extended spatial profile. We believe that our source produces spatial and spectral entanglement, in addition to polarization entanglement, thus providing additional degrees of freedom that can be used for quantum communication. Further work with this source is needed to experimentally demonstrate these additional forms of entanglement. If successful, we would then have a source that could be used to demonstrate fundamental quantum properties zeilinger  and in cryptographic protocols with improved security genovese .

This work was supported by a DoD Multidisciplinary University Research Initiative (MURI) program administered by the Army Research Office under grant DAAD-19-00-1-0177, the Defense Advanced Research Projects Agency (DARPA) and Air Force Research Laboratory under agreement F30602-01-2-0546, and the National Reconnaissance Office.


format: text

text:
# 1 Introduction

## 1 Introduction

In the theoretical studies of biomolecular systems, it is now indispensable to use molecular simulation techniques such as Monte Carlo (MC) and molecular dynamics (MD) methods. Most such simulations use the potential energy functions within the framework of classical mechanics consisting of certain energy terms with force-field parameters. Well-known force fields (or potential energy functions) are AMBER , CHARMM , OPLS , GROMOS , and ECEPP . These force fields have been parameterized to fit experimental data of small molecules and, for some terms, quantum chemistry calculations. In fact, it is rather difficult to test the validity of these force-field parameters, because simulations of biomoleuclar systems will be hampered by the multiple-minima problem; the simulations will get trapped in states of energy local minima and will yield unreliable results. Hence, it is essential to employ powerful simulation algorithms such as generalized-ensemble algorithms (for a recent review, see Ref. ) in order to have accurate comparisons of force fields. For instance, we have recently carried out detailed comparisons of three versions (parm94 , parm96 , and parm99 ) of AMBER, CHARMM , OPLS-AA/L , and GROMOS  by generalized-ensemble simulations of small peptides .

If any problems are found in the force fields, we have to improve the force-field parameters by some means. In this Letter, we propose a novel method for optimization of force-field parameters. Our method utilizes the protein structures in the Protein Data Bank (PDB) , in which the number of entries is rapidly increasing. Actually, there already exist many works of similar knowledge-based methods of force-field development . Our method is different from most of previous works in two points: we deal with all-atom force fields such as AMBER, and we use only the PDB data without introducing so-called decoys (or, misfolded structures).

In section 2 the details of the formulation of our method are described. In section 3 an example of the application of our method is presented. Conclusions and future prospects are presented in section 4.

## 2 Methods

The existing all-atom force fields for protein systems such as AMBER and CHARMM use essentially the same functional forms for the conformational potential energy $`E_{\mathrm{conf}}`$ except for minor differences. $`E_{\mathrm{conf}}`$ can be written as, for instance,

$$E_{\mathrm{conf}}=E_{\mathrm{BL}}+E_{\mathrm{BA}}+E_{\mathrm{torsion}}+E_{\mathrm{nonbond}},$$
(1)
$$E_{\mathrm{BL}}=\underset{\mathrm{bond}\mathrm{length}\mathrm{}}{}K_{\mathrm{}}(\mathrm{}\mathrm{}_{\mathrm{eq}})^2,$$
(2)
$$E_{\mathrm{BA}}=\underset{\mathrm{bond}\mathrm{angle}\theta }{}K_\theta (\theta \theta _{\mathrm{eq}})^2,$$
(3)
$$E_{\mathrm{torsion}}=\underset{\mathrm{dihedral}\mathrm{angle}\mathrm{\Phi }}{}\underset{n}{}\frac{V_n}{2}[1+\mathrm{cos}(n\mathrm{\Phi }\gamma _n)],$$
(4)
$$E_{\mathrm{nonbond}}=\underset{i<j}{}\left[\frac{A_{ij}}{r_{ij}^{12}}\frac{B_{ij}}{r_{ij}^6}+\frac{q_iq_j}{ϵr_{ij}}\right].$$
(5)
Here, $`E_{\mathrm{BL}}`$, $`E_{\mathrm{BA}}`$, and $`E_{\mathrm{torsion}}`$ represent the bond-stretching term, the bond-bending term, and the torsion-energy term, respectively. The bond-stretching and bond-bending energies are given by harmonic terms with the force constants $`K_{\mathrm{}}`$ and $`K_\theta `$, and the equilibrium positions, $`\mathrm{}_{\mathrm{eq}}`$ and $`\theta _{\mathrm{eq}}`$. The torsion energy is, on the other hand, described by the Fourier series in Eq. (4), where the sum is taken over all dihedral angles $`\mathrm{\Phi }`$, $`n`$ is the number of waves, $`\gamma _n`$ is the phase, and $`V_n`$ is the Fourier coefficient. The nonbonded energy in Eq. (5) is represented by the Lennard-Jones and Coulomb terms between pairs of atoms, $`i`$ and $`j`$, separated by the distance $`r_{ij}`$. The parameters $`A_{ij}`$ and $`B_{ij}`$ in Eq. (5) are the coefficients for the Lennard-Jones term, $`q_i`$ is the partial charge of the $`i`$-th atom in the electrostatic term, and $`ϵ`$ is the dielectric constant, where we set $`ϵ=1`$. Hence, we have five classes of force-field parameters, namely, those in the bond-stretching term ($`K_{\mathrm{}}`$ and $`\mathrm{}_{\mathrm{eq}}`$), those in the bond-bending term ($`K_\theta `$ and $`\theta _{\mathrm{eq}}`$), those in the torsion term ($`V_n`$ and $`\gamma _n`$), those in the Lennard-Jones term ($`A_{ij}`$ and $`B_{ij}`$), and those in the electrostatic term ($`q_i`$).

We now describe our new method for optimizing these force-field parameters. We first select $`N`$ molecules in PDB. We try to choose proteins from different folds (such as all $`\alpha `$-helix, all $`\beta `$-sheet, $`\alpha `$/$`\beta `$, etc.) and different homology classes as much as possible. If the force-field parameters are of ideal values, we expect that all the chosen native structures are stable without any force acting on each atom in the molecules. Hence, we expect

$$F=0,$$
(6)
where

$$F=\underset{m=1}{\overset{N}{}}\frac{1}{N_m}\underset{i_m=1}{\overset{N_m}{}}\left|\stackrel{}{f}_{i_m}\right|^2,$$
(7)
and

$$\stackrel{}{f}_{i_m}=\frac{E_{\mathrm{tot}}}{\stackrel{}{x}_{i_m}}.$$
(8)
Here, $`N_m`$ is the total number of atoms in molecule $`m`$, $`E_{\mathrm{tot}}`$ is the total potential energy, and $`\stackrel{}{f}_i`$ is the force acting on atom $`i`$. In reality, $`F0`$, and because $`F0`$, we expect that we can optimize the force-field parameters by minimizing $`F`$ with respect to these parameters. In practice, we perform a simulation in the force-field parameter space for this minimization.

Proteins are usually in aqueous solution, and hence we also have to incorporate some kind of solvent effects. Because the more the total number of proteins ($`N`$) is, the better the force-field parameter optimizations are expected to be, we want to minimize our efforts in the calculations of the solvent effects. Here, we employ the generalized-Born/surface area (GB/SA) terms for the solvent contributions . Hence, we use in Eq. (8)

$$E_{\mathrm{tot}}=E_{\mathrm{conf}}+E_{\mathrm{solv}},$$
(9)
where

$$E_{\mathrm{solv}}=E_{\mathrm{GB}}+E_{\mathrm{SA}},$$
(10)
$$E_{\mathrm{GB}}=166\left(1\frac{1}{ϵ_s}\right)\underset{i,j}{}\frac{q_iq_j}{\sqrt{r_{ij}^2+\alpha _{ij}^2e^{D_{ij}}}},$$
(11)
$$E_{\mathrm{SA}}=\underset{k}{}\sigma _kA_k.$$
(12)
Namely, in the GB/SA model, the total solvation free energy in Eq. (10) is given by the sum of a solute-solvent electrostatic polarization term, a solvent-solvent cavity term, and a solute-solvent van der Waals term. A solute-solvent electrostatic polarization term can be calculated by the generalized Born equation (11), where $`\alpha _{ij}=\sqrt{\alpha _i\alpha _j}`$, $`\alpha _i`$ is the so-called Born radius of atom $`i`$, $`D_{ij}=r_{ij}^2/(2\alpha _{ij})^2`$, and $`ϵ_s`$ is the dielectric constant of bulk water (we take $`ϵ_s=78.3`$). A solvent-solvent cavity term and a solute-solvent van der Waals term can be approximated by the term that is proportional to the solvent accessible surface area in Eq. (12). Here, $`A_k`$ is the total solvent-accessible surface area of atoms of type $`k`$ and $`\sigma _k`$ is an empirically determined proportionality constant .

The flowchart of our method for the optimization of force-field parameters is shown in Fig. 1.

In Step 1 we try to obtain as many structures as possible from PDB. The number is limited by the computer power that we have available in our laboratory. We want to choose proteins with different sizes (numbers of amino acids), different folds, and different homology classes as much as possible. We also want to use only those with high experimental resolutions. Note that only atomic coordinates of proteins are extracted from PDB (and coordinates from other molecules such as crystal water are neglected).

If we use data from X-ray experiments, hydrogen atoms are missing, and thus in Step 2 we have to add hydrogen coordinates. Many protein simulation software packages provide with routines that add hydrogen atoms to the PDB coordinates, and one can use one of such routines.

We now have $`N`$ protein coordinates ready, but usually such “raw data” result in very high total potential energy and strong forces will be acting on some of the atoms in the molecules. This is because the hydrogen coordinates we added as above are not based on experimental results and have rather large uncertainties. The coordinates of heavy atoms from PDB also have experimental uncertainties. We take the position that we leave the coordinates of heavy atoms as they are in PDB as much as possible, and adjust the hydrogen coordinates to adjust this mismatch. This is why we want to include as many PDB data as possible with high experimental resolution (so that the effects of experimental errors in PDB may be minimal). We thus minimize the total potential energy $`E_{\mathrm{tot}}=E_{\mathrm{conf}}+E_{\mathrm{solv}}+E_{\mathrm{constr}}`$ with respect to the coordinates for each protein conformation, where $`E_{\mathrm{constr}}`$ is the constraint energy term that is imposed on the heavy atoms in PDB (it is referred to as the “predefined constraints” in Steps 3 and 5 in Fig. 1):

$$E_{\mathrm{constr}}=\underset{\mathrm{heavy}\mathrm{atom}}{}K_x(\stackrel{}{x}\stackrel{}{x}_0)^2.$$
(13)
Here, $`K_x`$ is the force constants of the restriction (in the present work, we use the value of 100 kcal/mol/Å2), and $`\stackrel{}{x}_0`$ are the original coordinate vectors of heavy atoms in PDB. Because we are searching for the nearest local-minimum states, usual minimization routines such as conjugate-gradient method and Newton-Raphson method can be employed here. As one can see from Eq. (13), the coordinates of hydrogen atoms will be mainly adjusted, but unnatural heavy-atom coordinates will also be modified. We perform this minimization for all $`N`$ protein structures separately, and obtain $`N`$ refined structures.

Given $`N`$ set of “ideal” reference coordinates in Step 3, we now optimize the first set of force-field parameters in Step 4. In Eq. (1) we have five classes of force-field parameters as mentioned above. Namely, the force-field parameters are those in the bond-stretching term ($`K_{\mathrm{}}`$ and $`\mathrm{}_{\mathrm{eq}}`$), those in the bond-bending term ($`K_\theta `$ and $`\theta _{\mathrm{eq}}`$), those in the torsion term ($`V_n`$ and $`\gamma _n`$), those in the Lennard-Jones term ($`A_{ij}`$ and $`B_{ij}`$), and those in the electrostatic term ($`q_i`$). Because they are of very different nature, we believe that it is better to optimize these classes of force-field parameters separately (as in Steps 4, 6, and so on in Fig. 1). For each set of force-field parameters, the optimization is carried out by minimizing $`F`$ in Eq. (7) with respect to these parameters. Here, $`E_{\mathrm{tot}}`$ in Eq. (8) is given by Eq. (9) For this purpose usual minimization routines such as conjugate-gradient method are not adequate, because we need a global optimization. One can employ more powerful methods such as simulated annealing  and generalized-ensemble algorithms . We perform this minimization simulation in the above parameter space to obtain the parameter values that give the global minimum of $`F`$.

These processes are repeated until the optimized force-field parameters converge.

We can, in principle, optimize all the force-field parameters following the flowchart in Fig. 1. In the example given in the next section, however, we just optimize two classes of the force-field parameters; namely, the partial charges ($`q_i`$) and the backbone torsion energy parameters ($`V_1`$, $`V_2`$, and $`V_3`$). Here, we fix the phases ($`\gamma _n`$) and the main-chain torsion energy term is written as (see Eq. (4))

$$E_{\mathrm{\Phi }=\varphi ,\psi }=\frac{V_1}{2}\left[1+\mathrm{cos}(\mathrm{\Phi })\right]+\frac{V_2}{2}\left[1+\mathrm{cos}(2\mathrm{\Phi }\pi )\right]+\frac{V_3}{2}\left[1+\mathrm{cos}(3\mathrm{\Phi })\right],$$
(14)
for each backbone dihedral angle $`\varphi `$ and $`\psi `$. As for the partial charge optimization, we impose a condition that the total charge of each amino acid remains constant, which is the usual assumption adopted by force fields of Eqs. (2)–(5) based on classical mechanics.

We believe that these two classes of parameters have the most uncertainty among all the force-field parameters. This is because partial charges are usually obtained by quantum chemistry calculations of an isolated amino acid in vacuum separately, which is a very different condition from that in amino acids of proteins in aqueous solution and because torsion-energy term is the most problematic (for instance, parm94, parm96, and parm99 versions of AMBER differ mainly in torsion-energy parameters).

## 3 Results and discussion

We now present an example of the application of our force-field optimization scheme presented in the previous section. The force field that we optimized is the AMBER parm96 version . We have optimized two sets of parameters. The first set is the partial charge parameters ($`q_i`$ in Eq. (5)) and there are 602 such parameters altogether in AMBER. The second set is the backbone torsion-energy parameters ($`V_1`$, $`V_2`$, and $`V_3`$ in Eq. (14)) and there are six such parameters (three each for $`\varphi `$ and $`\psi `$). We used the program package TINKER  for all the calculations in the present work.

In Step 1 of the flowchart of Fig. 1, we chose 100 PDB files ($`N=100`$) with resolution 1.8 Å or better and with less than 200 residues (the average number of resiudes is 120.4) from PISCES . Their PDB code names are 2LIS, 1EP0, 1TIF, 1EB6, 1C1L, 1CCW, 2PTH, 1I6W, 1DBF, 1KPF, 1LRI, 1AAP, 1C75, 1CC8, 1FK5, 1KQR, 1K1E, 1CZP, 1GP0, 1KOI, 1IQZ, 3EBX, 1I40, 1EJG, 1AMM, 1I07, 1GK8, 1GVP, 1M4I, 1EYV, 1E29, 1I2T, 1VCC, 1FM0, 1EXR, 1GUT, 1H4X, 1GBS, 1B0B, 119L, 1IFC, 1DLW, 1EAJ, 1GGZ, 1JR8, 1RB9, 1VAP, 1JZG, 1M55, 1EN2, 1C9O, 2ERL, 1EMV, 1F41, 1EW6, 2TNF, 1IFR, 1JSE, 1KAF, 1HZT, 1HQK, 1FXL, 1BKR, 1ID0, 1LQV, 1G2R, 1KR7, 1QTN, 1D4O, 1EAZ, 2CY3, 1UGI, 1IJV, 3VUB, 1BZP, 1JYR, 1DZK, 1QFT, 1UTG, 2CPG, 1I6W, 1C7K, 1I8O, 1LO7, 1LNI, 1EQO, 1NDD, 1HD2, 3PYP, 1FD3, 1DK8, 1WHI, 1FAZ, 4FGF, 2MHR, 1JB3, 2MCM, 1IGD, 1C5E, and 1JIG.

The minimization for the coordinate refinement in Step 3 of the flowchart was done with the constraint of Eq. (13). In Fig. 2(a), a part of the original PDB structure of one of the molecules (2LIS) and the corresponding part in the refined structure are superposed. As expected, we see good coincidence of heavy-atom coordinates and small deviations in the hydrogen coordinates. RMSD (Root Mean Square Distance) (for heavy atoms) of these two conformations of the entire molecule (2LIS) is indeed small: 0.027 Å (the average of all 100 proteins is 0.03 Å).

As also expected, the structure of the raw PDB data (that obtained in Step 2 of the flowchart) is not stable in the sense that the force acting on some atoms in the molecule will be large, while the refined structure (that obtained in Step 3 of the flowchart) will not have abnormally large forces acting. This is illustrated in Fig. 2(b). The average absolute value of the components of the force before minimization and after minimization is 9.0 kcal/mol/Å and 1.5 kcal/mol/Å, respectively.

In Step 4 of the flowchart, we performed the optimization of partial-charge parameters by MC simulated annealing. Namely, we minimized $`F`$ in Eq. (7) by MC simulated annealing simulations of these parameters (the parameters are updated and the updates are accepted or rejected according to the Metropolis criterion). For this we introduced an effective “temperature” for the parameter space. Each simulation run consisted of 50,000 MC sweeps with the temperature decreased exponentially from 200 to 0.01. The simulation was repeated 10 times with different initial random numbers. The time series of $`F`$ from one of the simulations is shown in Fig. 3(a). We see that $`F`$ decreases quickly in the beginning until about 5,000 MC sweeps and then it decreases very slowly; the total number of MC sweeps (50,000) seems sufficient. The optimized partial charges are taken from those that resulted in the lowest $`F`$ value.

Each term contributing to $`F`$ (i.e., each component of the force acting on the atoms) before and after the optimization is compared in Fig. 2(c). We see that many terms decreased in magnitude as a whole.

In Tables 1, 2, and 3, three examples (alanine, glutamic acid, and tyrosine) of the obtained partial charges together with the original AMBER values are listed. The magnitudes of the charges seem to decrease a little in general, although there are exceptions (see, for instance, CG and CE of tyrosine). Although the sign of the partial charges remains the same for those with large magnitude, charges with small magnitude sometimes change their sign (see, for example, CG of glutamic acid and CA, CB, and HB of tyrosine).

In Step 5 of the flowchart, the original coordinates obtained in Step 2 were again minimized with the constraint of Eq. (13), but this time the optimized parameters were used. The average RMSD of 100 proteins is 0.03 Å, and the coordinates of heavy atoms have little changed.

In Step 6 of the flowchart, we carried out the optimization of the six torsion-energy parameters ($`V_1`$, $`V_2`$, and $`V_3`$ in Eq. (14) for both $`\varphi `$ and $`\psi `$) by minimizing $`F`$ in Eq. (7) with MC simulated annealing simulations in this parameter space. Each simulation run consisted of 10,000 MC sweeps with the temperature decreasing from 1,000 to 1.0. The simulation was repeated several times with different random numbers, but the six optimized parameters all agreed with one another. The time series of $`F`$ from one of the simulations is shown in Fig. 3(b), and the obtained torsion-energy parameters are listed in Table 4. Each term contributing to $`F`$ before and after the optimization is compared in Fig. 2(d). The term correspoding to the backbone atoms that are relevant to the torsion-energy parameters decreased, as expected. Non-zero contributions are presumably to be reduced by optimization of other force-field parameters.

In the present example, we stopped our process in Step 6 of the flowchart and did not iterate the optimizations.

We tested the effectiveness of the new force field by applying it to the folding simulations of two peptides, C-peptide of ribonuclease A and the C-terminal fragment of the B1 domain of streptococcal protein G, which is sometimes referred to as G-peptide . The C-peptide has 13 residues and its amino-acid sequence is Lys-Glu-Thr-Ala-Ala-Ala-Lys-Phe-Glu-Arg-Gln-His-Met. This peptide has been extensively studied by experiments and is known to form an $`\alpha `$-helix structure . Because the charges at peptide termini are known to affect helix stability , we blocked the termini by a neutral COCH3\- group and a neutral -NH2 group. The G-peptide has 16 residues and its amino-acid sequence is Gly-Glu-Trp-Thr-Tyr-Asp-Asp-Ala-Thr-Lys-Thr-Phe-Thr-Val-Thr-Glu. The termini were kept as the usual zwitter ionic states, following the experimental conditions . This peptide is known to form a $`\beta `$-hairpin structure by experiments .

Simulated annealing MD simulations of 1 ns were performed for both peptides from extended initial conformations. The unit time step was set to 1.0 fs (hence, each run consists of 1,000,000 MD steps). The temperature during MD simulations was controlled by Berendsen’s method . For each simulation the temperature was decreased from 2,000 K to 200 K. As for solvent effects, in order to carry out quantitative test, we want to include accurate contributions such as explicit water molecules. Here, for simplicity, however, we used the GB/SA model in Eqs. (10)–(12) was used. For both peptides, these folding simulations were repeated 16 times with different initial conformations. As for force fields, both the original AMBER parm96 force field and the optimized force field were used.

In Fig. 4 the lowest-energy conformations of C-peptide obtained from each of the folding simulations are shown. They are numbered in increasing order of energy. Almost all conformations are $`\alpha `$-helix structure in the case of the original AMBER parm96 (see Fig. 4(a)), although a short $`\beta `$-hairpin structure is also found in two conformations (Nos. 12 and 16). RMSD (only Cα atoms are taken into account) from the corresponding part of the native X-ray structure of the entire ribonulcease A was the lowest for conformations Nos. 3 (2.47 Å), 14 (2.59 Å), and 15 (2.81 Å). For the optimized force field (see Fig. 4(b)), five out of 16 conformations (Nos. 2, 7, 8, 9, and 11) are $`\alpha `$-helix structure, and $`\beta `$-hairpin structure also appeared in Nos. 12, 14, and 15. Though both $`\alpha `$-helix and $`\beta `$-hairpin structures appear with the optimized force field, $`\alpha `$-helix is more favored than $`\beta `$-hairpin in the sense that the energy is lower for conformations with $`\alpha `$-helix than those with $`\beta `$-hairpin. RMSD was the lowest for conformations Nos. 13 (3.17 Å), 16 (3.50 Å), and 11 (3.51 Å). Although the original AMBER parm96 gave more conformations closer to the native structure (smaller RMSD) than the optimized force field, it seems to favor $`\alpha `$-helix too much. The experiments imply only 30 % $`\alpha `$-helix formations around $`0^{}`$ C , which is more consistent with the results for the optimized force field.

In Fig. 5 the lowest-energy conformations of G-peptide obtained from the folding simulations are shown. The results for the original AMBER parm96 (see Fig. 5(a)), are again almost all $`\alpha `$-helical. This is in contradiction with the experimental results  which imply the formation of a $`\beta `$-hairpin. RMSD from the corresponding part of the native X-ray structure of the entire protein G was the lowest for conformations Nos. 14 (4.52 Å), 13 (5.30 Å), and 8 (5.53 Å). For the optimized force field (see Fig. 5(b)), as many as eight out of 16 conformations (Nos. 3, 4, 5, 8, 10, 11, 14, and 15) yielded $`\beta `$-hairpin structure, and $`\alpha `$-helix structure also appeared in Nos. 5 and 14. RMSD was the lowest for conformations Nos. 7 (3.74 Å), 4 (3.84 Å), and 11 (4.00 Å). Hence, for the case of G-peptide, the improvement of the force field by our optimization procedure is great (both in secondary structure content and in RMSD).

## 4 Conclusions

In this Letter, we proposed a novel method to optimize the force-field parameters for protein systems. Our method is knowledge-based in the sense that we utilize the protein structures in the Protein Data Bank. With this method one can optimize not only the existing force-field parameters such as AMBER, CHARMM, OPLS, and so on, but also any force field as long as its functional form is given.

We have presented an example of the application of the present method to the original AMBER parm96 force field by optimizing the partial-charge and backbone torsion-energy parameters. There is still a lot of testing to be done. For instance, other parameter sets in AMBER parm96 as well as other force fields should also be tried. We used only 100 protein structures, and this number should be increased. As for solvent effects in the optimization procedure, we employed the generalized-Born/surface area model, but the validity of this approximation should be tested. Namely, we have to check how insensitive or sensitive the final optimized parameters are to the solvation theory used. Finally, because the parameter space we optimize may be large, we want to use a more powerful optimization algorithm than simulated annealing. Generalized-ensemble algorithm  may be one choice of such optimization methods.

Acknowledgements:
The computations were performed on the computers at the Research Center for Computational Science, Okazaki National Research Institutes and at ITBL, Japan Atomic Energy Research Institute. This work was supported, in part, by the NAREGI Nanoscience Project, Ministry of Education, Culture, Sports, Science and Technology, Japan.


text:
# Dynamo-Driven Outflows in Pre-Planetary Nebulae

## 1. Observations Suggest a Role for Magneto-Rotational Outflows

Powerful anisotropic outflows appear in a wide range of astrophysical sources and magneto-rotationally driven models are often a favored explanation. Maybe asymmetric outflows (e.g. Balick & Frank 2003) in planetary nebulae (PPNe) and pre-planetary nebulae (PPNe) also involve dynamically important magnetic fields (Aller 1958). Magneto-rotational driving is particularly appealing for PPNe winds which are more powerful than PNe winds: The observationally inferred rates of mechanical momentum and energy deposition for PNe are (Kwok 2000) $`\dot{\mathrm{\Pi }}10^{27}(M_{pn}/0.04M_{})(v_{pn}/40\mathrm{k}\mathrm{m}/\mathrm{s})(\mathrm{\Delta }t_{PN}/10^4\mathrm{yr})^1`$ erg/cm, and $`\dot{E}10^{34}(\dot{\mathrm{\Pi }}/10^{27}\mathrm{erg}/\mathrm{cm})(v_{pn}/40\mathrm{k}\mathrm{m}/\mathrm{s})`$ erg/s, where $`M_{pn}`$ is the swept up mass in the PN, $`v_{pn}`$ is the speed of this mass, and the acceleration time $`\mathrm{\Delta }t_{PN}`$ age of the PN. But observations of PPNe (Bujarrabal et al. 2001) reveal values up to $`\dot{\mathrm{\Pi }}\text{ }<3\times 10^{29}(M_{ppn}/0.5M_{})(v_{ppn}/200\mathrm{k}\mathrm{m}/\mathrm{s})(\mathrm{\Delta }t_{ppn}/10^3\mathrm{yr})^1`$ erg/cm, and $`\dot{E}\text{ }<10^{37}(\mathrm{\Pi }/10^{40}\mathrm{erg}/\mathrm{cm})(v_{ppn}/200\mathrm{k}\mathrm{m}/\mathrm{s})\mathrm{erg}/\mathrm{s}`$ where the acceleration time $`\mathrm{\Delta }t_{ppn}<<`$ age of the PPN. The PPNe wind momenta are often $`10^3`$ times that carried by radiation even when multiple scattering is considered.

A large scale magnetic field can act as a drive belt that extracts rotational energy to drive a wind from the engine where the field lines are anchored (e.g. Spruit 1996). When the engine is a star, the extraction spins the star down. When the engine is an accretion disk, a fraction of the accretion luminosity gets redirected into the wind. In the Blandford & Payne (1982) type model, mass outflow is driven by the toroidal magnetic field pressure gradient along poloidal field lines, though in the rotator’s co-rotating frame near the footpoints, the processes can be thought of as centrifugal launching along favorably inclined poloidal field lines. Consideration of magnetically mediated outflows raises natural questions in the PPNe/PNe context: Where do the required fields come from? How do these mechanisms work for single star or binary/disk systems? Are the outflow powers and momenta consistent with the observations?

## 2.  Why Would Dynamos be Needed to Produce the Fields?

Convection driven and magneto-rotationally driven MHD turbulence (Balbus & Hawley 1998) are fully developed in in AGB stars and in accretion disks (the two relevant environments) respectively. Large scale magnetic fields can then experience turbulent diffusion so intense that the field is not flux-frozen on the relevant dynamical or advection time scales (Lubow et al. 1994; Blackman & Tan 2003 (BT)). This highlights the need for dynamo amplification. The fields also be large enough to avoid diffusing before escaping to the corona where they are needed. In coronae, the fields can further open up. Coronal loops which carry magneto-centrifugal winds should be at least as large as the radial scale of the anchoring rotator so that material can be accelerated to super-Alfvénic velocities before reaching the loop-tops and avoid sliding back down.

## 3. Some Recent Developments in Large Scale Dynamo Theory

Dynamos amplify or sustain magnetic fields by draining energy from random and shear motion. Nonhelical (direct) dynamos amplify fields at wave numbers $`kk_f`$ whereas helical (inverse) dynamos can also amplify fields at $`k<k_f`$, where $`k_f`$ is the smallest wavenumber at which the flow is turbulent. Here “helical” refers to whether the driving turbulence possesses a non-vanishing pseudoscalar helicity such as kinetic helicity $`𝐯\times 𝐯`$, where $`𝐯`$ is the turbulent velocity, and the brackets imply averaging within a hemisphere. The nonhelical dynamo increases magnetic energy as the turbulence stretches magnetic field lines in a nearly random walk, but it is the helical dynamo (e.g. “$`\alpha _d\mathrm{\Omega }`$,” see Moffatt 1978) which produce the large scale fields needed for outflows. Helical turbulence can be supplied by convection in a star, or by the magneto-rotational instability in a stratified disk (Brandenburg et al. 1995; Brandenburg 1998; BT).

While 3-D MHD simulations are essential for understanding dynamos, they are not always practical for modeling. For the latter, mean field theory (e.g. Moffatt 1978) is used. Here the field is broken into mean and fluctuating components. The helical dynamo growth coefficients become correlations of fluctuating quantities, and the mean field evolution is solved for. Traditional treatments do not dynamically determine the saturation strength of the large scale field; they either take the growth to be steady (=kinematic theory), or presume a form of quenching by the growing magnetic field without conserving magnetic helicity, a fundamental invariant of ideal MHD. A recent theory which alleviates these limitations by including the dynamical evolution of the velocity and magnetic helicity (Blackman & Field 2002) is illustrated below.

Fig. 1 below (from Blackman & Brandenburg 2003) shows the helical $`\alpha _d\mathrm{\Omega }`$ dynamo in a rotating northern hemisphere. The kinematic theory is shown in (a) and (b) while (c), (d) and (e) show the dynamic theory. (a): To conserve angular momentum, rising (falling) blobs in a stratified medium twist in the opposite (same) sense to the underlying rotation. This implies a non-vanishing $`𝐯\times 𝐯`$, averaged in a given hemisphere. This “$`\alpha _d`$” effect implies that a large scale toroidal loop in the northern hemisphere incurs a right-handed writhe, and a radial field component as it rises. (b): Differential rotation (the “$`\mathrm{\Omega }`$” effect) at the base of the loop shears the radial component, amplifying the toroidal field. Dissipation or ballooning/opening of the top part of the loop allows for a net flux gain through the rectangle. (c): Same as (a) but now with the field represented as a ribbon. The right-handed writhe of the large scale loop is accompanied by a left-handed twist along the tube, incorporating magnetic helicity conservation. (d): Same as (b) but with field represented as ribbon/tube. (e): Top view of the combined twist and writhe The backreaction force which resists bending is the small scale twist. Diffusing the top part of the loops allows a net flux generation in the rectangles of (a)-(d), and helps alleviate the backreaction.

The dynamic role of magnetic helicity is seen in simulations (Brandenburg 2001) and modeled analytically (Field & Blackman 2002; Blackman & Field 2002). Although the small scale magnetic helicity growth drives the helical dynamo to a steady state, the large scale field grows large enough for most astrophysical applications. Application of these principles to disks and stars with boundaries is work in progress. In addition, buoyancy rather than kinetic helicity may be the initial driver of the helical $`\alpha _d`$ effect in disks (Brandenburg 1998). This affects the initial sign of $`\alpha _d`$ in each hemisphere, but magnetic helicity evolution remains a key to understanding the saturation.

## 4. Single Star Scenario: AGB Dynamo and Magnetic Explosion

This scenario (Blackman et al. 2001 (BFMTV)) begins with a main sequence (MS) star $`3M_{}`$ of rotational speed $`200`$ km/s. At the end the MS lifetime, the core collapses and the envelope expands, and angular momentum is assumed to be conserved on spherical mass shells. The rotation profile satisfies $`\mathrm{\Omega }1/r^2`$ which is then combined with a numerical AGB stellar model (SDK1401, Steve Kawalar, personal communication). In SDK1401, the convection zone extends from the outer stellar radius $`5\times 10^{12}`$cm down to the core radius $`r_c=9\times 10^{10}`$cm. At the “interface” between the core and the envelope, $`\mathrm{\Omega }2\times 10^5`$ cm/s.

The combination of a strong differential rotation layer and overlying convection zone may lead to an “interface dynamo” (e.g. Markiel & Thomas 1999). Here the poloidal field is generated by a helical $`\alpha _d`$-effect layer in the convection zone. Diffusive transport (or turbulent pumping; Tobias et al. 1998) bleeds the poloidal field to the shear layer beneath the convection zone where it is further amplified by linear stretching (“$`\mathrm{\Omega }`$” effect). BFMTV applied the code of Markiel & Thomas (1999) to the AGB star, (The $`\alpha _d`$ formula used was consistent with magnetic helicity conservation.) At $`rr_c`$ the saturated field is $`B5\times 10^4`$ G.

The scenario requires this dynamo to be steady for $`10^6`$ yr until the envelope of the AGB star is radiatively driven away in the late AGB stages. Eventually, the density falls such that the field can no longer be stored: Turbulent pumping can be gauged by assuming it to be proportional to a fixed convective energy density transfer rate $`\rho v^3`$. As $`\rho `$ decreases, $`v`$ increases and the turbulent pressure, $`\rho v^2`$, decreases with decreasing $`\rho `$, eventually dropping below $`B^2/8\pi `$ at the core. At this point, a “magnetic explosion” associated with toroidal pressure could ensue and drive the PPNe outflow. The magnetic explosion is driven by Poynting flux, and surface integration reveals a maximum mechanical luminosity (Blackman, Frank, Welch 2001 (BFW)) $`L_m10^{37}\left(\frac{B}{5\times 10^4\mathrm{G}}\right)^2\left(\frac{\mathrm{\Omega }_c}{10^5/\mathrm{s}}\right)\left(\frac{r_c}{r_{}}\right)^3\mathrm{erg}/\mathrm{s},`$ where $`r_c`$ and $`\mathrm{\Omega }_c`$ are the core radius and angular speed respectively. This high $`L_m`$ lasts for a magnetic spin down time of the star ($`100`$ yr, starting after the envelope is lost) and then decays exponentially. The asymptotic outflow speed should approach $`v_{\mathrm{}}\mathrm{\Omega }_0r_A`$, where $`r_A`$ is the Alfvén radius, which can be less than or comparable to the escape speed $`v_{esc}520\left(r_c/10^{11}\mathrm{cm}\right)^{1/2}\left(M_{}/M_{}\right)^{1/2}`$ km/s, depending on $`r_A`$. Matt et al. (2003) (these proceedings) have simulated a magnetic explosion. The above results are consistent with the requirements of PPNe (see Sec. 1). (Imai et al. (2002) have observed a collimated AGB jet.)

Unfortunately, the present interface dynamo models, while including a saturation of the $`\alpha _d`$ effect, do not include saturation of the $`\mathrm{\Omega }`$ effect. In reality, amplification of the magnetic field at the shear layer drains shear energy, and transfers angular momentum from the core to the envelope, slowing down the core within $`100`$ yr, unless the differential rotation is re-seeded throughout the $`10^6`$ yr lifetime of the AGB phase. This is important because matching the PPNe observations requires the field to be anchored in a rapidly rotating core.

The so called “Lambda effect” (c.f. Rüdiger 1989) might be able to re-seed the differential rotation. This mechanism involves the interplay between the overall rotation and asymmetric convection to produce a steady-state rotation profile. The helical interface dynamo has to be solved with a dynamically coupled Lambda effect. This awaits attention, but note that there is enough energy in the fusion driven turbulent convection to sustain the differential rotation: The turbulent energy density can be calculated from from the aforementioned AGB stellar model SDK1401 and the result is (J. Nordhaus personal comm.) $`5\times 10^{43}`$ erg. This can be compared with the magnetic field energy in the shear layer from BFMTV, which is $`E_m10^{41}`$ erg, which in turn is about 1/8 the shear energy in the interface layer. Thus only $`2`$% of the energy from the convection needs to be steadily drained into the desired differential rotation.

## 5. Binary Scenarios: Disk Dynamos and/or Primary Core Spin-up

Binary scenarios do not suffer from the sustained shear problem just discussed because they deliver the angular momentum near the end of the AGB phase (e.g. Soker 1998), just when needed to drive the outflows.

Such angular momentum transfer could occur in the late AGB stage by common envelope (CE) evolution (e.g. Iben & Livio 1993), as the secondary star penetrates the distended AGB envelope and transfers orbital angular momentum to it. The AGB star can shed $`80\%`$ of its mass within years of the onset of CE, depending on the binary parameters. During the shedding, the companion spirals nearer to the core. While accretion of envelope material can occur onto the secondary (e.g. Morris 1987; Soker & Livio 1994; Mastrodemos & Morris 1998), much larger (needed) accretion rates arise when an accretion disk forms around the primary via disruption of the secondary (Reyes-Ruiz & López 1999). Favored systems for this mode of accretion involve an evolved AGB star with mass $`2.6M_{}/M_{}3.6`$, a secondary with mass ($`0.1M_{}`$), and the initial binary separation $`200R_{}`$. Reyes-Ruiz & López (1999) find that the disk accretion rate onto the residual primary core evolves in time (after an initial $`1`$yr viscous adjustment period) as $`\dot{M}_d=\dot{M}_{do}(t/\mathrm{yr})^{5/4}\mathrm{M}_{}\mathrm{yr}^1.`$ Typical values of the initial accretion rate $`\dot{M}_{do}10^3\dot{M}`$/yr.

Once the disk forms, a helical disk dynamo (e.g. Pudritz 1981; Brandenburg et al. 1995; Reyes-Ruiz & Stepinksi 1995; Campbell 2000; BT) provides the mean field strengths from which one can determine the magnetic luminosity $`L_m`$, available to drive an outflow. This can be estimated as the integrated Poynting flux flowing through the disk surface (or, the Poynting flux through the Alfvén surface (BFW; Frank & Blackman 2003 (FB)). Tan & Blackman (2003) find $`L_m\overline{B}_p\overline{B}_\varphi \mathrm{\Omega }_ir_i^3\frac{\alpha _{ss}^{\frac{1}{2}}G\dot{M}_dM_{}}{r_i}10^{37}\left(\frac{\alpha _{ss}}{0.1}\right)^{\frac{1}{2}}\left(\frac{M_{}}{M_{}}\right)\left(\frac{r_i}{10^{11}\mathrm{cm}}\right)^1\left(\frac{t}{1\mathrm{y}\mathrm{r}}\right)^{\frac{5}{4}}\frac{\mathrm{erg}}{\mathrm{s}}`$ where $`\alpha _{ss}\text{ }<0.1`$ is the assumed disk viscosity parameter (Shakura & Sunyaev 1979), $`r_i`$ is the inner radius, and $`c_s`$ and $`\mathrm{\Omega }_i`$ are the sound and rotation speeds there. The $`\dot{M}_d`$ enters because $`\rho \dot{M}_d`$, and the mean surface poloidal and toroidal fields satisfy $`\overline{B}_p\rho ^{1/2}\alpha _{ss}c_s`$ and $`\overline{B}_\varphi \rho ^{1/2}\alpha _{ss}^{1/2}c_s`$. A calculation of acceleration along a poloidal field line gives the asymptotic velocity $`v_{\mathrm{}}\mathrm{\Omega }_ir_A2^{1/2}\mathrm{\Omega }_ir_i730\left(M_{}/M_{}\right)^{\frac{1}{2}}\left(r_i/10^{11}\mathrm{cm}\right)^{\frac{1}{2}}`$ km/s (FB). For a rotationally supported disk, the asymptotic wind velocity is always near or greater than the escape speed. The mechanical luminosity and momentum delivered are consistent with the PPNe requirements of Sec. 1.

Rather than form a disk as just described, an alternative is that the companion could spin up both a slowly rotating core and envelope, leading to significant differential rotation between the two (assuming a turbulent viscosity in the envelope). This could rejuvenate the interface dynamo of Sec. 4 during CE in the late AGB stage just when needed. The time scales of spin-orbit synchronization (Zahn 1977), envelope ejection, and viscous dissipation (Iben & Livio 1993) must be compared, but this alternative might mean that a wider range of binary systems could produce magneto-rotational outflows.

## 6. Open Questions

Asymmetric outflows in PPNe and PNe signature the transport of angular momentum, and magnetically mediated outflows may produce the required power and momenta. Maser (Miranda et al. 2001; Vlemmings et al. 2003) and core X-ray (Kastner et al. 2003) observations loosely support this general paradigm.

However, key unresolved issues remain: (1) Is there a need for collimated outflows from both the AGB star and disk? Are both self-collimated or does the latter collimate the former? What might misaligned nested winds produce (BFW)? (2) Do all PPN/PNe produce collimated jets when the sources are young? (3) Large scale nonlinear dynamo theory is in its infancy for realistic boundary terms. How do helical dynamo fields open up and dynamically relax in the coronae? (4) How far is momentum carried as Poynting flux? (e.g. hydromagnetic “fling” models vs. magnetically dominated “spring” models, or nested hybrid models? Ustyugova et al. 2000). (5) MHD jet simulations have been separate from MHD dynamo simulations. Simulations capturing both together have not been done, though using a dynamo produced field and then simulating the outflow provides encouraging results (von Rekowski et al. 2003). Non-steady calculations are important in this context. (6) The single AGB wind scenario will only work for PPNe if convection + rotation can steadily re-seed the differential rotation to ensure a large enough Poynting flux at the end of the AGB phase. If not, binaries may be absolutely required for asymmetric PNe and PPNe (e.g. Soker & Livio 1994; Soker 2001). Stellar evolution models which include magnetic fields and rotation are needed. (7) To what extent are the strong outflows in PPNe correlated with the presence of binaries? Binaries with brown dwarfs or planets are not easily detected, but are favored for accretion around the AGB core. (8) How restrictive are the initial binary parameters for which common envelope evolution simulations (Demarco et al. 2003) would produce a suitably accreting disk to power PPNe? Is this relaxed by the scenario at the end of Sec. 5? (9) Consequences of internal dynamo produced fields for surface X-ray emission in AGB stars and disks need more study. (10) Could winds from the primary vs. secondary disks be distinguished by abundance differences (FB)?

## References

Aller L.H. 1958, AJ, 63, 47

Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1

Balick, B. & Frank, A. 2002, ARA&A, 40, 439

Blackman, E.G.& Brandenburg, A. 2003, ApJL, 584, L99

Blackman, E.G.& Field, G.B. 2002, Phys. Rev. Lett., 89, 265007

Blackman EG, Frank A, Welch C. 2001, ApJ, 546, 288 (BFW)

Blackman E.G., Frank A., Markiel J.A., Thomas J.H., Van Horn H.M. 2001, Nat. 409, 485 (BFMTV)

Blackman, E.G., & Tan J.C. 2003, in Proc. of the workshop on “Magnetic Fields and Star Formation” Madrid, 2003 (Dodrehcht: Kluwer) eds. A.I. Gómez de Castro et al., astro-ph/0307455 (BT)

Blandford, R.D. & Payne, D.G. 1982, MNRAS, 199, 883.

Brandenburg, A., Nordlund, A., Stein, R.F., & Torkelsson, U. 1995, ApJ, 446, 741

Brandenburg, A, 1998, in Theory of Black Hole Accretion Disks eds. M.A. Abramowicz, G. Bjornsson, & J.E. Pringle (Cambridge: CUP) p.61

Brandenburg, A. 2001, ApJ, 550, 824

Bujarrabal V, Castro-Carrizo A, Alcolea J, Sánchez Contreras C. 2001, A & A 377, 868

Campbell, C.G. 2000, MNRAS, 317, 501

Field, G.B. & Blackman, E.G. 2002, ApJ, 572, 685

Frank, A. & Blackman E.G. 2003, in press ApJ (FB)

Demarco O., Sandquist E.L., Mac Low M.-M., Herwig F., Taam R.E. 2003, in press Rev. Mex Ast.

Iben, I, Livio M. 1993, PASP, 105, 1373.

Imai, H., Obara, K., Diamond, P. J., Omodaka, T., & Sasao, T. 2002, Nat., 417, 829

Kastner, J. H., Balick, B., Blackman, E. G., Frank, A., Soker, N., Vrtílek, S. D., & Li, J. 2003, ApJ, 591, L37

Knapp, G.R. 1986, ApJ 311, 731

Kwok, S. 2000, Origin and Evolution of Planetary Nebulae (Cambridge: CUP)

Lubow, S.H., Papaloizou, J.C.B., & Pringle, J.E. 1994, MNRAS, 267, 235

Lynden-Bell, D. 1996, MNRAS, 279, 389

Markiel, J.A. & Thomas, J.H. 1999, ApJ, 523, 827

Mastrodemos N, Morris M. 1998, ApJ, 497, 303

Miranda, L.F., Gómez, Y., Anglada, G., & Torrelles, J.M. 2001, Nat., 414, 284

Moffatt, H. K. 1978, Magnetic Field Generation in Electronically Conducting Fluids (Cambridge, UK: CUP)

Pelletier, G. & Pudritz, R.E. 1992, ApJ, 394, 117

Pudritz, R.E. 1981, MNRAS, 195, 897

Reyes-Ruiz M. & Stepinksi T.F. 1995, ApJ, 438 750.

Reyes-Ruiz M, López J.A. 1999, ApJ 524, 952

Rüdiger, G. 1989, Differential rotation and stellar convection  (Berlin: Verlag)

Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337

Soker, N. 1998, ApJ, 496, 833

Soker N. 2001, ApJ, 558, 157

Soker N, Livio M. 1994, ApJ, 421, 219

Soker N, Rappaport S. 2001, ApJ, 557, 256

Soker, N. & Zoabi, E. 2002, MNRAS, 329, 204

Spruit, H. C. 1996, NATO ASIC Proc. 477, 249

Tan J.C. & Blackman E.G. 2003, sub. to ApJ, astro-ph/0307455

Tobias, S.M., Brummell, N.H., Clune, T.L., & Toomre, J. 1998, ApJL, 502, L177

Ustyugova, G.V., Lovelace, R.V.E., Romanova, M.M., Li, H., & Colgate, S.A. 2000, ApJ, 541, L21

Vlemmings W.H.T., Diamond P.J. van Langevelde H.J., in to appear in the PASP proc. of the VLBA 10th Anniversary Meeting astro-ph/0309185

von Rekowski, B., Brandenburg, A., Dobler, W., & Shukurov, A. 2003, A&A, 398, 825

Zahn J.-P. 1977, A & A, 57, 383


text:
# Fluxon dynamics by microwave surface resistance measurements in MgB2

## 1 Introduction

Since the discovery of superconductivity in the MgB2 compound, several authors have investigated the properties of such superconductor in the presence of magnetic fields. It has been shown that the electric, magnetic and transport properties of MgB2 are strongly affected by the magnetic field . Such enhanced field dependence of the properties of MgB2 has been ascribed to the peculiar two-gap structure, with a larger gap, associated with the two-dimensional $`\sigma `$ band, and a smaller gap, associated with the three-dimensional $`\pi `$ band, which is rapidly suppressed by the applied magnetic field .

Measurements of the microwave (mw) surface resistance, $`R_s`$, allow investigating the processes responsible for mw energy losses in superconductors as well as determining specific properties of the samples. It is well known that the mw surface resistance of superconductors in the mixed state depends on the magnetic field through several mechanisms . The different vortex states, in the different regions of the H-T plane, determine the temperature and field dependencies of $`R_s`$. Therefore, measurements of $`R_s(H,T)`$ may provide important information on the fluxon dynamics in the different regimes of fluxon motion.

It has been shown that mw losses in MgB2 superconductor are strongly affected by the magnetic field in the whole range of temperatures below $`T_c`$, even for relatively low field values . In particular, an unusually enhanced field dependence of $`R_s`$ has been observed at $`TT_c`$, where pinning effects should hinder the energy losses.

In this paper, we report experimental results of the field-induced variations of the mw surface resistance in high-density ceramic MgB2. Below $`T=0.95T_c`$, we have observed a magnetic hysteresis in the $`R_s(H)`$ curves, whose peculiarities cannot straightforwardly be justified by models reported in the literature. On the contrary, the experimental results obtained near $`T_c`$ are well accounted for in the framework of the Coffey and Clem model , with fluxons moving in the flux-flow regime, provided that the anisotropy of the upper critical field is taken into due account.

## 2 Experimental

The field-induced variations of $`R_s`$ have been studied in a bulk ceramic MgB2 sample of approximate dimensions 2$`\times `$3$`\times `$0.3 $`mm^3`$, with $`T_c`$ 39 K. The sample has been extracted from a high-density pellet (2.4 g/cm3) obtained by reactive infiltration of liquid Mg on a powdered B preform . After the reaction in a sealed stainless steel container, lined with a Nb foil, a thermal treatment has been performed for two hours at $`T`$ 900 C.

The mw surface resistance has been measured using the cavity perturbation technique. The cavity, of cylindrical shape with golden-plated walls, is tuned in the TE011 mode resonant at 9.6 GHz ($`Q`$-factor $``$ 40,000 at $`T=4.2`$ K). The sample is put in the center of the cavity by a sapphire rod, in the region of maximum mw magnetic field. The cavity is placed between the poles of an electromagnet, which generates DC magnetic fields up to $``$ 10 kOe. Two additional coils allow reducing to zero the residual field and working at low magnetic fields. The quality factor of the cavity has been measured by an hp$`8719`$D Network Analyzer. $`R_s`$ has been investigated as a function of the DC magnetic field, $`H_0`$, in the range 0 $`÷`$ 10 kOe, at fixed values of temperature. Before each measurement was performed, the sample was zero-field cooled (ZFC) to the desired value of temperature; the external field was then increased up to 10 kOe and successively decreased to zero, at constant temperature. The sample and field geometry is shown in Fig.1 a); in this geometry the mw current induces a tilt motion of all the fluxons present in the sample, as shown in Fig.1 b).

Fig.2 shows the field dependence of $`R_s`$, normalized to the normal-state surface resistance, $`R_n`$, measured at $`T`$ = 40 K, for different values of the temperature. Open and solid symbols refer to results obtained on increasing and decreasing $`H_0`$, respectively. As one can see, up to $``$ 4 K below $`T_c`$ the $`R_s(H)`$ curves show a magnetic hysteresis. The inset shows the temperature dependence of the height of the hysteresis loop at $`H_0`$ = 0, $`\mathrm{}R_s(0)/R_n`$, i.e. the variation of $`R_s/R_n`$ after a complete cycle of the DC magnetic field from 0 to 10 kOe and back. On increasing the temperature up to $`T`$ 34 K, $`\mathrm{}R_s(0)/R_n`$ takes a constant value; on further increasing the temperature it quickly decreases; the hysteretic behavior vanishes at temperatures very close to $`T_c`$.

Both the increasing and decreasing-field branches of the $`R_s(H)`$ curves show a plateau in the field ranges 0 $`÷`$ $`H_p`$ and 0 $`÷`$ $`H^{}`$, respectively. Since the sample was ZFC, $`H_p`$ marks the first-penetration field. Fig.3 a) shows $`H_p`$ as function of temperature; the line has been obtained by fitting the experimental data with the law $`H_p(T)=H_p(0)[1(T/T_c)^\beta ]`$; we have obtained, as best-fit parameters, $`H_p(0)`$ = 390 $`\pm `$ 20 Oe and $`\beta `$ = 2.4 $`\pm `$ 0.2. The $`\beta `$ value gives a temperature dependence of $`H_p`$ consistent with that of $`H_{c1}(T)`$. $`H_p(0)`$ is slightly larger than the values of $`H_{c1}(0)`$ reported in the literature for MgB2, suggesting that in our sample surface-barrier effects are weak.

Fig.3 b) shows the temperature dependence of $`H^{}`$, though the meaning of this characteristic parameter is not clear, one can note that its temperature dependence is very similar to that of $`H_p`$.

## 3 Discussion

Microwave losses induced by DC magnetic fields in superconductors in the mixed state are proportional to the complex penetration depth of the mw field, $`\stackrel{~}{\lambda }`$; they are influenced by the fluxon motion and the very presence of vortices which bring along normal fluid in their core . $`\stackrel{~}{\lambda }`$ can be expressed in terms of the normal-fluid density, flux-flow resistivity and $`\nu /\nu _0`$ ratio, where $`\nu `$ is the working frequency and $`\nu _0=\alpha _L/\eta `$ the depinning frequency, with $`\alpha _L`$ the Labush parameter and $`\eta `$ the viscous-drag coefficient. For $`\nu \nu _0`$ the induced mw current makes fluxons moving in the flux-flow regime. Otherwise, the fluxon motion is ruled by the strength of the restoring-pinning forces.

Our results show that for $`T<`$ 35 K the $`R_s(H)`$ curves exhibit an hysteretic behavior, suggesting that pinning effects are important in this range of temperatures. On the contrary, for $`T>`$ 35 K the hysteresis vanishes and fluxons should move in the flux-flow regime. Actually, the results near $`T_c`$ can be accounted for by using the Coffey and Clem theory , with fluxons moving in the flux-flow regime, provided that the anisotropy of the upper critical field is taken into due account. Following the same procedure of Ref., based on the assumption that the anisotropic Ginzburg-Landau theory is valid for MgB2 in this range of temperatures, we have fitted the experimental data by taking the anisotropy factor, $`\gamma `$, as fitting parameter. The best-fit curves, along with the experimental data, are shown in Fig.4. For all the curves, the best-fit value of the anisotropy factor is $`\gamma `$ = 3, showing that in our sample the anisotropy of $`H_{c2}`$ is constant in a range of temperatures of about 3 K below $`T_c`$.

The hysteresis of the $`R_s(H)`$ curves, observed for $`T<`$ 35 K, should be related to the different value of the induction field at increasing and decreasing fields. Both the normal fluid density and the flux-flow resistivity depend on the number of fluxons moving under the action of the mw current. In the field geometry we used the mw current induces a tilt motion of all the fluxons present in the sample (see Fig.1). So, the hysteresis of $`R_s(H)`$ should be directly related to the hysteretic behavior of the magnetization curve due to pinning effects. However, several anomalies can be highlighted in the $`R_s(H)`$ curves at low temperatures. Firstly, from Fig.2 one can see that an unusually enhanced field variation of $`R_s`$ is observed, even at the lowest temperatures; a magnetic field of 10 kOe induces an $`R_s`$ variation of about 1/3 of the normal-state value. We suggest that such enhanced field-induced energy losses are ascribable to the effects of the two-gap structure of MgB2. Indeed, it has been shown that the structure of vortices in MgB2 is characterized by two coherence lengths, associated with each band . On increasing the applied field, the giant cores start to overlap at field values much smaller than $`H_{c2}`$. Therefore, modest fields can suppress superconductivity between the vortices, resulting in an enhancement of the normal-fluid density . Furthermore, this particular vortex structure can modify the interaction among vortices and pinning centers, giving rise to an unconventional shape of the pinning potential well.

Another anomaly concerns the shape of the decreasing-field branch of the $`R_s(H)`$ hysteresis loop. By supposing fluxons in the critical state, after the external field starts decreasing the variation of $`B(H)`$ should be initially slow and then faster, no matter the field dependence of the critical current. So, a downward curvature of the decreasing-field branch of the $`R_s(H)`$ curve is expected, at variance with our results. The origin of this anomalous behavior, as well as the presence of the wide plateau in the decreasing-field $`R_s(H)`$ curves, is not yet understood and needs further investigation.

In conclusion, we have reported a detailed study of the field-induced variations of the mw surface resistance in high-density ceramic MgB2. An unusual behavior of the $`R_s(H)`$ curves has been highlighted at low temperatures, which cannot be explained in the framework of models reported in the literature. We suggest that it is ascribable to the unconventional vortex structure related to the presence of the two gaps. On the contrary, the results obtained at temperatures close to $`T_c`$ have been accounted for in the framework of the Coffey and Clem model, with fluxons moving in the flux-flow regime, by taking into account the anisotropy of upper critical field.