text:
# The Anisotropy of MHD Alfvénic Turbulence

## 1 Introduction

Many astrophysical plasmas, including the interstellar medium and the solar wind, often show magnetic fields whose energy density is greater than or equal to the local kinetic energy density. In these plasmas the magnetic fields play a dominant dynamical role, mediated by magnetohydrodynamic (MHD) waves. In the incompressible limit, there are only two types of linear modes: shear Alfvén waves and pseudo Alfvén waves. While these two modes have different polarization directions, they have the same dispersion relation and propagate along the magnetic field lines at the Alfvén speed. Therefore, the nonlinear interactions of wave packets moving along the magnetic field lines at the Alfvén speed determine the dynamics of incompressible magnetized plasmas with a strong background field. In this paper, we study the anisotropy of the MHD turbulence in this regime. We will refer to this turbulence as incompressible Alfvénic turbulence.

Nonlinear processes and the corresponding energy spectrum of incompressible Alfvénic turbulence are still among the most controversial problems in MHD. Since the pioneering works of Iroshnikov (1963) and Kraichnan (1965), the Iroshnikov-Kraichnan (IK) theory has been widely accepted as a model for incompressible, highly conducting MHD turbulence. The IK theory predicts $`E_M(k)E_K(k)k^{3/2}`$ from a Kolmogorov-like dimensional analysis. Here, $`E_M(k)`$ and $`E_K(k)`$ are the magnetic and kinetic energy spectra respectively. In this framework, two counter-traveling eddies (i.e. Alfvén wave packets) interact and transfer energy to smaller spatial scales only when they collide, as they move in opposite directions along the magnetic field lines. Since the duration of such a collision is shorter than the conventional eddy turnover time by a factor of $`t_v(l)/t_A(l)`$, this collisional process is inefficient and the spectral energy transfer time as a function of scale $`l`$ ($`=t_{cas}(l)`$) increases by the same factor compared to the eddy turnover time ($`l/v_l`$) in ordinary hydrodynamic turbulence. Here $`t_v(l)=l/v_l`$ and $`t_A(l)=l/V_A`$ are eddy turnover time and Alfvén time respectively, $`V_AB/\sqrt{4\pi \rho }`$, and $`B`$ is the rms magnetic field strength. When the external field is strong, as assumed in IK theory, this quantity is usually set to $`B_0`$, the strength of the uniform background field. If the spectral energy cascade rate

$$ϵ\frac{v_l^2}{t_{cas}(l)}\frac{v_l^3}{l}\frac{t_A(l)}{t_v(l)}$$
(1)
is scale-independent and $`E_M(k)E_K(k)`$, then we obtain the IK energy spectra.

The IK theory assumes an isotropic distribution of energy in $`𝐤`$-space. However, many researchers have argued that anisotropy is an important characteristic in MHD turbulence (for example, Shebalin et al 1983, Montgomery and Matthaeus 1995). This anisotropy results from the resonant conditions for 3-wave interactions (or 4-wave interactions, when 3-wave interactions are null). The resonant conditions for the 3-wave interactions are

$`𝐤_1+𝐤_2`$ $`=`$ $`𝐤_3,`$ (2)
$`\omega _1+\omega _2`$ $`=`$ $`\omega _3,`$ (3)
where $`𝐤`$’s are wavevectors and $`\omega `$’s are wave frequencies. The first condition can be viewed as momentum conservation and the second as energy conservation. Alfvén waves satisfy the dispersion relation: $`\omega =V_A|k_{}|`$, where $`k_{}`$ is the component of wavevector parallel to the background magnetic field. Since only opposite-traveling wave packets interact, $`𝐤_1`$ and $`𝐤_2`$ must have opposite signs. Then from equations (2) and (3), either $`k_{,1}`$ or $`k_{,2}`$ must be equal to 0. That is, zero frequency modes are essential for energy transfer. If $`k_{,2}=0`$, we have

$`k_{,1}`$ $`=`$ $`k_{,3},`$ (4)
$`k_{,2}`$ $`=`$ $`0`$ (5)
(Shebalin et al 1983). Therefore, in the wavevector space, 3-wave interactions make energy cascade in directions perpendicular to the mean magnetic field. Since the energy cascade is strictly perpendicular to the mean magnetic field, the actives modes in wavevector space have a slab-like geometry with a constant width. The implication is that the nonlinear cascade of energy works against isotropy in $`𝐤`$ space. Furthermore, it is important to note that equations (2) and (3) are true only when wave amplitudes are constant. In reality, nonlinear interactions provide a natural broadening mechanism, following the uncertainty relation, $`\mathrm{\Delta }t\mathrm{\Delta }\omega 1`$. In particular, if a wave has a frequency less than or comparable to the nonlinear interaction rate, it is effectively a zero frequency mode.

Goldreich and Sridhar (1995, 1997) showed that in the strong incompressible shear Alfvénic turbulence regime (i.e. $`\tau _{NL}^1\omega `$), these arguments lead to a new scaling law with a scale-dependent anisotropy. In this model smaller eddies are more elongated. Their arguments are based on the assumption of a critically balanced cascade, $`k_{}V_Ak_{}v_l`$, where $`k_{}`$ and $`k_{}`$ are wave numbers perpendicular and parallel to the external dc field111Later in this paper, we will use $`\stackrel{~}{k}_{}`$ and $`\stackrel{~}{k}_{}`$, instead of $`k_{}`$ and $`k_{}`$, to represent their scaling relation in a slightly different, yet we believe equivalent, viewpoint. Here, $`\stackrel{~}{k}_{}`$ and $`\stackrel{~}{k}_{}`$ are wavenumbers with respect to the direction of the local magnetic field, not the external field.. The argument given above for 3-wave interactions makes it clear that $`k_{}`$ will tend to increase until it becomes important in the plasma dynamics. The assumption of strong nonlinearity implies that wave packets lose their identity after they travel one wavelength along the field lines. Consequently the eddy turnover time ($`(k_{}v_l)^1`$) is actually the same as Alfvénic time ($`(k_{}V_A)^1`$). In this model, the cascade time, $`t_{cas}(l)`$ can be determined without ambiguity: $`t_{cas}(k_{}v_l)^1(k_{}V_A)^1`$. Since the cascade time is comparable to the period of the Alfvén wave, the 3-wave resonant condition can be violated according to the uncertainty relation $`\mathrm{\Delta }\omega \mathrm{\Delta }tV_Ak_{}t_{cas}1`$. The quantity $`k_{}`$ is the width of the active region in wavevector space. Finally the assumption of a scale-independent cascade rate $`ϵv_l^2/t_{cas}(l)E_{waves}V_A/L`$ gives

$`k_{}k_{}^{2/3}L^{1/3}\left({\displaystyle \frac{E_{waves}}{V_A^2}}\right)^{1/3},`$ (6)
$`v_lV_A(k_{}L)^{1/3}\left({\displaystyle \frac{E_{waves}}{V_A^2}}\right)^{1/3},`$ (7)
where $`E_{waves}`$ is the wave energy per mass. These formulae assume that all scales, from $`L`$ on down, are within the inertial range of MHD turbulence. Here the typical $`k_{}`$ should be interpreted as the size of the range of parallel wavevectors, corresponding to a given $`k_{}`$, that contain significant energy.

Matthaeus et al. (1998) recently tested this model numerically, and showed that the anisotropy of low frequency MHD turbulence scales linearly with the ratio of perturbed and total magnetic field strength ($`b/B`$), a result which seems inconsistent with Goldreich and Sridhar’s model. To explain this scaling relation, they suggested that the region of Fourier space where the energy transfer takes place actively is given by

$$|𝐤\frac{𝐁_0}{\sqrt{4\pi \rho }}|<\frac{1}{\tau _{NL}},$$
(8)
where $`\tau _{NL}`$ is the eddy turnover time of the energy containing length $`L`$. Consequently, “the region of the wave number space where spectral transfer is most vigorous” has a slab-like geometry with a constant width proportional to $`1/(\tau _{NL}B_0)`$.

