text:
# 1 Introduction

## 1 Introduction

One of the main goals in the commissioning and operation of modern high energy accelerators and storage rings is the achievement of higher and higher beam currents and charge densities. That is why, of particular importance are the effects of intense self-fields due to space charge and current, influencing the beam propagation, its stability and transport properties. In general, a complete description of collective processes in intense charged particle beams is provided by the Vlasov-Maxwell equations for the self-consistent evolution of the beam distribution function and the electromagnetic fields. Although the analytical basis (as mentioned above) for modelling the dynamics and behaviour of space-charge dominated beams is well established, a thorough and satisfactory understanding of collective processes, detailed equilibrium and formation of patterns and coherent structures is far from being complete.

The present paper may be regarded as the first (to our knowledge) attempt to take a view at the description of the evolution and the collective behaviour of intense charged particle beams from an entirely different perspective, as compared to the ones available in the literature. We will be mainly interested in describing the slow evolution of some coarse-grained quantities that are easily measurable, such as the amplitudes of density waves. Due to the nonlinear wave interaction contingent on the nonlinear coupling between the Vlasov and Maxwell equations, one can expect a formation of nontrivial coherent structure that might be fairly stable in space and time. Here, we show that solitary wave patterns in the beam density distribution are an irrevocable feature, characteristic of intense beams. Moreover, density condensates in the special case where a parametric resonance in terms of wave frequency between a particular mode of the fundamental density waves and the external focusing occurs, can be formed.

The paper is organized as follows. It was previously shown  that in the case of a sheet beam with constant phase-space density the Vlasov-Maxwell equations are fully equivalent to a hydrodynamic model with zero heat flow and triple-adiabatic equation-of-state. For consistency, in section 2, we repeat the derivation from a slightly different perspective and then use the hydrodynamic equations as a fundament for the subsequent analysis. In section 3, we consider the case of a smooth focusing where the time variation of the $`\beta `$-function (respectively, the time variation of the density envelope function) can be neglected. Further, we demonstrate that starting from the hydrodynamic equations, and using the renormalization group (RG) technique  a system of coupled nonlinear Schrodinger equations for the slowly varying amplitudes of density waves can be derived. The purpose of section 4 is twofold. First, we study the case where a parametric resonance between a particular mode of density waves and a resonant Fourier harmonic of the external focusing (a resonant harmonic from the $`\sqrt{\beta }`$-function Fourier decomposition) occurs. It is shown that under certain conditions the evolution of the resonant amplitudes of the forward and backward density waves is governed by a system of coupled Gross-Pitaevskii equations (see e.g. the review  and the references therein). Secondly, it is demonstrated that in the non-resonant case the results obtained in the smooth focusing approximation coincide with the ones for a periodic focusing up to second order in the formal perturbation parameter. In section 5, the reduction of the infinite system of coupled nonlinear Schrodinger equations to a system of two coupled nonlinear Schrodinger equations is discussed. Such reduction is possible if one neglects the effect of all other modes and takes into account only the self-interaction of a single mode with a particular wave number. In the non-resonant case it is shown that the renormalized solution for the beam density describes the process of formation of holes (cavitons) in intense particle beams. The case where a parametric resonance is present is more interesting. In the cold-beam limit it is demonstrated that the evolution of the forward and the backward resonant wave amplitudes can be well described by a system of two coupled Gross-Pitaevskii equations. Finally, conclusions are drawn in section 6.

## 2 Derivation of the Hydrodynamic Model

We begin with the Hamiltonian describing the one-dimensional betatron motion in the presence of space-charge field

$$\widehat{H}=\frac{R}{2}\widehat{p}^2+\frac{G(\theta )\widehat{x}^2}{2R}+V(\widehat{x};\theta )+\frac{eR}{E_s\beta _s^2}\phi _{sc}(\widehat{x};\theta ),$$
(2.1)
where $`e`$ is the electron charge, $`R`$ is the mean radius of the accelerator, $`G(\theta )`$ is the focusing strength of the linear machine lattice, $`V(\widehat{x};\theta )`$ is a contribution coming from nonlinear lattice elements (sextupoles, octupoles, etc.) and $`E_s`$ and $`\beta _s`$ are the energy and the relative velocity of the synchronous particle respectively. In addition, $`\phi _{sc}(\widehat{x};\theta )`$ is the self-field potential due to space-charge, which satisfies the one-dimensional Poisson equation

$$\frac{^2\phi _{sc}}{\widehat{x}^2}=\frac{en}{ϵ_0}d\widehat{p}f(\widehat{x},\widehat{p};\theta ).$$
(2.2)
Here $`n=N_p/V_t`$ is the density of beam particles ($`N_p`$ is the number of particles in the beam, while $`V_t`$ is the volume of the area occupied by the beam in the transverse direction), $`ϵ_0`$ is the dielectric susceptibility and $`f(\widehat{x},\widehat{p};\theta )`$ is the distribution function in phase space. Equations (2.1) and (2.2) can be written in an alternate form as

$$\widehat{H}=\frac{R}{2}\widehat{p}^2+\frac{G(\theta )\widehat{x}^2}{2R}+V(\widehat{x};\theta )+\lambda \stackrel{~}{\phi }(\widehat{x};\theta ),$$
(2.3)
$$\frac{^2\stackrel{~}{\phi }}{\widehat{x}^2}=d\widehat{p}f(\widehat{x},\widehat{p};\theta ),$$
(2.4)
where

$$\lambda =\frac{Re^2n}{ϵ_0E_s\beta _s^2},$$
(2.5)
is the beam perveance, and $`\phi _{sc}=en\stackrel{~}{\phi }/ϵ_0`$. Note that the beam perveance $`\lambda `$ is a dimensionless quantity.

Next, we perform a canonical transformation

$$\widehat{x}=x\sqrt{\beta },\widehat{p}=\frac{1}{\sqrt{\beta }}(p\alpha x),$$
(2.6)
defined by the generating function

$$F_2(\widehat{x},p;\theta )=\frac{\widehat{x}p}{\sqrt{\beta }}\frac{\alpha \widehat{x}^2}{2\beta },$$
(2.7)
where $`\alpha `$ and $`\beta `$ (and $`\gamma `$) are the well-known Twiss parameters satisfying the equations

$$\frac{\mathrm{d}\alpha }{\mathrm{d}\theta }=\frac{G\beta }{R}R\gamma ,\frac{\mathrm{d}\beta }{\mathrm{d}\theta }=2R\alpha ,\beta \gamma \alpha ^2=1.$$
(2.8)
As a result, we obtain the new Hamiltonian

$$H=\frac{\dot{\chi }}{2}\left(p^2+x^2\right)+V(x;\theta )+\lambda \sqrt{\beta }U(x;\theta ),$$
(2.9)
where the self-field potential $`U(x;\theta )`$ $`\left(\stackrel{~}{\phi }=\sqrt{\beta }U\right)`$ satisfies the equation

$$\frac{^2U}{x^2}=dpf(x,p;\theta ),$$
(2.10)
and

$$\dot{\chi }=\frac{\mathrm{d}\chi }{\mathrm{d}\theta }=\frac{R}{\beta },$$
(2.11)
is the derivative of the phase advance with respect to $`\theta `$.

We are now ready to write the Vlasov equation for the one-particle distribution function $`f(x,p;\theta )`$ in the two-dimensional phase space $`(x,p)`$. It reads as

$$\frac{f}{\theta }+\dot{\chi }p\frac{f}{x}\left(\dot{\chi }x+\frac{V}{x}+\lambda \sqrt{\beta }\frac{U}{x}\right)\frac{f}{p}=0,$$
(2.12)
and should be solved self-consistently with the Poisson equation (2.10). Following Davidson et al. , we consider the case where the distribution function $`f(x,p;\theta )`$ is constant (independent of $`x`$, $`p`$ and $`\theta `$) inside a region in phase space confined by the simply connected boundary curves $`p_{(+)}(x;\theta )`$ and $`p_{()}(x;\theta )`$, and zero outside. In other words,

$$f(x,p;\theta )=𝒞,\mathrm{for}p_{()}(x;\theta )<p<p_{(+)}(x;\theta ),$$
(2.13)
and $`f(x,p;\theta )=0`$ otherwise. The proof that such solution to the Vlasov equation (2.12) exists can be found in appendix A.

The equations for the boundary curves $`p_{(\pm )}(x;\theta )`$ can be derived as follows. We substitute the explicit form

$$f(x,p;\theta )=𝒞\left[\left(pp_{()}(x;\theta )\right)\left(pp_{(+)}(x;\theta )\right)\right]$$
(2.14)
of the distribution function into the Vlasov equation (2.12). Here $`(z)`$ is the well-known Heaviside function. Using the fact that the derivative of the Heaviside function with respect to its argument equals the Dirac $`\delta `$-function, we obtain

$`\delta \left(pp_{()}\right){\displaystyle \frac{p_{()}}{\theta }}+\delta \left(pp_{(+)}\right){\displaystyle \frac{p_{(+)}}{\theta }}+\dot{\chi }p\left[\delta \left(pp_{()}\right){\displaystyle \frac{p_{()}}{x}}+\delta \left(pp_{(+)}\right){\displaystyle \frac{p_{(+)}}{x}}\right]`$
$$\left(\dot{\chi }x+\frac{V}{x}+\lambda \sqrt{\beta }\frac{U}{x}\right)\left[\delta \left(pp_{()}\right)\delta \left(pp_{(+)}\right)\right]=0.$$
(2.15)
Multiplying equation (2.15) first by $`1`$ and then by $`p`$, and integrating the result over $`p`$, we readily obtain

$$\frac{}{\theta }\left(p_{(+)}p_{()}\right)+\frac{\dot{\chi }}{2}\frac{}{x}\left(p_{(+)}^2p_{()}^2\right)=0,$$
(2.16)
$$\frac{1}{2}\frac{}{\theta }\left(p_{(+)}^2p_{()}^2\right)+\frac{\dot{\chi }}{3}\frac{}{x}\left(p_{(+)}^3p_{()}^3\right)=\left(p_{(+)}p_{()}\right)\left(\dot{\chi }x+\frac{V}{x}+\lambda \sqrt{\beta }\frac{U}{x}\right),$$
(2.17)
$$\frac{^2U}{x^2}=𝒞\left(p_{(+)}p_{()}\right),$$
(2.18)
It is convenient to recast equations (2.16) - (2.18) in a more familiar form, widely used in hydrodynamics. Let us define

$$\rho =\underset{\mathrm{}}{\overset{\mathrm{}}{}}dpf(x,p;\theta )=𝒞\left(p_{(+)}p_{()}\right),$$
(2.19)
$$\rho v=\underset{\mathrm{}}{\overset{\mathrm{}}{}}dppf(x,p;\theta )=\frac{𝒞}{2}\left(p_{(+)}^2p_{()}^2\right),$$
(2.20)
where $`\rho (x;\theta )`$ and $`v(x;\theta )`$ are the density and the current velocity, respectively. From the two definitions above, it follows that

$$v=\frac{1}{2}\left(p_{(+)}+p_{()}\right).$$
(2.21)
Defining further the pressure $`𝒫(x;\theta )`$ and the heat flow $`𝒬(x;\theta )`$ and using equation (2.21), we have

$$𝒫=\underset{\mathrm{}}{\overset{\mathrm{}}{}}dp(pv)^2f(x,p;\theta )=\frac{1}{12𝒞^2}\left[𝒞\left(p_{(+)}p_{()}\right)\right]^3,$$
(2.22)
$$𝒬=\underset{\mathrm{}}{\overset{\mathrm{}}{}}dp(pv)^3f(x,p;\theta )=0.$$
(2.23)
From equation (2.19), it follows that the pressure can be expressed as

$$𝒫=\frac{𝒫_0}{\widehat{\rho }_0^3}\rho ^3,\frac{𝒫_0}{\widehat{\rho }_0^3}=\frac{1}{12𝒞^2}.$$
(2.24)
In addition, it is straightforward to verify that

$$\frac{𝒞}{3}\left(p_{(+)}^3p_{()}^3\right)=𝒫+\rho v^2,$$
(2.25)
which provides a closure for the equations (2.16) and (2.17) governing the evolution of the moments (boundary curves). With all the above definitions and relations in hand, we can write the system of hydrodynamic equations (supplemented by the Poisson equation) in the form

$$\frac{\rho }{\theta }+\dot{\chi }\frac{}{x}(\rho v)=0,$$
(2.26)
$$\frac{}{\theta }(\rho v)+\dot{\chi }\frac{}{x}\left(𝒫+\rho v^2\right)=\rho \left(\dot{\chi }x+\frac{V}{x}+\lambda \sqrt{\beta }\frac{U}{x}\right),$$
(2.27)
$$\frac{^2U}{x^2}=\rho .$$
(2.28)
Using the continuity equation (2.26), we cast the system of equations (2.26) - (2.28) in its final form

$$\frac{\rho }{\tau }+\frac{}{x}(\rho v)=0,$$
(2.29)
$$\frac{v}{\tau }+v\frac{v}{x}+v_T^2\frac{}{x}\left(\rho ^2\right)=x\frac{\beta }{R}\frac{V}{x}\lambda \frac{\beta \sqrt{\beta }}{R}\frac{U}{x},$$
(2.30)
$$\frac{^2U}{x^2}=\rho ,$$
(2.31)
which will be the starting point for the subsequent analysis. Here

