text:
# Spin/Orbital Pattern-Dependent Polaron Absorption in Nd1-xSrxMnO3

\[

## Abstract

We investigated optical properties of Nd1-xSrxMnO3 (x= 0.40, 0.50, 0.55, and 0.65) single crystals. In the spin/orbital disordered state, their conductivity spectra look quite similar, and the strength of the mid-infrared absorption peak is proportional to x($`1`$x) consistent with the polaron picture. As temperature lowers, the Nd1-xSrxMnO3 samples enter into various spin/orbital ordered states, whose optical responses are quite different. These optical responses can be explained by the spin/orbital ordering pattern-dependent polaron hopping.

\]

Recently, lots of attention has been paid to the roles of orbital degree of freedom in the physical properties of manganites. Resonant x-ray scattering of a layered manganite provided evidences of orbital ordering. Polarization microscope images demonstrated the existence of optical anisotropy , which might be due to the orbital ordering. A new elementary excitation caused by symmetry breaking of orbital ordering, called ‘orbiton’, was proposed and claimed to be observed in LaMnO3 by Raman spectroscopy measurements. Since the corresponding Raman peaks were observed in the energy region of a few tenth of eV, similar excitations could be observed by infrared spectroscopy.

Indeed, an anomalous incoherent absorption peak has been observed in the mid-infrared (IR) region for numerous manganites. However, the origin of this mid-IR absorption peak is not fully understood yet. In one scenario, it was suggested that the mid-IR absorption should be due to orbital excitations. Even though the observed anisotropy might be explained by this model, there have been few experimental reports to compare with the theoretical predictions even on a qualitative level. In the other scenario, it has been suggested that the mid-IR absorption should arise from hopping of the Jahn-Teller (JT) polaron . The polaron binding energy obtained from the mid-IR peak position agrees with that from other measurements. However, this scenario has not been extended to explain the observed optical anisotropy. Till now, there is no consensus on these two scenarios. Therefore, it is important to understand the role of the orbital degree of freedom on the mid-IR absorption peaks.

Nd1-xSrxMnO3 (NSMO) shows various types of spin/orbital (S/O) ordering patterns at low temperature ($`T`$), depending on the value of the doping concentration (x). Optical investigations on NSMO in various S/O ordered states can provide us further insights on the role of S/O for the mid-IR absorption in manganites. In spite of such importance, few investigations have been reported due to the difficulties in growing high quality single-domain crystals.

In this Letter, we present optical responses of NSMO (x=0.40, 0.50, 0.55, and 0.65) multi-domain single crystals and a single-domain crystal with x=0.55. At room $`T`$, optical conductivity spectra $`\sigma (\omega )`$ of all the samples appear quite similar and can be explained within the JT polaron picture. At low $`T`$, by taking account of the effects of the S/O ordering pattern on the polaron hopping motion, the anisotropy in $`\sigma (\omega )`$ for the single-domain sample can be explained quite well. This S/O pattern-dependent hopping model was successfully extended to quantitatively explain the reflectivity spectra $`R(\omega )`$ of other NSMO multi-domain crystals. This clearly demonstrates that the polaron scenario is indeed valid, and that the polaron hopping motion should be strongly affected by the S/O ordering patterns.

The NSMO single crystals were grown by the floating zone method. After polishing and proper annealing, $`R(\omega )`$ were measured by near normal incident unpolarized light over a wide photon energy region 5 meV$``$30 eV. For the single-domain x=0.55 crystal, the crystal axes were determined using the neutron scattering measurement. Using polarized light, we obtained $`R(\omega )`$ for the ab-plane, i.e. $`R_{ab}(\omega )`$, and for the c-axis, i.e. $`R_c(\omega )`$. Details of the optical measurements were described in Ref. .

Figure 1(a) shows $`R(\omega )`$ for NSMO multi-domain single crystals at room $`T`$. Note that all the samples are insulators with disordered S/O states at room $`T`$. Clearly, $`R(\omega )`$ are similar for all the samples despite the different hole concentrations. On the other hand, Fig. 1(b) shows $`R(\omega )`$ for the NSMO multi-domain samples at low $`T`$. $`R(\omega )`$ show quite different characteristics for the samples with different doping level x. For the x=0.40 (i.e., $`F`$-type ordered) sample, $`R(\omega )`$ increases significantly below 1.2 eV and approach 1.0 in the dc limit, showing metallic behavior. For the x=0.50 (i.e., $`CE`$-type ordered) sample, $`R(\omega )`$ shows an insulator-like behavior with rather strong suppression below 0.2 eV. For the x=0.55 (i.e., $`A`$-type ordered) sample, $`R(\omega )`$ keeps increasing and approaches to about 0.6 as $`\omega `$ decreases. Finally, for the x=0.65 (i.e., $`C`$-type ordered) sample, $`R(\omega )`$ shows an insulating behavior. Such large differences in $`R(\omega )`$ at low $`T`$ are not easily understood by keeping in consideration of the similarity in $`R(\omega )`$ at room $`T`$ and the small differences in x. It is quite apparent that S/O ordering patterns could be playing an important role in optical properties of NSMO.

We have obtained $`\sigma (\omega )`$ using the Kramers-Kronig (KK) analysis of $`R(\omega )`$ at room $`T`$. Since the optical anisotropy in these disordered S/O states is expected to be small, the KK analysis can be applied at room $`T`$ even for the multi-domain crystals. As shown in Fig. 2(a), $`\sigma (\omega )`$ are quite similar for all the samples. In order to compare the polaron and the orbital excitation scenarios, we should obtain terms related to the coherent and the incoherent motions of free carriers from $`\sigma (\omega )`$.

In the perovskite manganites, there are numerous optical transitions. As shown in Fig. 2(b) for the x=0.55 sample, the room temperature $`\sigma (\omega )`$ below 4.0 eV was fitted with four peaks represented by the dotted lines. The first peak (Peak I) was due to either polaron hopping motions or orbital excitations. The second peak (Peak II) located around 1.5 eV was suggested to be due to the transition between JT split eg bands. \[A Monte Carlo simulation using a Hamiltonian with the double exchange and JT interactions also clearly predicted the existence of Peaks I and II.\] The third peak (Peak III) was suggested to be due to a transition between the Hund’s rule exchange ($`J`$ 2 eV) split bands. Finally, the fourth peak (Peak IV) was proposed to be due to the charge transfer transition between the O 2p and Mn 3d levels. Note that the position of each peak are quite consistent with the results of other measurements, and the validity of those peaks has been suggested by many groups and such analyses have already been performed quite successfully to explain the optical responses of (Bi,Ca)MnO3, Nd0.7Sr0.3MnO3, and (La,Ca)MnO3.

In the polaron scenario, Peak I should be due to the JT polaron hopping motion. Since it can be considered as the interatomic transition from the JT distorted Mn3+\[x ions per unit volume\] to Mn4+\[(1$``$x) ions per unit volume\], it should be proportional to x(1$``$x). The inset in Fig. 2(a) shows the x-dependent strength of Peak I ($`S_I`$). The dotted line in the inset shows, within our fitting error bars, that the x-dependence of $`S_I`$ is well fitted with x(1$``$x). A similar behavior was also observed in (La,Ca)MnO3, supporting the agreement with the polaron scenario. However, as far as we know, the x-dependence of Peak I in the orbital scenario, i.e., the incoherent motion of free carriers due to orbital excitations, has not been predicted yet.

To obtain further insights, we have investigated $`\sigma (\omega )`$ for the x=0.55 single-domain crystal at low $`T`$ in the $`A`$-type S/O ordered state. The open triangles in Figs. 3(a) and 3(b) show $`\sigma (\omega )`$ along the c-axis, $`\sigma _c(\omega )`$, and in the ab-plane, $`\sigma _{ab}(\omega )`$, respectively. Clearly, there exists a rather strong anisotropy of $`\sigma (\omega )`$ below 1.5 eV. The solid lines show the fitting curves based on the above-mentioned mode assignments. As shown in the dotted lines, the strengths of Peaks II and IV were assumed to be $`T`$-independent, and the strength of Peak III were assumed to be nearly zero. Since Peaks II and IV are due to the interband transitions between JT split eg bands and the charge transfer transition between the O 2p and Mn 3d bands, respectively, they should be nearly $`T`$-independent. And, since Peak III is due to a transition between the Hund’s rule exchange ($`J`$ 2 eV) split bands, it should vanish in the spin ordered states. On the other hand, the strength of Peak I is suppressed along the c-axis but increases in the ab-plane. The enhancement factor of Peak I in the $`A`$-type S/O ordered state along the c-axis, $`\eta _A^c`$, was estimated to be about 0.4, and that in the ab-plane, $`\eta _A^{ab}`$, about 1.8. The difference between $`\eta _A^c`$ and $`\eta _A^{ab}`$ for Peak I can explain the reported optical anisotropy in the S/O ordered state. It is important to note that all the parameters, except the enhancement factors, are fixed in the fitting process: namely, the solid lines in Figs. 3(a) and 3(b) are the results of single parameter fittings, keeping all other parameters the same as in the room temperature $`\sigma (\omega )`$ fitting.

How can we explain the difference between $`\eta _A^c`$ and $`\eta _A^{ab}`$? Recently, $`\sigma (\omega )`$ for the planar orbital ordered phase of the eg electron with x$`{}_{}{}^{2}`$y2 symmetry (the $`A`$-type orbital ordered state) were calculated based on the orbital excitation scenario. This 2-dimensional model predicted that the spectral weight of the incoherent absorption, $`S_I`$, should decrease with decreasing $`T`$. On the other hand, that of the coherent absorption, corresponding to the Drude peak, should increase. Contrary to the predictions, the experimental value of $`\eta _A^{ab}`$ becomes larger than 1.0. Moreover, Fig. 3(b) does not show any significant Drude-like absorption term. \[Note that the value of $`\sigma _{ab}(\omega )`$ at the dc limit in Fig. 3(b) was estimated at around 2$`\times `$102 $`\mathrm{\Omega }^1cm^1`$, which is consistent with the transport experiment value on a single-domain crystal \]. Therefore, the model based on the orbital excitation scenario cannot explain our experimental data.