All these theories (except the IK theory) share a common prediction for anisotropy: anisotropy should be more pronounced on smaller scales. Oughton et al. (1994) and Ghosh and Goldstein (1997) already reported this scale-dependent anisotropy through numerical simulations. The former extended Shebalin et al. (1983)’s 2-D calculations to 3-D cases. To measure anisotropy, they introduced the Shebalin angles, $`\theta _Q`$, defined by $`\mathrm{tan}^2\theta _Q=(k_{}^2|𝐐(𝐤,t)|^2)/(k_{}^2|𝐐(𝐤,t)|^2),`$ where $`𝐐`$ is vector potential A, magnetic field B, or current J, etc. Greater $`\theta _Q`$ means greater anisotropy. They found that $`\theta _A<\theta _B<\theta _J`$. If the energy spectrum of B scales as $`E_M(k)k^s`$, then the spectra of vector potential and current scale as $`E_A(k)k^{s2}`$ and $`E_J(k)k^{s+2}`$, respectively. The spectra of vector potential has the steepest slope among the three spectra. This means that the vector potential is least strongly dependent on small scales (and the current is most strongly dependent on small scales). Therefore they concluded that anisotropies are more pronounced at smaller scales222 The width of the active region in Fourier space, $`k_{}`$, is a function of $`k_{}`$. If
$$\left[\frac{k_{}}{k_{}}\right]_{\text{small }k_{}}>\left[\frac{k_{}}{k_{}}\right]_{\text{large }k_{}},$$
then we will have the ordering of Shebalin angles as observed by Shebalin et al. and Oughton et al. However, if $`k_{}k_{}`$, then we do not expect any ordering among the angles. In §3, we will find that, no matter what the true functional form $`k_{}=k_{}(k_{})`$ is, Fourier transformation smooths out the true relation and leads to a fake linear relationship between $`k_{}`$ and $`k_{}`$. This suggests that their results are in contradiction to our discussion in §3. However, the apparent linear relationship between $`k_{}`$ and $`k_{}`$ in §3 is not perfectly linear. Instead, we expect $`k_{}=c_1k_{}+c_2`$, where $`c_1`$ is a decreasing function of $`B_0`$ and $`c_2`$ depends on the $`k_{}`$ of forcing terms (or initial values of $`k_{}`$ for decaying turbulence). The presence of $`c_2`$ does not seem to be important in our Fig. 3 in §3. However, it does affect the calculation of Shebalin angles. That is, because of $`c_2`$, the ratio $`k_{}/k_{}=c_1+c_2/k_{}`$ becomes a function of $`k_{}`$. At the largest energy containing eddy scale, $`k_{}c_2`$ and hence $`k_{}/k_{}c_1+O(1)O(1)`$. But, at small scales, the ratio can be much smaller than unity. Therefore, because of $`c_2`$, we can obtain a scale-dependent anisotropy: $`(k_{}/k_{})`$ at small $`k_{}`$ is greater than that at large $`k_{}`$. Note that this is a result of the initial conditions, rather than a true scaling relation. This will lead the ordering of the angles as observed by Shebalin et al. and Oughton et al.. This interpretation is qualitatively consistant with their results. For example. the ratio $`k_{}/k_{}=c_1+c_2/k_{}^{peak}`$, therefore $`(\mathrm{tan}\theta _Q)^1`$, approaches to a constant value $`c_2/k_{}^{peak}`$ as $`B_0`$ becomes strong. Here $`k_{}^{peak}`$ is the wavenumber of the peak of the energy spectrum. It is also possible to explain the increased anisotropy at high magnetic Reynolds numbers. If the magnetic Reynolds number increases, then $`k_{}^{peak}`$ increases and, therefore, the ratio decreases. . They also found a similar ordering for velocity field (and vorticity, etc). On the other hand, Ghosh and Goldstein (1997) calculated the Shebalin angles as a function of wavenumber bin. They found that the Hall-MHD333Hall MHD includes the Hall term, which is important at ion-cyclotron scales. In this paper, we consider only the standard MHD equations. simulations show increased anisotropy at small scales (i.e. greater Shebalin angles at smaller scales). However their standard MHD simulations do not show increased anisotropies at small scales. We refrain from comparing their work with ours because they used different physics (the Hall effect) and their simulations are $`2\frac{1}{2}`$ dimensional. The simulations given in this paper are 3-D standard MHD simulations. We note that none of the previous papers quantitatively compared their results with particular theories of anisotropy.

In this paper, we examine the scaling law for Alfvénic MHD turbulence numerically and resolve the controversy concerning the anisotropic structure of the turbulence. Our results are consistent with Goldreich and Sridhar’s model. We describe our numerical methods in §2. In §3, we describe our results for anisotropy in wavevector space. In this section, we demonstrate that none of the scaling laws mentioned above agrees with our data, and explain why a straightforward evaluation of the distribution of spectral power does not correspond to a physically meaningful set of scaling laws. We attribute this to the systematic effects of large scale curvature of the magnetic fields. In §4, we determine the shape of individual eddies, avoiding the systematic error described above. We compare our results to Goldreich and Sridhar’s model and give our conclusions in §5.

## 2 Numerical Methods and General Results

We have employed a pseudospectral code to solve the incompressible MHD equations in a periodic box of size $`2\pi `$:

$$\frac{𝐕}{t}=(\times 𝐕)\times 𝐕V_{A0}^2(\times 𝐁)\times 𝐁+\nu ^2𝐕+𝐟+P^{},$$
(9)
$$\frac{𝐁}{t}=𝐁𝐕𝐕𝐁+\eta ^2𝐁,$$
(10)
where $`𝐟`$ is random driving force and $`P^{}P+𝐕𝐕/2`$. Other variables have their usual meaning. $`V_{A0}=B_0/(4\pi \rho )^{1/2}`$ is the Alfvén velocity of the background field, which is set to be order unity in our simulations. In pseudo spectral methods, we calculate the temporal evolution of the equations (9) and (10) in Fourier space. To obtain the Fourier components of nonlinear terms, we first calculate them in real space, and transform back into Fourier space. We use 21 forcing components with $`2k\sqrt{12}`$. Each forcing component has correlation time of one. The peak of energy injection occurs at $`k2.5`$. The amplitudes of the forcing components are tuned to ensure $`V1`$. Since we expect the perturbed magnetic field strength $`b`$ to be comparable to or less than $`V`$ (i.e. $`b1`$), the resulting turbulence can be regarded as strong MHD turbulence (i.e. $`VB_0B`$). The average helicity in these simulations is not zero. However, our results are insensitive to the value of helicity, possibly because the box size is only slightly bigger than the energy injection scale. We use an appropriate projection operator to calculate $`P^{}`$ term in Fourier space and also to enforce divergence-free condition ($`𝐕=𝐁=0`$). We use up to $`256^3`$ collocation points. We use integration factor technique for kinetic and magnetic dissipation terms and leap-frog method for nonlinear terms. We eliminate $`2\mathrm{\Delta }t`$ oscillation of the leap-frog method by using an appropriate average. At $`t=0`$, magnetic field has only uniform component and velocity has a support between $`2k4`$ in wavevector space.

We use hyperviscosity and hyperdiffusivity for dissipation terms to maximize the inertial range. The only exception is Run 256P-$`B_0`$1, where we use physical viscosity and diffusivity. The power of hyperviscosity is set to 8, so that the dissipation term in the above equation is replaced with

$$\nu _8(^2)^8𝐕,$$
(11)
where we set the value of $`\nu _8`$ from the condition $`\nu _h(N/2)^{2h}\mathrm{\Delta }t0.5`$ (see Borue and Orszag 1996). Here $`\mathrm{\Delta }t`$ is the time step and $`N`$ is the number of grids in each direction. We use exactly same expression for the magnetic dissipation term. We list parameters used for the simulations in Table 1. We use the notation 256X-Y, where X = H or P refers to hyper- or physical viscosity; Y = $`B_0`$1 or $`B_0`$0.5 refers to the strength of the external magnetic fields.

Fig. 1 shows the evolution of $`V^2`$ and $`B^2`$ from Runs 256H-$`B_0`$1 and REF2. Results of 256H-$`B_0`$1 are plotted as solid lines, and results of REF2 are plotted as dotted lines. Both runs have similar results for the overlapped time, which means that the values of $`V^2`$ and $`B^2`$ are not sensitive to the spatial resolution. The magnetic energy density grows fast at the beginning of the simulations, as a result of the stretching of external field lines by turbulent motion. Since the external field is strong, magnetic forces soon become strong enough to balance the stretching effect. This balance happens at $`t1`$ and, after a transient stage ($`1t5`$), the turbulence reaches the saturation stage ($`t5`$). In this case the saturation stage begins quite early, since the external uniform magnetic field is strong and the magnetic diffusivity is effectively 0. In general, when the external field is weaker and/or magnetic diffusivity is larger, the saturation stage occurs later. At the saturation stage, $`B^21.45`$ and $`V^20.6`$. Since $`b^2=B^2B_0^2(1)0.45`$, there is a rough energy equipartition between $`b`$ and $`V`$. Note that, since $`VB`$, the condition for strong MHD turbulence is met.

## 3 Fourier Analysis: A False Scaling Law?

### 3.1 Fourier Analysis

In this subsection, we investigate the spectral structure of Alfvénic turbulence. It is important to note that we lose some information about individual eddies when we perform a global transformation, such as the Fourier transformation. In particular, the scaling law given in this subsection may $`not`$ be true for individual eddies (See §3.2 for details).

We plot the 1-dimensional energy spectra $`E_K(k)`$ and $`E_M(k)`$ in Fig. 2. Both spectra peak at the same $`k`$, which is also the scale of the energy injection. This reflects that there is a rough energy equipartition between $`b`$ and $`V`$ at the largest energy containing scale. In general, $`E_K(k)`$ always peaks at the energy injection scale. However, $`E_M(k)`$ peaks at a larger k than $`E_K(k)`$ when the external field is weak. The inertial range exhibits of two different scaling ranges: a small $`k`$ range where the slopes are steep and a large $`k`$ range where the slopes are mild. We believe that the former reflects a true inertial range, whereas the latter is a result of the $`bottleneck`$ effect. The $`1/k`$ bottleneck effect is a common feature in numerical hydrodynamic simulations with hyperviscosity (see Borue and Orzag 1996, 1995). Interestingly enough, the slope of the bottleneck is actually steeper than $`1`$. This might mean that the bottleneck effect is less serious in the MHD case. The kinetic and magnetic energy spectra have slightly different powers in the inertial range. Although the power indexes of both spectra are close to $`5/3`$ in the (true) inertial range, the data are also compatible with IK theory, where the power index is $`3/2`$.

We plot normalized 3-dimensional energy spectra of Run 256H-$`B_0`$1 in Fig. 3. In the figure, we plot contours of same

$$E_3(k_{},k_{})/E_3(k_{},0)$$
(12)
in the $`(k_{},k_{})`$ plane. As the energy cascades in the directions perpendicular to the mean field, this ratio drops as we move away from the $`k_{}`$ axis. The figure implies most of the energy is confined in a region around the $`k_{}`$ axis. The thickness of the region depends on the strength of the external magnetic fields (Fig. 4). We use the contours of $`E_3(k_{},k_{})/E_3(k_{},0)=0.5`$ to measure the thickness. The angle $`\mathrm{\Theta }_{0.5}`$ is the angle formed by the contours and the horizontal axis. We can see that $`\mathrm{tan}\mathrm{\Theta }_{0.5}`$ ($`=k_{}/k_{}`$) is proportional to $`b/B_0`$. This result confirms the scaling relation found by Matthaeus et al. (1998). Both velocity and magnetic fields have very similar structures. From the fact that $`k_{}<k_{}`$ in the active region (see Fig. 3), one can conclude that eddies do have anisotropic structure: eddies are stretched along the direction of the mean field. Apparently, Fig. 3 suggests that $`k_{}k_{}`$ and, hence, that anisotropy is scale-independent. No theory mentioned in §1 agrees with our result. However, it is not clear from the figure whether or not the true anisotropy is a function of scale: although the contours show a linear relationship between $`k_{}`$ and $`k_{}`$, the relationship doesn’t mean that all individual eddies have the same major axis to minor axis ratio. As explained in next subsection, the $`rotation`$ of eddies by large-scale waves in the magnetic fields can distort the actual scaling relation and lead to the linear relationship shown in the figure.