$$v_T^2=\frac{3𝒫_0}{2\widehat{\rho }_0^3},$$
(2.32)
is the normalized thermal speed-squared, and $`\tau =\chi (\theta )+\tau _0`$ is the new independent ”time”-variable.

## 3 Renormalization Group Reduction of the Hydrodynamic Equations

Before proceeding further, we make an important remark. Let $`F(\theta )`$ be a periodic function of $`\theta `$ with period $`2\pi `$. Noting that the phase advance $`\chi (\theta )`$ can be represented in the form

$$\chi (\theta )=\nu \theta +\chi _p(\theta )\chi _p(\pi ),$$
(3.1)
where $`\nu `$ is the betatron tune, and $`\chi _p(\theta +2\pi )=\chi _p(\theta )`$, we can expand $`F(\theta )`$ regarded as a function of $`\chi `$ (respectively $`\tau `$) in a Fourier series in the new variable $`\chi `$ (respectively $`\tau `$) as follows

$$F(\theta )=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}𝒜_n\mathrm{exp}\left(in\frac{\chi }{\nu }\right).$$
(3.2)
The expansion coefficients $`𝒜_n`$ are given by the expressions

$$𝒜_n=\frac{1}{2\pi \nu }\underset{\pi \nu }{\overset{\pi \nu }{}}d\chi F(\theta )\mathrm{exp}\left(in\frac{\chi }{\nu }\right).$$
(3.3)
Using the definition of the phase advance (3.1), and choosing $`\tau _0=\chi _p(\pi )`$, we can rewrite the Fourier expansion (3.2) as

$$F(\theta )=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}_n\mathrm{exp}\left(in\frac{\tau }{\nu }\right),$$
(3.4)
where the expansion coefficients $`_n`$ are given by the expressions

$$_n=\frac{R}{2\pi \nu }\underset{\pi }{\overset{\pi }{}}d\theta \frac{F(\theta )}{\beta (\theta )}e^{in\theta }\mathrm{exp}\left[in\frac{\chi _p(\theta )}{\nu }\right].$$
(3.5)
In what follows, we consider the case where the external potential $`V`$ is zero, so that the set of equations to be analyzed acquires the form

$$\frac{\rho }{\tau }+\frac{}{x}(\rho v)=0,$$
(3.6)
$$\frac{v}{\tau }+v\frac{v}{x}+v_T^2\frac{}{x}\left(\rho ^2\right)=x\lambda g(\tau )\frac{U}{x},$$
(3.7)
$$\frac{^2U}{x^2}=\rho ,$$
(3.8)
where

$$g(\tau )=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}a_n\mathrm{exp}\left(in\frac{\tau }{\nu }\right),$$
(3.9)
$$a_n=\frac{1}{2\pi \nu }\underset{\pi }{\overset{\pi }{}}d\theta \sqrt{\beta (\theta )}e^{in\theta }\mathrm{exp}\left[in\frac{\chi _p(\theta )}{\nu }\right].$$
(3.10)
Further, we assume that there exists a nontrivial solution to equations (3.6) - (3.8) in the interval $`x(x_{()},x_{(+)})`$, and that the sheet beam density is zero $`\left(\varrho =0\right)`$ outside of the interval. We introduce also the ansatz

$$\rho (x;\tau )=\frac{1}{}+ϵR(x;\tau ),v(x;\tau )=\frac{x}{}\frac{\mathrm{d}}{\mathrm{d}\tau }+ϵu(x;\tau ),U(x;\tau )=\frac{x^2}{2}+ϵ𝒰(x;\tau ),$$
(3.11)
where the envelope function $`(\tau )`$ is a solution to the equation

$$\frac{\mathrm{d}^2}{\mathrm{d}\tau ^2}+=\lambda g(\tau ),$$
(3.12)
and therefore can be represented in explicit form as

$$(\tau )=\lambda a_0+\lambda \underset{n0}{}b_ne^{in\tau /\nu },b_n=\frac{a_n}{1n^2/\nu ^2}.$$
(3.13)
It enables us to rewrite equations (3.6) - (3.8) as follows

$$\frac{R}{\tau }+\frac{1}{}\frac{u}{x}+\frac{\dot{}}{}\frac{}{x}(xR)+ϵ\frac{}{x}\left(Ru\right)=0,$$
(3.14)
$$\frac{u}{\tau }+\frac{\dot{}}{}\frac{}{x}(xu)+\frac{2v_T^2}{}\frac{R}{x}+ϵu\frac{u}{x}+ϵv_T^2\frac{}{x}\left(R^2\right)=\lambda g\frac{𝒰}{x},$$
(3.15)
$`{\displaystyle \frac{^2𝒰}{x^2}}=R.`$
Next, we differentiate the first equation (3.14) of the above system with respect to $`\tau `$ and the second equation (3.15) with respect to $`x`$. After summing these up, we obtain

$`{\displaystyle \frac{}{\tau }}\left({\displaystyle \frac{R}{\tau }}\right)2v_T^2{\displaystyle \frac{^2R}{x^2}}+\lambda gR`$
$$=\dot{}\frac{^2}{x^2}(xu)\frac{}{\tau }\left[\dot{}\frac{}{x}(xR)\right]+ϵ\frac{^2}{x^2}\left(\frac{u^2}{2}+v_T^2R^2\right)ϵ\frac{}{\tau }\left[\frac{}{x}(Ru)\right].$$
(3.16)
Equation (3.16) supplemented with equation (3.14) rewritten as

$$\frac{R}{\tau }+\frac{u}{x}+\dot{}\frac{}{x}(xR)+ϵ\frac{}{x}\left(Ru\right)=0,$$
(3.17)
comprises the basic system of equations for the analysis in the subsequent exposition.

In what follows, we consider the case of smooth focusing, where the time variation of $`g(\tau )`$ \[$`(\tau )`$, respectively\] can be neglected. In the next section, we will show that when the frequencies of all fundamental modes are sufficiently far from a parametric resonance, this assumption holds true to second order in the perturbation parameter $`ϵ`$ even in the case when such time variation is present. Thus, the basic equations (3.16) and (3.17) acquire the form

$$_0^2\frac{^2R}{\tau ^2}2v_T^2\frac{^2R}{x^2}+_0^2R=ϵ_0\frac{^2}{x^2}\left(\frac{u^2}{2}+v_T^2R^2\right)ϵ_0^2\frac{^2}{\tau x}(Ru),$$
(3.18)
$$_0\frac{R}{\tau }+\frac{u}{x}+ϵ_0\frac{}{x}\left(Ru\right)=0,$$
(3.19)
where

$$_0=\lambda a_0.$$
(3.20)
Let us further assume that the actual dependence of $`R`$ and $`u`$ on the independent variables is given by

$$R=R(x,\xi ;\tau ),u=u(x,\xi ;\tau ),$$
(3.21)
where

$$\xi =ϵx$$
(3.22)
is a slow spatial variable. Therefore, the basic equations (3.18) and (3.19) can be rewritten as

$`_0^2{\displaystyle \frac{^2R}{\tau ^2}}2v_T^2\left({\displaystyle \frac{^2}{x^2}}+2ϵ{\displaystyle \frac{^2}{x\xi }}+ϵ^2{\displaystyle \frac{^2}{\xi ^2}}\right)R+_0^2R`$
$$=ϵ_0\left(\frac{^2}{x^2}+2ϵ\frac{^2}{x\xi }+ϵ^2\frac{^2}{\xi ^2}\right)\left(\frac{u^2}{2}+v_T^2R^2\right)ϵ_0^2\left(\frac{^2}{\tau x}+ϵ\frac{^2}{\tau \xi }\right)(Ru),$$
(3.23)
$$_0\frac{R}{\tau }+\frac{u}{x}+ϵ\frac{u}{\xi }+ϵ_0\left(\frac{}{x}+ϵ\frac{}{\xi }\right)\left(Ru\right)=0,$$
(3.24)
Following the basic idea of the RG method, we represent the solution to equations (3.23) and (3.24) in the form of a standard perturbation expansion  in the formal small parameter $`ϵ`$ as

$$R=\underset{k=0}{\overset{\mathrm{}}{}}ϵ^kR_k,u=\underset{k=0}{\overset{\mathrm{}}{}}ϵ^ku_k,$$
(3.25)
The zeroth-order equations (3.23) and (3.24) read as

$$\lambda ^2a_0^2\frac{^2R_0}{\tau ^2}2v_T^2\frac{^2R_0}{x^2}+\lambda ^2a_0^2R_0=0,$$
(3.26)
$$\lambda a_0\frac{R_0}{\tau }+\frac{u_0}{x}=0.$$
(3.27)
Their solutions can be found in a straightforward manner to be

$$R_0(x,\xi ;\tau )=\underset{m0}{}\left[A_m(\xi )e^{iz_m^{(+)}}+B_m(\xi )e^{iz_m^{()}}\right],$$
(3.28)
$$u_0(x,\xi ;\tau )=\frac{2\pi }{\sigma ^2}\underset{m0}{}\frac{\omega _m}{m}\left[A_m(\xi )e^{iz_m^{(+)}}B_m(\xi )e^{iz_m^{()}}\right].$$
(3.29)
Here $`A_m`$ and $`B_m`$ are constant complex amplitudes (to this end dependent on the slow variable $`\xi `$ only) of the backward and the forward wave solution to equation (3.26), respectively. These will be the subject of renormalization at the final step of the renormalization group procedure resulting in RG equations governing their slow evolution. Furthermore,

$$z_m^{(\pm )}(x;\tau )=\omega _m\tau \pm m\sigma x,\sigma =\frac{2\pi }{x_{(+)}+x_{()}},\lambda a_0=\frac{2\pi }{\sigma }.$$
(3.30)
The discrete mode frequencies $`\omega _m`$ are determined from the dispersion relation

$$\omega _m^2=1+K^2m^2,K=\frac{v_T^2\sigma ^4}{2\pi ^2}.$$
(3.31)
It can be easily verified that the above choice of parameters leads to

$$\underset{x_{()}}{\overset{x_{(+)}}{}}dxR_0(x;\theta )=0,$$
(3.32)
which means that linear perturbation to the uniform density $`\rho _0=_0^1`$ average to zero and do not affect the normalization properties on the interval $`(x^{()},x^{(+)})`$. In addition, the following conventions and notations

$$\omega _m=\omega _m,A_m=A_m^{},B_m=B_m^{}$$
(3.33)
have been introduced.

The first order perturbation equation for $`R_1`$ reads as

$`{\displaystyle \frac{^2R_1}{\tau ^2}}{\displaystyle \frac{K^2}{\sigma ^2}}{\displaystyle \frac{^2R_1}{x^2}}+R_1={\displaystyle \frac{2iK^2}{\sigma }}{\displaystyle \underset{m0}{}}m\left({\displaystyle \frac{A_m}{\xi }}e^{iz_m^{(+)}}{\displaystyle \frac{B_m}{\xi }}e^{iz_m^{()}}\right)`$
$$\frac{\pi }{\sigma }\underset{m,n0}{}\left[\gamma _{mn}^{(+)}A_mA_ne^{i\left(z_m^{(+)}+z_n^{(+)}\right)}+2\gamma _{mn}^{()}A_mB_ne^{i\left(z_m^{(+)}+z_n^{()}\right)}+\gamma _{mn}^{(+)}B_mB_ne^{i\left(z_m^{()}+z_n^{()}\right)}\right],$$
(3.34)
where

$$\gamma _{mn}^{(\pm )}=(m\pm n)\left[\left(\omega _m+\omega _n\right)\left(\frac{\omega _m}{m}\pm \frac{\omega _n}{n}\right)+(m\pm n)\left(K^2\pm \frac{\omega _m\omega _n}{mn}\right)\right].$$
(3.35)
The solution for $`R_1`$ is readily obtained to be

$`R_1(x,\xi ;\tau )={\displaystyle \frac{K^2\tau }{\sigma }}{\displaystyle \underset{m0}{}}{\displaystyle \frac{m}{\omega _m}}\left({\displaystyle \frac{A_m}{\xi }}e^{iz_m^{(+)}}{\displaystyle \frac{B_m}{\xi }}e^{iz_m^{()}}\right)`$
$$\frac{\pi }{\sigma }\underset{m,n0}{}\left[\alpha _{mn}^{(+)}A_mA_ne^{i\left(z_m^{(+)}+z_n^{(+)}\right)}+2\alpha _{mn}^{()}A_mB_ne^{i\left(z_m^{(+)}+z_n^{()}\right)}+\alpha _{mn}^{(+)}B_mB_ne^{i\left(z_m^{()}+z_n^{()}\right)}\right],$$
(3.36)
where

$$\alpha _{mn}^{(\pm )}=\frac{\gamma _{mn}^{(\pm )}}{𝒟_{mn}^{(\pm )}},$$
(3.37)
$$𝒟_{mn}^{(\pm )}=1\left(\omega _m+\omega _n\right)^2+K^2(m\pm n)^2.$$
(3.38)
Note that the (infinite dimensional) matrices $`\widehat{\alpha }^{(\pm )}`$ are symmetric, i.e.