As far as we know, there are no explicit predictions on how the S/O ordering pattern will affect the polaron absorption. However, as shown in Fig. 3(c), it is quite clear that the polaron hopping motion along the c-axis should be suppressed since the antiparallel spin alignment in the $`A`$-type S/O ordered state will cause a large energy cost of $`J`$. Due to the x$`{}_{}{}^{2}`$y2 eg orbital ordering within the ab-plane, the overlap of two Mn eg orbitals via an oxygen pz orbital along the c-axis becomes small, thereby suppressing $`S_I`$ even further. On the contrary, $`S_I`$ in the ab-plane will become enhanced due to the ferromagnetic spin alignment and the increase of orbital overlap in the x$`{}_{}{}^{2}`$y2 orbital ordered state.

In the spin/orbital pattern dependent polaron (SOPDP) model, the changes in $`S_I`$ in the S/O ordered state should be calculated by considering the orbital contributions in the conventional polaron models, which usually include the double exchange and the JT polaron terms. However, this calculation might be quite difficult. In this paper, we made simple and crude estimations of the enhancement factors using the matrix element terms in the dipole transition. First, we considered the JT polaron absorption as a transition from a lower eg level of a JT distorted Mn3+ site to one of doubly degenerate eg levels of Mn4+ sites. Second, for the $`A`$-type ordered state, we evaluated the sum of possible dipole matrix elements considering the neighboring S/O configurations. For the Mn 3d eg and the O 2p orbital wave functions, we used the wave functions of isolated ions located at corresponding sites. Third, for the disordered state, we evaluated the sum of possible dipole matrix elements, assuming that there is no preference of S/O for the initial state. Finally, we divided the summed values to obtain enhancement factors of $`S_I`$ in the ordered state. Details of this estimation procedure will be published elsewhere. We found that $`\eta _A^c`$ 0 and $`\eta _A^{ab}`$ 1.7. The estimated value of $`\eta _A^{ab}`$ agrees quite well with the experimental value, suggesting the validity of the SOPDP model. The suppression of $`S_I`$ along the c-axis can be explained at least qualitatively. However, the difference between the experimental and the estimated values of $`\eta _A^c`$ might be due to the experimental errors and/or the over-simplification of our estimation.

For the other multi-domain single crystal samples, the polarized optical microscope images showed that they consist of multi-domains whose typical size is a little larger than 10 $`\mu m`$. Thus, the KK analysis based on $`R(\omega )`$ for multi-domain crystals may lead to errors in the positions and strengths of Peaks I and II . To avoid such errors, we decided to analyze our $`R(\omega )`$ directly. The optical responses of inhomogeneous media have been investigated extensively. In the long wavelength limit, where the wavelength $`\lambda `$ of the incoming light is larger than the domain size $`L`$, the medium should look homogeneous with an effective conductivity. On the other hand, if $`\lambda <L`$, the light will see the responses of individual domains, so we will obtain averaged optical responses. If the c-axes of the domains are randomly oriented along any one of the crystallographic axes, the measured $`R(\omega )`$ is given by \[$`\frac{2}{3}R_{ab}(\omega )`$+$`\frac{1}{3}R_c(\omega )`$\].

To test the above argument, we have calculated $`R(\omega )`$ for the multi-domain case of $`A`$-type (i.e., x=0.55) ordering, whose polarization-dependent $`\sigma (\omega )`$ were obtained experimentally. As mentioned above, the predictions of the SOPDP model are quite well matching to the real $`\sigma (\omega )`$ of a single-domain crystal. From $`\sigma _{ab}(\omega )`$ and $`\sigma _c(\omega )`$ of the SOPDP model, as shown in Figs. 3(a) and 3(b) (solid lines), respectively, corresponding reflectance spectra for given polarization directions were obtained. Using the predicted $`R_{ab}(\omega )`$ and $`R_c(\omega )`$, $`R(\omega )`$ were evaluated and is shown as the open triangles in Fig. 1(b). The excellent agreement between the experimental data and the SOPDP model predictions suggests that the method of obtaining the $`R(\omega )`$ for a multi-domain sample using the SOPDP model is accurate for $`\lambda <L10`$ $`\mu m`$.

For the $`F`$-type (i.e., x=0.40) sample, which has isotropic S/O structure even at low $`T`$, the enhancement of the mid-IR region and the metallic behavior can be explained well with the existing JT polaron model, which considers the JT distortion and double exchange interaction, as observed in La0.7Ca0.3MnO3. For other multi-domain samples, similar analyses based on the SOPDP model were also performed. For the $`C`$-type ordered (i.e., x=0.65) sample, $`\eta _C^c`$ $``$ 2.9 and $`\eta _C^{ab}`$ $``$ 0 can provide a good description to our experimental data. \[Here we used the same method with $`\eta _A`$ using 3z2-r2 orbital for the initial state instead of x2-y2 at low T.\] For the $`CE`$-type (i.e., x=0.50) sample , this approach can provide a reasonable fit to the experimental $`R(\omega )`$ when Peak I becomes asymmetric. It might be due to the opening of the charge gap at the low frequency region in the S/O ordered state. As shown in Fig. 1(b), most experimental $`R(\omega )`$ can be explained by the theoretical predictions using the SOPDP model with a single fitting parameter ($`\eta `$). The good agreements suggest that the optical response observed in various spin/orbital ordered states of NSMO can be explained quite well within the SOPDP model.

In conclusion, we have presented the temperature-dependent spectral weight changes of Nd1-xSrxMnO3 single crystals, which strongly depend on the pattern of the spin/orbital ordering at low temperature. These optical responses can be described quantitatively by taking account of spin/orbital pattern-dependent polaron hopping.

This work was financially supported by the Ministry of Science and Technology through the Creative Research Initiative Program. The work by Y.M. was supported by a Grant-In-Aid for Scientific Research from the Ministry of Education, Science, Sports, and Culture, and from PRESTO, JST. KHK is also supported by the BK-21 Project of the Ministry of Education.


format: text

text:
# 1 One loop diagrams contributing to the matrix elements of the twist-two operators. The self-energy insertions into the external line, which are not displayed explicitly, have to be added.

1. The fundamental problem of hadronic physics is the understanding of the structure of the nucleon, — the building block of our planet and its inhabitants, — in terms of quark and gluonic degrees of freedom. The nucleon’s mass and momentum content is one of the issues being under study for decades. Higgs mechanism giving quarks their masses generates only a tiny fraction of the total nucleon mass, the bulk of which comes entirely from the complicated dynamics of colored glue. This idea was supported by deeply inelastic scattering experiments which revealed that quarks and gluons are equally responsible for the energy-momentum structure of the nucleon with both of them carrying about $`50\%`$ of its momentum. A similar study has to be addressed for the spin content of hadrons. The nucleon being a composite particle builds its spin from the angular momenta111For a gauge particle the decomposition of its angular momentum into its spin and orbital components is not possible. The hand-waiving argument goes as follows: In quantum mechanics, the wave function of a spin-$`s`$ particle is a symmetric rank-$`2s`$ spinor having $`(2s+1)`$ components which transform into each other under the rotation of the coordinate system. The orbital wave function is related to the coordinate dependence of wave functions and is given by the spherical harmonics of order $`l`$ for the angular momentum $`l`$ of the system. Therefore, in order to distinguish clearly between the spin and orbital momentum, the spin and coordinate properties must be independent. As it is obvious, this condition is not fulfilled for gauge particles whose description in terms of field operators inevitably involves a gauge condition. For instance, in the Coulomb gauge the gluon wave function is given by the three-potential $`𝑨^a(x)`$ equivalent to the second rank spinor which, however, is a subject for the gauge condition $`\mathrm{div}𝑨^a(x)=0`$. As a result the coordinate dependence of the vector cannot be independently defined for each of its components and leads to the inability to separate the spin and orbital degrees of freedom. of its constituents $`J_{q,g}`$, i.e., in the case of quarks from their spin $`\mathrm{\Delta }\mathrm{\Sigma }/2`$ and orbital motion $`L_q`$,

$$\frac{1}{2}=\frac{\mathrm{\Delta }\mathrm{\Sigma }}{2}+L_q+J_g.$$
(1)
A next natural question arises about the numerical apportion between these contributions. The naive quark model (QM) attributes the whole of the nucleon spin to the sum of the ones of quarks $`\mathrm{\Delta }\mathrm{\Sigma }_{\mathrm{QM}}=1`$ . However, the experimental data tells the opposite: only a small fraction of the total spin is carried by quarks $`\mathrm{\Delta }\mathrm{\Sigma }_{\mathrm{exp}}=0.16\pm 0.08`$ . Thus in order to test the spin sum rule (1) one has to have access to the other sources of the nucleon spin $`L_q`$ and $`J_g`$. The nucleon’s angular momentum is expressed in QCD by a hadronic matrix element of the coordinate-space moment of the quark and gluon momentum densities given by specific components of the energy-momentum tensor. At this point recall that the magnetic moment of a particle is the space moment of the electric current flow and can be measured as a matrix element of the electric current at non-zero momentum transfer, i.e., in elastic electron-nucleon scattering by accessing electromagnetic form factors. The same methodology is applicable in the current circumstances so that the missing orbital momentum parts can be deduced from off-forward matrix elements of the QCD energy-momentum tensor. Thus, in order to understand the nucleon spin content one has to study nucleon gravitational form factors. These functions, which arise in the nucleon scattering off a weak gravitational field, are measurable experimentally in exclusive processes such as deeply virtual Compton scattering  and hard exclusive meson production . Moreover, they can be simulated on a lattice. As we just said, apart from the interest in their own right, gravitational form factors at zero recoil determine the quark and gluon contributions to the angular momentum. In practical determination of the latter either in laboratory or lattice experiments one needs in all cases to perform an extrapolation to the limit of zero momentum transfer $`\mathrm{\Delta }^2=0`$ to get the net effect. The experimental kinematics is restricted to low $`t`$-channel momenta of order $`0.1\mathrm{GeV}^2`$  and, thus, favors the applicability of the chiral perturbation theory ($`\chi `$PT)  to unravel the small-momentum transfer behavior of the form factors in question. On a lattice, the minimal achieved momentum is limited by the size of the lattice. The first Monte Carlo simulations  were done with lattices $`16^2\times 32`$ where $`\mathrm{\Delta }_{\mathrm{min}}^20.4\mathrm{GeV}^2`$. Thus in order to enter the regime of validity of $`\chi `$PT one needs repeat the simulations with bigger lattices. In the present Letter we compute one-loop contributions to off-forward matrix elements of local quark and gluon composite operators and thus determine their momentum transfer dependence.