Fig. 5 shows $`t_{phase}`$ as a function of wavenumber $`k`$. This time scale is defined as the average correlation time $`\mathrm{\Delta }t`$ such that Fourier components at $`𝐤=k`$ have a phase shift of $`60^o`$ with respect to the original phase. We plot the result for zero frequency modes (i.e. $`k_{}`$=0 modes). We see that $`t_{phase}k^1`$ (that is, $`k_{}^1`$). How can we interpret this result? If a turbulent structure, characterized by a wavenumber $`k`$, moves a distance $`l`$ then the phase is shifted by roughly $`kl`$, even if the eddy is unaffected by the motion. This implies that a large scale velocity $`V`$ will change the phase at a rate $`kV`$, so that $`t_{phase}k^1`$. This implies that our calculation of $`t_{phase}`$ is dominated by large scale motions rather than by the local cascade of energy, as long as the nonlinear cascade rate is proportional to $`k`$ to some power less than one, which is generally the case. This is an example of how large-scale fluctuations can complicate attempts to find physically meaningful scaling relations.

### 3.2 Rotation Effect

In the previous subsection, we showed that anisotropy appears to be scale-independent. In this subsection, we will show that this apparent scaling relation is an artifact caused by the Fourier transformation. That is, we will demonstrate that large-scale modes in the magnetic field can distort the actual scaling law at smaller scales when we perform a Fourier transformation. In this manner, a straightforward evaluation of anisotropy in Fourier space is strongly contaminated by the curvature of the large-scale magnetic fields and does not reflect the actual local anisotropy when anisotropy is more pronounced at smaller scales. Consequently, figures 3 and 4 are compatible with any scaling law that predicts that smaller eddies are more elongated.

Fig. 3 implies that eddies have anisotropic shapes: on average, eddies are stretched along the direction of $`𝐁_0`$. However, not all eddies are aligned along the large scale field. The elongation of an eddy is determined by its interaction with the $`local`$ magnetic field, not the background field. Since the large-scale magnetic field lines wander with respect to $`𝐁_0`$, all smaller scale eddies have similar angular distributions around $`𝐁_0`$. We illustrate this effect in Fig. 6.

In this way, we can explain the results of Goldreich and Sridhar (1995) and Matthaeus et al. (1998) simultaneously. Suppose that eddies are oriented according to the local field lines (Fig. 6a). Goldreich and Sridhar’s result implies smaller eddies are relatively more elongated. When we perform a Fourier transformation and measure the ratio of $`k_{}/k_{}`$, what we actually measure is not the ratio of the minor axis to the major axis of individual eddies. Instead, because the direction of the local magnetic field varies according to location and Fourier transformation is none other than a (weighted) averaging process, we actually measure the ratio averaged over all possible orientation of the eddies (Fig. 6b). The Fourier transformation sees that $`L_1`$ (and $`l_1`$) is the minor axis and $`L_2`$ (and $`l_2`$) is the major axis of an eddy. The ratio of $`L_1/L_2`$ ($`=l_1/l_2`$) is determined by the degree of the wandering of the large scale magnetic field lines with respect to $`𝐁_0`$ and, therefore, the ratio is nearly constant for all eddies, regardless of their sizes and shapes. In fact, the ratio will depend on the tangent of the angle between $`𝐁_0`$ and $`𝐁`$. Since $`k_{}1/L_1`$ and $`k_{}1/L_2`$, the measured value of the ratio $`k_{}/k_{}`$ is nearly scale-independent. We expect the angle $`\theta `$ ($`90^o\theta _w`$) between $`𝐁_0`$ and $`𝐁`$ to be

$$\mathrm{tan}\theta =b/B_0.$$
(13)
Therefore, we have

$$\mathrm{sin}\theta =b/B=\mathrm{cos}\theta _w,$$
(14)
which is exactly the scaling relation found by Matthaeus et al.

## 4 Measuring Eddy Shapes

In this section, we analyze the shape of eddies in real space. As explained in the previous section, we assume that elongated eddies are aligned in the direction of the local magnetic fields. If we want to visualize the shape of eddies, we need to first identify the direction of the local magnetic fields. It is important to note that, although we use the term ‘local magnetic fields’ for simplicity, the fields are different from $`𝐁(𝐫)`$. Let us consider an eddy of size $`l`$. The ‘local magnetic field’ of the eddy must act as the ‘large-scale magnetic field’ for the eddy. Therefore, the ‘local magnetic fields’ for eddies of size $`l`$, must be smoothly varying vector fields whose characteristic length of variation is $`>l`$. Note also that another eddy of size $`l^{}`$ ($`l`$) at the same location can have a slightly different direction for the ‘local magnetic field.’ The ‘local magnetic fields’ are functions of position ($`𝐫`$) and eddy size ($`l`$). Then, how can we define the direction of the local magnetic fields? We implement 2 independent numerical algorithms to calculate the direction of the $`local`$ magnetic fields.

In the first method, the local magnetic fields are calculated by

$$𝐁_l=(𝐁(𝐫_1)+𝐁(𝐫_2))/2$$
(15)
and the second order structure functions for $`v`$ and $`b`$ are given by

$`F_2^v(R,z)=<|𝐕(𝐫_1)𝐕(𝐫_2)|^2>,`$ (16)
$`F_2^b(R,z)=<|𝐛(𝐫_1)𝐛(𝐫_2)|^2>,`$ (17)
where $`R=|\widehat{𝐳}\times (𝐫_2𝐫_1)|,`$ $`z=\widehat{𝐳}(𝐫_2𝐫_1)`$ and $`\widehat{𝐳}=𝐁_l/|𝐁_l|.`$ That is, $`R`$ and $`z`$ are coordinates in a cylindrical coordinate system in which the z-axis is parallel to $`𝐁_l`$ (Fig. 7). $`𝐁(𝐫)`$ and $`𝐛(𝐫)`$ are the total and perturbed magnetic fields at a point $`𝐫`$. As usual, brackets denotes a spatial average.

In Fig. 8 we plot the second order structure functions in z-R plane. The horizontal axis (z-axis) is parallel to $`𝐁_l`$. The contours reflect the shapes of the eddies. Eddies are clearly elongated along the local field lines, and smaller eddies are more elongated. The velocity and magnetic fields show different structure at large scales. However, their small scale structure is quite similar.

In Fig. 9, we plot R-intercepts and z-intercepts of the contours. The R-intercept and z-intercept of a given contour can be regarded as a measure of $`\stackrel{~}{k}_{}`$ and $`\stackrel{~}{k}_{}`$ for the corresponding eddy scale444We interpret $`\stackrel{~}{k}_{}`$ as the inverse of the major axis of eddies and $`\stackrel{~}{k}_{}`$ as that of the minor axis. Since the major axis is assumed to be parallel to the local magnetic fields, $`\stackrel{~}{k}_{}`$ is the parallel wavenumber with respect to the local magnetic field direction. On the other hand, in §3.1, $`k_{}`$ is parallel to $`𝐁_0`$. Therefore, $`\stackrel{~}{k}_{}`$ and $`\stackrel{~}{k}_{}`$ in this section have different meanings from $`k_{}`$ and $`k_{}`$ in §3.1.. Fitting the results for velocity fields gives