$$\alpha _{mn}^{(\pm )}=\alpha _{nm}^{(\pm )},\alpha _{m,m}^{(\pm )}=0.$$
(3.39)
Having determined $`R_1(x;\tau )`$, the first-order current velocity $`u_1(x;\tau )`$ can be found in a straightforward manner. The result is

$`u_1(x,\xi ;\tau )={\displaystyle \frac{2\pi K^2\tau }{\sigma ^3}}{\displaystyle \underset{m0}{}}\left({\displaystyle \frac{A_m}{\xi }}e^{iz_m^{(+)}}+{\displaystyle \frac{B_m}{\xi }}e^{iz_m^{()}}\right)`$
$`{\displaystyle \frac{2\pi i}{\sigma ^3}}{\displaystyle \underset{m0}{}}{\displaystyle \frac{1}{m^2\omega _m}}\left({\displaystyle \frac{A_m}{\xi }}e^{iz_m^{(+)}}+{\displaystyle \frac{B_m}{\xi }}e^{iz_m^{()}}\right)`$
$$+\frac{2\pi ^2}{\sigma ^3}\underset{m,n0}{}\left[\beta _{mn}^{(+)}A_mA_ne^{i\left(z_m^{(+)}+z_n^{(+)}\right)}+2\beta _{mn}^{()}A_mB_ne^{i\left(z_m^{(+)}+z_n^{()}\right)}\beta _{mn}^{(+)}B_mB_ne^{i\left(z_m^{()}+z_n^{()}\right)}\right],$$
(3.40)
where

$$\beta _{mn}^{(\pm )}=\frac{\omega _m}{m}\pm \frac{\omega _n}{n}+\alpha _{mn}^{(\pm )}\frac{\omega _m+\omega _n}{m\pm n},\beta _{m,m}^{(+)}=0.$$
(3.41)
Note that the (infinite dimensional) matrix $`\widehat{\beta }^{(+)}`$ possesses the same symmetry properties as those displayed by equation (3.39) possessed by the matrix $`\widehat{\alpha }^{(+)}`$, while $`\widehat{\beta }^{()}`$ is antisymmetric (and evidently $`\beta _{mm}^{()}=0`$).

In second order, the equation for $`R_2(x,\xi ;\tau )`$ acquires a form similar to that of equation (3.34). It is important to emphasize that this assertion holds true in every subsequent order. Each entry on the right-hand-sides of the corresponding equations can be calculated explicitly utilizing the already determined quantities from the previous orders. The right-hand-side of the equation for $`R_2(x,\xi ;\tau )`$ contains terms which yield oscillating contributions with constant amplitudes to the solution for $`R_2(x,\xi ;\tau )`$. Apart from these, there are resonant terms (proportional to $`\mathrm{e}^{iz_m^{(\pm )}(x;\tau )}`$) leading to a secular contribution. To complete the renormalization group reduction of the hydrodynamic equations, we select these particular resonant second-order terms on the right-hand-side of the equation determining $`R_2(x,\xi ;\tau )`$. The latter can be written as

$`{\displaystyle \frac{^2R_2}{\tau ^2}}{\displaystyle \frac{K^2}{\sigma ^2}}{\displaystyle \frac{^2R_2}{x^2}}+R_2={\displaystyle \frac{K^2}{\sigma ^2}}{\displaystyle \underset{m0}{}}\left({\displaystyle \frac{^2A_m}{\xi ^2}}e^{iz_m^{(+)}}+{\displaystyle \frac{^2B_m}{\xi ^2}}e^{iz_m^{()}}\right)`$
$`+{\displaystyle \frac{2iK^4\tau }{\sigma ^2}}{\displaystyle \underset{m0}{}}{\displaystyle \frac{m^2}{\omega _m}}\left({\displaystyle \frac{^2A_m}{\xi ^2}}e^{iz_m^{(+)}}+{\displaystyle \frac{^2B_m}{\xi ^2}}e^{iz_m^{()}}\right)`$
$$+\frac{4\pi ^2}{\sigma ^2}\underset{m,n0}{}\left[\left(\mathrm{\Gamma }_{mn}^{(+)}\left|A_n\right|^2+\mathrm{\Gamma }_{mn}^{()}\left|B_n\right|^2\right)A_me^{iz_m^{(+)}}+\left(\mathrm{\Gamma }_{mn}^{()}\left|A_n\right|^2+\mathrm{\Gamma }_{mn}^{(+)}\left|B_n\right|^2\right)B_me^{iz_m^{()}}\right],$$
(3.42)
where

$$\mathrm{\Gamma }_{mn}^{(+)}=m^2\left[\beta _{mn}^{(+)}\left(\frac{\omega _m}{m}+\frac{\omega _n}{n}\right)+\alpha _{mn}^{(+)}\left(K^2+\frac{\omega _m\omega _n}{mn}\right)\right]\left(1\frac{\delta _{mn}}{2}\right),$$
(3.43)
$$\mathrm{\Gamma }_{mn}^{()}=m^2\left[\beta _{mn}^{()}\left(\frac{\omega _m}{m}\frac{\omega _n}{n}\right)+\alpha _{mn}^{()}\left(K^2\frac{\omega _m\omega _n}{mn}\right)\right],$$
(3.44)
Some straightforward algebra yields the solution for $`R_2(x,\xi ;\tau )`$ in the form

$`R_2(x,\xi ;\tau )={\displaystyle \frac{K^4\tau ^2}{2\sigma ^2}}{\displaystyle \underset{m0}{}}{\displaystyle \frac{m^2}{\omega _m^2}}\left({\displaystyle \frac{^2A_m}{\xi ^2}}e^{iz_m^{(+)}}+{\displaystyle \frac{^2B_m}{\xi ^2}}e^{iz_m^{()}}\right)`$
$`+{\displaystyle \frac{K^2\tau }{2i\sigma ^2}}{\displaystyle \underset{m0}{}}{\displaystyle \frac{1}{\omega _m^3}}\left({\displaystyle \frac{^2A_m}{\xi ^2}}e^{iz_m^{(+)}}+{\displaystyle \frac{^2B_m}{\xi ^2}}e^{iz_m^{()}}\right)`$
$$+\frac{2\pi ^2\tau }{i\sigma ^2}\underset{m,n0}{}\frac{1}{\omega _m}\left[\left(\mathrm{\Gamma }_{mn}^{(+)}\left|A_n\right|^2+\mathrm{\Gamma }_{mn}^{()}\left|B_n\right|^2\right)A_me^{iz_m^{(+)}}+\left(\mathrm{\Gamma }_{mn}^{()}\left|A_n\right|^2+\mathrm{\Gamma }_{mn}^{(+)}\left|B_n\right|^2\right)B_me^{iz_m^{()}}\right],$$
(3.45)
where non-secular oscillating terms have not been written out in full explicitly. The final step is to collect the terms proportional to the fundamental modes $`e^{iz_m^{(+)}}`$ and $`e^{iz_m^{()}}`$ in all orders and renormalize the amplitudes $`A_m`$ and $`B_m`$. As a result one obtains the following RG equations

$$2i\omega _m\frac{A_m}{\tau }2im\frac{K^2}{\sigma }\frac{A_m}{x}=\frac{K^2}{\sigma ^2\omega _m^2}\frac{^2A_m}{x^2}+\frac{4\pi ^2}{\sigma ^2}A_m\underset{n0}{}\left(\mathrm{\Gamma }_{mn}^{(+)}\left|A_n\right|^2+\mathrm{\Gamma }_{mn}^{()}\left|B_n\right|^2\right),$$
(3.46)
$$2i\omega _m\frac{B_m}{\tau }+2im\frac{K^2}{\sigma }\frac{B_m}{x}=\frac{K^2}{\sigma ^2\omega _m^2}\frac{^2B_m}{x^2}+\frac{4\pi ^2}{\sigma ^2}B_m\underset{n0}{}\left(\mathrm{\Gamma }_{mn}^{()}\left|A_n\right|^2+\mathrm{\Gamma }_{mn}^{(+)}\left|B_n\right|^2\right),$$
(3.47)
## 4 The Parametric Resonance

Let us now address the system of equations (3.16) and (3.17). Without loss of generality, we assume that the time variation of $`g(\tau )`$ can be treated as a second-order perturbation (which is usually the case), that is

$$g(\tau )=a_0+ϵ^2\underset{n0}{}a_ne^{in\tau /\nu }.$$
(4.1)
The same holds true for the envelope function $`(\tau )`$

$$(\tau )=\lambda a_0+ϵ^2\lambda \underset{n0}{}b_ne^{in\tau /\nu },b_n=\frac{a_n}{1n^2/\nu ^2}.$$
(4.2)
Thus, the basic equations (3.16) and (3.17) can be rewritten in the form

$`{\displaystyle \frac{}{\tau }}\left({\displaystyle \frac{R}{\tau }}\right)2v_T^2{\displaystyle \frac{^2R}{x^2}}+\lambda gR`$
$$=ϵ^2\dot{}_2\frac{^2}{x^2}(xu)ϵ^2\frac{}{\tau }\left[\dot{}_2\frac{}{x}(xR)\right]+ϵ\frac{^2}{x^2}\left(\frac{u^2}{2}+v_T^2R^2\right)ϵ\frac{}{\tau }\left[\frac{}{x}(Ru)\right].$$
(4.3)
$$\frac{R}{\tau }+\frac{u}{x}+ϵ^2\dot{}_2\frac{}{x}(xR)+ϵ\frac{}{x}\left(Ru\right)=0,$$
(4.4)
Let us reiterate that the assumption concerning the second order of magnitude in $`ϵ`$ of the time variation of $`g(\tau )`$ and $`(\tau )`$ does not restrict the generality of the subsequent analysis. If this variation were of first order in $`ϵ`$, the proper perturbation parameter to use would be $`\sqrt{ϵ}`$ instead of $`ϵ`$. In addition, the variables $`R`$ and $`u`$ have to be rescaled accordingly

$$R\frac{R}{\sqrt{ϵ}},u\frac{u}{\sqrt{ϵ}}.$$
(4.5)
In terms of the ansatz (3.21) and (3.22), equations (4.3) and (4.4) become

$`{\displaystyle \frac{}{\tau }}\left({\displaystyle \frac{R}{\tau }}\right)2v_T^2\left({\displaystyle \frac{^2}{x^2}}+2ϵ{\displaystyle \frac{^2}{x\xi }}+ϵ^2{\displaystyle \frac{^2}{\xi ^2}}\right)R+\lambda gR`$
$`=ϵ\dot{}_2\left({\displaystyle \frac{^2}{x^2}}+2ϵ{\displaystyle \frac{^2}{x\xi }}+ϵ^2{\displaystyle \frac{^2}{\xi ^2}}\right)\left(\xi u\right)ϵ{\displaystyle \frac{}{\tau }}\left[\dot{}_2\left({\displaystyle \frac{}{x}}+ϵ{\displaystyle \frac{}{\xi }}\right)\left(\xi R\right)\right]`$
$$+ϵ\left(\frac{^2}{x^2}+2ϵ\frac{^2}{x\xi }+ϵ^2\frac{^2}{\xi ^2}\right)\left(\frac{u^2}{2}+v_T^2R^2\right)ϵ\frac{}{\tau }\left[\left(\frac{}{x}+ϵ\frac{}{\xi }\right)\left(Ru\right)\right],$$
(4.6)
$$\frac{R}{\tau }+\frac{u}{x}+ϵ\frac{u}{\xi }+ϵ\dot{}_2\left(\frac{}{x}+ϵ\frac{}{\xi }\right)\left(\xi R\right)+ϵ\left(\frac{}{x}+ϵ\frac{}{\xi }\right)\left(Ru\right)=0.$$
(4.7)
Although the general case can be in principle treated through more labour-intensive manipulations, for the sake of simplicity in what follows, we select a particular mode with mode number $`m`$. The zeroth-order solution can be written in the form

$$R_0(x,\xi ;\tau )=A\left(\xi \right)e^{iz^{(+)}(x,\tau )}+B\left(\xi \right)e^{iz^{()}(x,\tau )}+c.c.,$$
(4.8)
$$u_0(x,\xi ;\tau )=\frac{2\pi \omega }{m\sigma ^2}\left[A\left(\xi \right)e^{iz^{(+)}(x,\tau )}B\left(\xi \right)e^{iz^{()}(x,\tau )}\right]+c.c.,$$
(4.9)
where again

$`z^{(\pm )}(x,\tau )=\omega \tau \pm m\sigma x,\omega ^2=1+K^2m^2.`$
The particular mode is chosen such that an exact parametric resonance (if possible)

$$2\omega =\frac{n_0}{\nu }$$
(4.10)
occurs for some integer $`n_0`$. The first order perturbation equation for $`R_1`$ can be written in the form