2. The spin and momentum structure of the nucleon is carried by off-forward matrix elements of quark and gluon composite operators. In our consequent presentation we will be concerned with the parity-even sector which can shed some light on the spin structure of the nucleon. Under chiral $`SU(2)_LSU(2)_R`$ group the operators transform either as isovector $`(3,1)(1,3)`$ or isoscalar $`(1,1)`$ representations and correspond to flavor non-singlet $`^{q,a}`$ and singlet $`^{q,0}`$ quark operators, respectively,

$$_{\mu _1\mu _2\mathrm{}\mu _j}^{q,A}=\underset{\mu _1\mu _2\mathrm{}\mu _j}{\text{S}}\overline{\psi }\tau ^A\gamma _{\mu _1}i\underset{\mu _1}{\overset{}{𝒟}}i\underset{\mu _2}{\overset{}{𝒟}}\mathrm{}i\underset{\mu _j}{\overset{}{𝒟}}\psi ,$$
(2)
where $`\tau ^A=(\text{1}\text{l},\tau ^a)`$ with $`\text{1}\text{l}1_{[2]\times [2]}`$. The operator of symmetrization and trace subtraction is defined as follows, e.g., $`\text{S}_{\mu _1\mu _2}t_{\mu _1\mu _2}=\frac{1}{2!}(t_{\mu _1\mu _2}+t_{\mu _2\mu _1}\frac{1}{2}g_{\mu _1\mu _2}t_{\mu \mu })`$. The isoscalar combination $`^{q,0}`$ mixes under the renormalization group evolution with gluonic operators,

$$_{\mu _1\mu _2\mathrm{}\mu _j}^g=\underset{\mu _1\mu _2\mathrm{}\mu _j}{\text{S}}G_{\nu \mu _1}i\underset{\mu _1}{\overset{}{𝒟}}i\underset{\mu _2}{\overset{}{𝒟}}\mathrm{}i\underset{\mu _{j1}}{\overset{}{𝒟}}G_{\mu _j\nu }.$$
(3)
The nucleon matrix element of these operators are parametrized via ‘form factors’ as follows

$`p_2|_{\mu _1\mathrm{}\mu _j}^A|p_1`$ $`=`$ $`\underset{\mu _1\mathrm{}\mu _j}{\text{S}}\overline{u}(p_2)\tau ^A\gamma _{\mu _1}u(p_1)\left\{p_{\mu _2}\mathrm{}p_{\mu _j}A_{j,j}(\mathrm{\Delta }^2)+\mathrm{}+\mathrm{\Delta }_{\mu _2}\mathrm{}\mathrm{\Delta }_{\mu _j}A_{j,1}(\mathrm{\Delta }^2)\right\}`$ (4)
$`+`$ $`\underset{\mu _1\mathrm{}\mu _j}{\text{S}}\overline{u}(p_2)\tau ^A{\displaystyle \frac{i\sigma _{\mu _1\nu }\mathrm{\Delta }_\nu }{2M}}u(p_1)\left\{p_{\mu _2}\mathrm{}p_{\mu _j}B_{j,j}(\mathrm{\Delta }^2)+\mathrm{}+\mathrm{\Delta }_{\mu _2}\mathrm{}\mathrm{\Delta }_{\mu _j}B_{j,1}(\mathrm{\Delta }^2)\right\}`$
$`+`$ $`\underset{\mu _1\mathrm{}\mu _j}{\text{S}}{\displaystyle \frac{\overline{u}(p_2)\tau ^Au(p_1)}{2M}}\mathrm{\Delta }_{\mu _1}\mathrm{}\mathrm{\Delta }_{\mu _j}C_j(\mathrm{\Delta }^2),`$
where we have introduced the vectors $`pp_1+p_2`$ and $`\mathrm{\Delta }=p_2p_1`$ and assumed the normalization for hadronic states $`p|p=2p_0V`$ and bispinors $`\overline{u}u=2M`$. We have reduced in the above equation all tensor structures of the scalar Dirac bilinear containing at least one vector $`p_\mu `$ by means of the Gordon identity to the vector and tensor Dirac bilinears of the first two lines. There is an important constraints on the matrix elements (4) imposed by the time-reversal invariance, namely, vanishing of Lorentz structures with odd powers of the momentum transfer $`\mathrm{\Delta }`$ . We will keep this condition in mind without implementing it explicitly in our subsequent considerations.

The reduced matrix elements $`A`$, $`B`$ and $`C`$ are related to the moments of parity-even quark or gluon generalized parton distributions  measurable in deeply virtual Compton scattering via

$`{\displaystyle _1^1}𝑑xx^{j1}H(x,\eta ,\mathrm{\Delta }^2)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{j1}{}}}\eta ^kA_{j,jk}(\mathrm{\Delta }^2)+\eta ^jC_j(\mathrm{\Delta }^2),`$
$`{\displaystyle _1^1}𝑑xx^{j1}E(x,\eta ,\mathrm{\Delta }^2)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{j1}{}}}\eta ^kB_{j,jk}(\mathrm{\Delta }^2)\eta ^jC_j(\mathrm{\Delta }^2).`$ (5)
The angular momentum sum rule (1) involves the form factors at zero recoil entering their second moment

$$J=\frac{1}{2}\underset{\mathrm{\Delta }^20}{lim}\left\{A_{2,2}(\mathrm{\Delta }^2)+B_{2,2}(\mathrm{\Delta }^2)\right\},$$
(6)
which are attributed to quarks or gluons depending on the contributing composite operators.

We consider single-nucleon systems and, in order to have a consistent power counting, use the formalism of the heavy baryon $`\chi `$PT which treats the nucleon as a non-relativistic infinitely heavy particle . For the purpose of a unique translation of the Lorentz-covariant matrix elements in Eq. (4) into non-relativistic ones its is convenient to use the Breit reference frame . It is defined by the condition $`𝒑=0`$ which leads in turn to $`\mathrm{\Delta }_0=0`$, $`p_0=2M\sqrt{1\mathrm{\Delta }^2/(4M^2)}2M\delta `$, and $`𝒑_2=𝒑_1=𝚫/2`$. In the Breit frame, the Dirac bilinears arising in Eq. (4) are reduced via equations

$`\overline{u}(p_2)\gamma _\mu u(p_1)=v_\mu \overline{u}_v(p_2)u_v(p_1)+{\displaystyle \frac{1}{M}}\overline{u}_v(p_2)[S_\mu ,S\mathrm{\Delta }]u_v(p_1),`$
$`\overline{u}(p_2){\displaystyle \frac{i\sigma _{\mu \nu }\mathrm{\Delta }_\nu }{2M}}u(p_1)=v_\mu {\displaystyle \frac{\mathrm{\Delta }^2}{4M^2}}\overline{u}_v(p_2)u_v(p_1)+{\displaystyle \frac{1}{M}}\overline{u}_v(p_2)[S_\mu ,S\mathrm{\Delta }]u_v(p_1),`$
to the ones constructed from the projected large components of the heavy baryon bispinor $`u_v(p_i)(1+\overline{)}v)u(p_i)/\sqrt{2+2(vp_i)/M}`$ normalized as $`\overline{u}_vu_v=2M`$. The velocity $`v`$ arises in the heavy-mass decomposition of the incoming and outgoing nucleon momenta $`p_{1,2}=Mv+k\mathrm{\Delta }/2`$ and the residual momentum is defined as $`k=(M(\delta 1),\mathrm{𝟎})`$ when $`v=(1,\mathrm{𝟎})`$. The spin vector $`S_\mu \frac{i}{2}\sigma _{\mu \nu }\gamma _5v_\nu `$ reduces to the spin three-vector $`S=(0,\mathrm{\Sigma }/2)`$ in the nucleon rest frame. Exploiting these results, the decomposition (4) transforms in the Breit frame to the equation

$`p_2|_{\mu _1\mathrm{}\mu _j}^A|p_1=\underset{\mu _1\mathrm{}\mu _j}{\text{S}}\overline{u}_v(p_2)\tau ^Au_v(p_1)v_{\mu _1}\left\{(2M\delta )^{j1}E_{j,j}(\mathrm{\Delta }^2)v_{\mu _2}\mathrm{}v_{\mu _j}+\mathrm{}+E_{j,1}(\mathrm{\Delta }^2)\mathrm{\Delta }_{\mu _2}\mathrm{}\mathrm{\Delta }_{\mu _j}\right\}`$
$`+\underset{\mu _1\mathrm{}\mu _j}{\text{S}}{\displaystyle \frac{1}{M}}\overline{u}_v(p_2)\tau ^A[S_{\mu _1},S\mathrm{\Delta }]u_v(p_1)\left\{(2M\delta )^{j1}M_{j,j}(\mathrm{\Delta }^2)v_{\mu _2}\mathrm{}v_{\mu _j}+\mathrm{}+M_{j,1}(\mathrm{\Delta }^2)\mathrm{\Delta }_{\mu _2}\mathrm{}\mathrm{\Delta }_{\mu _j}\right\}`$
$`+\underset{\mu _1\mathrm{}\mu _j}{\text{S}}{\displaystyle \frac{\overline{u}_v(p_2)\tau ^Au_v(p_1)}{2M\delta ^1}}\mathrm{\Delta }_{\mu _1}\mathrm{}\mathrm{\Delta }_{\mu _j}C_j,`$ (7)
where we have introduced a short-hand notation for the generalized electric- and magnetic-like form factors

$$E_{j,k}A_{j,k}+\frac{\mathrm{\Delta }^2}{4M^2}B_{j,k},M_{j,k}A_{j,k}+B_{j,k}.$$
(8)
Now we are in a position to address the calculation of the small-$`\mathrm{\Delta }`$ dependence of these form factors. As it is obvious from the tensor decomposition, one cannot compute form factors accompanying Lorentz tensors involving several $`\mathrm{\Delta }_\mu `$, i.e., $`k<j2`$, while restricting oneself to next-to-leading order in the chiral expansion. Also note that $`E_{1,1}`$ and $`M_{1,1}`$ correspond to the standard Sachs nucleon electromagnetic form factors extensively discussed in the literature  and, therefore, their evaluation will not be repeated presently.