$$v:\stackrel{~}{k}_{}\{\begin{array}{cc}\stackrel{~}{k}_{}^{0.69},\hfill & (256HB_00.5)\hfill \\ \stackrel{~}{k}_{}^{0.70},\hfill & (256HB_01)\hfill \\ \stackrel{~}{k}_{}^{0.73},\hfill & (256PB_01)\hfill \end{array}$$
(18)
On the other hand, for the magnetic fields we find

$$b:\stackrel{~}{k}_{}\{\begin{array}{cc}\stackrel{~}{k}_{}^{0.64},\hfill & (256HB_00.5)\hfill \\ \stackrel{~}{k}_{}^{0.50},\hfill & (256HB_01)\hfill \\ \stackrel{~}{k}_{}^{0.53},\hfill & (256PB_01)\hfill \end{array}$$
(19)
The velocity fields show good agreement with the scaling relation, $`\stackrel{~}{k}_{}\stackrel{~}{k}_{}^{2/3}`$, proposed by Goldreich and Sridhar (1995). The power indices are insensitive to the strength of $`B_0`$ or the form of viscosity (and diffusivity). However, the magnetic field shows different scaling behavior. When the external field is moderately strong (256H-$`B_00.5`$), the power index is very close to 2/3. On the other hand, for stronger external fields (256H-$`B_01`$ and 256P-$`B_01`$), the power indices are smaller than 2/3. The results are insensitive to the choice of viscosity (and diffusivity). It is not clear whether the existence of a separate scaling law for the magnetic field represents a physical effect not included in Goldreich and Sridhar’s model, e.g. the first signs of small scale intermittency in the magnetic field distribution, or merely the failure of the numerical models used here to fully resolve the inertial range of strong MHD turbulence.

In the second method, we employ a completely different approach. We obtain the local large scale magnetic fields by filtering out the small scale magnetic fields:

$$𝐁_\sigma (𝐫)=\underset{𝐫^{}}{}𝐁(𝐫^{})\varphi (|𝐫𝐫^{}|),$$
(20)
where $`\varphi (r)\mathrm{exp}(r^2/\sigma _r^2)`$ is a gaussian function. To determine the shape of small scale eddies (i.e. eddies smaller than the filter size, $`\sigma _r`$), we consider the following quantities:

$`𝐯_\sigma (𝐫)=𝐕(𝐫)𝐕_\sigma (𝐫),`$ (21)
$`𝐛_\sigma (𝐫)=𝐁(𝐫)𝐁_\sigma (𝐫),`$ (22)
where $`𝐕_\sigma `$ and $`𝐁_\sigma `$ are filtered fields (cf. eq. (20)). The fields $`𝐯_\sigma (𝐫)`$ and $`𝐛_\sigma (𝐫)`$ represent small scale fluctuation of velocity and magnetic fields, the shape of which can be regarded as an adequate approximation of small scale eddies. From these two fields we calculate the structure functions

$`F^v(R,z)=|𝐯_\sigma (𝐫_2)𝐯_\sigma (𝐫_1)|,`$ (23)
$`F^b(R,z)=|𝐛_\sigma (𝐫_2)𝐛_\sigma (𝐫_1)|,`$ (24)
where $`R`$ and $`z`$ are similarly defined as in the first method with $`\widehat{𝐳}𝐁_\sigma (𝐫_1)`$ (Fig. 10).

In Fig. 11, we plot the results of the second method. We can clearly observe the flattening effect: When the filter size is large, eddies are less anisotropic. When filter size is small, eddies show highly anisotropic structure. In the figure, we plot only the results for magnetic fields. Velocity fields show similar trends. The thick contours represent $`F^b(R,z)=<|𝐛_\sigma |^2>^{1/2}`$. The values next to the thick contours are $`<|𝐛_\sigma |^2>^{1/2}`$.

This second method is useful for visualizing small scale eddies, but may not be as useful for quantitative analysis. The difficulty comes from the fact that there is no well defined eddy scale associated with filtered fields $`𝐯_\sigma (𝐫)`$ and $`𝐛_\sigma (𝐫)`$.

In this paper, we will not pursue quantitative analysis using the method 2. Instead, we just wish to point out that both methods describe the same scaling law: smaller eddies are relatively more stretched along the local magnetic field lines. In particular, the results from the first method are consistent with the relation $`\stackrel{~}{k}_{}\stackrel{~}{k}_{}^{2/3}`$ proposed by Goldreich and Sridhar (1995).

## 5 Discussion and Conclusions

Here we rederive the scaling law $`\stackrel{~}{k}_{}\stackrel{~}{k}_{}^{2/3}`$ in the framework of 3-wave interactions. Except for the use of the uncertainty principle, the work in this section is independent of Goldreich and Sridhar’s derivation. As noted by Goldreich and Sridhar (1995), 3-wave interactions are an adequate proxy for wave-wave interactions of all orders in a strong MHD turbulence. As long as we assume the locality of interactions, it is pointless to distinguish $`k_{}`$ and $`k_{}`$ from $`\stackrel{~}{k}_{}`$ and $`\stackrel{~}{k}_{}`$. Hence, for simplicity, we use $`k_{}`$ and $`k_{}`$ instead of $`\stackrel{~}{k}_{}`$ and $`\stackrel{~}{k}_{}`$ during the derivation.

Suppose the 3-dimensional energy spectrum is given by

$$E_3(k_{},k_{})|\widehat{𝐕}(𝐤)|^2k_{}^{2\alpha }f,$$
(25)
where $`\widehat{𝐕}(𝐤)`$ is the amplitude of the mode whose wavevector is $`𝐤`$ and $`f(u)`$ is a positive, symmetric function of $`u`$ (cf. equation (7) of Goldreich and Sridhar (1995)) which describes the power distribution as a function of eddy shape. If the $`width`$ (or, $`thickness`$) of the energy spectrum in the direction of $`k_{}`$ is $`k_{}^\beta `$, then we can write

$$E_3(k_{},k_{})k_{}^{2\alpha }f(k_{}/k_{}^\beta ).$$
(26)
If the $`width`$ is caused by the uncertainty principle ($`\mathrm{\Delta }t\mathrm{\Delta }\omega 1`$ with $`\mathrm{\Delta }tt_{cas}(l)`$ and $`\mathrm{\Delta }\omega k_{}`$), then

$$t_{cas}(l)k_{}^\beta .$$
(27)
Suppose the energy cascade rate $`ϵv_l^2/t_{cas}(l)`$, where $`l=2\pi /k_{}`$, is scale-independent. Because $`v_l^2k_{}^{2\alpha }k_{}^2k_{}^\beta `$ ($`k_{}E(k_{})`$, $`E(k_{})=`$ 1-dimensional spectrum) and $`t_{cas}(l)k_{}^\beta `$, we have

$$k_{}^{2\alpha +2+2\beta }=const.$$
(28)
Therefore,

$$1\alpha +\beta =0.$$
(29)
Now, let’s pick up a mode at a wavevector $`𝐩`$ and consider nonlinear interactions with other wave modes. First, the strength of the interaction with another mode at $`𝐪`$ is $`p|\widehat{𝐕}(𝐩)||\widehat{𝐕}(𝐪)|`$. (This comes from the $`(\times 𝐕)\times 𝐕`$ term in equation (9)) Hereafter we assume $`p|𝐩|p_{}`$. Second, the number of interactions is $`p^2p^\beta `$. This is the number of modes $`near`$ $`𝐩`$. Here we use locality of 3-wave interactions. If the interactions are random, the net change of amplitude per unit time will be the strength of the interaction times the square root of the number of interactions, or

$$|\mathrm{\Delta }\widehat{𝐕}(𝐩)|p|\widehat{𝐕}(𝐩)||\widehat{𝐕}(𝐪)|(p^2p^\beta )^{1/2}.$$
(30)
Therefore, we have

$$t_{cas}|\widehat{𝐕}(𝐩)|/|\mathrm{\Delta }\widehat{𝐕}(𝐩)|p^2p^\alpha p^{\beta /2},$$
(31)
where we assumed $`pq`$. Equating this with $`t_{cas}p^\beta `$, we can write

$$\alpha 2=\beta /2.$$
(32)
From equations (29) and (32), we have

$$\alpha =5/3,\beta =2/3,$$
(33)
which is just the result of Goldreich and Sridhar (1995):

$$k_{}k_{}^{2/3}.$$
(34)
As a consequence, the 3-D energy spectrum becomes

$$E_3(k_{},k_{})k_{}^{10/3}f(k_{}/k_{}^{2/3})$$
(35)
and the corresponding 1-D spectrum is given by

$$E(k)k^{5/3}.$$
(36)
In summary, we have shown that the anisotropy of Alfvénic turbulence depends on the spatial scales of eddies. In particular, our results confirm the claim by Goldreich and Sridhar (1995, 1997) that smaller eddies are relatively more elongated along the direction of the local magnetic field lines than larger ones. Quantitative measurements of the anisotropy using the velocity fields show good agreement with their proposed scaling law, $`k_{}k_{}^{2/3}`$ as long we interpret these wavenumbers as referring to the local magnetic field direction. However, when the external magnetic field is very strong, magnetic fields scale somewhat differently, showing a slightly more rapid increase in anisotropy at smaller scales.

It is important to note that the correct scaling laws depend on comparing the eddy shape to the local magnetic field direction.

As a final note, we wish to stress that our results are not in agreement with the IK theory. The IK theory is based on the assumption of isotropy in wavenumber space, which may be true when the external magnetic field is very weak or zero. However, even in these cases, the turbulence is globally isotropic, but locally very anisotropic. In this paper, we showed that eddies do show anisotropy and that the anisotropy is scale-dependent when there is a strong large scale field (which should apply to very small scales within any MHD turbulence cascade). On the other hand, our results are consistent with Goldreich and Sridhar’s theory of strong MHD turbulence. More precisely, if we consider the ratio of hydrodynamic to Alfvénic rates, that is $`(kv_k/\stackrel{~}{k}_{}V_A)`$, we find from equations (18) and (19) that

$$\frac{kv_k}{\stackrel{~}{k}_{}V_A}\stackrel{~}{k}_{}^{0.30.5}v_k,$$
(37)
where $`k\stackrel{~}{k}_{}`$. From Fig. 2 we see that for the inertial range this implies a ratio which is either constant or increasing with wavenumber. A constant ratio is predicted by Goldreich and Sridhar’s model. The IK model predicts a slow decline.

This work was partially supported by National Computational Science Alliance under CTS980010N and utilized the NCSA SGI/CRAY Origin2000.


format: text

text:
# Realistic Calculation of the ℎ𝑒𝑝 Astrophysical Factor

## Abstract

The astrophysical factor for the proton weak capture on 3He is calculated with correlated-hyperspherical-harmonics bound and continuum wave functions corresponding to a realistic Hamiltonian consisting of the Argonne $`v_{18}`$ two-nucleon and Urbana-IX three-nucleon interactions. The nuclear weak charge and current operators have vector and axial-vector components, that include one- and many-body terms. All possible multipole transitions connecting any of the $`p^3`$He S- and P-wave channels to the 4He bound state are considered. The $`S`$-factor at a $`p^3`$He center-of-mass energy of $`10`$ keV, close to the Gamow-peak energy, is predicted to be $`10.1\times 10^{20}`$ keV b, a factor of five larger than the standard-solar-model value. The P-wave transitions are found to be important, contributing about 40 % of the calculated $`S`$-factor.

preprint: JLAB-THY-00-09

Recently, there has been a revival of interest in the reaction 3He($`p`$,$`e^+\nu _e`$)4He . This interest has been spurred by the Super-Kamiokande collaboration measurements of the energy spectrum of electrons recoiling from scattering with solar neutrinos . At energies larger than 14 MeV more recoil electrons have been observed than expected on the basis of standard-solar-model (SSM) predictions . The $`hep`$ process, as the proton weak capture on 3He is known, is the only source of solar neutrinos with energies larger than 15 MeV–their end-point energy is about 19 MeV. This fact has naturally led to questions about the reliability of the currently accepted SSM value for the astrophysical factor at zero energy, $`2.3\times 10^{20}`$ keV b . In particular, Bahcall and Krastev  have shown that a large enhancement, by a factor in the range 20–30, of the SSM $`S`$-factor value given above would essentially fit the observed excess  of recoiling electrons.

The theoretical description of the $`hep`$ process, as well as that of the neutron and proton radiative captures on deuteron and 3He, constitute a challenging problem from the standpoint of nuclear few-body theory. Its difficulty can be appreciated by comparing the measured values for the cross section of thermal neutron radiative capture on 1H, 2H, and 3He. Their respective cross sections are: $`334.2\pm 0.5`$ mb , $`0.508\pm 0.015`$ mb , and $`0.055\pm 0.003`$ mb . Thus, in going from $`A`$=2 to 4 the cross section has dropped by almost four orders of magnitude. These processes are induced by magnetic-dipole transitions between the initial two-cluster state in relative S-wave and the final bound state. In fact, the inhibition of the $`A`$=3 and 4 captures has been understood for a long time : the 3H and 4He states are approximate eigenstates of the magnetic dipole operator $`𝝁`$, and consequently matrix elements of $`\mu _z`$ between $`nd`$ ($`n^3`$He) and 3H (4He) vanish (approximately) due to orthogonality. This orthogonality argument fails in the case of the deuteron, since then $`\mu _z`$ can connect the large S-wave component of the deuteron to an isospin $`T`$=1 1S0 $`np`$ state.

This quasi-orthogonality, while again invalid in the case of the proton weak capture on protons , is also responsible for inhibiting the $`hep`$ process. Both these reactions are induced by the Gamow-Teller operator, which differs from the (leading) isovector spin part of the magnetic dipole operator essentially by an isospin rotation. As a result, the $`hep`$ weak capture and $`nd`$, $`pd`$ and $`n^3`$He radiative captures are extremely sensitive to: i) the D-state admixtures generated by tensor interactions, and ii) many-body terms in the electroweak current operator. For example, many-body current contributions provide, respectively, 50 % and over 90 % of the calculated $`pd`$  and $`n^3`$He  cross sections at very low energies.