$`{\displaystyle \frac{^2R_1}{\tau ^2}}{\displaystyle \frac{K^2}{\sigma ^2}}{\displaystyle \frac{^2R_1}{x^2}}+R_1={\displaystyle \frac{2imK^2}{\sigma }}\left({\displaystyle \frac{A}{\xi }}e^{iz^{(+)}}{\displaystyle \frac{B}{\xi }}e^{iz^{()}}\right)`$
$`+{\displaystyle \frac{im\sigma ^2}{2\pi }}\xi \lambda {\displaystyle \underset{n0}{}^{}}{\displaystyle \frac{n}{\nu }}\left(2\omega +{\displaystyle \frac{n}{\nu }}\right)b_ne^{in\tau /\nu }\left(Ae^{iz^{(+)}}Be^{iz^{()}}\right)`$
$$\frac{4\pi }{\sigma }\left(4\omega ^21\right)\left(A^2e^{2iz^{(+)}}+B^2e^{2iz^{()}}\right)+\frac{8\pi }{\sigma }AB^{}e^{i\left(z^{(+)}z^{()}\right)}+c.c..$$
(4.11)
Here $`^{}`$ implies exclusion of the harmonic with $`n=n_0`$ from the sum. It is important to note that terms giving rise to secular contribution due to the parametric resonance vanish identically in the first order. The solution for $`R_1`$ can be written explicitly as

$`R_1(x,\xi ;\tau )={\displaystyle \frac{mK^2}{\sigma \omega }}\tau \left({\displaystyle \frac{A}{\xi }}e^{iz^{(+)}}{\displaystyle \frac{B}{\xi }}e^{iz^{()}}\right){\displaystyle \frac{im\sigma ^2}{2\pi }}\xi \lambda {\displaystyle \underset{n0}{}^{}}b_ne^{in\tau /\nu }\left(Ae^{iz^{(+)}}Be^{iz^{()}}\right)`$
$$+\frac{4\pi }{3\sigma }\left(4\omega ^21\right)\left(A^2e^{2iz^{(+)}}+B^2e^{2iz^{()}}\right)+\frac{8\pi }{\sigma \left(4\omega ^23\right)}AB^{}e^{i\left(z^{(+)}z^{()}\right)}+c.c..$$
(4.12)
Straightforward calculations yield the solution for $`u_1`$

$`u_1(x,\xi ;\tau )={\displaystyle \frac{2\pi i}{m^2\sigma ^3\omega }}\left({\displaystyle \frac{A}{\xi }}e^{iz^{(+)}}+{\displaystyle \frac{B}{\xi }}e^{iz^{()}}\right){\displaystyle \frac{2\pi K^2}{\sigma ^3}}\tau \left({\displaystyle \frac{A}{\xi }}e^{iz^{(+)}}+{\displaystyle \frac{B}{\xi }}e^{iz^{()}}\right)`$
$`{\displaystyle \frac{i\lambda n_0}{\nu }}\xi \left(b_{n_0}e^{in_0\tau /\nu }b_{n_0}^{}e^{in_0\tau /\nu }\right)\left(Ae^{iz^{(+)}}+Be^{iz^{()}}\right)`$
$$+i\omega \xi \lambda \underset{n0}{}^{}b_ne^{in\tau /\nu }\left(Ae^{iz^{(+)}}+Be^{iz^{()}}\right)\frac{8\pi ^2}{3m\sigma ^3}\omega \left(4\omega ^25/2\right)\left(A^2e^{2iz^{(+)}}B^2e^{2iz^{()}}\right)+c.c..$$
(4.13)
In second order, we take into account terms providing secular contribution to the solution only. Thus, we write the equation for $`R_2`$ as

$`{\displaystyle \frac{^2R_2}{\tau ^2}}{\displaystyle \frac{K^2}{\sigma ^2}}{\displaystyle \frac{^2R_2}{x^2}}+R_2={\displaystyle \frac{2iK^4m^2}{\omega \sigma ^2}}\tau \left({\displaystyle \frac{^2A}{\xi ^2}}e^{iz^{(+)}}+{\displaystyle \frac{^2B}{\xi ^2}}e^{iz^{()}}\right)`$
$`+{\displaystyle \frac{\lambda \sigma }{\pi }}b_{n_0}\left(1\omega ^2\right)\left[\left(\xi {\displaystyle \frac{}{\xi }}1\right)B^{}e^{iz^{(+)}}+\left(\xi {\displaystyle \frac{}{\xi }}1\right)A^{}e^{iz^{()}}\right]`$
$`+{\displaystyle \frac{K^2}{\sigma ^2}}\left({\displaystyle \frac{^2A}{\xi ^2}}e^{iz^{(+)}}+{\displaystyle \frac{^2B}{\xi ^2}}e^{iz^{()}}\right)+{\displaystyle \frac{2\lambda ^2\sigma ^2}{\pi ^2}}\omega ^2m^2\sigma ^2\left|b_{n_0}\right|^2\xi ^2\left(Ae^{iz^{(+)}}+Be^{iz^{()}}\right)`$
$`{\displaystyle \frac{8\pi ^2}{3\sigma ^2}}\left(16\omega ^411\omega ^2+1\right)\left(\left|A\right|^2Ae^{iz^{(+)}}+\left|B\right|^2Be^{iz^{()}}\right)`$
$$+\frac{16\pi ^2}{\sigma ^2\left(4\omega ^23\right)}\left(\left|B\right|^2Ae^{iz^{(+)}}+\left|A\right|^2Be^{iz^{()}}\right)+c.c..$$
(4.14)
Similar to equation (3.45), we obtain in a straightforward manner

$`R_2(x,\xi ;\tau )={\displaystyle \frac{K^4m^2}{2\omega ^2\sigma ^2}}\tau ^2\left({\displaystyle \frac{^2A}{\xi ^2}}e^{iz^{(+)}}+{\displaystyle \frac{^2B}{\xi ^2}}e^{iz^{()}}\right)+{\displaystyle \frac{K^2\tau }{2i\omega ^3\sigma ^2}}\left({\displaystyle \frac{^2A}{\xi ^2}}e^{iz^{(+)}}+{\displaystyle \frac{^2B}{\xi ^2}}e^{iz^{()}}\right)`$
$`+{\displaystyle \frac{\lambda \sigma \tau }{2\pi i\omega }}b_{n_0}\left(1\omega ^2\right)\left[\left(\xi {\displaystyle \frac{}{\xi }}1\right)B^{}e^{iz^{(+)}}+\left(\xi {\displaystyle \frac{}{\xi }}1\right)A^{}e^{iz^{()}}\right]`$
$`+{\displaystyle \frac{2\lambda ^2\sigma ^2\tau }{2i\pi ^2\omega }}\omega ^2m^2\sigma ^2\left|b_{n_0}\right|^2\xi ^2\left(Ae^{iz^{(+)}}+Be^{iz^{()}}\right)`$
$`{\displaystyle \frac{\tau }{2i\omega }}{\displaystyle \frac{8\pi ^2}{3\sigma ^2}}\left(16\omega ^411\omega ^2+1\right)\left(\left|A\right|^2Ae^{iz^{(+)}}+\left|B\right|^2Be^{iz^{()}}\right)`$
$$+\frac{\tau }{2i\omega }\frac{16\pi ^2}{\sigma ^2\left(4\omega ^23\right)}\left(\left|B\right|^2Ae^{iz^{(+)}}+\left|A\right|^2Be^{iz^{()}}\right)+c.c..$$
(4.15)
Collecting once again the zeroth-order term together with all secular terms in higher orders and renormalizing the amplitudes $`A`$ and $`B`$, we obtain the RG equations

$`2i\omega {\displaystyle \frac{A}{\tau }}2im{\displaystyle \frac{K^2}{\sigma }}{\displaystyle \frac{A}{x}}={\displaystyle \frac{K^2}{\sigma ^2\omega ^2}}{\displaystyle \frac{^2A}{x^2}}+{\displaystyle \frac{\lambda \sigma }{\pi }}b_{n_0}\left(1\omega ^2\right)\left(x{\displaystyle \frac{}{x}}1\right)B^{}`$
$$+\frac{2\lambda ^2\sigma ^4}{\pi ^2}m^2\omega ^2\left|b_{n_0}\right|^2x^2A\frac{8\pi ^2}{3\sigma ^2}\left(16\omega ^411\omega ^2+1\right)\left|A\right|^2A+\frac{16\pi ^2}{\sigma ^2\left(4\omega ^23\right)}\left|B\right|^2A,$$
(4.16)
$`2i\omega {\displaystyle \frac{B}{\tau }}+2im{\displaystyle \frac{K^2}{\sigma }}{\displaystyle \frac{B}{x}}={\displaystyle \frac{K^2}{\sigma ^2\omega ^2}}{\displaystyle \frac{^2B}{x^2}}+{\displaystyle \frac{\lambda \sigma }{\pi }}b_{n_0}\left(1\omega ^2\right)\left(x{\displaystyle \frac{}{x}}1\right)A^{}`$
$$+\frac{2\lambda ^2\sigma ^4}{\pi ^2}m^2\omega ^2\left|b_{n_0}\right|^2x^2B\frac{8\pi ^2}{3\sigma ^2}\left(16\omega ^411\omega ^2+1\right)\left|B\right|^2B+\frac{16\pi ^2}{\sigma ^2\left(4\omega ^23\right)}\left|A\right|^2B.$$
(4.17)
## 5 The Nonlinear Schrodinger Equation for a Single Mode

Equations (3.46) and (3.47) represent a system of coupled nonlinear Schrodinger equations for the mode amplitudes $`A_m`$ and $`B_m`$. Neglecting the contribution from modes with $`n\pm m`$ and introducing a new amplitude $`\stackrel{~}{B}_m(x;\tau )=B_m(x;\tau )`$, for the single mode amplitudes $`A_m`$ and $`B_m`$, we obtain the equations

$$i\frac{A}{\tau }\frac{^2A}{\zeta ^2}\frac{8\pi ^2}{\omega \sigma ^2\left(4\omega ^23\right)}\left(G\left|A\right|^2+\left|B\right|^2\right)A=0,$$
(5.1)
$$i\frac{B}{\tau }\frac{^2B}{\zeta ^2}\frac{8\pi ^2}{\omega \sigma ^2\left(4\omega ^23\right)}\left(\left|A\right|^2G\left|B\right|^2\right)B=0,$$
(5.2)
where

$$\zeta =\sqrt{2\omega }\left(mK\tau +\frac{\omega \sigma }{K}\right),G=\frac{4\omega ^23}{6}\left(16\omega ^411\omega ^2+1\right),$$
(5.3)
the index $`m`$ and the tilde sign over the new amplitude $`\stackrel{~}{B}`$ has been omitted. Clearly enough, equations (5.1) and (5.2) follow directly from the system (4.16) - (4.17) if $`b_{n_0}`$ is set to zero. This implies that in the case where a parametric resonance does not occur in the original system, the smooth approximation is valid to second order in the formal perturbation parameter. Moreover, it is straightforward to verify that $`G`$ is always positive. The system of two coupled nonlinear Schrodinger equations (5.1) and (5.2) is in general non integrable. Its integrability is proven by Manakov  only in the simplest case of $`G=1`$.

If one of the amplitudes ($`B`$ or $`A`$) in its capacity of being a particular solution to the corresponding nonlinear Schrodinger equation is identically zero, the equation for the other amplitude (say $`A`$) becomes

$$i\frac{A}{\tau }\frac{^2A}{\zeta ^2}+\mathrm{\Gamma }\left|A\right|^2A=0,$$
(5.4)
where

$$\mathrm{\Gamma }=\frac{4\pi ^2}{3\omega \sigma ^2}\left(16\omega ^411\omega ^2+1\right).$$
(5.5)
In nonlinear optics equation (5.4) is known to describe the formation and evolution of the so-called dark solitons . In the case of charged particle beams these correspond to the formation of holes or cavitons in the beam. The solution to equation (5.4) can be written as

$$A(\zeta ;\tau )=\frac{r\left(\xi \right)}{\sqrt{\mathrm{\Gamma }}}e^{i\left[n\tau +\theta \left(\xi \right)\right]},\xi =\zeta c\tau ,$$
(5.6)
where

$$r\left(\xi \right)=\sqrt{n2a^2\mathrm{sech}^2\left(a\xi \right)},\theta \left(\xi \right)=\mathrm{arctan}\left[\frac{2a}{c}\mathrm{tanh}\left(a\xi \right)\right],$$
(5.7)
for all $`c`$ and

$$a=\frac{1}{2}\sqrt{2nc^2},$$
(5.8)
provided $`n>c^2/2`$.

Let us now examine equations (4.16) and (4.17). In the cold-beam limit $`v_T0`$ (or equivalently $`\omega 1`$) the second term on their right-hand-sides can be neglected as compared to the other terms. Therefore, we can write

$`2i\omega {\displaystyle \frac{\stackrel{~}{A}}{\tau }}={\displaystyle \frac{K^2}{\sigma ^2\omega ^2}}{\displaystyle \frac{^2\stackrel{~}{A}}{x^2}}+{\displaystyle \frac{2\lambda ^2\sigma ^4}{\pi ^2}}m^2\omega ^2\left|b_{n_0}\right|^2x^2\stackrel{~}{A}`$
$$\frac{8\pi ^2}{3\sigma ^2}\left(16\omega ^411\omega ^2+1\right)\left|\stackrel{~}{A}\right|^2\stackrel{~}{A}+\frac{16\pi ^2}{\sigma ^2\left(4\omega ^23\right)}\left|\stackrel{~}{B}\right|^2\stackrel{~}{A},$$
(5.9)
$`2i\omega {\displaystyle \frac{\stackrel{~}{B}}{\tau }}={\displaystyle \frac{K^2}{\sigma ^2\omega ^2}}{\displaystyle \frac{^2\stackrel{~}{B}}{x^2}}+{\displaystyle \frac{2\lambda ^2\sigma ^4}{\pi ^2}}m^2\omega ^2\left|b_{n_0}\right|^2x^2\stackrel{~}{B}`$
$$\frac{8\pi ^2}{3\sigma ^2}\left(16\omega ^411\omega ^2+1\right)\left|\stackrel{~}{B}\right|^2\stackrel{~}{B}+\frac{16\pi ^2}{\sigma ^2\left(4\omega ^23\right)}\left|\stackrel{~}{A}\right|^2\stackrel{~}{B},$$
(5.10)
where the local gauge transformation