3. The calculational procedure is rather straightforward. First, one constructs composite operators from the nucleon and pion fields which match the quantum numbers of the ones on the quark-gluon level (2,3) and adds them to the effective Lagrangian . Then one uses the heavy-baryon chiral perturbation theory , see also the review , for the computation of tree and one-loop contributions. For the case at hand, one needs the leading order pion-nucleon Lagrangian

$$=\overline{N}_v\left\{iv𝒟+2g_AS𝒜\right\}N_v+\frac{f_\pi ^2}{4}\mathrm{tr}\left\{_\mu \mathrm{\Sigma }_\mu \mathrm{\Sigma }^{}\right\}+\lambda \mathrm{tr}\left\{M_q\left(\mathrm{\Sigma }+\mathrm{\Sigma }^{}\right)\right\}.$$
(9)
The axial-vector and vector (in the covariant derivative $`𝒟_\mu _\mu +𝒱_\mu `$) pion potentials are

$$𝒱_\mu \frac{1}{2}\left(\xi ^{}_\mu \xi +\xi _\mu \xi ^{}\right),𝒜_\mu \frac{i}{2}\left(\xi ^{}_\mu \xi \xi _\mu \xi ^{}\right).$$
(10)
Here $`\mathrm{\Sigma }`$ is constructed from the multiplet of the pion fields

$$\mathrm{\Sigma }\xi ^2=\mathrm{exp}\left(\frac{i}{f_\pi }\stackrel{}{\pi }\stackrel{}{\tau }\right),$$
(11)
and $`f_\pi =93\mathrm{MeV}`$ the pion decay constant. Expanding to the first non-trivial order in the pion fields we get for the axial and vector potentials: $`𝒜_\mu =\frac{1}{2f_\pi }_\mu \stackrel{}{\pi }\stackrel{}{\tau }+\mathrm{}`$ and $`𝒱_\mu =\frac{i}{(2f_\pi )^2}\epsilon _{abc}\tau ^c\pi ^a_\mu \pi ^b+\mathrm{}`$, respectively. In the mass term, the coefficient $`\lambda `$ is related to the quark condensate via $`\lambda =\frac{1}{2}\overline{\psi }\psi `$ and the quark mass matrix is $`M_q=\mathrm{diag}(m_u,m_d)`$.

Ultraviolet divergences generated by one-loop diagrams have to be absorbed into coefficients of next-to-leading order counterterms in the chiral expansion.

Now we are in a position to discuss the construction of a basis of hadronic twist-two operators corresponding to the quark and gluon ones in Eqs. (2) and (3). There are two types of local hadronic operators: built either from the pions alone or bilinear in the nucleon fields. The nucleon operators have the form which mimics the tensor decomposition of the off-forward matrix elements. The leading operators contributing to isoscalar and isovector combinations of $`E_{j,j}`$ and $`M_{j,j}`$ are

$`𝒪_{\mu _1\mu _2\mathrm{}\mu _j}^{N,A}`$ $`=`$ $`a_j^N(2M)^{j1}\underset{\mu _1\mathrm{}\mu _j}{\text{S}}v_{\mu _1}\mathrm{}v_{\mu _j}\overline{N}_v\tau _{\xi +}^AN_v`$ (12)
$`+`$ $`b_j^N(2M)^{j1}(i_\nu )\underset{\mu _1\mathrm{}\mu _j}{\text{S}}v_{\mu _1}\mathrm{}v_{\mu _{j1}}\overline{N}_v\tau _{\xi +}^A{\displaystyle \frac{[S_{\mu _j},S_\nu ]}{M}}N_v+\mathrm{},`$
where

$$\tau _{\xi +}^A\frac{1}{2}\left(\xi \tau ^A\xi ^{}+\xi ^{}\tau ^A\xi \right).$$
(13)
The matching coefficients $`a^N`$ and $`b^N`$ from the partonic to the hadronic level are unknown and have to be determined from experimental data. Subleading operators are deduced from this expression by a mere replacement of $`v`$s in the Lorentz structure by the derivatives. However, due to the time-reversal symmetry restrictions one has to replace an even number of velocities. Here and everywhere later, the nucleon operators are normalized as follows $`N_v(0)|p=u_v(p)|0`$.

The pion sector has to be considered separately for the isoscalar and isovector components. As we will see momentarily, one can essentially reduce the number of operators at each level in the derivative expansion provided one uses all underlying symmetries. A number of relations arises on the basis of the following ‘magic’ property of the Pauli matrices $`\tau ^2\stackrel{}{\tau }\tau ^2=\stackrel{}{\tau }^T`$ which results into a useful conjugation property for the nonlinear pion field $`\mathrm{\Sigma }`$,

$$\tau ^2\mathrm{\Sigma }^{}\tau ^2=\mathrm{\Sigma }^T.$$
(14)
4. Let us start from the isoscalar sector which transforms as $`(1,1)`$ under chiral $`SU(2)_LSU(2)_R`$. The leading nucleon operators are displayed in Eq. (12) for $`A=0`$ so that $`\tau _{\xi +}^0=\text{1}\text{l}`$. Next, one can immediately conclude that the lowest local operator in Eq. (2) does not correspond to any pion operator as a consequence of the unitarity of $`\mathrm{\Sigma }`$ and Eq. (14),

$$\mathrm{tr}\left\{\mathrm{\Sigma }_\mu \mathrm{\Sigma }^{}\right\}=0.$$
(15)
Its generalization for an arbitrary odd $`j=2k+1`$ reads

$$\underset{\mu _1\mu _2\mathrm{}\mu _j}{\text{S}}\mathrm{tr}\left\{_{\mu _1}\mathrm{}_{\mu _k}\mathrm{\Sigma }_{\mu _{k+1}}\mathrm{}_{\mu _j}\mathrm{\Sigma }^{}\right\}=\frac{1}{2}\underset{\mu _1\mu _2\mathrm{}\mu _j}{\text{S}}_{\mu _1}\mathrm{tr}\left\{_{\mu _2}\mathrm{}_{\mu _{k+1}}\mathrm{\Sigma }_{\mu _{k+2}}\mathrm{}_{\mu _j}\mathrm{\Sigma }^{}\right\}.$$
(16)
Therefore, all pion operators with an odd number of derivatives $`j=2k+1`$ can be reduced to the ones containing an odd number of total derivatives, so that we have

$`\underset{\mu _1\mathrm{}\mu _j}{\text{S}}\left\{{\displaystyle \underset{l=0}{\overset{k1}{}}}a_{j,2l+1}^\pi (i)^{2l+1}\mathrm{tr}\left\{(i)^{kl}\mathrm{\Sigma }(i)^{kl}\mathrm{\Sigma }^{}\right\}\right\}_{\mu _1\mathrm{}\mu _j}.`$
Here, on the right-hand side of the equation, we imply that the Lorentz indices are attached to the partial derivatives $``$ in a particular order, e.g., $`\{^3\}_{\mu _1\mu _2\mu _3}=_{\mu _1}_{\mu _2}_{\mu _3}`$. Note, however, that matrix elements of these operators violate the time reversal properties of the off-forward matrix elements in question (4) and thus cannot contribute.

Using the symmetry properties established above we can easily write composite pion operators for even $`j=2k`$,

$$𝒪_{\mu _1\mu _2\mathrm{}\mu _j}^\pi =f_\pi ^2\underset{\mu _1\mathrm{}\mu _j}{\text{S}}\left\{\underset{l=0}{\overset{k1}{}}a_{j,2l}^\pi (i)^{2l}\mathrm{tr}\left\{(i)^{kl}\mathrm{\Sigma }(i)^{kl}\mathrm{\Sigma }^{}\right\}\right\}_{\mu _1\mathrm{}\mu _j}.$$
(17)
Since quantum numbers of hadronic operators match both quark and gluon ones in Eqs. (2) and (3), the difference between the latter on the hadronic level arises only in the value of the matching coefficients which appear in two species, $`a_{j,k}^{\pi ,q}`$ and $`a_{j,k}^{\pi ,g}`$. Below, the index which attributes hadronic operators either to the quark or gluon sector will not be displayed explicitly for brevity. Note that $`a_{j,0}^\pi a_j^\pi `$ survive in the forward matrix elements and, thus, are measurable in conventional deeply inelastic scattering experiments. For instance, $`a_2^\pi `$ is related to the momentum fraction of the pion carried by quarks and gluons forming it. For an on-shell pion at low normalization point $`\mu ^20.3\mathrm{GeV}^2`$ one has  $`a_2^{\pi ,q}=x_q_\pi 0.7`$ and respectively $`a_2^{\pi ,g}=x_g_\pi 0.3`$ due to the momentum conservation. There is a number of pion operators with more than two $`\mathrm{\Sigma }`$ fields since their insertion does not alter the twist of the composite operator, however, they do not contribute to the one-loop matrix elements we are interested in and will be totally omitted.

To summarize, the isoscalar nucleon gravitational form factors, — the second moment of the generalized parton distributions $`H`$ and $`E`$, — receive contributions from the following pion and nucleon operators:

$`𝒪_{\mu _1\mu _2}^\pi `$ $`=`$ $`f_\pi ^2a_2^\pi \underset{\mu _1\mu _2}{\text{S}}\mathrm{tr}\left\{_{\mu _1}\mathrm{\Sigma }_{\mu _2}\mathrm{\Sigma }^{}\right\},`$ (18)
$`𝒪_{\mu _1\mu _2}^N`$ $`=`$ $`a_2^N(2M)\underset{\mu _1\mu _2}{\text{S}}v_{\mu _1}v_{\mu _2}\overline{N}_vN_v+b_2^N(2M)(i_\nu )\underset{\mu _1\mu _2}{\text{S}}v_{\mu _1}\overline{N}_v{\displaystyle \frac{[S_{\mu _2},S_\nu ]}{M}}N_v`$ (19)
$`+`$ $`c_2^N{\displaystyle \frac{1}{2M}}\underset{\mu _1\mu _2}{\text{S}}(i_{\mu _1})(i_{\mu _2})\overline{N}_vN_v,`$
and we have also included an operator with two total derivatives which generates the structure $`C_2(\mathrm{\Delta }^2)`$.