In this respect, the $`hep`$ weak capture is a particularly delicate reaction, for two additional reasons: firstly and most importantly, the one- and many-body current contributions are comparable in magnitude, but of opposite sign ; secondly, many-body axial currents, specifically those arising from excitation of $`\mathrm{\Delta }`$ isobars which give the dominant contribution, are model dependent . This destructive interference between one- and many-body currents also occurs in the $`n^3`$He (“$`hen`$”) radiative capture , with the difference that there the leading components of the many-body currents are model independent, and give a much larger contribution than that associated with the one-body current.

The cancellation in the $`hep`$ process between the one- and two-body matrix elements has the effect of enhancing the importance of P-wave capture channels. Indeed, one of the results of the present work is that these channels give about 40 % of the $`S`$-factor calculated value. That the $`hep`$ process could proceed as easily through P- as S-wave capture was not not sufficiently appreciated  in all earlier studies of this reaction we are aware of, with the exception of Ref. , in which Horowitz suggested, on the basis of a very simple one-body reaction model, that the 3P0 channel may be important.

Most of the earlier studies  had attempted to relate the matrix element of the axial current occurring in the $`hep`$ capture to that of the electromagnetic current in the $`hen`$ capture, exploiting (approximate) isospin symmetry. This approach led, however, to $`S`$-factor values ranging from $`3.7`$ to $`57`$, in units of $`10^{20}`$ keV b. In an attempt to reduce the uncertainties in the predicted values for both the radiative and weak capture rates, ab initio microscopic calculations of these reactions were performed in the early nineties , using variational wave functions corresponding to a realistic Hamiltonian, and a nuclear electroweak current consisting of one- and many-body components. These studies showed that inferring the $`hep`$ $`S`$-factor from the measured $`hen`$ cross section can be misleading, because of different initial-state interactions in the $`n^3`$He and $`p^3`$He channels, and because of the large contributions associated with the two-body components of the electroweak current operator, and their destructive interference with the one-body current contributions.

The significant progress made in the last few years in the modeling of two- and three-nucleon interactions and the nuclear weak current, and the description of the bound and continuum four-nucleon wave functions, have prompted us to re-examine the $`hep`$ reaction. In the present work we briefly summarize the salient points in the calculation, and report our results for the $`S`$-factor in the energy range 0–10 keV. An exhaustive account of this study , however, will be published elsewhere.

The cross section for the 3He($`p`$,$`e^+\nu _e`$)4He reaction at a c.m. energy $`E`$ is written as

$`\sigma (E)={\displaystyle }`$ $`2\pi \delta \left(\mathrm{\Delta }m+E{\displaystyle \frac{q^2}{2m_4}}E_eE_\nu \right){\displaystyle \frac{1}{v_{\mathrm{rel}}}}`$ (2)
$`\times {\displaystyle \frac{1}{4}}{\displaystyle \underset{s_es_\nu }{}}{\displaystyle \underset{s_1s_3}{}}|f|H_W|i|^2{\displaystyle \frac{d𝐩_e}{(2\pi )^3}}{\displaystyle \frac{d𝐩_\nu }{(2\pi )^3}},`$
where $`\mathrm{\Delta }m=m+m_3m_4`$ = 19.8 MeV ($`m`$, $`m_3`$, and $`m_4`$ are the proton, 3He, and 4He rest masses, respectively), $`v_{\mathrm{rel}}`$ is the $`p^3`$He relative velocity, and the transition amplitude is given by

$$f|H_W|i=\frac{G_V}{\sqrt{2}}l^\sigma 𝐪;^4\mathrm{He}|j_\sigma ^{}(𝐪)|𝐩;p^3\mathrm{He}.$$
(3)
Here $`G_V`$ is the Fermi constant, $`𝐪=𝐩_e+𝐩_\nu `$, $`|𝐩;p^3\mathrm{He}`$ and $`|𝐪;^4\mathrm{He}`$ represent, respectively, the $`p^3`$He scattering state with relative momentum $`𝐩`$ and 4He bound state recoiling with momentum $`𝐪`$, $`l_\sigma `$ is the leptonic weak current, $`l_\sigma =\overline{u}_\nu \gamma _\sigma (1\gamma _5)v_e`$ (the lepton spinors are normalized as $`v_e^{}v_e=u_\nu ^{}u_\nu =1`$), and $`j^\sigma (𝐪)`$ is the nuclear weak current, $`j^\sigma (𝐪)=(\rho (𝐪),𝐣(𝐪))`$. The dependence of the amplitude upon the spin projections of the leptons, proton and 3He has been omitted for ease of presentation. Since the energies of interest are of the order of 10 keV or less–the Gamow-peak energy is 10.7 keV for the $`hep`$ reaction–it is convenient to expand the $`p^3`$He scattering state into partial waves, and perform a multipole decomposition of the nuclear weak charge, $`\rho (𝐪)`$, and current, $`𝐣(𝐪)`$, operators. Standard manipulations lead to

$$\frac{1}{4}\underset{s_es_\nu }{}\underset{s_1s_3}{}|f|H_W|i|^2=(2\pi )^2G_V^2L_{\sigma \tau }N^{\sigma \tau },$$
(4)
where the lepton tensor $`L^{\sigma \tau }`$ is written in terms of electron and neutrino four-velocities as $`L^{\sigma \tau }=\mathrm{v}_e^\sigma \mathrm{v}_\nu ^\tau +\mathrm{v}_\nu ^\sigma \mathrm{v}_e^\tau g^{\sigma \tau }\mathrm{v}_e\mathrm{v}_\nu +\mathrm{i}ϵ^{\sigma \alpha \tau \beta }\mathrm{v}_{e,\alpha }\mathrm{v}_{\nu ,\beta }`$, while the nuclear tensor is defined as

$$N^{\sigma \tau }\underset{s_1s_3}{}W^\sigma (𝐪;s_1s_3)W^\tau (𝐪;s_1s_3),$$
(5)
with

$`W^{\sigma =0,3}(𝐪;s_1s_3)`$ $`=`$ $`{\displaystyle \underset{LSJ}{}}X_0^{LSJ}(\widehat{𝐪};s_1s_3)T_J^{LSJ}(q),`$ (6)
$`W^{\sigma =\lambda }(𝐪;s_1s_3)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle \underset{LSJ}{}}X_\lambda ^{LSJ}(\widehat{𝐪};s_1s_3)`$ (8)
$`\left[\lambda M_J^{LSJ}(q)+E_J^{LSJ}(q)\right],`$
where $`\lambda =\pm 1`$ denote spherical components. The functions $`X_{\lambda =0,\pm 1}`$ depend upon the direction $`\widehat{𝐪}`$, and the proton and 3He spin projections $`s_1`$ and $`s_3`$  (note that the quantization axis for the hadronic states is taken along $`\widehat{𝐩}`$, the direction of the $`p^3`$He relative momentum), while $`T_J^{LSJ}`$=$`C_J^{LSJ}`$ or $`L_J^{LSJ}`$ for $`\sigma `$=0 or 3. The quantities $`C_J^{LSJ}`$, $`L_J^{LSJ}`$, $`E_J^{LSJ}`$ and $`M_J^{LSJ}`$ are the reduced matrix elements (RMEs) of the Coulomb $`(C_\mathrm{}\mathrm{}_z)`$, longitudinal $`(L_\mathrm{}\mathrm{}_z)`$, transverse electric $`(E_\mathrm{}\mathrm{}_z)`$, and transverse magnetic $`(M_\mathrm{}\mathrm{}_z)`$ multipole operators between the initial $`p^3`$He state with orbital angular momentum $`L`$, channel spin $`S`$ ($`S`$=0,1), and total angular momentum $`J`$, and final 4He state. The present study includes S- and P-wave capture channels, i.e. the 1S0, 3S1, 3P0, 1P1, 3P1, and 3P2 states, and retains all contributing multipoles connecting these states to the $`J^\pi `$=0+ ground state of 4He.

The bound- and scattering-state wave functions are obtained variationally with the correlated-hyperspherical-harmonics (CHH) method, developed in Refs. . The nuclear Hamiltonian consists of the Argonne $`v_{18}`$ two-nucleon  and Urbana-IX three-nucleon  interactions. This realistic Hamiltonian, denoted as AV18/UIX, reproduces the experimental binding energies and charge radii of the trinucleons and 4He in “exact”Green’s function Monte Carlo (GFMC) calculations . The binding energy of 4He calculated with the CHH method  is within 1 % of that obtained with the GFMC method. The accuracy of the CHH method to calculate scattering states has been successfully verified in the case of three-nucleon systems, by comparing results for a variety of $`Nd`$ scattering observables obtained by a number of groups using different techniques . Studies along similar lines  to assess the accuracy of the CHH solutions for the four-nucleon continuum have already begun.