$$\stackrel{~}{A}=A\mathrm{exp}\left[im\omega \left(\frac{mK^2}{2}\tau +\sigma \omega x\right)\right],\stackrel{~}{B}=B\mathrm{exp}\left[im\omega \left(\frac{mK^2}{2}\tau \sigma \omega x\right)\right],$$
(5.11)
has been performed on the amplitudes $`A`$ and $`B`$. Equations (5.9) and (5.10) represent a system of two coupled Gross-Pitaevskii equations, known to govern the formation of condensate structures in atomic gases confined in magnetic traps. In addition, superfluidity in helium as well as patterns in the gas of paraexcitons in semiconductors are considered as a possible manifestation of Bose-Einstein condensation. It is remarkable that under certain conditions similar phenomenon can be observed in space-charge dominated beams.

## 6 Concluding Remarks

The analysis performed in the present paper is based on the Vlasov-Maxwell equations for the self-consistent evolution of the beam distribution function and the electromagnetic fields. We considered the propagation of an intense beam through a periodic focusing lattice. It has been proven that a special class of solutions to the Vlasov-Maxwell equations exists, which defines a uniform phase-space density inside of simply connected boundary curves. These solutions provide an exact closure to the hierarchy of moments resulting in a set of hydrodynamic equations for the beam density (continuity equation) and for the current velocity. The latter is characterized by a triple adiabatic pressure law. The major drawback of the model discussed here is that due to the constancy of the distribution function inside of the evolving region of phase space, it fails to take into account an important feature such as the well-known Landau damping. Inclusion of the Landau damping in the description can be achieved correctly only by direct analysis of the Vlasov-Maxwell system .

Further, we studied first the case when the smooth focusing approximation applies. Based on the RG method, a system of coupled nonlinear Schrodinger equations has been derived for the slowly varying amplitudes of interacting beam-density waves. Under the approximation of an isolated mode neglecting the effect of the rest of the modes, this system reduces to two coupled nonlinear Schrodinger equations for the amplitudes of the forward and of the backward wave. The particular solution to the latter system asserting that the amplitude of either of the waves (the forward or the backward) can vanish identically, leads to a single nonlinear Schrodinger equation for the other amplitude. The latter is characterized by a repulsive nonlinearity and therefore, it describes the evolution of holes in intense charged particle beams.

The analysis of periodic focusing clearly showed that the results obtained in the case of smooth focusing remain unchanged up to second order in the perturbation parameter if a parametric resonance between a particular mode and an appropriate Fourier harmonic of the external focusing does not occur. If however, an exact parametric resonance takes place, it was shown that the evolution of the wave amplitudes of the resonant mode in the cold-beam limit is described by a system of coupled Gross-Pitaevskii equations. Quite remarkably, it was found that there exist a possibility of formation of density condensates in space-charge dominated beams.

### Acknowledgments

It is a pleasure to thank Y. Oono and S.-I. Goto for careful reading of the manuscript and for making valuable suggestions. I am also greatly indebted to B. Baizakov and F. Illuminati for many helpful discussion concerning the subject of the present paper and for their valuable comments on its early version. This research was supported by the INFM grant number OA03002735.

## Appendix

In order to prove that the expression (2.14) represents the solution of the Vlasov equation (2.12), we write the equations for the characteristics of the latter in the form

$$\frac{\mathrm{d}\theta }{1}=\frac{\mathrm{d}x}{\dot{\chi }p}=\mathrm{d}p\left(\dot{\chi }x+\frac{V}{x}+\lambda \sqrt{\beta }\frac{U}{x}\right)^1=\frac{\mathrm{d}f}{0}.$$
(A.1)
Let us also assume that the distribution function $`f(x,p;\theta _0)`$ at some initial time $`\theta _0`$ is given by

$$f(x,p;\theta _0)=f_0(x,p)=𝒞\left[\left(pp_{()}^{(0)}\left(x\right)\right)\left(pp_{(+)}^{(0)}\left(x\right)\right)\right].$$
(A.2)
Suppose now that we have been able to solve the equations for the characteristics (A.1), and we have expressed the solution subsequently according to

$$p=P(x;\theta ),$$
(A.3)
where $`P(x;\theta )`$ is an appropriate function of its arguments. At each instant of time $`\theta `$ equation (A.3) defines two curves $`p_{(+)}(x;\theta )`$ and $`p_{()}(x;\theta )`$ in the phase space $`(x,p)`$, such that

$$p_{(+)}(x;\theta _0)=p_{(+)}^{(0)}\left(x\right),p_{()}(x;\theta _0)=p_{()}^{(0)}\left(x\right).$$
(A.4)
Therefore, if initially the distribution function is given by equation (A.2), then the solution of the Vlasov equation (2.12) at every subsequent instant of time $`\theta `$ is represented by the expression (2.14).

It is interesting to note that the chain of equations (excluding the last equation) for the characteristics (A.1) of the Vlasov equation formally coincide with those for the equation

$$\frac{P}{\theta }+\dot{\chi }P\frac{P}{x}=\dot{\chi }x\frac{V}{x}\lambda \sqrt{\beta }\frac{U}{x}.$$
(A.5)
As a matter of fact, equation (A.5) is equivalent to the two equations

$$\frac{p_{(\pm )}}{\theta }+\dot{\chi }p_{(\pm )}\frac{p_{(\pm )}}{x}=\dot{\chi }x\frac{V}{x}\lambda \sqrt{\beta }\frac{U}{x},$$
(A.6)
with the initial conditions (A.4), provided $`U`$ is defined as a solution to the equation (2.18).


format: text

text:
# On computational irreducibility and the predictability of complex physical systems

## Abstract

Using elementary cellular automata (CA) as an example, we show how to coarse-grain CA in all classes of Wolfram’s classification. We find that computationally irreducible (CIR) physical processes can be predictable and even computationally reducible at a coarse-grained level of description. The resulting coarse-grained CA which we construct emulate the large-scale behavior of the original systems without accounting for small-scale details. At least one of the CA that can be coarse-grained is irreducible and known to be a universal Turing machine.

Can one predict the future evolution of a physical process which is described or modeled by a computationally irreducible (CIR) mathematical algorithm? For such systems, in order to know the system’s state after (e.g.) one million time steps, there is no faster algorithm than to solve the equation of motion a million time steps into the future. Wolfram has suggested that the existence of CIR systems in nature is at the root of our apparent inability to model and understand complex systems Wolfram (1984a, 1985, 2002); Ilachinski (2001).

Complex physical systems that are CIR might therefore seem to be inherently unpredictable. It is tempting to conclude from this that the enterprise of physics itself is doomed from the outset; rather than attempting to construct solvable mathematical models of physical processes, computational models should be built, explored and empirically analyzed. This argument, however, assumes that infinite precision is required for the prediction of future evolution. Usually coarse-grained or even statistical information is sufficient: indeed, a physical model is usually correct only to a certain level of resolution, so that there is little interest in predictions from such a model on a scale outside its regime of validity.

In this Letter, we report on experiments with nearest neighbour one-dimensional cellular automata, which show that because in practice one only seeks coarse-grained information, complex physical systems can be predictable and even computationally reducible at some level of description. The implication of these results is that, at least for systems whose complexity is the outcome of very simple rules, useful approximations can be made that enable predictions about future behavior.

Cellular automata (CA) are dynamical systems composed of a lattice of cells. Each cell in the lattice can assume a value from a given finite alphabet. The system evolves in time according to an update rule that gives a cell’s new state as a function of values in its finite neighborhood. CA were originally introduced by von Neumann and Ulam von Neumann (1966) in the 1940’s as a possible way of simulating self reproduction in biological systems. Since then, CA have attracted a great deal of interest in physics Wolfram (1983); phy ; Ilachinski (2001); Wolfram (1994) because they capture two basic ingredients of many physical systems: 1) they evolve according to a local uniform rule. 2) CA can exhibit rich behavior even with very simple update rules. For similar and other reasons, CA have also attracted attention in computer science Mitchel (1998); Sarkar (2000), biology Ermentrout and Edelstein-Keshet (1993), material science Raabe (2002) and other fields.

In early work Wolfram (1984b, a, 1985, 2002), Wolfram proposed that CA can be grouped into four classes of complexity. Class 1 consists of CA whose dynamics reaches a steady state regardless of the initial conditions. Class 2 consists of CA whose long time evolution produces periodic or nested structures. CA from both of these classes are simple in the sense that their long time evolution can be deduced from running the system a small number of time steps. On the other hand, class 3 and class 4 consist of “complex” CA. Class 3 CA produce structures that seem random. Class 4 CA produce a mixture of random structures and periodic behavior. This transition in CA complexity was later regarded as a phase transitionLangton (1990) in the CA rule space. For a review on CA classification see Refs. Mitchel, 1998; Ilachinski, 2001; Wolfram, 2002. There is no generally agreed upon algorithm for classifying a given CA. The assignment of CA to these four classes is somewhat subjective and, we will argue, may need to be refined. Based on numerical experiments, Wolfram hypothesized Wolfram (1984b, a, 2002) that most CA from class 3 and 4 are CIR.

There is no unique way to define coarse-graining, but here we will mean that our information about the CA is locally coarse-grained in the sense of being stroboscopic in time, but that nearby cells are grouped into a supercell according to some specified rule (as is frequently done in statistical physics). A system which can be coarse-grained is compactable since it is possible to calculate its future time evolution (or some coarse aspects of it) using a more compact algorithm than its native description. Note that our use of the term compactable refers to the phase space reduction associated with coarse-graining, and is agnostic as to whether or not the coarse-grained system is reducible or irreducible. Accordingly, we define predictable to mean that a system is computationally reducible or has a computationally reducible coarse-graining. Thus, it is possible to calculate the future time evolution of a predictable system (or some coarse aspects of it) using an algorithm which is more compact than both the native and coarse-grained descriptions.

In order to quantify the implications of looking at coarse-grained information only, we have systematically attempted to coarse-grain the 256 nearest neighbor one-dimensional binary CA that were the subject of Wolfram’s investigationsWolfram (2002, 1983). The outcome, described in detail below, is surprising: we found that many CA can be coarse-grained and that in some cases CIR CA are coarse-grained by computationally reducible ones. In other words, even though microscopically a given system might be CIR, its coarse-grained dynamics can be compactable and predictable.

We start by defining a simple procedure for coarse-graining a CA. Other constructions are undoubtedly possible. For simplicity we limit our treatment to one-dimensional systems with nearest neighbor interactions. Generalizations to higher dimensions and different interaction radii are straightforward. Let $`A=(a\left(t\right),S_A,f_A)`$ be a cellular automaton defined on an array of cells $`a\left(t\right)=\left\{a_n\left(t\right)\right\}|_{n=\mathrm{}}^{\mathrm{}}`$. Each cell accepts an alphabet of $`S_A`$ symbols, namely $`a_n\left(t\right)\left\{0\mathrm{}S_A1\right\}`$. The values of the cells evolve in time according to the update rule $`a_n\left(t+1\right)=f_A[a_{n1}\left(t\right),a_n\left(t\right),a_{n+1}\left(t\right)]`$, where $`f_A:\left\{S_A\right\}^3\left\{S_A\right\}`$ is the transition function. The update rule is applied simultaneously to all the cells and we denote this application by $`a\left(t+1\right)=f_Aa\left(t\right)`$.

Our goal is to find a modified CA $`B=(b\left(t\right),S_B,f_B)`$ and an irreversible coarse-graining function $`b=C\left(a\right)`$, which are capable of a coarse-grained emulation of $`A`$. For every initial condition $`a(0)`$, $`B`$ and $`C`$ must satisfy:

$$C\left(f_A^{Tt}a(0)\right)=f_B^tC\left(a(0)\right).$$
(1)
Namely, running the original CA for $`Tt`$ time steps and then coarse-graining is equivalent to coarse-graining the initial condition and then running the modified CA $`t`$ time steps. The constant $`T`$ is a time scale associated with the coarse-graining.

To search for explicit coarse-graining rules, we define the $`N`$’th block version $`A^N=(a^N,S_{A^N},f_{A^N})`$ of $`A`$. $`S_{A^N}=\left(S_A\right)^N`$ and each cell in $`A^N`$ represents a block of $`N`$ cells in $`A`$. Cell values are translated between $`A`$ and $`A^N`$ according to the base $`S_A`$ value of $`N`$ cells in $`A`$. The transition function $`f_{A^N}:\left\{S_{A^N}\right\}^3\left\{S_{A^N}\right\}`$ is computed by running $`A`$ for $`N`$ time steps on all possible initial conditions of length $`3N`$. In this way $`A^N`$ computes in one time step $`N`$ time steps of $`A`$. Note that $`A^N`$ is not a coarse-graining of $`A`$ because no information was lost in the cell translation.