A calculation of the diagrams (a) and (c), (b) and (d) do not contribute for isoscalar operators, gives

$`E_{2,2}(\mathrm{\Delta }^2)`$ $`=`$ $`a_2^N+3a_2^\pi {\displaystyle \frac{g_A^2}{64\pi f_\pi ^2M}}\left\{\left(2m^2\mathrm{\Delta }^2\right){\displaystyle _0^1}𝑑x\sqrt{m^2(x)}+{\displaystyle \frac{4}{3}}m^3\right\},`$
$`M_{2,2}(\mathrm{\Delta }^2)`$ $`=`$ $`b_2^N+3{\displaystyle \frac{g_A^2}{(4\pi f_\pi )^2}}{\displaystyle _0^1}𝑑x\left(a_2^\pi m^2(x)b_2^Nm^2\delta (x)\right)\left\{\mathrm{ln}{\displaystyle \frac{m^2(x)}{\mathrm{\Lambda }_\chi ^2}}1\right\},`$
$`C_2(\mathrm{\Delta }^2)`$ $`=`$ $`c_2^N+3a_2^\pi {\displaystyle \frac{g_A^2}{16\pi f_\pi ^2}}M\left(2m^2\mathrm{\Delta }^2\right){\displaystyle _0^1}𝑑x{\displaystyle \frac{x(1x)}{\sqrt{m^2(x)}}},`$ (20)
where $`m^2(x)m^2x(1x)\mathrm{\Delta }^2`$. Note the absence of one-loop contributions to the $`E`$ and $`C`$ form factors from the nucleon operators. It is a mere consequence of a cancellation of the one-loop contribution (c) by self-energy insertions into the external lines. We have minimally absorbed ultraviolet divergences generated in $`M_{2,2}`$ into higher order counterterms that results into a redefinition of their coefficients, with which they enter in the effective Lagrangian $`r_ir_i+\mathrm{coeff}L_\epsilon `$, with $`L_\epsilon \frac{1}{\epsilon }\gamma _E+\mathrm{ln}4\pi `$. So that the right-hand side of $`M_{2,2}`$ has extra analytic terms

$`{\displaystyle \frac{m^2}{(4\pi f_\pi )^2}}r_1+{\displaystyle \frac{\mathrm{\Delta }^2}{(4\pi f_\pi )^2}}r_2.`$
Note, that one can also completely absorb analytical contributions present in (S0.Ex11) into counterterms. This simply sets the renormalization scheme prescription which has to be used for all observables once $`r_i`$ are fitted to experimental data in one.

Another comment concerns $`E_{2,2}`$. Pion operators generate $`1/M`$-suppressed contributions to the structure of interest. Therefore, by power counting we have to add $`𝒪(1/M)`$ bilinear nucleon operators constructed from one large and one small components of the nucleon field as well as analogous terms stemming from the chiral Lagrangian. Lorentz invariance fixes unambiguously their coefficients and no new low energy constants arise. However, the effect of these contributions in form factors is $`\mathrm{\Delta }^2`$-independent and analytic in the pion mass. Thus they will not computed by us presently, although can be anticipated to generate an addendum $`5a_2^Ng_A^2m^3/\left(32\pi f_\pi ^2M\right)`$ on the right-hand side of $`E_{2,2}`$ in Eq. (S0.Ex11) required by momentum sum rules discussed below. Also, there arises a Foldy-like term $`(b_j^Na_j^N)\mathrm{\Delta }^2/(4M^2)`$ in $`E_{j,j}`$, similarly to the form factor case , from contributions of the small nucleon field components in the heavy-mass expansion of relativistic nucleon operators.

Sum rules for the total momentum and spin of the nucleon (pion) impose constraints on the coefficients $`a_2^N`$ and $`b_2^N`$ ($`a_2^\pi `$),

$$a_2^{N,q}+a_2^{N,g}=1,b_2^{N,q}+b_2^{N,g}=1,a_2^{\pi ,q}+a_2^{\pi ,g}=1.$$
(21)
The latter two equations imply that the total gravitomagnetic moment of the nucleon vanishes, i.e., $`B_{2,2}^q(0)+B_{2,2}^g(0)=0`$ .

The leading $`(j,j)`$-structures of the higher $`j`$-moments ($`j>2`$) do not receive non-analytic contributions in the momentum transfer at next-to-leading order in the chiral expansion due to the absence of relevant pion operators. Thus one gets

$$E_{j,j}(\mathrm{\Delta }^2)=a_j^N+\mathrm{},M_{j,j}(\mathrm{\Delta }^2)=b_j^N\left(13\frac{g_A^2m^2}{(4\pi f_\pi )^2}\mathrm{ln}\frac{m^2}{\mathrm{\Lambda }_\chi ^2}+\mathrm{}\right),$$
(22)
where ellipsis stand for, at least, $`m^2/(4\pi f_\pi )^2`$-suppressed analytic contributions from one-loop diagrams and counterterms.

5. In the isovector sector, the use of Eq. (14) results into

$$\underset{\mu _1\mu _2}{\text{S}}\mathrm{tr}\left\{\tau ^a_{\mu _1}\mathrm{\Sigma }_{\mu _2}\mathrm{\Sigma }^{}\right\}=0.$$
(23)
Its generalizations for even $`j=2k`$ are straightforward and read

$$\underset{\mu _1\mu _2\mathrm{}\mu _j}{\text{S}}\mathrm{tr}\left\{\tau ^a_{\mu _1}\mathrm{}_{\mu _k}\mathrm{\Sigma }_{\mu _{k+1}}\mathrm{}_{\mu _j}\mathrm{\Sigma }^{}\right\}=0.$$
(24)
Using this property it is easy to convince oneself that pion operators for even $`j=2k`$ are reduced to the ones involving an odd number of total derivatives, translated in the momentum space to odd powers of the momentum transfer,

$`\underset{\mu _1\mathrm{}\mu _j}{\text{S}}\left\{{\displaystyle \underset{l=0}{\overset{k1}{}}}a_{2l}^{\pi ;(j)}(i)^{2l+1}\mathrm{tr}\tau ^a\left\{(i)^{kl}\mathrm{\Sigma }(i)^{kl+1}\mathrm{\Sigma }^{}+(i)^{kl}\mathrm{\Sigma }^{}(i)^{kl+1}\mathrm{\Sigma }\right\}\right\}_{\mu _1\mathrm{}\mu _j},`$
which, therefore, cannot contribute to the matrix elements in question due to the time-reversal condition alluded to above.

For odd $`j=2k+1`$, an elementary consideration along the same line as above leads to the following structure of the non-singlet pion operators

$`𝒪_{\mu _1\mu _2\mathrm{}\mu _j}^{\pi ,a}`$ $`=`$ $`{\displaystyle \frac{f_\pi ^2}{2}}\underset{\mu _1\mathrm{}\mu _j}{\text{S}}`$
$`\times `$ $`\left\{{\displaystyle \underset{l=0}{\overset{k}{}}}a_{2l+1}^{\pi ;(j)}(i)^{2l}\mathrm{tr}\tau ^a\left\{(i)^{kl}\mathrm{\Sigma }(i)^{kl+1}\mathrm{\Sigma }^{}+(i)^{kl}\mathrm{\Sigma }^{}(i)^{kl+1}\mathrm{\Sigma }\right\}\right\}_{\mu _1\mathrm{}\mu _j}.`$
The leading nucleon operators are given in Eq. (12) with $`A=a`$.

Thus, for the isovector gravitational form factors there no contributions from coupling to pions while the nucleon operators read

$`𝒪_{\mu _1\mu _2}^{N,a}`$ $`=`$ $`a_2^N(2M)\underset{\mu _1\mu _2}{\text{S}}v_{\mu _1}v_{\mu _2}\overline{N}\tau _{\xi +}^aN+b_2^N(2M)(i_\nu )\underset{\mu _1\mu _2}{\text{S}}v_{\mu _1}\overline{N}\tau _{\xi +}^a{\displaystyle \frac{[S_{\mu _2},S_\nu ]}{M}}N`$ (26)
$`+`$ $`c_2^N{\displaystyle \frac{1}{2M}}\underset{\mu _1\mu _2}{\text{S}}(i_{\mu _1})(i_{\mu _2})\overline{N}\tau _{\xi +}^aN,`$
with $`\tau _{\xi +}^a\frac{1}{2}\left(\xi \tau ^a\xi ^{}+\xi ^{}\tau ^a\xi \right)`$.

Due to the absence of the pion cloud contribution at this order of $`\chi `$PT, i.e., diagrams (a) and (b), no non-analytic dependence on the momentum transfer arises. However, the diagrams (c) and (d) develop chiral logarithms in the pion mass of the form

$`E_{2,2}(\mathrm{\Delta }^2)`$ $`=`$ $`a_2^N\left\{1{\displaystyle \frac{m^2}{(4\pi f_\pi )^2}}\left((3g_A^2+1)\mathrm{ln}{\displaystyle \frac{m^2}{\mathrm{\Lambda }_\chi ^2}}g_A^21\right)\right\},`$
$`M_{2,2}(\mathrm{\Delta }^2)`$ $`=`$ $`b_2^N\left\{1{\displaystyle \frac{m^2}{(4\pi f_\pi )^2}}\left((2g_A^2+1)\mathrm{ln}{\displaystyle \frac{m^2}{\mathrm{\Lambda }_\chi ^2}}1\right)\right\},`$
$`C_2(\mathrm{\Delta }^2)`$ $`=`$ $`c_2^N\left\{1{\displaystyle \frac{m^2}{(4\pi f_\pi )^2}}\left((3g_A^2+1)\mathrm{ln}{\displaystyle \frac{m^2}{\mathrm{\Lambda }_\chi ^2}}g_A^21\right)\right\}.`$ (27)
We imply that one adds counterterms to the right-hand side of these equations linear in $`m^2`$ whose (unknown) coefficients absorb minimally the ultraviolet divergences stemming from loops, $`\mathrm{coeff}L_\epsilon `$, and the Foldy-like term as discussed after Eq. (S0.Ex11). The leading structures of the higher moments $`M_{j,j}`$ and $`E_{j,j}`$, apart from the change of an overall normalization $`b_2^Nb_j^N`$ and $`a_2^Na_j^N`$, have the same dependence on the chiral logarithms as in Eq. (27).