The CHH predictions  for the $`n^3`$H total elastic cross section and coherent scattering length have been found to be in excellent agreement with the corresponding experimental values. The $`n^3`$H cross section is known over a rather wide energy range, and its extrapolation to zero energy is not problematic . The situation is different for the $`p^3`$He channel, for which the singlet and triplet scattering lengths $`a_\mathrm{s}`$ and $`a_\mathrm{t}`$ have been determined from effective range extrapolations of data taken above 1 MeV, and are therefore somewhat uncertain, $`a_\mathrm{s}=(10.8\pm 2.6)`$ fm  and $`a_\mathrm{t}=(8.1\pm 0.5)`$ fm  or $`(10.2\pm 1.5)`$ fm . Nevertheless, the CHH results are close to the experimental values above: the AV18/UIX Hamiltonian predicts  $`a_\mathrm{s}=11.5`$ fm and $`a_\mathrm{t}=9.13`$ fm. At low energies (below 4 MeV) $`p^3`$He elastic scattering proceeds mostly through S- and P-wave channels, and the CHH predictions, based on the AV18/UIX model, for the differential cross section  are in good agreement with the experimental data.

The nuclear weak current has vector and axial-vector parts, with corresponding one- and many-body components. The one-body components have the standard expressions obtained from a non-relativistic reduction of the covariant single-nucleon vector and axial-vector currents, including terms proportional to $`1/m^2`$. The two-body weak vector currents are constructed from the isovector two-body electromagnetic currents in accordance with the conserved-vector-current hypothesis, and consist  of “model-independent”and “model-dependent”terms. The model-independent terms are obtained from the nucleon-nucleon interaction, and by construction satisfy current conservation with it. The leading two-body weak vector current is the “$`\pi `$-like”operator, obtained from the isospin-dependent spin-spin and tensor nucleon-nucleon interactions. The latter also generate an isovector “$`\rho `$-like”current, while additional isovector two-body currents arise from the isospin-independent and isospin-dependent central and momentum-dependent interactions. These currents are short-ranged, and numerically far less important than the $`\pi `$-like current. With the exception of the $`\rho `$-like current, they have been neglected in the present work. The model-dependent currents are purely transverse, and therefore cannot be directly linked to the underlying two-nucleon interaction. The present calculation includes the currents associated with excitation of $`\mathrm{\Delta }`$ isobars which, however, are found to give a rather small contribution in weak-vector transitions, as compared to that due to the $`\pi `$-like current. The $`\pi `$-like and $`\rho `$-like contributions to the weak vector charge operator  have also been retained in the present study.

The leading many-body terms in the axial current due to $`\mathrm{\Delta }`$-isobar excitation are treated non-perturbatively in the transition-correlation-operator (TCO) scheme, originally developed in Ref.  and further extended in Ref. . In the TCO scheme–essentially, a scaled-down approach to a full $`N`$$`+`$$`\mathrm{\Delta }`$ coupled-channel treatment–the $`\mathrm{\Delta }`$ degrees of freedom are explicitly included in the nuclear wave functions. The axial charge operator includes, in addition to $`\mathrm{\Delta }`$-excitation terms (which, however, are found to be unimportant ), the long-range pion-exchange term , required by low-energy theorems and the partially-conserved-axial-current relation, as well as the (expected) leading short-range terms constructed from the central and spin-orbit components of the nucleon-nucleon interaction, following a prescription due to Riska and collaborators .

The largest model dependence is in the weak axial current. To minimize it, the poorly known $`N`$$`\mathrm{\Delta }`$ transition axial coupling constant $`g_A^{}`$ has been adjusted to reproduce the experimental value of the Gamow-Teller matrix element in tritium $`\beta `$-decay . While this procedure is model dependent, its actual model dependence is in fact very weak, as has been shown in Ref. .

The calculation proceeds in two steps : first, the matrix elements of $`\rho (𝐪)`$ and $`𝐣(𝐪)`$ between the initial $`p^3`$He $`LSJJ_z`$ states and final 4He are calculated with Monte Carlo integration techniques; second, the contributing RMEs are extracted from these matrix elements, and the cross section is calculated by performing the integrations over the electron and neutrino momenta in Eq. (2) numerically, using Gauss points.

The results for the $`S`$-factor, defined as $`S(E)=E\sigma (E)\mathrm{exp}(4\pi \alpha /v_{\mathrm{rel}})`$ ($`\alpha `$ is the fine structure constant), at $`p^3`$He c.m. energies of $`0`$, $`5`$, and $`10`$ keV are reported in Table I. In the table, the column labelled S includes both the 1S0 and 3S1 channel contributions, although the former are at the level of a few parts in $`10^3`$. The energy dependence is rather weak: the value at $`10`$ keV is only about 4 % larger than that at $`0`$ keV. The P-wave capture states are found to be important, contributing about 40 % of the calculated $`S`$-factor. However, the contributions from D-wave channels are expected to be very small. We have verified explicitly that they are indeed small in 3D1 capture. The many-body axial currents associated with $`\mathrm{\Delta }`$ excitation play a crucial role in the (dominant) 3S1 capture, where they reduce the $`S`$-factor by more than a factor of four. Thus the destructive interference between the one- and many-body current contributions, first obtained in Ref. , is confirmed in the present study. The (suppressed) one-body contribution comes mostly from transitions involving the D-state components of the 3He and 4He wave functions, while the many-body contributions are predominantly due to transitions connecting the S-state in 3He to the D-state in 4He, or viceversa.

It is important to stress the differences between the present and all previous studies. Apart from ignoring, or at least underestimating, the contribution due to P-waves, the latter only considered the long-wavelength form of the weak multipole operators, namely, their $`q`$=$`0`$ limit. In 3P0 capture, for example, only the $`C_0`$-multipole, associated with the weak axial charge, survives in this limit, and the corresponding $`S`$-factor is calculated to be $`2.2\times 10^{20}`$ keV b, including two-body contributions. However, when the transition induced by the longitudinal component of the axial current (via the $`L_0`$-multipole, which vanishes at $`q`$=$`0`$) is also taken into account, the $`S`$-factor becomes $`0.82\times 10^{20}`$ keV b, because of a destructive interference between the $`C_0`$ and $`L_0`$ RMEs. Thus use of the long-wavelength approximation in the calculation of the $`hep`$ $`S`$-factor leads to inaccurate results.

Finally, besides the differences listed above, the present calculation also improves that of Ref.  in a number of other important respects: firstly, it uses accurate CHH wave functions, corresponding to the last generation of realistic interactions; secondly, the model for the nuclear weak current has been extended to include the axial charge as well as the vector charge and current operators. Thirdly, the one-body operators now take into account the $`1/m^2`$ relativistic corrections, which were previously neglected. In 3S1 capture, for example, these terms increase by 25 % the dominant (but suppressed) $`L_1`$ and $`E_1`$ RMEs calculated with the (lowest order) Gamow-Teller operator. These improvements in the treatment of the one-body axial current indirectly affect also the contributions of the $`\mathrm{\Delta }`$-excitation currents , because of the procedure used to determine the coupling constant $`g_A^{}`$.

To conclude, we have carried out a realistic calculation of the $`hep`$ reaction, predicting a value for the $`S`$-factor five times larger than that used in the SSM. This enhancement, while very significant, is far smaller than that required by fits to the Super-Kamiokande data. Although the present result is inherently model dependent, it is unlikely the model dependence be so large to accomodate a drastic increase in the prediction obtained here.

The authors wish to thank J.F. Beacom, J. Carlson, V.R. Pandharipande, D.O. Riska, P. Vogel, and R.B. Wiringa for useful discussions. The support of the U.S. Department of Energy under contract number DE-AC05-84ER40150 is gratefully acknowledged by L.E.M. and R.S.


text:
# 1 Introduction

## 1 Introduction

Stable quasiperiodic (QP) planforms have been first discovered as equilibrium patterns in metallic alloys . The simplest and most important types of the QP planforms are the ten-, eight-, and twelvefold ones (note an argument according to which only these types of the two-dimensional quasicrystals may occur in physical systems ). A $`2n`$-fold pattern is a superposition of $`n4`$ spatial harmonics based on a set of equal-length wave vectors $`𝐤_j`$ with the equal angles $`\pi /n`$ between $`𝐤_j`$ and $`𝐤_{j+1}`$ (more general aperiodic patterns with unequal wave vector lengths and/or unequal angles between the vectors may also exist, but they have probably never been considered, except for the case of unequilateral hexagons with $`n=3`$ ).

Then, it was predicted that the $`2n`$-fold structures of different particular types may exist as dynamical planforms of a stationary form in nonequilibrium systems, such as thermal convection and the like . Later, the predicted pattern (a twelvefold one, in particular) was indeed experimentally observed, first, as a stable dynamical nonequilibrium planform in the Faraday ripples  (a large-aspect-ratio liquid layer subject to high-frequency shaking in the vertical direction), and recently in an optical cell filled with Na vapors . Nevertheless, experimental generation of dynamical QP patterns is far from being straightforward; in particular, in the work  it was necessary to shake the liquid layer with a two-frequency quasiperiodic force at a specially selected ratio between the two frequencies. These experimental difficulties are related to the general fact that, in terms of the corresponding coupled Ginzburg-Landau (GL) equations for amplitudes of the spatial modes, a superposition of which gives rise to the pattern, a formal solution for the QP planforms always exists, but it may be stable only in a relatively narrow parametric region .