Next we attempt to generate the coarse CA $`B`$ by projecting the alphabet of $`A^N`$ on a subset of $`\left\{0\mathrm{}S_{A^N}1\right\}`$. This is the key step where information is being lost, a manipulation which distinguishes between coarse-graining and emulation blocking transformations Wolfram (2002); Ilachinski (2001). The transition function $`f_B`$ is constructed from $`f_{A^N}`$ by projecting its arguments and outcome:

$$f_B[\overline{x_1},\overline{x_2},\overline{x_3}]=\overline{f_{A^N}[\overline{x_1},\overline{x_2},\overline{x_3}]}.$$
(2)
Here $`\overline{x}P\left(x\right)`$ denotes the projection operation. This construction is possible only if

$$\overline{f_{A^N}[x_1,x_2,x_3]}=\overline{f_{A^N}[y_1,y_2,y_3]},\left(x,y|\overline{x_i}=\overline{y_i}\right).$$
(3)
Otherwise, $`f_B`$ is multi-valued and our coarse-graining attempt fails for the specific choice of $`N`$ and $`P`$.

In cases where the above conditions are satisfied, the resulting CA $`B`$ is a coarse-graining of $`A^N`$ with a time scale $`T=1`$. For every step $`a_n^N(t+1)=f_{A^N}[a_{n1}^N(t),a_n^N(t),a_{n+1}^N(t)]`$ of $`A^N`$, $`B`$ makes the move

$`b_n(t+1)`$ $`=`$ $`f_B[b_{n1}(t),b_n(t),b_{n+1}(t)]`$
$`=`$ $`\overline{f_{A^N}[\overline{a_{n1}^N(t)},\overline{a_n^N(t)},\overline{a_{n+1}^N(t)}]}=\overline{a_n^N(t+1)},`$
where we have used Eq. (3) in the last step. $`B`$ therefore satisfies Eq. (1) with $`P`$ as the coarse-graining function. Since a single time step of $`A^N`$ computes $`N`$ time steps of $`A`$, $`B`$ is also a coarse-graining of $`A`$ with a coarse-grained time scale $`T=N`$. The coarse-graining function $`C`$ in this case is composed of the translation from $`A`$ to $`A^N`$ followed by the projection operator $`P`$. Analogies of these operators have been used in attempts to reduce the computational complexity of certain stochastic partial differential equations Hou et al. (2001); Degenhard and Rodriguez-Laguna (2002). Similar ideas have been used to calculate critical exponents in probabilistic CA de Oliveira and Satulovsky (1997); Monetti and Satulovsky (1998).

It is interesting to notice that the above coarse-graining procedure can lose two very different types of dynamic information. To see this, consider Eq. (3). This equation can be satisfied in two ways. In the first case

$$f_{A^N}[x_1,x_2,x_3]=f_{A^N}[y_1,y_2,y_3],\left(x,y|\overline{x_i}=\overline{y_i}\right),$$
(5)
which necessarily leads to Eq. (3). $`f_{A^N}`$ in this case is insensitive to the projection of its arguments. The distinction between two variables which are identical under projection is therefore irrelevant to the dynamics of $`A^N`$, and by construction to the long time dynamics of $`A`$. By eliminating irrelevant degrees of freedom (DOF), coarse-graining of this type removes information which is redundant on the microscopic scale. The coarse CA in this case accounts for all possible long time trajectories of the original CA and the complexity classification of the two CA is therefore the same.

In the second case Eq. (3) is satisfied even though Eq. (5) is violated. Here the distinction between two variables which are identical under projection is relevant to the dynamics of $`A`$. Replacing $`x`$ by $`y`$ in the initial condition may give rise to a difference in the dynamics of $`A`$. Moreover, the difference can be (and in many occasions is) unbounded in space and time. Coarse-graining in this case is possible because the difference is constrained in the symbol space by the projection operator. Namely, projection of all such different dynamics results in the same coarse-grained behavior. Note that the coarse CA in this case cannot account for all possible long time trajectories of the original one. It is therefore possible for the original and coarse CA to fall into different complexity classifications.

Coarse-graining by elimination of relevant DOF removes information which is not redundant with respect to the original system. The information becomes redundant only when moving to the coarse scale. In fact, “redundant” becomes a subjective qualifier here since it depends on our choice of coarse description. In other words, it depends on what aspects of the microscopic dynamics we want the coarse CA to capture. In a sense, this is analogous to the subtleties encountered in constructing renormalization group transformations for the critical behavior of antiferromagnets Goldenfeld (1992); van Leeuwen (1975).

We now give specific examples of coarse-graining CA. In the sequel, CA rules are numbered using Wolfram’s notationWolfram (1983, 2002). Figure 1 (a) and (b) shows a coarse-graining of rule 146 by rule 128. Rule 146 produces a complex, seemingly random behavior which falls into the class 3 group. We use a super-cell size $`N=3`$, and project the alphabet $`\{0\mathrm{}7\}`$ of the super-cells back to the $`\{0,1\}`$ alphabet with $`P(7)=1`$ and $`P(x7)=0`$pro . A triplet of cells in rule 146 are therefore coarse-grained to a single cell and the value of the coarse cell is $`1`$ only when the triplet are all $`1`$. Using this projection operator we construct the transition function of the coarse CA which is found to be rule 128, a class 1 CA. This choice of coarse-graining eliminates the small scale details of rule 146. Only structures of lateral size of three or more cells are accounted for. The decay of such structure in rule 146 is accurately described by rule 128.

Note that a class 3 CA was coarse-grained to a class 1 CA in the above example. Our gain was therefore two-fold. In addition to the phase space reduction associated with coarse-graining we have also achieved a reduction in complexity. Our procedure found predictable coarse-grained aspects of the dynamics even though the small scale behavior of rule 146 is complex, potentially CIR. As we explained earlier, this type of simplification can be achieved only by eliminating relevant DOF.

As a second example we show a transition between rules with a comparable complexity. Fig. 1 (c) and (d) shows a coarse-graining of rule 105 by rule 150. $`N=2`$ in this example and $`P(x)=1`$ only when $`x=0,3`$pro .

The coarse-graining procedure we described above is not constructive, but instead is a self-consistency condition on a putative coarse-graining rule with a specific block size $`N`$ and projection operator $`P`$. In many cases the coarse-graining fails and one must try other choices of $`N`$ and $`P`$. Can all CA be coarse-grained? If not, which CA can be coarse-grained and which cannot?

To answer these questions we tried systematically to coarse-grain Wolfram’s 256 elementary rules. We applied the coarse-graining procedure to each rule and scaned the $`N`$,$`P`$ space for valid solutions. In this way we were able to coarse-grain 240 out of the 256 CA dif . These 240 coarsen-able rules include members of all four classes. Many elementary CA can be coarse-grained by other elementary CA. Figure 2 shows a map of the coarse-graining transitions that we found within the family of elementary rules. Only coarse-grainings with $`N4`$ are shown due to limited computing power. Other transitions may exist with larger $`N`$. We observe that rule complexity never increases along the map’s transitions, i.e., coarse-graining introduces (at least here) partial order among CA rules.

As mentioned above, we were unable to coarse-grain 16 elementary rules. 12 out of the 16 are the class 3 rules 30,45,106 and their symmetries. The other four are the class 2 rule 154 and it’s symmetries. We don’t know if our inability to coarse-grain these 16 rules comes from limited computing power or from something deeper.

It is worth noticing that a subset of the fixed points in the transition map is composed of all elementary additive rules (see page 952 in Ref. Wolfram, 2002) and their symmetries. This result is not limited to elementary rules. All additive CA whose alphabet sizes $`S`$ are prime numbers coarse-grain themselves with $`N=S`$ and $`P`$ the modulo $`S`$ sum of $`S`$ cells. We conjecture that there are situations where reducible fixed points exist for a wide range of systems, analogous to the emergence of amplitude equations in the vicinity of bifurcations points.

Coarse-graining transitions can also exit the elementary CA family. This happens whenever the alphabet of the coarse CA consists of more than two symbols. One such example which is of special importance is rule 110. Rule 110 is interesting because this class 4 CA is universal Wolfram (2002) in the Turing sense Herken (1995) and is therefore CIR Wolfram (1984b). It is capable of emulating all computations done by other computing devices in general and CA in particular.

We found several ways to coarse-grain rule 110. Using $`N=6`$, it is possible to project the 64 possible states onto an alphabet of 63 symbols. A more impressive reduction in the alphabet size is obtained by going to larger values of $`N`$. For $`N=7,8,9,10,11`$ we found an alphabet reduction of $`6/128`$, $`22/256`$, $`67/512`$, $`181/1024`$ and $`463/2048`$ respectively. We expect this behavior to persist for larger values of $`N`$.

Another interesting coarse-graining of rule 110 that we found is the transition to rule 0. Rule 0 has the trivial dynamics where all initial states evolve to the null configuration in a single time step. The transition to rule 0 is possible because the cell combination “01010” is not generated by rule 110 and can only appear in the initial state. Coarse-graining by rule 0 is achieved using $`N=5`$ and projecting “01010” to 1 and all other five cell combinations to 0. This example is important because it shows that even though rule 110 is CIR it has a predictable coarse-grained dynamics (however trivial). To our knowledge rule 110 is the only proven CIR elementary CA and therefore this is the only example of irreducible to reducible transition between elementary rules that we found.

We did find other complex, potentially CIR rules that can be coarse-grained by reducible CA. Rules 18, 54, 126 and their symmetries are coarse-grained by rule 0. As we showed above, rule 146 and it’s symmetries can be coarse-grained by rule 128 in a non-trivial way. We don’t know if these rules are CIR for lack of proof. Nevertheless, non-trivial irreducible to reducible transitions can in principle exist. Consider for example the CIR CA generated by the product of rules 110 and 128. Let the CA alphabet be $`\{a,b\}\times \{0,1\}`$, where the letters evolve according to rule 110 and the digits according to 128. We can recover the reducible, coarse-grain-able rule 128 by projecting the alphabet onto $`\{0,1\}`$.

The fact that CIR rules can be coarse-grained and that they has predictable coarse-grained dynamics shows that CIR is not a good measure of physical complexity. As in the case of rule 110, a CIR system may still yield an efficient predictable theory, provided that we are willing to ask coarse-grained questions. It seems that a better classification of physical complexity is related to what classes of projection operator are required to coarse-grain the system: local, real space projections or more complex non-geometric projections?

In summary, we have found that many CA, including CIR ones can be locally coarse-grained in space and time. In some cases CIR systems are predictable, if coarse-grained information only is required.

NG wishes to thank Stephen Wolfram for useful discussions and correspondence. This work was partially supported by the National Science Foundation through grant NSF-DMR-99-70690 (NG) and by the National Aeronautics and Space Administration through grant NAG8-1657.


text:
# The 3-Dimensional Structure of HH 32 from GMOS IFU SpectroscopyBased on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Particle Physics and Astronomy Research Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), CNPq (Brazil), and CONICET (Argentina).

## 1 Introduction

HH 32 is one of the brightest sources in the original catalog of Herbig-Haro objects (Herbig 1974). Spectroscopic studies have shown that HH 32 has a high excitation spectrum with significant emission from H, He and a wide range of metals in the 3700 to 10800 Å region (Dopita 1978; Brugel, Böhm & Mannery 1981; Raga, Böhm & Cantó 1996). HH 32 also has an unusually large extinction (E(B-V)$`0.7`$; Brugel, Böhm & Mannery 1981), hence it is surprising that this object was detected in the ultraviolet with IUE (Böhm & Böhm-Vitense 1984; Lee et al. 1988; Moro-Martín et al. 1996). The UV detection results from of the strong emission in the high exitation \[C III\] 1909 Å and \[C IV\] 1550 Å ultraviolet lines. Interestingly, this high excitation object also shows H2 v = 1-0 S(1) emission in the infrared which is excited in shocks as the flow encounters ambient cloud material. (Zealey et al. 1986; Zinnecker et al. 1989; Davis, Eislöffel & Smith 1996).

The HH 32 outflow is clearly associated with the T Tauri star AS 353A, with condensations A, B and D at $`20^{\prime \prime }`$ to the west, and condensation C at $`10^{\prime \prime }`$ to the east of the star (see Figure 1 from Curiel et al. 1997). AS 353A has a rich emission line spectrum (Eislöffel, Solf & Böhm 1990), but shows only two forbidden emission lines: \[O I\] 5577, 6300 Å. High resolution optical spectroscopy of the HH condensations has shown that HH 32 A, B and D have very broad, redshifted line profiles, with full widths of $`400`$ km s-1 (Herbig & Jones 1983; Hartigan, Mundt & Stocke 1986; Solf, Böhm & Raga 1986; Hartigan, Raymond & Hartmann 1987). The faint knot C to the east of the exciting source is blueshifted, indicating that the HH 32 emission knots are part of a bipolar outflow from AS 353A. The very large line widths, the shape of the line profiles, and the position-velocity diagrams have been interpreted successfully in terms of “3/2-D” bow shock models (Solf et al. 1986; Raga, Böhm & Solf 1986; Raga & Böhm 1986; Hartigan et al. 1987). In such models, the bow shock is modeled with a parametrized shape (i.e., Hartigan et al. 1987) and the emission is computed by assuming that the post-bow shock recombination region resembles the structure of a 1D stationary shock. Comparisons with the models show a $`v_{bs}300350`$ km s-1 velocity for the bow shock (Solf et al. 1986; Hartigan et al. 1987).