6. In this Letter we used the framework of the heavy baryon chiral perturbation theory in one-loop approximation in order to predict the small momentum transfer dependence of the moments of generalized parton distributions. At next-to-leading order in the chiral expansion only the first two moments develop a non-analytic $`\mathrm{\Delta }^2`$ behavior: isovector electromagnetic form factors and isosinglet gravitational form factors. Higher moments can have at most a linear $`\mathrm{\Delta }^2`$-dependence at this order of $`\chi `$PT.

Although, the momentum and spin sum rules suggest that the total quark and gluon anomalous gravitomagnetic moment is zero , the gravitomagnetic radius of the nucleon does not vanish and reads

$$r^26\frac{dM_{2,2}(\mathrm{\Delta }^2)}{d\mathrm{\Delta }^2}|_{\mathrm{\Delta }^2=0}=3a_2^\pi \frac{g_A^2}{(4\pi f_\pi )^2}\left\{\mathrm{ln}\frac{m^2}{\mathrm{\Lambda }_\chi ^2}+1\right\}+6\frac{r_2}{(4\pi f_\pi )^2}.$$
(28)
A pragmatic application of our present results consists of their use for extrapolation of experimental measurements of the angular momentum sum rule (6) to zero momentum transfer  and as a fit formula for lattice simulations of the parton angular momentum by means of the nucleon gravitational form factors . Our considerations can be extended to other single-nucleon observables such as parity-odd generalized parton distributions as well as off-forward quark transversity and tensor gluon form factors.

We would like to thank J.-W. Chen and T.D. Cohen for useful discussions. This work was supported by the US Department of Energy under contract DE-FG02-93ER40762.


text:
# Pulsar bow-shock nebulae

## 1 Introduction

Pulsars moving supersonically with respect to the ambient medium are expected to give rise to a bow shock. In the case of interaction with the interstellar medium (ISM) the H$`\alpha `$ emission from the nebula may be detected. Although pulsars have high typical velocities, and there are more than one thousand known radio pulsars, only four bow shocks have been discovered so far: PSR 1957+20 (Kulkarni & Hester kulkarni88 (1988)), PSR 2224+65 (Cordes et al. cordes93 (1993)), PSR J0437+4715 (Bell bell95 (1995)), and PSR 0740–28 (Jones et al. jones01 (2001)); while it is not clear whether the cometary-like nebula associated with the isolated neutron star RX J185635–3754 (van Kerkwijk & Kulkarni vankerkwijk01 (2001)) is a bow shock or results from photo-ionization. In this kind of object the relativistic particles in the pulsar wind produce a synchrotron emission too weak to be detectable in radio wavelengths (Gaensler et al. gaensler00 (2000)). These nebulae may be detected instead in optical Balmer lines (mostly H$`\alpha `$ ) as a signature of a non-radiative shock moving through a partially neutral medium (Chevalier & Raymond chevalier78 (1980)). The Balmer lines are due to the de-excitations of neutral H atoms, following collisional excitations or exciting charge-exchange processes. When recombination times are long these are the dominant processes giving emission in the optical band, as it has been widely studied on various supernova remnants (Smith et al. smith91 (1991)).

A puzzling point is that two out of the four known nebulae have a shape closely matching that of a classical bow shock, while the others show a more peculiar shape with a conical tail: the nebula associated with PSR 2224+65 in fact shows a remarkably conical tip followed by a bubble (which justifies why it has been named the “Guitar Nebula”); (Cordes cordes96 (1996)). Various hypothesis have been put forward to justify its shape: a peculiar ISM density distribution in the surroundings of the PSR 2224+65 (Cordes et al. cordes93 (1993)); effects like mass loading due to neutral atoms which might penetrate in the external layers of shocked material (Bucciantini & Bandiera bucciantini01 (2001), hereafter Paper I) . The fourth nebula, which has recently been discovered (Jones et al. jones01 (2001)) near PSR 0740–28, seems to represent an intermediate case, with a “standard” head and a conical tail.

A pulsar wind interacting with a homogeneous ISM is in principle a much easier problem to model than that with a surrounding supernova remnant, since in the latter case the ambient medium might be quite inhomogeneous. Thus our analysis will be restricted to the former case. The ISM is seen by the pulsar as a plane-parallel flow, likely with a constant density (at least on the typical length scale of the bow shock), and lasting long enough to produce a steady-state regime. However, modeling is made more complicated by the fact that the pulsar wind is relativistic and magnetized. Moreover the ISM is typically partly neutral and the neutral atoms may have collisional mean free paths comparable to the scale length of the system, then invalidating a purely fluid treatment. However as demonstrated in a previous paper (Paper I), if we suppose that the H atoms can be ionized only via collisions with electrons and protons of the shocked plasma, for a large number of pulsars the presence of a neutral component in the ISM can be neglected as far as the fluid dynamics is concerned, or at least it can be taken into account as a small perturbative effect.

The interaction of the relativistic magnetized pulsar wind with the ISM ionized component is not direct, because the relativistic particles themselves have very small cross sections, but is mediated by the magnetic field advected by the pulsar wind and compressed on the head of the nebula. The effective mean free path of particles is then their gyroradius, which is typically much smaller than the typical dimension of the nebula, scaling with the distance of the bow shock stagnation point from the pulsar:

$$d_\mathrm{o}=\sqrt{/4\pi c\rho _\mathrm{o}V_\mathrm{o}^2},$$
(1)
where $``$ and $`V_\mathrm{o}`$ are the pulsar luminosity and peculiar velocity, and $`\rho _\mathrm{o}`$ the density of the dragged component of ambient medium (Paper I). For this reason we have assumed little mixing between the two media, which allows the use of a fluid approach to investigate the structure and dynamics of these objects.

In Section 2 we present the numerical code and the parameters of the various simulations we have performed; Sect. 3 discusses the results of these simulations and compares them to what is expected from an analytic “two thin layer” model (Comeron & Kaper comeron98 (1998)); in Sect. 4 we evaluate the penetration thickness of the external layer for a neutral hydrogen atom of the ISM.

## 2 Simulation Setup and Initial Conditions

To simulate a spherically symmetric wind from a moving star, interacting with the ambient medium, we have used the multi-D code CLAWPACK (LeVeque leveque94 (1994)): a second order Godunov (Godunov godunov59 (1959)) with piecewise linear reconstruction and transverse flux corrections. The problem has an intrinsic cylindrical symmetry which allows us to use a two-dimensional grid with Strang Splitting (Strang strang68 (1968)) for the geometrical terms. The boundary conditions we have imposed are: 1. At every time step the values of all relevant quantities inside a sphere of given radius centered on the star are reset to the free flowing solution; 2. A uniform wind enters the computational box from the right hand side. We have put reflecting conditions on the axis of symmetry. On the left hand side of the box there is a free outflow boundary condition: all values of variables in the ghost zones are equal to the values in the corresponding active zones (extrapolating the flow beyond the boundary), taking care to avoid the possible inflow of material due to roll-up features (Stone & Norman stone92 (1992), LeVeque & al. lev97 (1997)). This boundary condition, despite its simplicity, works quite well, and is actually used by various codes. We know from solar-wind simulations that the choice of such boundary conditions may influence the evolution in the case of subsonic outflow, but even in more sophisticated codes, the only way to solve the problem seems to take boundaries very far from the zone of interest (Pauls & al. pauls95 (1995)). We have tested the code with our simple buondary, and we have found that there are no significant problems with reflected spurious waves, and even following the evolution of our simulations we do not find spurious wave. Both the stellar and the external winds are taken with high Mach numbers (hypersonic limit). In fact the pulsar wind has a thermal pressure much lower than its ram pressure and the same holds true for the ISM, because pulsars have typical peculiar velocities of 200 km s-1(Cordes & Chernoff cordes98 (1998)) while the typical temperature of the ISM is of the order of $`10^410^5`$ K.

The pulsar wind is highly relativistic, with a Lorentz factor up to $`10^6`$ (Rees & Gunn rees73 (1973); Kennel & Coroniti kennel84 (1984)), thus in order to handle it properly a relativistic code is required; however the geometry of such nebulae is essentially determined by the momentum flux of the winds (Wilkin wilkin96 (1996)), so even a classical model can give a reasonable approximation. We decided to use a different adiabatic coefficient for the pulsar wind material and the ambient medium, so we have modified the Roe Riemann Solver (Roe roe81 (1981)) used by the code, adding a tracer $`\varphi `$ to separate the ambient medium (with adiabatic coefficient $`\gamma _1=5/3`$) from the relativistic (with $`\gamma _2=4/3`$) fluid coming from the star (Shyue shyue98 (1998), Karni karni94 (1994)). We then integrate the set of equations:

$`{\displaystyle \frac{}{t}}\left(\begin{array}{c}\rho \\ \rho v_r\\ \rho v_z\\ E\\ \rho \varphi \end{array}\right)+{\displaystyle \frac{}{z}}\left(\begin{array}{c}\rho v_z\\ \rho v_rv_z\\ \rho v_z^2+P\\ \rho v_zH\\ \rho \varphi v_z\end{array}\right)+`$ (2)
$`{\displaystyle \frac{}{r}}\left(\begin{array}{c}\rho v_r\\ \rho v_r^2+P\\ \rho v_zv_r\\ \rho v_rH\\ \rho \varphi v_r\end{array}\right)+{\displaystyle \frac{1}{r}}\left(\begin{array}{c}\rho v_r\\ \rho v_r^2\\ \rho v_zv_r\\ \rho v_rH\\ \rho \varphi v_r\end{array}\right)=0,`$
where $`\rho `$ is the density, $`v_z`$ and $`v_r`$ are the components of the velocity, $`P`$ $`E`$ and $`H`$ are the pressure, energy density and enthalpy of the fluid. $`\varphi `$ takes the value $`\varphi _1`$ for the ambient medium and $`\varphi _2`$ for the relativistic material. So, the value of the adiabatic coefficient $`\gamma `$ in each cell is determined by the relation:

$$\frac{1}{\gamma 1}=\frac{1}{\gamma _11}\frac{\varphi _2\varphi }{\varphi _2\varphi _1}+\frac{1}{\gamma _21}\frac{\varphi _1\varphi }{\varphi _1\varphi _2}$$
(3)
The physical conditions of such a bow-shock guarantee conservation of energy, and little mixing between the two media, so that we can use a fluid treatment. However we added a little diffusion to avoid the growth of numerical instabilities (“carbuncle”) due to the solver, and associated with grid aligned shocks. In fact these instabilities develop on the axis and affect both the bow-shock and the contact discontinuity and then move towards the tail of the nebula, leading to a complete disruption of its structure. The introduction of a little diffusion, which avoids the growth of numerical instabilities, but may also damp the physical ones, should not be considered only as a numerical trick; in fact the four known nebulae show a well-defined shape, so there must be some process (perhaps connected to the magnetic field) which tends to stabilize the flow.

### 2.1 Choice of initial parameters

We have performed two different types of simulations: the first one limited to the head of the nebula, using a finer grid, and the second one extended to the tail up to about 4 times the distance from the stagnation point to the star, with a coarser resolution.

The unit length we have used, Ul, has the value $`10^{16}`$cm, chosen to reproduce typical pulsar bow-shock nebula dimensions (Eq. 1). The simulations of the head were performed on a 0.40Ul$`\times `$0.58Ul box, using a grid of 480$`\times `$700 cells with the star located in (100,0). For the simulations of the tail, we used a 1.45Ul$`\times `$1.00Ul box with a grid of 580$`\times `$400 cells and the star in (420,0). The stellar wind has a velocity of 1000 km s-1, a density of 0.04 cm-3 and a Mach number of about 13 upstream of the termination shock along the bow shock axis. We have performed various simulations using different velocities for the external flow (400, 300, 200, 150 km s-1) and densities varying accordingly, in order to keep the ram pressure constant. The following figures refer to the simulation of 400 km s-1(density of 0.25 cm-3). The time needed to reach the steady state condition is of the order of 4-5 times the external flow crossing time in the computational box.

## 3 Comparison with the analytic model

### 3.1 Shape

The geometry and internal structure of the bow shock nebula resulting from our simulations are shown in Fig. 1 and Fig. 2, while a grey-scale contour plot of the density is shown in Fig. 3. The positions of the termination shock, bow shock and contact discontinuity do not change among the mentioned above simulations, and also the position of the sonic point remains essentially the same. This is because the solution, in the hypersonic limit, depends only on the ratio between ram pressures of the stellar and ambient medium, but not separately on their densities or velocities: therefore we are confident that we have approximated well the hypersonic limit. The dashed lines represent three analytic solutions (Wilkin wilkin96 (1996)), one scaled to match the contact discontinuity near the head, one to match the bow-shock near the head, and the long-dashed one, which is evaluated for the correct wind/environment ram pressure ratio. It can be seen that the analytic model reproduces quite well the shape of the nebula, although an uncertainty remains in the size. Fig. 1 shows that the region of unperturbed stellar wind extends in the tail up to a “spherical shock”, whose position is strictly connected to the value of the thermal pressure in the ambient medium. In fact in the tail the pressure tends to the values of the ambient medium $`P_\mathrm{o}`$ so the distance of the spherical shock can be easily derived as:

$$d_{back}\sqrt{/4\pi cP_\mathrm{o}}=d_\mathrm{o}M_\mathrm{o}\gamma _1^{1/2},$$
(4)
where $`M_\mathrm{o}`$ is the external flow Mach number.

We remark that the outflow boundaries we use may introduce spurious waves, and this plays a role in the structure of the subsonic region (region A) especially as it may change the position of the triple point (where the spherical shock meets the termination shock). However the head of the simulated nebula, and the supersonic region of the external shocked layer are not affected, so that the following comparison with the analytic model is still valid.

Looking at the head (Fig. 2) we notice a bump in the contact discontinuity, due to an instability that has not been damped completely by diffusion. The contact discontinuity in the tail (Fig. 1 and Fig. 3) tends to reach a limiting distance from the axis (i.e. the shocked pulsar wind is confined to a cylinder), with radius of the same order of the stagnation point distance; instead, the external bow-shock tends to reach the Mach cone.

An upper limit to the radius of the contact discontinuity can be evaluated using an asymptotic treatment. Asymptotically in the tail the bow shock degenerates in a Mach cone, so the pressure in the external layer is essentially the same as the thermal pressure of the ambient medium $`P_\mathrm{o}`$. Inside the contact discontinuity there are two different regions: the innermost (labeled as region A) corresponds to the material shocked in the spherical shock, which moves subsonically; the other region (region B), which surrounds the previous one, is filled with the material which crossed the termination shock at distances from the pulsar less than the spherical shock. As shown in Fig. 1 and Fig. 2, the material in the latter region moves supersonically. In a steady state model the pressure in these two regions must asymptotically match the conditions in the ambient medium. So the pressure inside the contact discontinuity is also $`P_\mathrm{o}`$ and this justifies Eq. (4) for the position of the spherical shock.

If we consider a fluid particle moving in region B, from a position near the stagnation point (labeled with subscript ini) towards the asymptotic region in the tail (labeled with subscript asy), the entropy conservation requires:

$$P_{asy}/\rho _{asy}^{\gamma _2}=P_{ini}/\rho _{ini}^{\gamma _2}$$
(5)
On the other hand we have: $`P_{asy}=P_\mathrm{o}`$, and:

$$P_{ini}=A(\gamma _1)\rho _\mathrm{o}V_\mathrm{o}^2,\rho _{ini}=\rho _{ups}B(\gamma _2),$$
(6)
where $`\rho _{ups}`$ is the value of the density upstream of the wind termination shock along the axis, and :

$`A(\gamma )=`$ $`\left({\displaystyle \frac{\gamma +1}{2}}\right)^{\frac{\gamma +1}{\gamma 1}}\left({\displaystyle \frac{1}{\gamma }}\right)^{\frac{\gamma }{\gamma 1}},`$
$`B(\gamma )=`$ $`{\displaystyle \frac{\gamma +1}{\gamma 1}}\left({\displaystyle \frac{\gamma +1}{2}}\right)^{\frac{2}{\gamma 1}}\left({\displaystyle \frac{1}{\gamma }}\right)^{\frac{1}{\gamma 1}}.`$ (7)
If we consider an external wind with high Mach number, the position of the spherical shock is far from the star so we can assume that all the material of the stellar wind is confined to region B and that the distance $`Z`$ of the contact discontinuity from the axis is much larger than the radius of internal region A. The matter flux conservation law allows us to write:

$$4\pi d_\mathrm{o}^2\rho _{ups}V_{wind}V_{asy}\rho _{asy}\pi Z^2$$
(8)
From a two thin layer model we found that $`V_{asy}`$, taken to be constant over the whole transverse section, is close to $`V_{wind}`$, the stellar wind velocity. This can be easily understood because, if we assume that the tail extends to infinity (this constraint will be released below), the velocity scale for the stellar wind material, in the absence of mixing with the ambient medium, is $`V_{wind}`$ (in the hypersonic limit the dependence on the external flow Mach number may be neglected). More specifically from the Bernoulli equation:

$$\frac{V_{asy}^2}{2}+\frac{\gamma _2}{\gamma _21}\frac{P_{asy}}{\rho _{asy}}=\frac{\gamma _2}{\gamma _21}\frac{P_{ini}}{\rho _{ini}},$$
(9)
together with Eq. 5, it follows:

$$V_{asy}=V_{wind}\sqrt{1\left\{\gamma _1M_\mathrm{o}^2/A(\gamma _1)\right\}^{\frac{1\gamma _2}{\gamma _2}}}$$
(10)
Substituting $`\rho _{asy}`$ and $`V_{asy}`$ in Eq. (8), using respectively Eq. (5) and Eq. (10), we obtain:

$`Z^2`$ $`4d_\mathrm{o}^2\left\{\gamma _1M_\mathrm{o}^2A(\gamma _1)\right\}^{\frac{1}{\gamma _2}}/`$ (11)
$`B(\gamma _2)\sqrt{1\left\{\gamma _1M_\mathrm{o}^2/A(\gamma _1)\right\}^{\frac{1\gamma _2}{\gamma _2}}}`$
This can be considered as an upper limit to the radius of the contact discontinuity. For the stellar wind material to expand sideways beyond $`Z`$, some process is required that reduces the velocity in the layer or increases its pressure. This can be achieved via processes like mass accretion, or if some further shock develops.

The distance of the bow-shock from the star given by the analytic model is 0.175 Ul (long-dashed line) while the two limiting analytic solutions of Fig. 1 and Fig. 2 (dashed lines) are at 0.22 and 0.29 Ul. So the analytic model for the correct wind/environment ram pressure ratio gives a value which corresponds to the termination shock, while the thickness of the external layer is 0.025 Ul, in good agreement with what was expected (Paper I). Even though the shocked layers of the nebula are not thin at all, we see that the analytic solution reproduces quite well the morphology, at least of the external layer, which corresponds to the H$`\alpha `$ emitting region. In the tail simulation we see that the contact discontinuity has not reached the asymptotic regime, its distance from the axis being less than the maximum value given by Eq. (11).

The behaviour of the fluid in regions A and B may affect the synchrotron emission from the relativistic electrons and positrons coming from the pulsar and shocked in the termination shock. As we can see from Fig. 1, material flowing in region B becomes supersonic soon in the head. In a relativistic case, the Lorentz factor associated with the bulk motion of the fluid in the tail becomes greater than the Lorentz factor associated with the thermal motion of relativistic particles. If particles have relativistic motions their synchrotron emission can be detected only if their pitch angle is close to the observer direction. In the tail the relativistic bulk motion produces a beaming of the emission in the backward direction, so that the synchrotron emission from relativistic particles, in the case of a direction of observation perpendicular to the pulsar velocity direction, can be detected only in the head, while it fades in the tail. This holds true for the less energetic particles responsible for the radio emission and may even be true for those responsible for the X-ray emission. In region A the flow is subsonic, so the radio and X-ray emission can be detected. In this case we expect a cylindrical emitting region. If the spherical shock is far from the pulsar it may occur that the number of high energy particles confined in region A is so small that their emission might be under the threshold of detection, so that this region can be observed only in the radio band.