An objective of this work is to put forward a much easier possibility of generating eight- and twelvefold QP planforms of a finite but large size (a stripe) in “normal” systems, where direct generation of the quasicrystallic patterns does not seem feasible. To this end, we notice that the set of the wave vectors on which the eight- or twelvefold QP is based may be regarded, in an obvious way, as a superposition of two half-sets of $`2`$ or $`3`$ vectors, see Fig. 1 below. Each half-set, in turn, may give rise to a usual periodic pattern consisting of square cells (SC) or hexagonal cells, respectively. Thus, a natural way to generate a stripe filled with the QP pattern in a generic system (e.g., a convective layer), in which a stable QP pattern is not available, but stable periodic SC and/or hexagonal patterns do exist, is to produce it as a transient layer (“domain wall” ) between two large domains filled with the periodic cells, the (half-) sets of the two or three wave vectors in the two domains being oriented under the angle, respectively, $`45`$ or $`30`$ degrees relative to each other. The pattern in the transient layer will then be a superposition of the two periodic patterns, having a full set of the four or six wave vectors, respectively. The so generated stripe will feature a QP pattern, provided that its width is essentially larger than the wavelength (a size of the elementary cell in the corresponding periodic pattern). It is known that the latter condition can be achieved, under special but not unrealistic conditions, for domain walls in the nonequilibrium systems .

In fact, it may be simpler to generate, following this way, a twelvefold QP pattern between two hexagonal domains, as it is usually much easier to find a stable hexagonal structure than a stable SC one . However, the theoretical analysis of the transient layer between two SC domains is much simpler, therefore in this work we concentrate on the latter case. In any case, it should be stressed that, although stable square-cell patterns are rarer than their counterparts in the form of hexagons or rolls, examples of stable square cellular planforms are known in gas flames , thermal convection , and in some optical systems . The interest to the SC planforms has been recently revived by the discovery of this pattern in a double-layer Marangoni convection (see the e-print  and references therein).

The rest of the paper is organized as follows. In section 2, we formulate the model. Some analytical results, which predict, in particular, that the transient layer between the two SC patterns has a complex structure, consisting of three layers, two broad ones and a narrow sublayer sandwiched between them, are obtained in section 3. In section 4, we display results of a direct numerical solution of the stationary real GL equations, which comply with the analytical predictions. The paper is concluded by section 5.

## 2 The Model

In this work, we assume a simplest configuration, shown in Fig. 1, that is going to give rise to the QP transient layer. The layer is parallel to the $`y`$ axis, and the (half) sets of the two wave vectors are chosen so that in the left domain the vectors $`𝐤_{1,2}`$ are parallel to the $`x`$ and $`y`$ axes, while in the right domain both vectors $`𝐤_{3,4}`$ have the angle $`45`$ degrees relative to the axes. Accordingly, the complex wave field describing the spatial distribution of physical variables is assumed to be

$$u(𝐫,t)=\underset{j=1}{\overset{4}{}}B_j(𝐫,t)\mathrm{exp}\left(i𝐤_j𝐫\right),$$
(1)
where $`𝐫`$ is the two-dimensional coordinate, and $`B_j(𝐫,t)`$ are slowly varying amplitude functions.

In an experiments, the chosen configuration can be created, at least in a part of the system, by means of specially selected boundary conditions, which, in turn, are imposed by the sidewalls of the experimental cell. Although we do not analyze the sidewall boundary conditions in this work, it is obvious that, to support the configuration that we consider, one will need to have a large-aspect-ratio cell with two sidewalls forming an angle $`45`$ degrees, which is quite possible. More general configurations, with different orientations of the wave vectors relative to the layer between the two domains, can be considered similarly to what is done below, but their technical treatment will be more cumbersome.

A usual approach to the description of spatially nonuniform patterns of the domain-wall type is based on a system of coupled real Ginzburg-Landau (GL) equations for the slowly varying amplitudes, assuming that they do not depend on the coordinate $`y`$ running along the stripe (wall), which is a reasonable approximation if the stripe is long enough . Then, as it was shown in the mentioned works, the effective diffusion coefficient (the one in front of the term $`(B_j)_{xx}`$) in each GL equation is

$$D_j=\left(k_j^{(x)}\right)^2,$$
(2)
where $`k_j^{(x)}`$ is the $`x`$-component of the vector. For the configuration shown in Fig. 1, the system of the real GL equations can be easily cast into a final form

$`A_t`$ $`=`$ $`(1/2)A_{xx}+AA^32\left(g_2A^2+g_1B_1^2+g_1B_2^2\right)A,`$ (3)
$`\left(B_1\right)_t`$ $`=`$ $`\left(B_1\right)_{xx}+B_1B_1^32\left(2g_1A^2+g_2B_2^2\right)B_1,`$ (4)
$`\left(B_2\right)_t`$ $`=`$ $`B_2B_2^32\left(2g_1A^2+g_2B_1^2\right)B_2,`$ (5)
where $`B_{1,2}`$ are the amplitudes corresponding to the wave vectors $`𝐤_{1,2}`$ in the left domain (Fig. 1), and, using the obvious symmetry of the configuration shown in Fig. 1, we have set $`B_3=B_4A`$.

Note that the diffusion coefficient in Eq. (5) is zero due to Eq. (2). Generally speaking, in this case one should take into account a higher-order derivative term $`\left(B_2\right)_{xxxx}`$ . However, this is not necessary while the diffusion terms do not vanish in the two other equations (see details below).

In Eqs. (3) through (5), the linear gain coefficient and the coefficient of the nonlinear self-interaction of the spatial mode are normalized to be $`1`$, $`g_1`$ and $`g_2`$ being the coefficients of the nonlinear interaction between the modes with the angles $`\alpha =45`$ and $`90`$ degrees between their carrier wave vectors. Normally, the nonlinear interaction coefficient decreases with the increase of $`\alpha `$, so that

$$g_2<g_1<1.$$
(6)
The necessary and sufficient stability conditions for the SC pattern are well known ,

$$g_2<\frac{1}{2},g_1>\frac{1}{4}\left(1+2g_2\right),$$
(7)
while the conditions providing for stability of the eightfold QP planform are

$$g_2<\frac{1}{2},g_1<\frac{1}{4}\left(1+2g_2\right).$$
(8)
An obvious feature of the two sets of the stability conditions is their incompatibility, i.e., SC and QP can never be stable simultaneously.

In fact, the inequality $`g_2<1/2`$ in the set of the conditions (7) is a cause for the relative rarity of stable SC planforms, as, despite the general property (6), the actual dependence $`g(\alpha )`$ is usually weak, so that $`g_2`$ is not essentially smaller than $`1`$. Nevertheless, the full set of the SC stability conditions (7) can be satisfied in the above-mentioned physical systems.

A consequence of the conditions (8) necessary for the stability of the QP pattern is $`g_1<1/2`$, which, with regard to Eq. (6), makes the stability of the QP pattern still less feasible than that of the SC one. An objective of this work is to propose a way to produce a transient QP pattern between two stable SC domains with different orientations (Fig. 1) in the case when the uniform QP planform is unstable. Note that a similar approach is known as a way to generate of a broad stripe of a SC pattern between two domains of orthogonally oriented rolls in the case when the uniform SC pattern is unstable, while the rolls are stable .

In the case when the SC stability conditions (7) are met, it is quite reasonable to assume that the coefficients $`g_{1,2}`$ are close to $`1/2`$ (because it is physically implausible to have $`g_{1,2}`$ much smaller than $`1/2`$), which suggests to present them as

$$g_{1,2}\frac{1}{2}\left(1\mu _{1,2}\right),|\mu _{1,2}|1.$$
(9)
As it will be seen below, the smallness of $`\mu _{1,2}`$ naturally provides for the strip of the QP pattern to be broad, which is exactly the condition justifying the consideration of the transient-layer patterns. Note that $`\mu _2`$ must be positive according to Eq. (7), while $`\mu _1`$ may formally have either sign. However, it will be shown below that the necessary solution does not exist if $`\mu _1<0`$.

Below, we will also use a parameter

$$m2\mu _1/\mu _2,$$
(10)
which is, generally, $`1`$. In terms of $`m`$, the second stability condition (7) takes a very simple form, $`m<1`$.

Eqs. (3), (4), and (5) must be supplemented by boundary conditions (b.c.) to guarantee that, at $`x\pm \mathrm{}`$ (recall the system is formally assumed to be infinitely large), the pattern considered asymptotically coincides with either of the two SC planforms composed of the modes $`B_1`$ and $`B_2`$ or $`B_3=B_4A`$. Obviously, this implies

$`\underset{x+\mathrm{}}{lim}A(x)`$ $`=`$ $`1/\sqrt{2\mu _2}A_{lim},\underset{x+\mathrm{}}{lim}B_{1,2}(x)=0,`$ (11)
$`\underset{x\mathrm{}}{lim}A(x)`$ $`=`$ $`0,\underset{x\mathrm{}}{lim}B_{1,2}(x)=A_{lim}.`$ (12)
Formally, this set may seem overdetermined, as we add six b.c. to the fourth-order system of Eqs. (3) and (4) (Eq. (5) contains no $`x`$-derivatives). However, it is easy to check that the seemingly superfluous b.c. for $`B_2`$ are nothing else but direct corollaries of the four legitimate b.c. for $`A`$ and $`B_1`$.

## 3 Analytical Results

A stationary version ($`/t=0`$) of Eqs. (3), (4) and (5) can be essentially simplified, as in this case Eq. (5) becomes just an algebraic relation, that has two solutions: $`B_2=0`$, or

$$B_2^2=1\left(2m\mu _2\right)A^2\left(1\mu _2\right)B_1^2,$$
(13)
$`m`$ being the parameter defined by Eq. (10). First, we consider the case when the expression (13) holds (obviously, it may hold as long as it yields $`B_2^2>0`$). Substituting it into the stationary versions of Eqs. (3) and (4), we obtain

$`\mu _2^1A^{\prime \prime }+mA\left[\left(2\left(2m1\right)m^2\mu _2\right)A^2+\left(2m\mu _2\right)B_1^2\right]A`$ $`=`$ $`0,`$ (14)
$`\mu _2^1B_1^{\prime \prime }+B_1\left[\left(2\mu _2\right)B_1^2+\left(2m\mu _2\right)A^2\right]B_1`$ $`=`$ $`0,`$ (15)
the prime standing for $`d/dx`$. If, instead of Eq. (13), we take $`B_2=0`$, the stationary equations take the form