By comparing the proper motions with the radial velocities we know that HH 32A is moving away from the observer at a $`\varphi 70^{}`$ angle with respect to the plane of the sky (Herbig & Jones 1983). This angle is consistent with the orientation necessary to model the emission line profiles (Solf et al. 1986; Hartigan et al. 1987). The detailed morphological structure and the proper motions of the sub-condensations of HH 32A described by Curiel et al. (1997) have been interpreted in terms of 3/2-D bow shock models to obtain the same $`\varphi =70^{}`$ and $`v_{bs}300`$ km s-1 values that are deduced from the emission line profiles (Raga et al. 1997).

In this paper we present and discuss new spectroscopic observations with 2D spatial coverage of a field centered on the HH 32A condensation, which give new information on the kinematics and structure of the emission knots. The previous spatially resolved spectroscopy of Solf et al. (1986) covered all of HH 32 with several different long-slit positions centered on AS 353A, but did not have a comparable angular resolution nor as high a signal-to-noise ratio as the spectra we present here. Also, the HST images of Curiel et al. (1997) show the high resolution spatial structure of HH 32, but lack any kinematic information. In $`\mathrm{\S }`$2, we present the observations and discuss briefly the method of data reduction. We show the H$`\alpha `$ line profiles and velocity channel maps obtained from the spectra, and the velocity channel maps of different line ratios in $`\mathrm{\S }`$3 and $`\mathrm{\S }`$4. In $`\mathrm{\S }`$5 we interpret the observations in terms of a 3/2-D bow shock model and discuss its successes and failures in modeling the data. Finally, we present our conclusions and a discussion of the implications in $`\mathrm{\S }`$6.

## 2 The observations

The data for this project were obtained on 2002 Aug 03 with the Gemini Multi-Object Spectrograph (GMOS; Hook et al. 2002) Integral Field Unit (IFU; Allington-Smith et al. 2002) at the Gemini North Fredrick C. Gillett Telescope. These data were taken as System Verification for dithered observations on an emission line source, and are available to the public at: http://www.gemini.edu/sciops/instruments/gmos/SVobs/gmosSVobs73.html under Gemini program ID GN-2002B-SV-73. The observations were made in cloudy conditions (0.5 to 1 magnitude of extinction) with $``$0.′′5 average seeing. The source star of the HH 32 outflow, AS 353A, was used to provide guiding and tip-tilt corrections using the GMOS on-instrument wavefront sensor (OIWFS).

These data were obtained with the IFU in 1-slit mode, which provides a spatial field of $`3\stackrel{}{\mathrm{.}}5`$ x 5′′ for each resulting science datacube. With this observing configuration, the GMOS IFU is comprised of 750 fibers; each spans a $`0\stackrel{}{\mathrm{.}}2`$ hexagonal region on the sky. 500 fibers make up the $`3\stackrel{}{\mathrm{.}}5`$ x 5′′ science field of view, and 250 fibers make up a smaller, dedicated sky field which is fixed at a 1 distance from the science position (Allington-Smith et al. 2002). A total of 12 900 second frames were obtained; two steps through a 3 x 2 dither pattern resulting in total integration time of 30 minutes at each of the 6 dither positions. We used the R831 grating in GMOS, which has a 0.034 nm/pixel scale and results in $`15`$ km s-1 spectral resolution. We were able to obtain spectra over the 4920 - 7040 Å region, which includes the \[O III\] 4969/5007, \[O I\] 6300/63, \[N II\] 6548/83, H$`\alpha `$, and \[S II\] 6716/31 lines.

We used the GMOS IFU scripts available in the Gemini IRAF package for reduction of the data. The gfreduce task was used to prepare the frames for reduction, subtract off the bias and overscan levels, reject cosmic rays by calling the gscrrej routine, and call the additional reduction routines gfextract and gfskysub. gfextract traces and extracts each IFU fiber to 1-D spectral output in the flat field image, applies this trace to the science data, extracts the science spectra using a $`\pm `$2.5 pixel aperture on each side of the traces, and applies the flat field correction. gfskysub subtracts off a sky spectrum from the science data using averaged spectra from the dedicated sky field inherent in the IFU. The gswavelength task then establishes the wavelength calibration for the IFU arc lamp images, which is thus applied to the science frames using the gftransform routine. These 2-Dimensional data images were reformatted into a 3-D datacube format (x, y, $`\lambda `$) and resampled to square pixels with a 0.′′05 spatial resolution using the gfcube routine. The dithered data were mosaiced together to form a larger cube in the spatial dimension by manually offsetting and coadding the individual cubes using IDL. The final datacube has a spatial extent of 8′′.7 by 5.′′85. At the forementioned webpage we have made available the raw data, the final IRAF datacubes, the mosaiced cube constructed in IDL, an IRAF cl-reduction script, and a README file that describes our data reduction process.

In Figure 1 we present a 30 second exposure r-band image of the HH 32 system that was obtained for setup on the IFU field; the HH 32 A1, A2, B and D condensations are labeled (using the notation of Curiel et al. 1997). Overplotted on the image is the position and orientation in the spatial dimension of the final mosaiced datacube, and a vector indicating the direction to the outflow source, AS 353A. In Figure 2 we show a spectrum of the HH 32 A2 knot averaged over the inner $`1\stackrel{}{\mathrm{.}}5`$ in the spatial dimension; the \[O I\] 6300/63, \[N II\] 6548/83, H$`\alpha `$ and \[S II\] 6716/31 lines are identified. All of the lines show a pronounced double-peaked structure, consistent with previous observations (Hartigan, Mundt & Stocke 1986).

## 3 H$`\alpha `$ line profiles and channel maps

Figure 3 shows a contour map of the H$`\alpha `$ emission constructed by summing the line over its wavelength extent in the final mosaiced datacube. Overplotted in Figure 3 are the H$`\alpha `$ line profiles obtained by integrating the emission over boxes of $`10\times 10`$ pixels (corresponding to $`0\stackrel{}{\mathrm{.}}5\times 0\stackrel{}{\mathrm{.}}5`$). Many of the observed line profiles are double peaked, with a stronger peak at a $`v_r90`$ km s-1 heliocentric radial velocity and a weaker high velocity peak at $`v_r240276`$ km s-1. Several of the line profiles in the emission joining the maxima of the A1 and A2 condensations show a dominant high velocity peak, and a few of the profiles in the periphery of A1 show three components. In many of the line profiles along the edges of the observed field, only the low velocity component is visible.

In Figure 4 we present the H$`\alpha `$ radial velocity channel maps derived from the IFU observations. The H$`\alpha `$ emission is detected from heliocentric radial velocities ranging from $`v_r20`$ to $`+400`$ km s-1. The surface brightness of the condensations grows with increasing $`v_r`$, reaching a maximum emission level in the $`+72`$ km s-1 map. It then decreases, reaching a minimum at $`+165`$ km s-1, and increases again reaching a second maximum at $`v_r=+241`$ km s-1. For the larger radial velocities, the surface brightness decreases monotonically.

At $`v_r=10`$ km s-1, there are two condensations with a separation of $`2\stackrel{}{\mathrm{.}}40`$. The western condensation (A2) has a structure of two superimposed arcs and the eastern condensation (A1) also has an arc-like structure. For both condensations, this basic morphology is preserved to velocities of $`v_r100`$ km s-1. We find a pronounced arc of emission that connects the A1 and A2 condensations in the $`v_r25`$ km s-1 to $`v_r87`$ km s-1 maps. A1 and A2 no longer have their arc-like shapes and have a more compact morphology for $`v_r>118`$ km s-1. They approach each other in the spatial direction for increasing values of $`v_r`$, and have a separation of only $`1\stackrel{}{\mathrm{.}}15`$ in the $`v_r=304`$ km s-1 map. Condensation A1 dominates the total intensity for $`v_r>300`$ km s-1, and again develops a compact, arc-like shape for higher values of $`v_r`$.

Although not shown, the channel maps for the \[O I\], \[N II\], and \[S II\] lines show qualitatively similar morphologies as those for H$`\alpha `$ . However, the channel maps for the \[O III\] line show significant differences. The surface brightness of the two condensations in \[O III\] grows with increasing radial velocity, reaching a maximum level at +168 km s-1, and then the surface brightness decreases for larger velocities. No second maximum is detected. At the lowest radial velocities, subcondensation A2 shows an elongated structure that, despite its low signal-to-noise, resembles the arc-shaped structure also visible in the H$`\alpha `$ maps. At higher radial velocities, subcondensation A1 has a more compact morphology. Subcondensation A1 is detected at $`v_r`$ from +26 to +350 km s-1, and it dominates the total intensity for $`v_r`$ $`>`$ +188 km s-1. Subcondensation A1 shows a compact morphology at these radial velocities; no arc-like structure is detected in this subcondensation in the \[O III\] emission. As reported for the H$`\alpha `$ maps, both subcondenstions approach each other in the spatial direction as the radial velocity increases.

## 4 The line ratios

The \[S II\] 6716/6731 ratio depends only weakly on temperature, so it is a direct diagnostic of the electron density of the emitting plasma. Unfortunately, the other ratios between the observed lines do not have such a direct interpretation. However, ratios such as \[S II\] (6716+6731)/H$`\alpha `$ and \[O III\] 5007/H$`\alpha `$ do show clear trends as a function of shock velocity in models of HH jets (see, e. g., Hartigan et al. 1987), and can therefore be used to estimate the velocities of the shocks which dominate the emission along a given line of sight. Our spectra are not flux calibrated, so we are not able to give direct estimates of the values of the shock velocities. Therefore, our results can only be used to obtain a qualitative picture of the variations of the shock properties across the emitting region of HH 32.

From the position-velocity (PV) cubes of the \[O I\] 6300, \[N II\] 6583, \[S II\] 6716, 6731, \[O III\] 5007 and the H$`\alpha `$ lines, we have generated seven channel maps centered at heliocentric radial velocities of $`v_r=5`$, 48, 110, 172, 234, 296 and 358 km s-1 (each map has a velocity width of 47 km s-1). We have binned the maps spatially over $`5\times 5`$ pixels of the original spectrum to decrease the effective pixel size to $`0\stackrel{}{\mathrm{.}}25`$. Figure 5 shows the \[O I\] 6300 and \[N II\] 6583 channel maps divided by the corresponding H$`\alpha `$ maps, Figure 6 shows the \[S II\] (6716+6731)/H$`\alpha `$ and the \[S II\] 6716/6731 line ratio channel maps and Figure 7 shows the \[O III\] 5007 Å channel maps divided by the corresponding H$`\alpha `$ channel maps. The overall spatial size of the line ratio maps is the same as the H$`\alpha `$ channel maps, but the pixels have been binned by 5 in the x and y directions. The line ratio maps appear to cover a wider area than the H$`\alpha `$ emission maps because there is non-negligible flux in the lines at distances spatially extended from the condensations seen in H$`\alpha `$. As a function of radial velocity, the behavior of the line ratio maps can be described as:

* 6583/H$`\alpha `$ : at the position of knot A2, this line ratio has low values of $`0.4`$ at low radial velocities, and the ratio grows with increasing $`v_r`$. At the position of knot A1, the maximum values of $`0.7`$ for the line ratio occur at $`v_r=48`$ km s-1, and the ratios decrease for increasing $`v_r`$, reaching values of $`0.4`$ for $`v_r=296`$ km s-1. Interestingly, the region with the highest line ratio lies to the South of knot A1, with values of $`1`$ for $`v_r=48`$ km s-1.
* 6300/H$`\alpha `$ : knot A2 shows an initial increase in this line ratio, from $`0.2`$ for $`v_r=5`$ km s-1 to $`0.4`$ for $`v_r=110`$ km s-1. The line ratio has a $`0.3`$ value in the $`v_r=172`$ km s-1 velocity channel and then it increases for higher radial velocities, reaching a value of $`0.5`$ for $`v_r=296`$ km s-1. Knot A1 shows high values of $`0.7`$ at $`v_r=48`$ km s-1, decreasing to values of $`0.1`$ at $`v_r=296`$ km s-1. The highest values in this line ratio are found to the East of knot A1 in the $`v_r=110`$ km s-1 channel map.
* (6716+6731)/H$`\alpha `$ : this line ratio shows a general trend of decreasing values with increasing radial velocity for $`v_r=5172`$ km s-1, and then has similar values for all higher radial velocity channel maps. Knot A2 has line ratios of $`1`$ for $`v_r=5`$ km s-1, dropping to $`0.4`$ for $`v_r=110`$ km s-1, and then keeping this value for the higher radial velocity channels. Knot A1 has a more drastic drop in the line ratio, with values of $`2`$ at low velocities, decreasing to $`0.2`$ for $`v_r=234`$ km s-1.
* 6716/6731 : for both A1 and A2, this line ratio shows values of $`0.7`$ (corresponding to $`n_e230`$ cm-3) in the $`v_r=5`$ km s-1 map. The line ratio then decreases with increasing radial velocities, with values of $`0.6`$ ($`n_e420`$ cm-3) for $`v_r=48`$ km s-1 and $`0.55`$ ($`n_e650`$ cm-3) for $`v_r=110`$ km s-1. The line ratio at the positions of knots A1 and A2 remains close to the high density limit for all of the radial velocity channel maps with $`v_r>110`$ km s-1. It is clear in the lower heliocentric radial velocity channel maps (with $`v_r110`$ km s-1) that knots A1 and A2 are surrounded by a diffuse emitting region (filling most of the observed field in the $`v_r=5`$ km s-1 map) with line ratios of $`1`$ ($`n_e150`$ cm-3).
* 5007/H$`\alpha `$: at the position of condensation A1, this emission line ratio grows continually with radial velocities from -5 km s-1 up to +300 km s-1. At the condensation A2, this line ratio is roughly constant for radial velocities from +110 to +238 km s-1 and increases for larger radial velocities. This line ratio is clearly larger at the position of condensation A1 than at the knot A2.