### 3.2 Internal magnetic field

Using the numerical hydrodynamic solutions described above we have also evaluated the behaviour of a toroidal magnetic field wound around the symmetry axis and passively advected by the fluid. We considered the case of a pulsar wind with the energy associated with the ordered magnetic field much smaller than the ram pressure of the relativistic particles, and with the ratio between the two being independent of the direction. The value of such a field increases in the head, due to the compression in the termination shock, ranging from 7 times its upstream value, right after the shock, up to a factor of 30-35 with respect to its upstream value, near the contact discontinuity. In the tail it tends to decrease with respect to its value in the head due to the subsequent expansion of the fluid until it approaches a constant value:

$`B_{asy}`$ $`{\displaystyle \frac{\pi }{2}}B_{ups}\left\{\gamma _1M_\mathrm{o}^2A(\gamma _1)\right\}^{\frac{1}{2\gamma _2}}B(\gamma _2)^{\frac{1}{2}}\times `$ (12)
$`\left\{1\left(\gamma _1M_\mathrm{o}^2/A(\gamma _1)\right)^{\frac{1\gamma _2}{\gamma _2}}\right\}^{1/4}`$
In fact the contact discontinuity tends to a limit section and all the magnetic field produced by the star is confined in such region. If the original magnetic field is not too far below equipartition, it can dominate the fluid dynamics both in the tail and in the head. So an hydrodynamical approach is asymptotically valid if:

$$\beta _{asy}=\frac{8\pi P_\mathrm{o}}{B_{asy}^2}1$$
(13)
Following a procedure analogous to the one used to derive the asymptotic axial distance of the contact discontinuity we are able to relate the asymptotic condition to the initial values in the head:

$`\beta _{asy}`$ $`\left({\displaystyle \frac{2}{\pi }}\right)^2\beta _{ini}\left\{\gamma _1M_\mathrm{o}^2A(\gamma _1)\right\}^{\frac{1\gamma _2}{\gamma _2}}B(\gamma _2)\times `$ (14)
$`\left\{1\left(\gamma _1M_\mathrm{o}^2/A(\gamma _1)\right)^{\frac{1\gamma _2}{\gamma _2}}\right\}^{1/2}`$
So, even if initially the magnetic field’s energy is below equipartition, it can reach equipartition if the external wind Mach number is high enough. Full MHD simulations are under development.

### 3.3 Surface density and tangential velocity in the external layer

To verify that the two thin layer model could be taken as a reasonable approximation of the true structure of the nebula, at least as far as the external layer is concerned, and to evaluate how far in the tail it could be used, we have extracted from our simulations the surface density and the tangential velocity in the external layer, between the bow shock and the contact discontinuity.

Fig. 4 compares the surface density in the external layer, extracted from the tail simulation, with that evaluated by a two thin layer model. The numerical values have been obtained by integrating in the direction perpendicular to the bow shock. The difficulty in determining how good the analytic model approximation is follows from the fact that the analytic model depends on just one parameter, while the real problem has intrinsically two length scales (the stagnation point distance and the thickness in the head). In fact we have found two limit solutions (see Fig. 1 and Fig. 2), so that we could put only an upper and a lower limit to the length scale that should be used in the analytic model. The analytic curve that best reproduces the numerical points gives an external density of 0.26 cm-3 (instead of the true value 0.25 cm-3) and a stagnation point distance of 0.235 Ul, compatible with the range of the two limit solutions. The fluctuations of the numerical points are caused by round-off errors and uncertainties in the evaluation of the proper orthogonal direction, due to small oscillations in the shape of the bow shock, caused by discretisation errors in the evaluation of bow-shock position, eventually enhanced by the interpolation routine we use to find the orthogonal direction. In the simulation of the head we have noticed that, near the axis, the density values are lower than expected, and this effect corresponds to the bump of the contact discontinuity, that we have previously discussed. If corrected for the presence of the bump also in this zone we find reasonable values for the density.

In the same way, we have evaluated the tangential velocity in the layer between the bow shock and the contact discontinuity. Fig. 5 shows the numerical values of the tangential velocity in the external layer, normalized to the external wind value, and compares it to the two thin layer model. The analytic curve that reproduces the points has been obtained using a value of 0.25 Ul for the distance of the stagnation point, a value in agreement with the two limit solutions. The fluctuations in the numerical points are much smaller than those seen in the surface density, which implies that the velocity pattern is more uniform in the layer with respect to the density.

The fact that the analytic curves that best reproduce the numerical values correspond to different values of the stagnation point distance tells us that the density and velocity cannot be reproduced exactly by the same thin layer model. It should be stressed, though, that, in spite of the non-thin structure of the external layer, the thin layer model is not a bad approximation. Deviations from the continuous curves are present only in the tail where a steady state regime may have not been achieved yet.

## 4 Penetration thickness for neutral hydrogen

As we have demonstrated in a previous paper (Paper I) the penetration thickness $`\tau `$ for neutral hydrogen in the external layer, depends on the ambient density and the pulsar velocity. For a pulsar moving in a partially neutral medium (the ISM ionization fraction is of order 0.1; Taylor & Cordes taylor93 (1993); Dickey & Lockman dickey90 (1990)), this parameter plays an essential role in determining the true size of the nebula and the H$`\alpha `$ photon flux. The penetration thickness depends on the process that could drag the atoms: essentially by charge-exchange and by ionization via the interaction with the protons confined in the external layer. The cross sections of such processes are strongly dependent on velocity (McClure mcclure66 (1966); Newman et al. newman82 (1982); Peek peek66 (1966); Ptak & Stoner ptak73 (1973)), and variations of the surface density in the external layer give rise to changes of the thickness moving away from the stagnation point. In our previous paper we have evaluated the thickness using the two thin layer model limited to a region near the axis. Our aim is to verify if such an approximation is acceptable and how far in the tail can it be extended.

The thickness of the external layer has been evaluated integrating the ionization and charge-exchange probabilities along the path of a neutral test hydrogen atom. So we have taken into account variations of density, relative velocity and temperature. We have verified that, for typical pulsar velocities, collisional ionization can be neglected, even in non equilibrium cases.

Fig. 6 shows the penetration thickness, normalized to the external density. The differences in the values for the various simulations are due to variations of the cross section with velocity, while the effect of temperature is less sensitive. The thickness decreases at higher velocities as expected: in fact it follows the behaviour of the cross section for charge-exchange, which is the dominant process. Concerning the variation of $`\tau `$ with the distance from the stagnation point, we can see that it has a homogeneous trend for different velocity values. Moving away from the stagnation point, the surface density increases and the same holds for the cross section of charge-exchange, due to the decrease of relative velocity between the neutral hydrogen atoms coming from ISM and the protons flowing in the external layer. This effect, that tends to increase the penetration thickness, is however reduced because the interaction rate, which depends on the relative velocity, decreases even if the cross section is higher, so that the mean free path is longer. The thickness, at least for the head, ranges only over a factor 2. Thus we can say that the model restricted to the stagnation point, which we have proposed (Paper I), can reproduce within the same factor the behaviour of the whole head. This can be quite important, because, if the nebula is thin to the penetration of neutral hydrogen, it can remain thin in the tail, even far from the pulsar.

## 5 Conclusion

We have performed hydrodynamical simulations of the interaction with an ambient homogeneous medium, of the wind of a pulsar moving supersonically. The use of a two component (the ambient medium and the pulsar wind) code allows us to verify directly if the two thin layer model gives a reasonable approximation. We have also extended our simulations to the tail. Our results show that the analytic model is a good approximation for the shape, the surface density and the tangential velocity in the external layer, and can be extended in the tail up to about 3 or 4 times the stagnation point distance from the pulsar. For the range of pulsar peculiar velocities that we have take into account, we have verified that the penetration thickness for neutral hydrogen in the head of the nebula has a rage of a factor of 2.

This tells us that models based on the analytic two thin layer solution could be considered a reasonable starting point, over which we might introduce effects of mass loading. So it will be possible to derive a simple solution for the case in which the neutral component of the ambient medium cannot be considered as a small perturbation. The analysis we have performed gives information on the asymptotic behaviour of such nebulae. As we have pointed out in the introduction, two out of the four known objects seem not to correspond to our expectations. This is evidence that in such cases we cannot assume a purely hydrodynamical steady state condition, and we require some processes in the tail (mass loading, shocks) to model the observations. However the observational data on these objects are not good enough to allow for a clear choice of the more adequate description. Our simulations show that the penetration thickness of the nebula changes slowly moving away from the head, thus, if we assume that ionizing processes other than collisions are active, mass loading may affect the fluid motions not only in the head but also in the tail. It is known that mass loading or charge-exchange from slowly moving particles in a supersonic fluid tends to reduce the velocity and to increase the pressure. The shocked material coming from the stellar wind becomes supersonic soon in the head. In the supersonic region, mass loading due to ionization of the neutral atoms, penetrating from the ambient medium, may be due to ionizing photons coming from the hot shocked ISM, while photoionization from the pulsar may be neglected (Paper I). So mass loading may justify the presence of a conical tail, which requires an internal thermal pressure greater than what we found from purely hydrodynamical simulations, and this suggests the necessity of ionizing processes other than collisions.

We have also pointed out that the magnetic field might play an important role both in the head and in the tail, leading to a faster collimation of the flow. The presence of a component of the magnetic field wounded around the axis of the nebula may reduce the effect of mass loading preventing the flow in region B from expanding sideways and collimating it. In fact in the case of the “Guitar Nebula” (the only for which we know the shape even far from the pulsar) we see that there is a conical tail near the star (which we suppose due to mass loading) which stops increasing and become a cylindrical tube. To justify this change in the shape we need some other process (that we believe could be the action of magnetic fields, which may reach equipartition) that collimates the tail. Full MHD simulations are under development, as well as relativistic HD simulations, and will be presented in a future paper.

###### Acknowledgements.

I want to acknowledge Luca Del Zanna for the help in modifying and testing the numerical code, and Elena Amato for the fruitful discussion and the reading of this article. This work has been supported by the Italian Ministry for University and Research (MIUR) under grant Cofin2001–02–10