$`A^{\prime \prime }+2A\left[2\left(2\mu _2\right)A^2+\left(2m\mu _2\right)B_1^2\right]A`$ $`=`$ $`0,`$ (16)
$`B_1^{\prime \prime }+B_1\left[B_1^2+\left(2m\mu _2\right)A^2\right]B_1`$ $`=`$ $`0.`$ (17)
It is noteworthy that, although the SC patterns may be stable at $`\mu _1<0`$, i.e., $`m<0`$, Eq. (14) cannot have a solution for $`A(x)`$ exponentially decaying at $`x\mathrm{}`$ (see the b.c. (12)) unless $`m>0`$, hence the present problem has no solution with $`m<0`$.

Obviously, a solution to Eqs. (14) and (15) can satisfy, with regard to Eq. (13), the b.c. (12). However, the same set of equations (14) and (15) cannot satisfy the b.c. (11): setting $`B_1=0`$ and $`A=\mathrm{const}`$, one obtains from Eq. (14)

$$A^2=A_0^2\frac{m}{2\left(2m1\right)m^2\mu _2},$$
(18)
which is obviously different from the necessary limit value $`A_{lim}^21/\left(2\mu _2\right)`$. Moreover, the value of $`B_2^2`$ corresponding to $`A_0^2`$ as per Eq. (13) is different from zero, which also violates the b.c. (11).

In fact, the stationary state corresponding to the asymptotic value (18) with $`B_20`$ is another uniform pattern, which is a superposition of three spatial harmonics. This pattern is periodic in one spatial direction and quasiperiodic in the other one. As it was demonstrated in , this pattern is always dynamically unstable, hence a solution having it as an asymptotic state is physically irrelevant.

However, a solution to Eqs. (14) and (15) makes sense as long as it provides for $`B_2^2>0`$ according to Eq. (13). Taking the asymptotic state $`A^2=A_0^2`$, $`B_1^2=0`$, one finds that it gives rise to negative $`B_2^2`$ exactly in the case $`m<1`$, which is considered in this work, as this is the case when the SC pattern is stable, see above. Thus, Eqs. (14) and (15) should be used in the region $`\mathrm{}<x<x_0`$, where, by definition, $`x=x_0`$ is a point at which $`B_2^2\left(x\right)`$, as given by Eq. (13), vanishes, i.e.,

$$\left(2m\mu _2\right)A^2(x_0)+\left(1\mu _2\right)B_1^2(x_0)=1.$$
(19)
At the point $`x=x_0`$, one must switch from Eqs. (14) and (15) to Eqs. (16) and (17), setting $`B_20`$ at $`xx_0`$. The continuity dictates to take the values $`A(xx_0^{})`$ and $`B_1(xx_0^{})`$ as the b.c. to Eqs. (16) and (17) at $`x=x_0`$, the b.c. at $`x=+\mathrm{}`$ being fixed by Eqs. (11).

The transition from Eqs. (14) and (15) to Eqs. (16) and (17) at $`x=x_0`$ provides for the continuity of all the functions $`A(x)`$ and $`B_{1,2}(x)`$, but their first derivatives suffer a jump at $`x=x_0`$. In fact, if the above-mentioned fourth-derivative term is added to Eq. (5), the jump of the derivative will be smoothed down in a narrow boundary layer. It should be stressed that the presence of the jump does not violate the applicability of the description in terms of the GL equations; the only problem is the absence of a detailed description of the narrow boundary layer in which the jump is smoothed by the fourth-order derivative term, but details of the inner structure of the narrow layer do not affect the global picture.

Thus, the transient layer between the two SC domains consists of two sublayers, described, respectively, by Eqs. (14) and (15), and by Eqs. (16) and (17), with the discontinuity of the derivative between them. The first (left) sublayer contains all the four spatial harmonics, hence is locally seems as a QP pattern. An important finding is that, in the case of small $`\mu _2`$ considered here, the width of this sublayer is large, scaling $`\mu _2^{1/2}`$, according to Eqs. (14) and (15). The possibility to produce a broad QP stripe justifies all the consideration of the transient layer. The second (right) sublayer contains only three spatial harmonics, as $`B_20`$ in it, hence it locally looks like a pattern quasiperiodic only in one direction ), and its width is $`1`$ (i.e., not specifically large), according to Eqs. (16) and (17).

This qualitative analysis of the transient layer’s structure will be corroborated and illustrated by direct numerical results displayed in the next section. However, before using the numerical methods, one can notice that the structure of the second (right) sublayer can be described in an approximate analytical form if we consider a special case when the value $`A(x=x_0)`$ at the internal boundary between the two sublayers is already close to the asymptotic value $`A_{lim}`$, and, accordingly, the value $`B_1^2`$ is small (in other words, this is case when $`0<1m1)`$. Then, we may set

$$A(x)A_{lim}a(x),$$
(20)
where $`a(x)`$ is positive but small and vanishes at $`x=+\mathrm{}`$. Using the smallness of $`a(x)`$ and $`B_1`$, one can simplify Eqs. (16) and (17):

$`(1/\sqrt{2})a^{\prime \prime }\left[\sqrt{2}(2\mu _2)aB_1^2\right]\mu _1B_1^2`$ $`=`$ $`0,`$ (21)
$`B_1^{\prime \prime }+\left[\sqrt{2}(2\mu _2)aB_1^2\right]B_1{\displaystyle \frac{1}{2}}\left(\mu _12\mu _2\right)B_1`$ $`=`$ $`0.`$ (22)
The b.c. for $`a(x)`$ at $`x=x_0`$ can be obtained from the expansion of the exact b.c. (19):

$$a(x_0)=\frac{B_1^2(x_0)}{\sqrt{2}\left(2\mu _2\right)}.$$
(23)
The lowest-order approximate analytical solution to Eqs. (21) and (22), satisfying the necessary boundary conditions, is very simple:

$$B_1(x)=B_1(x_0)\mathrm{exp}(\lambda x),a(x)=a(x_0)\mathrm{exp}(2\lambda x),$$
(24)
with $`\lambda =\sqrt{\left(\mu _22\mu _1\right)/2}`$. Note that the condition necessary for $`\lambda `$ to be real is again exactly tantamount to $`m<1`$.

## 4 Numerical Results

To check the qualitative predictions for the structure of the transient layer obtained in the previous section, we performed a two-stage numerical integration of Eqs. (14) and (15), and then of Eqs. (16) and (17). To this end, the first pair of the equations was solved with the b.c. (12), continuing the solution until it hits a point where it satisfies Eq. (19) (recall the latter condition is equivalent to the vanishing of $`B_2^2`$). Actually, the numerical integration was performed, instead of the original coordinate $`x`$, in terms of a variable $`\xi \mathrm{tanh}x`$. This transformation is convenient because it maps the semi-infinite intervals $`(\pm \mathrm{},0)`$ of the variable $`x`$ into finite ones $`(\pm 1,0)`$.

The values $`A(x_0)`$ and $`B_1(x_0)`$, obtained from the solution of Eqs. (14) and (15), where then used to find a solution to Eqs. (16) and (17) in the interval $`\mathrm{tanh}x_0<\xi <1`$, satisfying the b.c. (11) at $`\xi =1`$. In accord with what said above, a condition of the continuity of the first derivatives across the point $`\xi =\mathrm{tanh}x_0`$ was not imposed.

Three typical examples of the thus obtained numerical solutions are displayed vs. the coordinate $`x`$ in Fig. 2, for three characteristic values $`m=0.75`$, $`m=0.5`$, and $`m=0.25`$, the parameter $`\mu _2`$, that must be small enough, being fixed in all the three cases as $`\mu _2=0.1`$. A characteristic feature clearly seen in all the cases is that, in accord with the prediction of the above analysis, a width of the left sublayer is essentially larger than that of the right one. It is also noteworthy that the solutions changes very little with a large change in $`m`$, i.e., the transient layer between the two SC domains is expected to be quite robust.

## 5 Conclusion

In this work, we have proposed an approach that makes it possible to the generate a two-dimensional quasiperiodic pattern in nonlinear dissipative systems where a direct generation of stable uniform quasiperiodic planforms is not possible. An eightfold pattern can be created in the form of a broad transient stripe between two domains filled by square cells, which are oriented under the angle of $`45`$ degrees relative to each other. Using the symmetry of the configuration considered, the structure of the pattern was described in terms of a system of three coupled real stationary Ginzburg-Landau equations, which were analyzed by means of analytical and numerical methods. It was found that the transient quasiperiodic pattern exists exactly in a parametric region in which the uniform square-cell pattern is stable. Further, it was found that the transient layer consists of two different sublayers, with a derivative jump between them (that can be smoothed into an additional narrow boundary layer, if higher-order derivatives are added to the Ginzburg-Landau equations). The width of the sublayer that features the eightfold quasiperiodic pattern is found to be large, while the other sublayer (filled with a less interesting pattern, which is quasiperiodic only in one direction) has a width $`1`$. A broad stripe of a twelvefold QP pattern can be similarly generated as a transient layer between two domains of hexagonal cells oriented at the angle of $`30`$ degrees.

It still remains to perform simulations of the full time-dependent Ginzburg-Landau equations, in order to directly test the dynamical stability of the broad transient layers. However, the numerically found robustness of the layers against the variation of the crucial control parameter $`m`$ suggests that they have a good chance to be dynamically stable.

## Figure Captions

Fig. 1. The configuration giving rise to the transient layer filled with the eightfold quasi-periodic pattern between two domains of square cells oriented under the angle $`45`$ degrees relative to each other..

Fig. 2. Numerically found structure of the transient layer corresponding to the configuration shown in Fig. 1. The small parameter $`\mu _2`$ defined by Eq. (9) is fixed to be $`0.1`$, while the control parameter $`m`$, defined by Eq. (10), takes values $`0.75`$ (a), $`0.50`$ (b), and $`0.25`$ (c) (recall only the values $`0<m<1`$ make sense in the present context).