These results show that the line ratios of HH 32 have strong spatial variations at arcsecond scales, as well as variations as a function of radial velocity. The ratio maps for the \[S II\] lines show that the brighter A1 and A2 condensates are enclosed within a diffuse emitting region with lower electron density than the knots. This diffuse low electron density region surrounding knots A1 and A2 has low H$`\alpha `$ emission; it has significant values for both the high excitation \[N II\] 6583/H$`\alpha `$ and the low excitation \[O I\] 6300/H$`\alpha `$ line ratios (see Figure 5). Indeed, the peak \[N II\] 6583/H$`\alpha `$ and \[O I\] 6300/H$`\alpha `$ values are found in this diffuse region rather than centered on the positions of the A1 or A2 condensations.

Although the behavior of the line ratios is very complex, some of the observed trends are consistent with a simple, plane shock wave interpretation. For example, the \[O III\]/H$`\alpha `$ ratio in knot A1 grows monotonically with radial velocity, consistent with the fact that this line ratio increases with increasing shock velocities. Also, the \[S II\] (6716+6731)/H$`\alpha `$ ratio in the A1 condensation first decreases and then stabilizes as a function of increasing $`v_r`$, which qualitatively follows the line ratio vs. shock velocity trend predicted from plane-parallel shock models (see Hartigan et al. 1987). However, different regions of HH 32 show different trends in the line ratio versus radial velocity which cannot be easily interpreted with the present plane shock models. As an example, knot A1 shows a decrease in the \[O I\] 6300/H$`\alpha `$ ratio as a function of increasing radial velocity, while knot A2 has the opposite trend. Interestingly, for this oxygen line ratio the plane shock models predict a decrease in intensity as a function of shock velocities (from 20 to 90 km/s), followed by a strong increase for larger velocities (from 100 to 300 km/s, see Hartigan et al. 1987). The observed line ratios of \[O I\] 6300/H$`\alpha `$ in the two knots are not consistent with the present models.

The behavior as a function of $`v_r`$ of the line emission of knots A1 and A2 can be seen in a more quantitative way in Figure 8. In this figure, we have plotted the ratios between the H$`\alpha `$, \[O III\] 5007, \[O I\] 6300, \[N II\] 6583 and \[S II\] (6716+31) line emission of A1 and A2 as a function of $`v_r`$. From this graph we see that knot A1 is brighter than A2 for heliocentric radial velocities $`v_r<150`$ km s-1 in the \[O I\] 6300 and \[S II\] (6716+31) lines, and A2 is brighter than A1 for $`150<v_r<350`$ km s-1. For $`v_r>350`$ km s-1, knot A1 is again brighter than A2.

The relative intensities between A1 and A2 in the H$`\alpha `$ and \[N II\] 6583 lines show a similar behavior to \[O I\] and \[S II\], but have a much shallower dependence on radial velocity for $`v_r<300`$ km s-1. The H$`\alpha `$ fluxes of A1 and A2 differ by at most 20 % in this radial velocity range. In the H$`\alpha `$ and \[N II\] lines, condensation A1 brightens quite steeply relative to A2 for $`v_r>300`$ km s-1, becoming twice as bright for $`v_r=390`$ km s-1. The relative \[O III\] intensity between A1 and A2 shows a different behavior than the other emission lines plotted (the \[O III\] channel maps are shown in Figure 7). The \[O III\] emission seen in A1 brightens quite steeply relative to A2 for $`v_r`$ from 0 to +120 km s-1. The relative intensities remain more or less constant with a mean value of 1.75 for radial velocities in the range from +120 to +300 km s-1, and decrease again for larger velocities ($`v_r`$ $`>`$ 300 km s-1).

## 5 A 3/2-D bow shock model

In order to interpret the observations of HH 32, we consider a traditional “3/2-D” bow shock model (Hartigan et al. 1987). It is assumed that the shock velocity is equal to the component of the pre-bow shock flow velocity normal to the shock surface. Such models have been used to interpret line profiles of several HH objects (see, e. g., Hartigan et al. 1987), with 1D or 2D spatial resolution (Raga & Böhm 1985; 1986; Solf et al. 1991; Morse et al. 1992). To date, one of the most successful applications of the 3/2-D bow shock models comes from the line profile fits to the HH 32 outflow (Solf et al. 1986; Hartigan et al. 1987).

We consider a bow shock with a functional form:

$$\frac{z}{a}=\left(\frac{r}{a}\right)^p,$$
(1)
where $`z`$ is measured along the symmetry axis and $`r`$ is the cylindrical radius. The constants $`a`$ and $`p`$ are the free parameters of the model, with $`p`$ determining the shape of the bow shock, and $`a`$ its physical size. The bow shock is assumed to be moving at a velocity $`v_{bs}`$ with respect to a uniform pre-shock medium, directed at an angle $`\varphi `$ with respect to the plane of the sky. Without any loss of generality, we can set $`a=1`$, and scale the intensity maps predicted from the model in an arbitrary way to produce spatial scales comparable to the ones of HH 32A.

In order for the model to fit closely with the observations, we have assumed that the wings of the bow shock (with the shape given by equation 1) suddenly end at a distance $`z_{max}`$, as measured along the symmetry axis from the head of the bow shock. We then use the geometry of the system to calculate: 1) the shock velocity, which is the velocity component of the bow shock normal to the surface of the shock, 2) the velocity of the post-bow shock flow which is approximately equal to the normal velocity for the compressions found in radiative shocks, and 3) the intensity per unit area of the bow shock derived by interpolating between the line emissions predicted from a series of plane-parallel shock models with differing shock velocities. The necessary geometrical construction we have used is described by Raga & Böhm (1986) and in more detail by Hartigan et al. (1987).

We computed the emission using the tabulation of plane-parallel models with a “self consistent” pre-ionization and pre-shock density of $`n_{pre}=100`$ cm-3 as described by Hartigan et al. (1987). In the models, we fixed the bow shock velocity to a $`v_{bs}=350`$ km s-1 value necessary to obtain line widths comparable to the observed ones, and we adopt an orientation angle of $`\varphi =70^{}`$ moving away from the observer, as derived by the observed line profiles, radial velocities and proper motions (Raga et al. 1986; Hartigan et al. 1987; Curiel et al. 1997). We further assume that the pre-shock medium is at rest with respect to the outflow source.

At this point, the power law index $`p`$ is the only free parameter in the model. If the shock wings extend to large distances from the head of the bow shock, then we can set $`z_{max}\mathrm{}`$ (see equation 1). In this case, the best fits for the model are obtained for $`p23`$. However, if we allow a finite value for $`z_{max}`$, the models more closely resemble the observations. In particular, we choose a model with $`p=2`$ and $`z_{max}=1.5a`$. From this bow shock model, we have computed velocity channel maps which are used to compare with the observations of HH 32 (shown in Figure 9).

From comparison of the velocity channel maps shown in Figures 4 and 9, we find that the shock model predicts structures that are in qualitative agreement with the HH 32 observations. At low radial velocities ($`v_r=1590`$ km s-1), the model shows strong H$`\alpha `$ emission distributed in a single arc-like filament. The observations of HH 32 also show strong H$`\alpha `$ emission but it is distributed in two or three arc-like structures. At $`v_r=100200`$ km s-1, the bow shock model shows fainter H$`\alpha `$ emission with a circular structure, which becomes more compact for increasing radial velocities. In this velocity range, the HH 32 data also shows fainter, more compact emission at increasing radial velocities, but with a more complex morphology of knots (see Figure 4). For $`v_r>200`$ km s-1, both the model and HH 32 data show an increase in H$`\alpha `$ emission, which is then followed by an intensity decrease, with tighter knot-like structures at the higher radial velocities. There is a lack of quantitative agreement because the high velocity H$`\alpha `$ peak occurs at $`v_r=308`$ km s-1 in the model, but at $`v_r=227`$ km s-1 in HH 32. Also, while the model shows a single condensation for $`200<v_r<340`$ km s-1, HH 32 shows two or three compact condensations, and it does not converge to a single condensation until $`v_r350`$ km s-1. Our simple model considers only one bow shock observed at an angle of $`70^{}`$, comparison with the observations shows that a more complicated structure of 2 (or more) bow shocks would better fit the observations (see Raga et al. 2003).

We computed line ratio maps from our bow shock model and did not obtain results that qualitatively agree with the observations presented in Figures 5, 6 and 7. This discrepancy is not surprising because of the differences seen in the H$`\alpha `$ channel maps predicted by the model and observed in HH 32. The significant structure observed in the line ratio maps and the lack of consistency with our calculations is difficult to interpret uniquely, as a superposition of several shocks along the viewing geometry of HH 32 will complicate the models.

There are two main criticisms for the 3/2-D bow shock model we have utilized to describe the observations of HH 32. First, it is obvious that the simple, parametrized shape (equation 1) which we have assumed for the emitting shock structure is not appropriate for describing the more complex line of sight superposition of shocks in HH 32. Secondly, a 3/2-D shock model is based on the assumption that the post-bow shock cooling distance is small compared to the size of the bow shock; this is probably not the case in this scenario. As can be seen from the models of Hartigan et al. (1987), the cooling distance behind the head of a 250 km/s shock has a value $`d_c10^{16}`$ cm, which is comparable to the sizes of condensations A1 and A2. The model we have used is qualitatively successful at describing the observations of HH 32, but the assumptions on which the 3/2-D bow shock models are based are not ideal for describing this object. We present and discuss a more detailed 3D simulation of multiple working surfaces in a different paper (Raga et al. 2003)

## 6 Conclusions

We have obtained a high resolution spectrum of HH 32A with 2D spatial resolution which provides a wealth of information about the excitation and kinematics of the HH 32 outflow. From these spectra, we have constructed channel maps at velocities ranging from $`20`$ to $`+400`$ km s-1 for the H$`\alpha `$, \[O III\] 5007, \[O I\] 6300, \[N II\] 6583, \[S II\] 6716 and \[S II\] 6731 emission lines; the H$`\alpha `$ maps are presented in Figure 4 and \[O III\] is shown in the central panels of Figure 7. The channel maps of the emission lines of H$`\alpha `$, \[O I\], \[N II\], and \[S II\] have a qualitatively similar behavior. In the lower heliocentric radial velocity maps ($`v_r<100`$ km s-1), HH 32A has a morphology of two or three superimposed arc-like features. At higher radial velocities, the emission becomes progressively more concentrated, showing two or three compact condensations. Finally, for $`v_r>350`$ km s-1, a single condensation is seen.

Figures 5, 6 and 7 show complex structures in the line ratio maps. The local maxima and minima do not coincide necessarily with the observed knots in the H$`\alpha `$ images, and the ratios vary for each emission line over the velocity width of the feature. The complex morphological structures suggest that we are likely observing the emission from several compact knots, superimposed on a diffuse emission component, with the knots and the diffuse component having different spectral properties. This interpretation has been previously suggested to explain 2D spatial resolution spectra of HH 2 (Böhm & Solf 1992).

The observations of HH 32 can be interpreted in a qualitative way with a simple “3/2-D” bow shock model. We constructed model H$`\alpha `$ velocity channel maps that show a transition between a single arc-like feature at low radial velocities, to a more concentrated and centrally peaked condensation at high radial velocities (Figure 9). This basic change in morphology is also seen in the velocity channel maps of HH 32A (Figure 4), but the structures in the data are much more complex.

From comparison with the simple bow shock model, we conclude that HH 32A is likely a superposition of two or three bow shocks, corresponding to different working surfaces along the HH 32 outflow. HH 32A shows at least two condensations (A1 and A2) with spatial and kinematic properties which resemble the predictions from a single bow shock model. The fact that we see such a superposition is not surprising; the $`\varphi 70^{}`$ orientation of HH 32 with respect to the plane of the sky will lead naturally to a line of sight overlap of features along the outflow axis. HH 32 might be intrinsically similar to the collimated outflow observed in HH 34, we might be viewing two working surfaces catching up with each other as seen in HH 34S (Reipurth et al. 2002; Raga et al. 2002). The arc-like structures and multiple condensations observed in the data of HH 32A could then correspond to a line of sight superposition of the working surfaces. We present a study of this scenario based on more detailed 3D numerical simulations in Raga et al. (2003).

We thank Inger Jørgensen, GMOS-N instrument scientist, for approving this project for System Verification observations. A. Riera acknowledges the ICN-UNAM for support during her sabbatical. A. Riera was supported by the MCyT grant AYA2002-00205 (Spain). The work of A. C. Raga was supported by the CONACyT grant 36572-E.