# Fast two-beam collisions in a linear optical medium with weak cubic loss in
spatial dimension higher than 1

Avner Peleg1, Toan T. Huynh2,3, and Quan M. Nguyen4 1Department of
Mathematics, Azrieli College of Engineering, Jerusalem 9371207, Israel
2Department of Mathematics, University of Science, Vietnam National
University-HCMC, Ho Chi Minh City, Vietnam  3Department of Mathematics,
University of Medicine and Pharmacy at Ho Chi Minh City, Ho Chi Minh City,
Vietnam  4Department of Mathematics, International University, Vietnam
National University-HCMC, Ho Chi Minh City, Vietnam

###### Abstract

We study the dynamics of fast two-beam collisions in linear optical media with
weak cubic loss in spatial dimension higher than 1. For this purpose, we
extend the perturbation theory that was developed for analyzing two-pulse
collisions in spatial dimension 1 to spatial dimension 2. We use the extended
two-dimensional version of the perturbation theory to show that the collision
leads to a change in the beam shapes in the direction transverse to the
relative velocity vector. Furthermore, we show that in the important case of a
separable initial condition for both beams, the longitudinal part in the
expression for the amplitude shift is universal, while the transverse part is
not universal. Additionally, we demonstrate that the same behavior holds for
collisions between pulsed optical beams in spatial dimension 3. We check these
predictions of the perturbation theory along with other predictions concerning
the effects on the collision of partial beam overlap and anisotropy in the
initial condition by extensive numerical simulations with the weakly perturbed
linear propagation model in spatial dimensions 2 and 3. The agreement between
the perturbation theory and the simulations is very good. Therefore, our study
significantly extends and generalizes the results of previous works, which
were limited to spatial dimension 1.

Optical beam propagation, cubic loss, beam collisions, perturbation theory,
amplitude dynamics.

##  I Introduction

Linear evolution models have an important role in many areas of science.
Examples include the linear diffusion equation Van_Kampen , the linear wave
equation Whitham99 , the linear propagation equation Ishimaru2017 ; Siegman86
; Kogelnik66 , and the linear Schrödinger equation Merzbacher98 . In many
cases, the physical systems that are described by these linear evolution
models include weak nonlinear dissipation Van_Kampen ; Agrawal2007a . As a
result, the latter physical systems are more accurately described by perturbed
linear evolution models with weak nonlinear dissipation. The presence of the
nonlinear dissipation induces new physical effects, which do not exist in the
unperturbed linear physical systems. A major example is the change in the
pulse mass or energy during fast collisions between pulses of the linear
propagation model PNH2017B ; QMN2018 ; NHP2020 . Since the pulses of the
linear physical systems (and also of their weakly perturbed counterparts) are
not shape preserving, one does not expect to observe simple dynamics in these
collisions. As a result, one also does not expect to be able to make simple
general statements about the collision-induced effects.

In two previous works PNH2017B ; NHP2020 , we showed that the opposite is in
fact true for fast two-pulse collisions. The latter are collisions, in which
the collision length, i.e., the length of the interval where the two pulses
overlap, is much smaller than all the other length scales in the problem
fast_collisions . In Refs. PNH2017B ; NHP2020 , we showed that the amplitude
shifts in fast two-pulse collisions in linear physical systems, weakly
perturbed by nonlinear dissipation, exhibit simple soliton-like behavior. The
behavior was demonstrated for the following two central cases: (1) systems
described by the linear propagation equation with weak cubic loss; (2) systems
described by the linear diffusion-advection equation with weak quadratic loss.
We first developed a perturbation method for analyzing the fast two-pulse
collision dynamics in these weakly perturbed linear systems. We then used the
perturbation method to show that in both cases, the expressions for the
collision-induced amplitude shifts in the presence of weak nonlinear loss have
the same simple form as the expression for the amplitude shift in a fast
collision between two solitons of the nonlinear Schrödinger equation in the
presence of weak cubic loss. Furthermore, in Ref. NHP2020 , we showed that the
expressions for the amplitude shifts are universal in the sense that they are
independent of the details of the initial pulse shapes. In addition, we found
that within the leading order of the perturbation theory, the pulse shapes are
not changed by the collision. The perturbation theory predictions for the
collision-induced amplitude shifts were verified by extensive numerical
simulations with the two perturbed linear evolution models for a variety of
initial pulse shapes PNH2017B ; NHP2020 . Additionally, in Ref. QMN2018 , we
showed that the amplitude shift in a fast two-pulse collision in systems
described by the linear propagation model with high-order nonlinear loss can
be calculated by the same perturbation method that was developed in Ref.
PNH2017B .

The three studies in Refs. PNH2017B ; QMN2018 ; NHP2020  were limited to
spatial dimension 1 spatial_dimension . As a result, these studies did not
consider important collisional effects, which exist only is spatial dimension
higher than 1, such as anisotropy and partial pulse overlap. Additionally, it
is not clear if the simple (universal) dependence of the expressions for the
collision-induced amplitude shifts on the physical parameters that was found
in Refs. PNH2017B  and NHP2020  remains valid in the high-dimensional problem.
It is also unclear if the pulse shapes remain unchanged in fast collisions in
the presence of cubic (or quadratic) loss in the high-dimensional problem.
Thus, all the key aspects of the fast two-pulse collision problem, which are
associated with the collision dynamics in spatial dimension higher than 1,
were not addressed in previous studies.

In the current paper, we address the important aspects of the high-dimensional
fast two-pulse collision problem that were mentioned in the preceding
paragraph. For this purpose, we first develop a perturbation method, which
generalizes the perturbation method that was introduced in Refs. PNH2017B ;
NHP2020  for the one-dimensional problem in three major ways. (1) It extends
the perturbative calculation from spatial dimension 1 to spatial dimension 2,
and enables further extension of the calculation to a general spatial
dimension in a straightforward manner. (2) It provides a perturbative
calculation of the collision-induced dynamics of the beam shape both inside
and outside of the collision interval. In contrast, the perturbative
calculation of Refs. PNH2017B ; NHP2020  was limited to the collision interval
only. (3) It enables the discovery and analysis of several collision-induced
effects, which exist only in the high-dimensional problem.

We use the generalized version of the perturbation method to derive formulas
for the collision-induced changes in the beam shapes and amplitudes in spatial
dimension 2. We find that for a general initial condition, the collision leads
to a change in the beam shape in the direction transverse to the relative
velocity vector between the beam centers. Additionally, we find that in the
important case of an initial condition that is separable for both beams, the
beam shape in the longitudinal direction is not changed by the collision
within the leading order of the perturbation theory. Furthermore, we show that
for a separable initial condition, the longitudinal part in the expression for
the amplitude shift is universal, while the transverse part is not universal
and is proportional to the integral of the product of the beam intensities
with respect to the transverse coordinate. We also show that the same behavior
of the longitudinal and transverse parts in the expression for the collision-
induced amplitude shift exists in collisions between pulsed-beams in spatial
dimension 3.

We check these predictions of the perturbation theory together with other
predictions concerning the effects of partial beam overlap and anisotropy in
the initial condition by extensive numerical simulations with the perturbed
linear propagation model in spatial dimensions 2 and 3. The simulations in
spatial dimension 2 are carried out for four different two-beam collision
setups. These setups demonstrate the following four major effects and
properties of the collision that either exist only in spatial dimension higher
than 1, or are qualitatively different from their one-dimensional
counterparts. (1) The universality of the longitudinal part in the expression
for the collision-induced amplitude shift. (2) The effect of partial beam
overlap. (3) The effect of anisotropy in the initial condition. (4) The
collision-induced change in the beam shape in the transverse direction. The
prediction for universal behavior of the longitudinal part in the expression
for the amplitude shift is also checked in spatial dimension 3 by numerical
simulations of collisions between pulsed optical beams. In all the simulation
setups we obtain very good agreement between the perturbation theory and the
numerical simulations. Therefore, the simulations validate the theoretical
predictions for the four high-dimensional effects and properties of the
collision and show that the extended perturbation approach can indeed be used
for analyzing the effects of fast two-beam collisions in spatial dimension
higher than 1.

In a related work, we studied the dynamics of fast two-pulse collisions in
systems described by linear diffusion-advection models with weak quadratic
loss in spatial dimension higher than 1 PHN2020 . We first developed a
perturbation method for analyzing the collision dynamics, which is similar to
the one introduced in the current paper. Using this perturbation method and
numerical simulations, we showed that the collision-induced changes in pulse
shapes and amplitudes in these systems exhibit similar behavior to the one
found in the current paper PHN2020 . Thus, the current paper and the related
study of Ref. PHN2020  significantly extend and generalize the results of the
previous works in Refs. PNH2017B ; QMN2018 ; NHP2020 , which were limited to
spatial dimension 1. We also comment that detailed analytic results on
collisions between pulse solutions of linear or nonlinear evolution models in
the presence of nonlinear dissipation in spatial dimension higher than 1 are
quite scarce. Therefore, the current work and the work in Ref. PHN2020  also
significantly extend the understanding of the more general high-dimensional
problem of two-pulse collisions in the presence of nonlinear dissipation.

We choose to study two-beam collisions in the presence of cubic loss, since
cubic loss is important in many optical systems, and is therefore a major
example for nonlinear dissipative perturbations. The optical medium’s cubic
loss typically arises due to two-photon absorption (2PA) Agrawal2007a ;
Dekker2007 ; Borghi2017 ; Boyd2008 . Propagation of optical pulses and optical
beams in the presence of cubic loss has been studied in many earlier works,
both in weakly perturbed linear media PNH2017B ; NHP2020 ; Perry97 ; Liang2005
; Cohen2005b ; Cohen2004 , and in nonlinear media Malomed89 ; Stegeman89 ;
Silberberg90 ; Aceves92 ; Kivshar95 ; Tsoy2001 ; Silberberg2008 ; PCDN2009 ;
PNC2010 ; Gaeta2012 ; PC2018 ; PC2020 . The subject attracted renewed
attention in recent years due to the importance of 2PA in silicon
nanowaveguides, which are expected to play a major role in many applications
in optoelectronic devices Agrawal2007a ; Dekker2007 ; Borghi2017 ; Gaeta2008 ;
Soref2006 . In the current paper, we assume that the optical medium is weakly
nonlinear and neglect the effects of cubic (Kerr) nonlinearity. We emphasize
that this assumption was successfully used in previous experimental and
theoretical works, see, e.g., Refs. Perry97 ; Liang2005 ; Cohen2005b ;
Cohen2004 . For the same reason, we also neglect the effects of high-order
nonlinear loss on the collision. We remark that the latter effects can be
described by the same perturbation method that is introduced in the current
paper (see also Ref. QMN2018 , where the calculation was carried out for
spatial dimension 1).

The rest of the paper is organized as follows. In Section II, we present the
extended perturbation method for calculating the amplitude and beam shape
dynamics in fast collisions between beams of the linear propagation model in
spatial dimension 2. In Section III, we present the perturbation theory
predictions and the results of numerical simulations with the weakly perturbed
linear propagation model for four major collision setups. These setups
demonstrate four major effects and properties of the collision that exist only
in spatial dimension higher than 1. In Section IV, we present the main
predictions of the perturbation theory for collisions between pulsed-beams of
the linear propagation equation in spatial dimension 3. We also present a
comparison between the perturbation theory prediction for the collision-
induced amplitude shift and the results of numerical simulations with the
weakly perturbed linear propagation model. In Section V, we summarize our
conclusions. The five Appendixes contain calculations that support the
material in the main body of the paper.

##  II The perturbation theory for fast two-beam collisions in spatial
dimension 2

###  II.1 Introduction

We consider fast collisions between two optical beams in a three-dimensional
linear optical medium with weak cubic loss. We assume that the beams propagate
along the $z$ axis with beam-steering in the $xy$ plane, and that the
propagation is accurately described by the paraxial approximation Ishimaru2017
; Siegman86 ; Kogelnik66 . For each given value of $z$, the distribution of
the electric field is a function of $x$ and $y$. Therefore, we can think about
the $z$ coordinate as a dynamical coordinate, and about the $x$ and $y$
coordinates as the actual spatial coordinates, which help describe the
distribution of the electric field for each value of $z$. We refer to the
dimension of the space, in which the distribution of the electric field is
described (for a given $z$) as the spatial dimension. Thus, in the current
problem, the spatial dimension is 2 and the total dimension is 3. The
propagation is described by a $(2+1)$-dimensional propagation model, where the
2 in this terminology corresponds to the spatial dimension, and the 1 is the
dimension of the dynamical axis (the $z$ axis).

We take into account the effects of isotropic diffraction and weak cubic loss,
as well as the velocity difference between the beam centers, which is
controlled by beam-steering McManamon2009 ; Carlson88 ; Brandl2013 ; Gomez2015
; Oh2016 ; Kivshar2012 ; Kartashov2009 ; Lederer2002 . For simplicity and
without loss of generality, we assume that the relative velocity vector
between the beam centers lies along the $x$ axis. This choice enables one to
obtain closed formulas for the collision-induced changes in beam shapes and
amplitudes, and in this manner, enables a significantly deeper insight into
the collision dynamics. Furthermore, in Appendix E, we show that the choice of
the relative velocity vector along the $x$ axis does not change the value of
the collision-induced amplitude shift obtained by our perturbation approach.
That is, the latter value is invariant under rotations of the $x$ and $y$
axes. Thus, the dynamics of the two-beam collision is described by the
following weakly perturbed linear propagation model:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!i\partial_{z}\psi_{1}+\partial_{x}^{2}\psi_{1}+\partial_{y}^{2}\psi_{1}=-i\epsilon_{3}|\psi_{1}|^{2}\psi_{1}-2i\epsilon_{3}|\psi_{2}|^{2}\psi_{1},$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!i\partial_{z}\psi_{2}+id_{11}\partial_{x}\psi_{2}+\partial_{x}^{2}\psi_{2}+\partial_{y}^{2}\psi_{2}=-i\epsilon_{3}|\psi_{2}|^{2}\psi_{2}-2i\epsilon_{3}|\psi_{1}|^{2}\psi_{2}.\\!\\!\\!\\!\\!\\!\\!\\!$
(1)

In Eq. (1), $\psi_{j}$ with $j=1,2$ are proportional to the electric fields of
the beams, and $x$, $y$, and $z$ are the spatial coordinates Dimensions1 . In
addition, $d_{11}$ is the coefficient related to the velocity difference
between the beam centers (the beam-steering coefficient), and $\epsilon_{3}$
is the cubic loss coefficient, which satisfies $0<\epsilon_{3}\ll 1$. The
terms $\partial_{x}^{2}\psi_{j}$ and $\partial_{y}^{2}\psi_{j}$ on the left
hand side of Eq. (1) describe the effects of isotropic diffraction, while
$id_{11}\partial_{x}\psi_{2}$ is related to the velocity difference between
the beam centers. The first and second terms on the right hand side of Eq. (1)
describe intra-beam and inter-beam effects due to cubic loss. In the current
paper, we do not take into account the effects of linear loss, since these
effects do not change the form of the expressions for the collision-induced
changes in beam amplitudes and shapes. Furthermore, the simple effects of
linear loss on amplitude dynamics can be incorporated into the analysis in
exactly the same manner as was done in Refs. PNH2017B ; NHP2020  for spatial
dimension 1 (see also Appendix B). We remark that the same perturbed linear
propagation model (with some changes in the physical variables) also describes
the dynamics of a fast collision between two pulsed-beams in a two-dimensional
linear optical medium (e.g., a planar waveguide) with weak cubic loss. In this
case, the coordinate $x$ is replaced by the time variable $t$, the term
$id_{11}\partial_{t}\psi_{2}$ describes the effects of the group velocity
difference, and the terms $\partial_{t}^{2}\psi_{j}$ describe the effects of
second-order dispersion. The more general case of fast collisions between
pulsed-beams in a three-dimensional medium (i.e., in spatial dimension 3 and
total dimension 4) is studied in section IV.

We consider fast collisions between beams with generic initial shapes and with
tails that decay sufficiently fast, such that the values of the integrals
$\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dy|\psi_{j}(x,y,0)|^{2}$ are
finite. We assume that the initial beams can be characterized by the following
parameters. (1) The initial amplitudes $A_{j}(0)$. (2) The initial beam
widths, i.e., the widths of the maxima of $|\psi_{j}(x,y,0)|$, which can be
expressed in terms of the widths along the $x$ and $y$ axes, $W_{j0}^{(x)}$
and $W_{j0}^{(y)}$, respectively. (3) The initial positions of the beam
centers, i.e., the locations of the maxima of $|\psi_{j}(x,y,0)|$, which are
denoted by $(x_{j0},y_{j0})$. (4) The initial phases $\alpha_{j0}$. Therefore,
the initial electric fields of the optical beams can be expressed as:

$\displaystyle\psi_{j}(x,y,0)=A_{j}(0)h_{j}(x,y)\exp(i\alpha_{j0}),$ (2)

where $h_{j}(x,y)$ is real-valued. Note that for brevity of notation, we did
not write the dependence of the function $h_{j}(x,y)$ on the beam parameters
explicitly. We are also interested in the important case, where the initial
electric fields of both beams are separable, i.e., where each of the functions
$\psi_{j}(x,y,0)$ can be expressed as a product of a function of $x$ and a
function of $y$ laser_modes . In this case, the initial electric fields are
given by:

$\displaystyle\\!\\!\\!\\!\\!\\!\psi_{j}(x,y,0)=A_{j}(0)h_{j}^{(x)}[(x-x_{j0})/W_{j0}^{(x)}]h_{j}^{(y)}[(y-y_{j0})/W_{j0}^{(y)}]\exp(i\alpha_{j0}).$
(3)

In what follows, we will also consider cases where the initial electric field
is separable for one beam and nonseparable for the other beam.

In the current paper, we study the collision-induced dynamics of complete fast
collisions. The complete collision assumption means that the beams are well-
separated before and after the collision. More specifically, in these
collisions, the values of the $x$ coordinate of the beam centers at $z=0$ and
at the final propagation distance $z_{f}$, $x_{j0}$ and $x_{j}(z_{f})$,
satisfy $|x_{20}-x_{10}|\gg W_{10}^{(x)}+W_{20}^{(x)}$ and
$|x_{2}(z_{f})-x_{1}(z_{f})|\gg W_{1}^{(x)}(z_{f})+W_{2}^{(x)}(z_{f})$, where
$W_{j}^{(x)}(z_{f})$ are the beam widths in the $x$ direction at $z=z_{f}$. To
obtain the condition for a fast collision, we first define the collision
length $\Delta z_{c}$, as the distance along which the beam widths in the $x$
direction overlap. From this definition it follows that $\Delta
z_{c}=2(W_{10}^{(x)}+W_{20}^{(x)})/|d_{11}|$. For a fast collision, we require
that $\Delta z_{c}$ would be much smaller than the smallest diffraction length
in the problem. We note that the diffraction lengths of the $j$th beam in the
$x$ and $y$ directions are $z_{Dj}^{(x)}=W_{j0}^{(x)2}/2$ and
$z_{Dj}^{(y)}=W_{j0}^{(y)2}/2$, respectively. Thus, the smallest diffraction
length $z_{D}^{(min)}$ is
$z_{D}^{(min)}=\min\left\\{z_{D1}^{(x)},z_{D2}^{(x)},z_{D1}^{(y)},z_{D2}^{(y)}\right\\}$.
Requiring that $\Delta z_{c}\ll z_{D}^{(min)}$ and using the definition of
$\Delta z_{c}$, we obtain that the condition for a fast collision can be
expressed as $2(W_{10}^{(x)}+W_{20}^{(x)})\ll|d_{11}|z_{D}^{(min)}$.

###  II.2 Calculation of the collision-induced changes in the beam shape and
amplitude for a general initial condition

####  II.2.1 Introduction

The perturbation method that we present here generalizes the perturbation
method presented in Refs. PNH2017B ; NHP2020  in three major ways. First, it
extends the calculation from spatial dimension 1 to spatial dimension 2 nD .
Second, it provides a perturbative calculation and analytic expressions for
the collision-induced change in the beam shape both in the collision interval
and away from the collision interval, whereas the calculation of the change in
the beam shape in Refs. PNH2017B ; NHP2020  was limited to the collision
interval only. Third, it helps uncover several collision-induced effects,
which exist only in spatial dimension higher than 1. In the first step in the
perturbative calculation, we look for a solution of Eq. (1) in the form:

$\displaystyle\psi_{j}(x,y,z)=\psi_{j0}(x,y,z)+\phi_{j}(x,y,z),$ (4)

where $j=1,2$, $\psi_{j0}$ are the solutions of the weakly perturbed linear
propagation equations without the inter-beam interaction terms, and $\phi_{j}$
describe corrections to the $\psi_{j0}$ due to the effects of inter-beam
interaction on the collision. By their definition, the $\psi_{j0}$ satisfy the
following two equations:

$\displaystyle
i\partial_{z}\psi_{10}+\partial_{x}^{2}\psi_{10}+\partial_{y}^{2}\psi_{10}=-i\epsilon_{3}|\psi_{10}|^{2}\psi_{10},$
(5)

and

$\displaystyle
i\partial_{z}\psi_{20}+id_{11}\partial_{x}\psi_{20}+\partial_{x}^{2}\psi_{20}+\partial_{y}^{2}\psi_{20}=-i\epsilon_{3}|\psi_{20}|^{2}\psi_{20}.$
(6)

Substituting the ansatz (4) into Eq. (1) and using Eqs. (5) and (6), we obtain
equations for the $\phi_{j}$. We concentrate on the calculation of $\phi_{1}$,
since the calculation of $\phi_{2}$ is similar. The equation for $\phi_{1}$ in
the leading order of the perturbative calculation is

$\displaystyle
i\partial_{z}\phi_{1}+\partial_{x}^{2}\phi_{1}+\partial_{y}^{2}\phi_{1}=-2i\epsilon_{3}|\psi_{20}|^{2}\psi_{10}.$
(7)

Note that in writing Eq. (7), we neglected the high-order terms containing
$\phi_{j}$ on the right hand side of the equation.

In solving the equation for $\phi_{1}$, we distinguish between two intervals
along the $z$ axis, the collision interval and the post-collision interval.
These intervals are defined in terms of the collision distance $z_{c}$, which
is the distance at which the $x$ coordinates of the beam centers coincide,
i.e., $x_{1}(z_{c})=x_{2}(z_{c})$. The collision interval is the small
interval $z_{c}-\Delta z_{c}/2\leq z\leq z_{c}+\Delta z_{c}/2$ centered about
$z_{c}$, in which the two beams are overlapping. The post-collision interval
is the interval $z>z_{c}+\Delta z_{c}/2$, in which the beams are no longer
overlapping.

####  II.2.2 Collision-induced effects in the collision interval

We substitute $\psi_{j0}(x,y,z)=\Psi_{j0}(x,y,z)\exp[i\chi_{j0}(x,y,z)]$ and
$\phi_{1}(x,y,z)=\Phi_{1}(x,y,z)\exp[i\chi_{10}(x,y,z)]$, where $\Psi_{j0}$
and $\chi_{j0}$ are real-valued, into Eq. (7), and obtain the following
equation for $\Phi_{1}$:

$\displaystyle
i\partial_{z}\Phi_{1}-\left(\partial_{z}\chi_{10}\right)\Phi_{1}+\left[\partial_{x}^{2}\Phi_{1}+2i\left(\partial_{x}\chi_{10}\right)\partial_{x}\Phi_{1}\right.$
$\displaystyle\left.+i\left(\partial_{x}^{2}\chi_{10}\right)\Phi_{1}-\left(\partial_{x}\chi_{10}\right)^{2}\Phi_{1}\right]+\left[\partial_{y}^{2}\Phi_{1}+2i\left(\partial_{y}\chi_{10}\right)\partial_{y}\Phi_{1}\right.$
$\displaystyle\left.+i\left(\partial_{y}^{2}\chi_{10}\right)\Phi_{1}-\left(\partial_{y}\chi_{10}\right)^{2}\Phi_{1}\right]=-2i\epsilon_{3}\Psi_{20}^{2}\Psi_{10}.$
(8)

Since the collision length $\Delta z_{c}$ is of order $1/|d_{11}|$, the term
$i\partial_{z}\Phi_{1}$ is of order $|d_{11}|\times O(\Phi_{1})$.
Additionally, the term $-2i\epsilon_{3}\Psi_{20}^{2}\Psi_{10}$ is of order
$\epsilon_{3}$. Equating the orders of $i\partial_{z}\Phi_{1}$ and
$-2i\epsilon_{3}\Psi_{20}^{2}\Psi_{10}$, we find that $\Phi_{1}$ is of order
$\epsilon_{3}/|d_{11}|$. In addition, $\Phi_{1}$ does not contain any fast
dependence on $x$ and $y$, and $\chi_{10}$ does not contain any fast
dependence on $x$, $y$, and $z$. As a result, all the other terms in Eq. (8)
are of order $\epsilon_{3}/|d_{11}|$ or higher, and can therefore be
neglected. It follows that in the leading order of the perturbative
calculation, the equation for $\Phi_{1}$ is

$\displaystyle\partial_{z}\Phi_{1}=-2\epsilon_{3}\Psi_{20}^{2}\Psi_{10}.$ (9)

Equation (9) has the same form as the equation obtained for a fast collision
between two pulses of the linear propagation equation in the presence of weak
cubic loss in spatial dimension 1 PNH2017B ; NHP2020 . It also has the same
form as the equation obtained for a fast collision between two solitons of the
NLS equation in the presence of weak cubic loss in spatial dimension 1 PNC2010
.

We now introduce the following approximations to the solutions
$\psi_{j0}(x,y,z)$ of Eqs. (5) and (6):

$\displaystyle\psi_{j0}(x,y,z)\simeq A_{j}(z)\tilde{\psi}_{j0}(x,y,z),$ (10)

where $A_{j}(z)$ are the $z$ dependent beam amplitudes, and

$\displaystyle\tilde{\psi}_{j0}(x,y,z)=\tilde{\Psi}_{j0}(x,y,z)\exp[i\chi_{j0}(x,y,z)],$
(11)

are the solutions to the unperturbed linear propagation equation with the
initial condition (2) with unit amplitude. From Eqs. (10) and (11), it follows
that

$\displaystyle\Psi_{j0}(x,y,z)\simeq A_{j}(z)\tilde{\Psi}_{j0}(x,y,z).$ (12)

Using Eqs. (4) and (10), we find that the total electric fields of the beams
can be approximated by

$\displaystyle\psi_{j}(x,y,z)\simeq
A_{j}(z)\tilde{\psi}_{j0}(x,y,z)+\phi_{j}(x,y,z).$ (13)

Note that the approximate expressions (10), (12), and (13) are used both
inside and outside of the collision interval. In addition, the dynamics of the
$A_{j}(z)$ that is associated with single-beam propagation is described in
Appendix B.

The collision-induced amplitude shift of beam 1 is calculated from the
collision-induced change in $\Phi_{1}$ in the collision interval,
$\Delta\Phi_{1}(x,y,z_{c})=\Phi_{1}(x,y,z_{c}+\Delta
z_{c}/2)-\Phi_{1}(x,y,z_{c}-\Delta z_{c}/2)$. To calculate
$\Delta\Phi_{1}(x,y,z_{c})$, we substitute the approximations (12) for the
$\Psi_{j0}$ into Eq. (9), and integrate with respect to $z$ over the collision
interval. This calculation yields

$\displaystyle\\!\\!\\!\\!\\!\\!\Delta\Phi_{1}(x,y,z_{c})\\!=\\!-2\epsilon_{3}\\!\\!\int_{z_{c}-\Delta
z_{c}/2}^{z_{c}+\Delta
z_{c}/2}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dz^{\prime}A_{1}(z^{\prime})A_{2}^{2}(z^{\prime})$
$\displaystyle\times\tilde{\Psi}_{10}(x,y,z^{\prime})\tilde{\Psi}_{20}^{2}(x,y,z^{\prime}).$
(14)

Note that $\tilde{\Psi}_{20}$ is the only function in the integrand on the
right hand side of Eq. (14) that contains fast variations with respect to $z$,
which are of order 1. Therefore, we can approximate the other functions
$A_{1}(z)$, $A_{2}(z)$, and $\tilde{\Psi}_{10}(x,y,z)$ by $A_{1}(z_{c}^{-})$,
$A_{2}(z_{c}^{-})$, and $\tilde{\Psi}_{10}(x,y,z_{c})$, where
$A_{j}(z_{c}^{-})$ is the limit from the left of $A_{j}$ at $z_{c}$.
Furthermore, in calculating the integral, we can take into account in an exact
manner only the fast dependence of $\tilde{\Psi}_{20}$ on $z$, i.e., the
dependence on $z$ that is contained in the factors
$\tilde{x}=x-x_{20}-d_{11}z$, and replace $z$ by $z_{c}$ everywhere else in
the expression for $\tilde{\Psi}_{20}$. We denote this approximation of
$\tilde{\Psi}_{20}(x,y,z)$ by $\bar{\Psi}_{20}(\tilde{x},y,z_{c})$. Carrying
out all the aforementioned approximations in Eq. (14), we obtain:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\Delta\Phi_{1}(x,y,z_{c})\\!=\\!-2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})\tilde{\Psi}_{10}(x,y,z_{c})$
$\displaystyle\times\\!\\!\int_{z_{c}-\Delta z_{c}/2}^{z_{c}+\Delta
z_{c}/2}\\!\\!\\!\\!\\!\\!\\!\\!dz^{\prime}\bar{\Psi}_{20}^{2}(x-x_{20}-d_{11}z^{\prime},y,z_{c}).$
(15)

We assume that the integrand on the right hand side of Eq. (15) is sharply
peaked at a small interval around the collision distance $z_{c}$. Under this
assumption, we can extend the integral’s limits to $-\infty$ and $\infty$. We
also change the integration variable from $z^{\prime}$ to
$\tilde{x}=x-x_{20}-d_{11}z^{\prime}$ and obtain:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\Delta\Phi_{1}(x,y,z_{c})\\!=\\!-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}\tilde{\Psi}_{10}(x,y,z_{c})\\!\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!d\tilde{x}\bar{\Psi}_{20}^{2}(\tilde{x},y,z_{c}).$
(16)

We see that the $y$ dependence of beam 2 at $z=z_{c}$ affects the $y$
dependence of $\Delta\Phi_{1}(x,y,z_{c})$, while the $x$ dependence of beam 2
does not affect the $x$ dependence of $\Delta\Phi_{1}(x,y,z_{c})$. Thus,
inside the collision interval, the beam shape in the longitudinal direction is
preserved, while the beam shape in the transverse direction is changed by the
collision. We also point out that the collision-induced change in the beam
shape is an effect that exists only in spatial dimension higher than 1.
Indeed, it was shown in Refs. PNH2017B ; NHP2020  that in the one-dimensional
case, the beam shape is preserved by the collision within the leading order of
the perturbative calculation.

In Appendix A, we show that the collision-induced amplitude shift of beam 1
$\Delta A_{1}^{(c)}$ is related to $\Delta\Phi_{1}(x,y,z_{c})$ by:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\Delta
A_{1}^{(c)}=C_{p1}^{-1}\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\;\tilde{\Psi}_{10}(x,y,z_{c})\Delta\Phi_{1}(x,y,z_{c}),$
(17)

where

$\displaystyle
C_{p1}=\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\;\tilde{\Psi}_{10}^{2}(x,y,0).$
(18)

Substituting Eq. (16) into Eq. (17), we obtain the following expression for
the collision-induced amplitude shift of beam 1 for the general initial
condition (2):

$\displaystyle\\!\\!\\!\\!\Delta
A_{1}^{(c)}=-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{C_{p1}|d_{11}|}$
$\displaystyle\times\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\;\tilde{\Psi}_{10}^{2}(x,y,z_{c})\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!d\tilde{x}\;\bar{\Psi}_{20}^{2}(\tilde{x},y,z_{c}).$
(19)

####  II.2.3 Dynamics of $\phi_{1}(x,y,z)$ in the post-collision interval

In the post collision interval, i.e., for $z>z_{c}+\Delta z_{c}/2$, the two
beams are no longer overlapping. As a result, the inter-beam interaction terms
$-2i\epsilon_{3}|\psi_{2}|^{2}\psi_{1}$ and
$-2i\epsilon_{3}|\psi_{1}|^{2}\psi_{2}$ are negligible in this interval.
Therefore, in the leading order of the perturbative calculation, the equation
describing the dynamics of $\phi_{1}(x,y,z)$ in the post-collision interval is
the unperturbed linear propagation equation

$\displaystyle
i\partial_{z}\phi_{1}+\partial_{x}^{2}\phi_{1}+\partial_{y}^{2}\phi_{1}=0.$
(20)

To find the initial condition for Eq. (20), we first note that for
$|d_{11}|\gg 1$, $\Delta\Phi_{1}(x,y,z_{c})$ can be written as

$\displaystyle\Delta\Phi_{1}(x,y,z_{c})\simeq\Phi_{1}(x,y,z_{c}^{+})-\Phi_{1}(x,y,z_{c}^{-})\simeq\Phi_{1}(x,y,z_{c}^{+}),$
(21)

where $\Phi_{1}(x,y,z_{c}^{+})$ is the limit from the right of
$\Phi_{1}(x,y,z)$ at $z=z_{c}$. Thus, using the relation
$\phi_{1}(x,y,z)=\Phi_{1}(x,y,z)\exp[i\chi_{10}(x,y,z)]$, the initial
condition for Eq. (20) is:

$\displaystyle\phi_{1}(x,y,z_{c}^{+})=\Phi_{1}(x,y,z_{c}^{+})\exp[i\chi_{10}(x,y,z_{c})],$
(22)

where $\Phi_{1}(x,y,z_{c}^{+})$ is given by Eq. (16). The solution of Eq. (20)
with the IC (22) can be written as

$\displaystyle\phi_{1}(x,y,z)={\cal
F}^{-1}\left(\hat{\phi}_{1}(k_{1},k_{2},z_{c}^{+})\exp[-i(k_{1}^{2}+k_{2}^{2})(z-z_{c})]\right),$
(23)

where $\hat{\phi}_{1}(k_{1},k_{2},z_{c}^{+})={\cal
F}\left(\phi_{1}(x,y,z_{c}^{+})\right)$, and ${\cal F}$ and ${\cal F}^{-1}$
are the Fourier transform and the inverse Fourier transform with respect to
$x$ and $y$. In addition, the total electric field of beam 1 in the post-
collision interval is given by Eq. (13), where $A_{1}(z)$ is given by Eq.
(111) in Appendix B.

###  II.3 Calculation of the collision-induced changes in the beam shape and
amplitude for a separable initial condition

####  II.3.1 Introduction

We now describe the collision-induced dynamics in the important case, where
the initial condition for both beams is separable, i.e., it is given by Eq.
(3). This case is of special importance for two main reasons. First, this
initial condition describes the output electric field from many types of
lasers Siegman86 ; Kogelnik66 . Second, in this case, it is possible to
further simplify the expressions for the collision-induced changes of the beam
shape and amplitude, and in this manner, obtain deeper insight into the
collision dynamics.

It is straightforward to show that the solutions of the unperturbed linear
propagation equation with the separable initial condition (3) and with unit
amplitude can be written as

$\displaystyle\tilde{\psi}_{j0}(x,y,z)=g_{j}^{(x)}(x,z)g_{j}^{(y)}(y,z)\exp(i\alpha_{j0})=$
$\displaystyle
G_{j}^{(x)}(x,z)G_{j}^{(y)}(y,z)\exp\left\\{i\left[\chi_{j0}^{(x)}(x,z)+\chi_{j0}^{(y)}(y,z)+\alpha_{j0}\right]\right\\},$
(24)

where

$\displaystyle
g_{1}^{(x)}(x,z)=(2\pi)^{-1/2}\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!dk_{1}\hat{f}_{1}^{(x)}(k_{1})\exp[-ik_{1}^{2}z+ik_{1}x],$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!g_{2}^{(x)}(x,z)=(2\pi)^{-1/2}\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!dk_{1}\hat{f}_{2}^{(x)}(k_{1})\exp[-id_{11}k_{1}z-ik_{1}^{2}z+ik_{1}x],$
(25) $\displaystyle
g_{j}^{(y)}(y,z)=(2\pi)^{-1/2}\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!dk_{2}\hat{f}_{j}^{(y)}(k_{2})\exp[-ik_{2}^{2}z+ik_{2}y],$
(26)

and $G_{j}^{(x)}$, $G_{j}^{(y)}$, $\chi_{j0}^{(x)}$, and $\chi_{j0}^{(y)}$ are
real-valued. The functions $\hat{f}_{j}^{(x)}$ and $\hat{f}_{j}^{(y)}$ in Eqs.
(25) and (26) are defined by:

$\displaystyle\hat{f}_{j}^{(x)}(k_{1})={\cal
F}\left(h_{j}^{(x)}[(x-x_{j0})/W_{j0}^{(x)}]\right),$ (27)

and

$\displaystyle\hat{f}_{j}^{(y)}(k_{2})={\cal
F}\left(h_{j}^{(y)}[(y-y_{j0})/W_{j0}^{(y)}]\right).$ (28)

Using Eqs. (11) and (24), we obtain

$\displaystyle\tilde{\Psi}_{j0}(x,y,z)=G_{j}^{(x)}(x,z)G_{j}^{(y)}(y,z),$ (29)

and

$\displaystyle\chi_{j0}(x,y,z)=\chi_{j0}^{(x)}(x,z)+\chi_{j0}^{(y)}(y,z)+\alpha_{j0}.$
(30)

In addition, using the conservation of the total energy for the unperturbed
linear propagation equation, the definitions of $G_{j}^{(x)}$ and
$G_{j}^{(y)}$, and the initial condition (3), we obtain

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!dx\,G_{j}^{(x)2}(x,z)=\int_{-\infty}^{\infty}\\!\\!\\!\\!dx\,G_{j}^{(x)2}(x,0)=W_{j0}^{(x)}\int_{-\infty}^{\infty}\\!\\!\\!\\!ds\,h_{j}^{(x)2}(s)=W_{j0}^{(x)}c_{pj}^{(x)},$
(31)

and

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!dy\,G_{j}^{(y)2}(y,z)=\int_{-\infty}^{\infty}\\!\\!\\!\\!dy\,G_{j}^{(y)2}(y,0)=W_{j0}^{(y)}\int_{-\infty}^{\infty}\\!\\!\\!\\!ds\,h_{j}^{(y)2}(s)=W_{j0}^{(y)}c_{pj}^{(y)},$
(32)

where $c_{pj}^{(x)}$ and $c_{pj}^{(y)}$ are constants.

####  II.3.2 Collision-induced effects in the collision interval

We first obtain the general expression for $\Delta\Phi_{1}(x,y,z_{c})$ for an
initial condition that is separable for both beams. For this purpose, we note
that from the definition of $\bar{\Psi}_{20}(\tilde{x},y,z)$ it follows that
at $z=z_{c}$,
$\bar{\Psi}_{20}(\tilde{x},y,z_{c})=\tilde{\Psi}_{20}(x,y,z_{c})$. Using this
relation along with Eq. (29), we obtain
$\bar{\Psi}_{20}(\tilde{x},y,z_{c})=G_{2}^{(x)}(x,z_{c})G_{2}^{(y)}(y,z_{c})$.
It follows that:

$\displaystyle\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!d\tilde{x}\,\bar{\Psi}_{20}^{2}(\tilde{x},y,z_{c})=G_{2}^{(y)2}(y,z_{c})\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\,G_{2}^{(x)2}(x,z_{c}).$
(33)

Employing the conservation law (31) in Eq. (33), we obtain:

$\displaystyle\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!d\tilde{x}\,\bar{\Psi}_{20}^{2}(\tilde{x},y,z_{c})=c_{p2}^{(x)}W_{20}^{(x)}G_{2}^{(y)2}(y,z_{c}).$
(34)

Substitution of Eq. (34) into Eq. (16) yields the following expression for
$\Delta\Phi_{1}(x,y,z_{c})$, which is valid for an initial condition that is
separable for beam 2:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\Delta\Phi_{1}(x,y,z_{c})\\!=\\!-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}c_{p2}^{(x)}W_{20}^{(x)}G_{2}^{(y)2}(y,z_{c})\tilde{\Psi}_{10}(x,y,z_{c}).$
(35)

Equation (35) is valid for a general initial condition for beam 1. When the
initial condition for beam 1 is also separable, Eq. (35) takes the form

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\Delta\Phi_{1}(x,y,z_{c})\\!=\\!-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}$
$\displaystyle\times
c_{p2}^{(x)}W_{20}^{(x)}G_{1}^{(x)}(x,z_{c})G_{1}^{(y)}(y,z_{c})G_{2}^{(y)2}(y,z_{c}).$
(36)

We see that as in the case of the general initial condition (2), the shape of
the beam in the longitudinal direction does not change inside of the collision
interval. Moreover, it follows from Eq. (36) that for a separable initial
condition, the shape of the beam in the longitudinal direction is not changed
by the collision at all, i.e., for any $z>z_{c}$ (within the leading order of
the perturbative calculation). Indeed, for $|d_{11}|\gg 1$, Eq. (36) is also
the initial condition for the dynamics of $\phi_{1}(x,y,z)$ in the post-
collision interval. We observe that this initial condition is separable. In
addition, at $z=z_{c}$ the $x$ dependences of $\phi_{1}$ and
$\tilde{\psi}_{10}$ are identical. Since in the post-collision region
$\phi_{1}$ and $\tilde{\psi}_{10}$ satisfy the same linear propagation
equation with separable initial conditions, which have the same dependence on
$x$, the $x$ dependences of $\phi_{1}$ and $\tilde{\psi}_{10}$ remain
identical for any $z>z_{c}$. Thus, for a separable initial condition, the
shape of the beam in the longitudinal direction is not changed by the
collision at all.

We now turn to obtain the expression for $\Delta A_{1}^{(c)}$ for an initial
condition that is separable for both beams. Using the conservation of the
total energy and Eqs. (31) and (32), we obtain
$C_{p1}=c_{p1}^{(x)}c_{p1}^{(y)}W_{10}^{(x)}W_{10}^{(y)}$. In addition, using
Eqs. (29) and (34), we find

$\displaystyle\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\;\tilde{\Psi}_{10}^{2}(x,y,z_{c})\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!d\tilde{x}\;\bar{\Psi}_{20}^{2}(\tilde{x},y,z_{c})=$
$\displaystyle
c_{p1}^{(x)}c_{p2}^{(x)}W_{10}^{(x)}W_{20}^{(x)}\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\,G_{1}^{(y)2}(y,z_{c})G_{2}^{(y)2}(y,z_{c}).$
(37)

Substituting Eq. (37) and the expression for $C_{p1}$ into Eq. (19), we obtain
the following expression for the collision-induced amplitude shift for a
separable initial condition:

$\displaystyle\\!\\!\\!\\!\Delta
A_{1}^{(c)}=-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}$
$\displaystyle\times\frac{c_{p2}^{(x)}W_{20}^{(x)}}{c_{p1}^{(y)}W_{10}^{(y)}}\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\,G_{1}^{(y)2}(y,z_{c})G_{2}^{(y)2}(y,z_{c}).$
(38)

Note that the expression for $\Delta A_{1}^{(c)}$ has the form

$\displaystyle\\!\\!\\!\\!\Delta A_{1}^{(c)}=-(\mbox{overall
factor})\times(\mbox{longitudinal factor})\times(\mbox{transverse factor}),$
(39)

where the overall factor is equal to
$2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})/|d_{11}|$, and the
longitudinal factor is $c_{p2}^{(x)}W_{20}^{(x)}$. This form of the expression
for $\Delta A_{1}^{(c)}$ is expected to be valid for a general spatial
dimension, when the initial condition is separable in the longitudinal
direction for both beams. It is interesting to note that Eq. (38) is a
generalization of the equation obtained for a fast two-pulse collision in
spatial dimension 1. Indeed, using the notation of the current paper, the
latter equation, which is Eq. (19) in Ref. NHP2020 , can be written as:

$\displaystyle\Delta
A_{1}^{(c)(1D)}=-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}c_{p2}^{(x)}W_{20}^{(x)}.$
(40)

We observe that the overall and longitudinal factors in the equation for the
amplitude shift in spatial dimension 1 have the same form as the overall and
longitudinal factors in spatial dimension 2, while the transverse factor in
the one-dimensional case is equal to 1. We also observe that the longitudinal
factor $c_{p2}^{(x)}W_{20}^{(x)}$ is universal in the sense that it does not
depend on the exact details of the initial pulse shapes and on the collision
distance $z_{c}$. In contrast, the transverse factor is not universal since it
does depend on the details of the initial pulse shapes and on the collision
distance. Therefore, the universality of the expression for $\Delta
A_{1}^{(c)}$ in the one-dimensional case, which was first demonstrated in Ref.
NHP2020 , is extended to spatial dimension 2 (and to spatial dimension $n$),
but in a somewhat restricted manner. That is, in the two-dimensional (and the
$n$-dimensional) case, only the overall and longitudinal parts of the
expression for $\Delta A_{1}^{(c)}$ are universal, and this is true when the
initial condition is separable in the longitudinal direction for both beams.

####  II.3.3 Dynamics of $\phi_{1}(x,y,z)$ in the post-collision interval

We now turn to analyze the dynamics of $\phi_{1}(x,y,z)$ in the post-collision
interval. This analysis is especially interesting for two main reasons. First,
we showed in subsection II.2 that the collision induces a change of the beam
shape in the transverse direction. Even though this effect exists for a
general initial condition, its simplest and clearest demonstration is realized
in the case of an initial condition that is separable for both beams.
Furthermore, since in both experiments and simulations the change in the beam
shape is measured in the post-collision interval, we must analyze the
evolution of the beam shape in this interval. Second, we claimed in section
2.3.2 that for a separable initial condition, the shape of the beam in the
longitudinal direction does not change at all due to the collision. This claim
can be directly proved by analyzing the dynamics of $\phi_{1}(x,y,z)$ in the
post-collision interval.

In the post-collision interval, $\phi_{1}$ satisfies the unperturbed linear
propagation equation (20). Using Eqs. (22), (29), and (36), we find that the
initial condition for Eq. (20) is

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\phi_{1}(x,y,z_{c}^{+})=-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}$
$\displaystyle\times
c_{p2}^{(x)}W_{20}^{(x)}g_{1}^{(x)}(x,z_{c})g_{1}^{(y)}(y,z_{c})G_{2}^{(y)2}(y,z_{c})\exp(i\alpha_{10}).$
(41)

This initial condition can be written as:

$\displaystyle\phi_{1}(x,y,z_{c}^{+})=-\tilde{a}_{1}(z_{c}^{-})g_{1}^{(x)}(x,z_{c})g_{12}^{(y)}(y,z_{c})\exp(i\alpha_{10}),$
(42)

where

$\displaystyle\tilde{a}_{1}(z_{c}^{-})=2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})c_{p2}^{(x)}W_{20}^{(x)}/|d_{11}|,$
(43)

and

$\displaystyle
g_{12}^{(y)}(y,z_{c})=g_{1}^{(y)}(y,z_{c})G_{2}^{(y)2}(y,z_{c}).$ (44)

The Fourier transform of the initial condition (42) is

$\displaystyle\hat{\phi}_{1}(k_{1},k_{2},z_{c}^{+})=-\tilde{a}_{1}(z_{c}^{-})\hat{g}_{1}^{(x)}(k_{1},z_{c})\hat{g}_{12}^{(y)}(k_{2},z_{c})\exp(i\alpha_{10}),$
(45)

where $\hat{g}_{1}^{(x)}$ and $\hat{g}_{12}^{(y)}$ are the Fourier transforms
of $g_{1}^{(x)}$ and $g_{12}^{(y)}$ with respect to $x$ and $y$, respectively.
Substituting Eq. (45) into Eq. (23), we obtain:

$\displaystyle\phi_{1}(x,y,z)=-\tilde{a}_{1}(z_{c}^{-}){\cal
F}^{-1}\left(\hat{g}_{1}^{(x)}(k_{1},z_{c})\exp[-ik_{1}^{2}(z-z_{c})]\right)$
$\displaystyle\times{\cal
F}^{-1}\left(\hat{g}_{12}^{(y)}(k_{2},z_{c})\exp[-ik_{2}^{2}(z-z_{c})]\right)\exp(i\alpha_{10}).$
(46)

Note that when the initial condition for beam 1 is separable,
$\hat{g}_{1}^{(x)}(k_{1},z_{c})\exp[-ik_{1}^{2}(z-z_{c})]$ is equal to
$\hat{g}_{1}^{(x)}(k_{1},z)$:

$\displaystyle\hat{g}_{1}^{(x)}(k_{1},z_{c})\exp[-ik_{1}^{2}(z-z_{c})]=\hat{g}_{1}^{(x)}(k_{1},0)\exp(-ik_{1}^{2}z_{c})\exp[-ik_{1}^{2}(z-z_{c})]=$
$\displaystyle\hat{g}_{1}^{(x)}(k_{1},0)\exp(-ik_{1}^{2}z)=\hat{g}_{1}^{(x)}(k_{1},z).$
(47)

Substituting this relation into Eq. (46), we obtain the expression for
$\phi_{1}(x,y,z)$ in the post-collision interval for a separable initial
condition:

$\displaystyle\phi_{1}(x,y,z)=-\tilde{a}_{1}(z_{c}^{-})g_{1}^{(x)}(x,z)$
$\displaystyle\times{\cal
F}^{-1}\left(\hat{g}_{12}^{(y)}(k_{2},z_{c})\exp[-ik_{2}^{2}(z-z_{c})]\right)\exp(i\alpha_{10}).$
(48)

We see that when the initial condition is separable for both beams, the $x$
dependences of $\phi_{1}(x,y,z)$ and $\tilde{\psi}_{10}(x,y,z)$ are identical
for $z>z_{c}$. Therefore, as argued in subsection 2.3.2, the shape of the beam
in the longitudinal direction does not change at all by the collision.
Furthermore, the calculation of the modified beam shape in the transverse
direction amounts to the calculation of the inverse Fourier transform of
$\hat{g}_{12}^{(y)}(k_{2},z_{c})\exp[-ik_{2}^{2}(z-z_{c})]$.

##  III Perturbative calculation and numerical simulations for new collisional
effects in spatial dimension 2

###  III.1 Introduction

We now use the perturbation method of subsections II.2 and II.3 along with
numerical simulations with Eq. (1) to demonstrate four important effects and
properties of the collision, which either exist only in spatial dimension
higher than 1, or are qualitatively different from their one-dimensional
counterparts. These four effects and properties are: (1) universality of the
longitudinal part in the expression for the collision-induced amplitude shift,
(2) the effect of partial beam overlap, (3) the effect of anisotropy in the
initial condition, (4) the collision-induced change in the beam shape in the
transverse direction. For each effect or property, we first use the
perturbation theory to obtain explicit formulas, which demonstrate the
collisional effect or property. Since these formulas are only approximate
expressions, which are based on a number of simplifying assumptions of the
perturbative calculation, it is important to check their validity by numerical
simulations with the perturbed linear propagation equation (1). Therefore, in
the current subsection, we also take on this important numerical
investigation, by carrying out extensive numerical simulations with Eq. (1),
and by comparing the simulations results with the approximate predictions of
the perturbation theory for each of the four collisional effects and
properties. We solve Eq. (1) numerically by the split-step method with
periodic boundary conditions Agrawal2001 ; Yang2010 .

###  III.2 Universality of the longitudinal part in the expression for the
amplitude shift

In subsection II.3, we showed that for a separable initial condition, the
longitudinal factor in the expression for the collision-induced amplitude
shift is universal in the sense that it does not depend on the exact details
of the initial beam shapes. In contrast, the transverse factor is not
universal, since it does depend on the details of the beam shapes and on the
collision distance. Thus, according to the perturbative calculation, the
universality of the expression for $\Delta A_{1}^{(c)}$ is extended from
spatial dimension 1 to higher spatial dimensions, but in a somewhat restricted
manner.

In the current subsection, we demonstrate the universality of the longitudinal
part in the expression for the collision-induced amplitude shift. For this
purpose, we first obtain explicit expressions for $\Delta A_{1}^{(c)}$ for two
initial beam shapes that have widely different dependences on the $x$
coordinate. Moreover, we verify the validity of the expressions for $\Delta
A_{1}^{(c)}$ by extensive numerical simulations with Eq. (1). This numerical
investigation is especially important, since it shows that the approximations
used in the perturbative calculation are indeed valid for widely different
beam shapes. In this manner, the extensive numerical simulations with Eq. (1)
help validate the universal nature of the longitudinal part in the expression
for $\Delta A_{1}^{(c)}$. The initial $x$ dependence for the first beam type
that we consider is Gaussian, i.e., it is rapidly decreasing with increasing
distance from the beam center. In contrast, the initial $x$ dependence for the
second beam type that we consider is given by a Cauchy-Lorentz distribution,
i.e., it decreases slowly (as a power-law) with increasing distance from the
beam center. The initial beam profile in the transverse direction is taken as
Gaussian, as this choice enables the explicit calculation of the integral with
respect to $y$ on the right hand side of Eq. (38). Therefore, the two initial
conditions for the two-beam collision problem are

$\displaystyle\psi_{1}(x,y,0)=A_{1}(0)\exp\left[-\frac{x^{2}}{2W^{(x)2}_{10}}-\frac{y^{2}}{2W^{(y)2}_{10}}+i\alpha_{10}\right],$
$\displaystyle\psi_{2}(x,y,0)=A_{2}(0)\exp\left[-\frac{(x-x_{20})^{2}}{2W^{(x)2}_{20}}-\frac{y^{2}}{2W^{(y)2}_{20}}+i\alpha_{20}\right],$
(49)

for Gaussian beams, and

$\displaystyle\psi_{1}(x,y,0)=A_{1}(0)\left[1+\frac{2x^{4}}{W^{(x)4}_{10}}\right]^{-1}\exp\left[-\frac{y^{2}}{2W^{(y)2}_{10}}+i\alpha_{10}\right],$
$\displaystyle\psi_{2}(x,y,0)=A_{2}(0)\left[1+\frac{2(x-x_{20})^{4}}{W^{(x)4}_{20}}\right]^{-1}\exp\left[-\frac{y^{2}}{2W^{(y)2}_{20}}+i\alpha_{20}\right],$
(50)

for Cauchy-Lorentz-Gaussian beams.

Let us obtain the expression for the collision-induced amplitude shift for the
initial conditions (49) and (50). From Eq. (134) in Appendix D it follows that
for both initial conditions

$\displaystyle
G_{j}^{(y)}(y,z_{c})=\frac{W_{j0}^{(y)}}{(W_{j0}^{(y)4}+4z_{c}^{2})^{1/4}}\exp\left[-\frac{W_{j0}^{(y)2}y^{2}}{2(W_{j0}^{(y)4}+4z_{c}^{2})}\right],$
(51)

where $j=1,2$. In addition, $c_{p1}^{(y)}=\pi^{1/2}$ for both initial
conditions, $c_{p2}^{(x)}=\pi^{1/2}$ for Gaussian beams, and
$c_{p2}^{(x)}=3\pi/2^{11/4}$ for Cauchy-Lorentz-Gaussian beams. We now
substitute Eq. (51) along with the values of $c_{p1}^{(y)}$ and $c_{p2}^{(x)}$
into Eq. (38), and perform the integration with respect to $y$. This
calculation yields the following expression for $\Delta A_{1}^{(c)}$:

$\displaystyle\Delta
A_{1}^{(c)}=-\frac{2b\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}$
$\displaystyle\times\frac{W_{10}^{(y)}W_{20}^{(x)}W_{20}^{(y)2}}{(W_{10}^{(y)2}+W_{20}^{(y)2})^{1/2}(W_{10}^{(y)2}W_{20}^{(y)2}+4z_{c}^{2})^{1/2}},$
(52)

where $b=\pi^{1/2}$ for Gaussian beams, and $b=3\pi/2^{11/4}$ for Cauchy-
Lorentz-Gaussian beams. The longitudinal part in the expression for $\Delta
A_{1}^{(c)}$, $c_{p2}^{(x)}W_{20}^{(x)}=bW_{20}^{(x)}$, is clearly universal.
In contrast, the transverse part in the expression, which is given by:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{transverse
factor}=\frac{W_{10}^{(y)}W_{20}^{(y)2}}{(W_{10}^{(y)2}+W_{20}^{(y)2})^{1/2}(W_{10}^{(y)2}W_{20}^{(y)2}+4z_{c}^{2})^{1/2}},$
(53)

depends on $z_{c}$, and does not seem to have a simple universal form. One
aspect of the nonuniversal nature of the expression for $\Delta A_{1}^{(c)}$
in spatial dimension 2 is the deviation of the dependence on $|d_{11}|$ from
the $1/|d_{11}|$ scaling, which is observed in the one-dimensional case
PNH2017B ; NHP2020 , and also in fast collisions between NLS solitons in the
presence of nonlinear dissipation in spatial dimension 1 PNC2010 ; PC2012 .
Note that the collision distance $z_{c}$ satisfies
$z_{c}=(x_{10}-x_{20})/d_{11}$. Therefore, the deviation of the $|d_{11}|$
dependence of $\Delta A_{1}^{(c)}$ from the $1/|d_{11}|$ scaling is due to the
term $4(x_{10}-x_{20})^{2}/d_{11}^{2}$ inside the factor
$(W_{10}^{(y)2}W_{20}^{(y)2}+4(x_{10}-x_{20})^{2}/d_{11}^{2})^{1/2}$ on the
right hand side of Eq. (52). It is useful to define the quantity $\Delta
A_{1}^{(c)(s)}$, which is the approximate expression for the amplitude shift
that is obtained from the full expression by neglecting the
$(x_{10}-x_{20})^{2}/d_{11}^{2}$ term. Carrying out the latter approximation,
we find:

$\displaystyle\Delta
A_{1}^{(c)(s)}=-\frac{2b\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}\frac{W_{20}^{(x)}W_{20}^{(y)}}{(W_{10}^{(y)2}+W_{20}^{(y)2})^{1/2}}.$
(54)

Therefore, the difference $|\Delta A_{1}^{(c)}-\Delta A_{1}^{(c)(s)}|$ is a
measure for the deviation of the $d_{11}$ dependence of $\Delta A_{1}^{(c)}$
from the $1/|d_{11}|$ scaling observed in the one-dimensional case. Since in a
complete collision $|x_{20}-x_{10}|\gg 1$, the term
$4(x_{10}-x_{20})^{2}/d_{11}^{2}$ is not necessarily small for intermediate
values of $|d_{11}|$. As a result, the deviation from the $1/|d_{11}|$ scaling
might be significant even for intermediate $|d_{11}|$ values.

To check the perturbation theory predictions for universality of the
longitudinal part in the expression for $\Delta A_{1}^{(c)}$, we carry out
numerical simulations with Eq. (1) with the two initial conditions (49) and
(50), which possess widely different initial beam profiles in the longitudinal
direction. The extensive simulations with these initial conditions provide a
careful test for the validity of the perturbation theory approximations for
widely different beam shapes, and in this manner, help confirm the
universality of the longitudinal part in the expression for $\Delta
A_{1}^{(c)}$. We carry out the simulations for $d_{11}$ values in the
intervals $4\leq|d_{11}|\leq 60$. For concreteness, we present the results of
the simulations with $\epsilon_{3}=0.01$. The parameter values of the initial
conditions (49) and (50) are $A_{j}(0)=1$, $\alpha_{j0}=0$, $x_{20}=\pm 20$,
$W_{10}^{(x)}=3$, $W_{10}^{(y)}=5$, $W_{20}^{(x)}=4$, and $W_{20}^{(y)}=6$.
The final propagation distance is $z_{f}=2z_{c}=-2x_{20}/d_{11}$. The values
of $x_{20}$ and $z_{f}$ ensure that the centers of the two beams are well
separated at $z=0$ and at $z=z_{f}$. We point out that results similar to the
ones described in the current subsection are obtained in simulations with
other parameter values. For each initial condition, we compare the dependence
of $\Delta A_{1}^{(c)}$ on $d_{11}$ obtained in the simulations with the
perturbation theory prediction of Eq. (52), and with the more crude
approximation $\Delta A_{1}^{(c)(s)}$ of Eq. (54). We also discuss the
behavior of the relative errors in the approximation of $\Delta A_{1}^{(c)}$
(in percentage), which are defined by $E_{r}^{(1)}=|\Delta
A_{1}^{(c)(num)}-\Delta A_{1}^{(c)(th)}|\times 100/|\Delta A_{1}^{(c)(th)}|$
and $E_{r}^{(2)}=|\Delta A_{1}^{(c)(num)}-\Delta A_{1}^{(c)(s)}|\times
100/|\Delta A_{1}^{(c)(s)}|$, respectively.

We start by discussing the results of the simulations for Gaussian beams,
which represent beams with rapidly decaying tails. The initial beam shapes
$|\psi_{j}(x,y,0)|$, and the beam shapes $|\psi_{j}(x,y,z)|$ obtained in the
simulation with $d_{11}=10$ at the intermediate distance $z_{i}=2.4>z_{c}$
zi_values , and at the final distance $z_{f}=4$ are shown in Fig. 1. We
observe that the beams undergo broadening due to diffraction without
developing significant tails. In addition, the maximum values of
$|\psi_{j}(x,y,z)|$ decrease with increasing $z$, mainly due to diffraction.
The dependence of $\Delta A_{1}^{(c)}$ on $d_{11}$ that is obtained in the
simulations is shown in Fig. 2 along with the analytic prediction of Eq. (52)
and the more crude approximation $\Delta A_{1}^{(c)(s)}$ of Eq. (54). We
observe that despite the diffraction-induced broadening of the beams, the
agreement between the simulations result and the analytic prediction of Eq.
(52) is very good. In particular, the relative error $E_{r}^{(1)}$ is less
than $3.5\%$ for $10\leq|d_{11}|\leq 60$ and less than $5.1\%$ for
$4\leq|d_{11}|<10$. We also observe that the two theoretical curves for
$\Delta A_{1}^{(c)}$ and $\Delta A_{1}^{(c)(s)}$ are close to each other. As a
result, the relative error $E_{r}^{(2)}$ is only somewhat larger than
$E_{r}^{(1)}$. More specifically, the value of $E_{r}^{(2)}$ is smaller than
$4.4\%$ for $10\leq|d_{11}|\leq 60$ and smaller than $10.0\%$ for
$4\leq|d_{11}|<10$. This means that the deviation of the $d_{11}$ dependence
of $\Delta A_{1}^{(c)}$ from the $1/|d_{11}|$ scaling is not significant for
$|d_{11}|\geq 10$. However, it should be noted that the latter result is due
to the choice of the values of $W_{10}^{(y)}$ and $W_{20}^{(y)}$, i.e., for
smaller values of these parameters we can observe significantly larger
deviations from the $1/|d_{11}|$ scaling.

Figure 1: (Color online) Contour plots of the beam shapes $|\psi_{j}(x,y,z)|$
at $z=0$ (a), $z=z_{i}=2.4$ (b), and $z=z_{f}=4$ (c) in a fast collision
between two Gaussian beams with parameter values $\epsilon_{3}=0.01$ and
$d_{11}=10$. The plots represent the beam shapes obtained by numerical
solution of Eq. (1) with the initial condition (49). Figure 2: (Color online)
Dependence of the collision-induced amplitude shift of beam 1 $\Delta
A_{1}^{(c)}$ on the beam-steering coefficient $d_{11}$ in a fast collision
between two Gaussian beams for $\epsilon_{3}=0.01$. The red circles represent
the result obtained by numerical simulations with Eq. (1) with the initial
condition (49). The solid blue and dashed green curves represent the
theoretical predictions of Eqs. (52) and (54), respectively.

We now turn to discuss the simulations results for beams, whose tails exhibit
slow decay in the longitudinal direction. For such beams it is unclear if the
sharp-peak approximation, which is used in the derivation of Eq. (16) from Eq.
(14), is valid. Therefore, in this case, the numerical simulations of the two-
beam collision serve as an important check of the perturbation theory
predictions. We use the Cauchy-Lorentz-Gaussian beams of Eq. (50) as
prototypical examples for beams, whose tails exhibit slow (power-law) decay
with increasing distance from the beam center. Figure 3 shows the initial beam
shapes $|\psi_{j}(x,y,0)|$, and the beam shapes $|\psi_{j}(x,y,z)|$ obtained
in the simulation with $d_{11}=10$ at $z=z_{i}=2.4$, and at $z=z_{f}=4$. We
observe that the beams experience broadening and develop extended tails due to
diffraction. In addition, the maximum values of $|\psi_{j}(x,y,z)|$ decrease
with increasing $z$, mainly due to diffraction. The latter decrease is
especially noticeable for beam 1. This can be explained by noting that
$W_{10}^{(x)}<W_{20}^{(x)}$, and as a result, diffraction-induced beam
broadening and generation of extended tails are stronger for beam 1 compared
with beam 2. The $d_{11}$ dependence of $\Delta A_{1}^{(c)}$ obtained in the
simulations is shown in Fig. 4 together with the analytic predictions $\Delta
A_{1}^{(c)}$ and $\Delta A_{1}^{(c)(s)}$ of Eqs. (52) and (54). The agreement
between the simulations result for $\Delta A_{1}^{(c)}$ and the perturbation
theory prediction of Eq. (52) is very good despite the beam broadening and the
generation of extended beam tails. In particular, the value of the relative
error $E_{r}^{(1)}$ is smaller than $2.8\%$ for $10\leq|d_{11}|\leq 60$ and
smaller than $5.1\%$ for $4\leq|d_{11}|<10$. Note that these values are
comparable to the values of $E_{r}^{(1)}$ for collisions between Gaussian
beams. Thus, based on the results of our simulations, we conclude that the
longitudinal part in the expression for $\Delta A_{1}^{(c)}$,
$c_{p2}^{(x)}W_{20}^{(x)}$, is indeed universal in the sense that it is not
sensitive to the exact form of the initial beam shapes. We also note that the
value of $E_{r}^{(2)}$ is smaller than $4.3\%$ for $10\leq|d_{11}|\leq 60$ and
smaller than $9.2\%$ for $4\leq|d_{11}|<10$. Thus, the deviation of the
$d_{11}$ dependence of $\Delta A_{1}^{(c)}$ from the $1/|d_{11}|$ scaling is
not significant for $|d_{11}|\geq 10$ for the parameter values used in our
simulations.

Figure 3: (Color online) Contour plots of the beam shapes $|\psi_{j}(x,y,z)|$
at $z=0$ (a), $z=z_{i}=2.4$ (b), and $z=z_{f}=4$ (c) in a fast collision
between two Cauchy-Lorentz-Gaussian beams with parameter values
$\epsilon_{3}=0.01$ and $d_{11}=10$. The plots represent the beam shapes
obtained by numerical solution of Eq. (1) with the initial condition (50).
Figure 4: (Color online) Dependence of the collision-induced amplitude shift
of beam 1 $\Delta A_{1}^{(c)}$ on the beam-steering coefficient $d_{11}$ in a
fast collision between two Cauchy-Lorentz-Gaussian beams for
$\epsilon_{3}=0.01$. The red circles represent the result obtained by
numerical simulations with Eq. (1) with the initial condition (50). The solid
blue and dashed green curves represent the perturbation theory predictions of
Eqs. (52) and (54), respectively.

###  III.3 Fast collisions between partially overlapping beams

Another important property of a complete two-beam collision is related to the
relative location of the beam centers at the collision distance $z_{c}$. When
the $y$ coordinates of the beams are equal at $z_{c}$,
$y_{1}(z_{c})=y_{2}(z_{c})$, we say that the beams are fully overlapping at
$z_{c}$. In contrast, when the $y$ coordinates of the beams are not equal at
$z_{c}$, $y_{1}(z_{c})\neq y_{2}(z_{c})$, we say that the beams are only
partially overlapping at $z_{c}$. It is clear that complete collisions between
two partially overlapping beams exist only in spatial dimension higher than 1,
since in spatial dimension 1, the two beams are always fully overlapping at
$z_{c}$ in a complete fast collision. It is therefore interesting to employ
the perturbation theory for studying the effect of the partial overlap between
the colliding beams in 2D on the collision-induced amplitude shift. This
problem is of further interest, since we can use it for checking the
perturbation theory prediction for the transverse part in the expression for
the amplitude shift in a nontrivial setup. Thus, in the current subsection, we
investigate the dynamics of the amplitude shift in fast collisions between
partially overlapping beams both analytically and by numerical simulations
with Eq. (1).

To demonstrate the effects of partial overlap on the collision-induced
amplitude shift, we consider an initial condition in the form of two Gaussian
beams with different initial values of the $y$ coordinates of the beam
centers, $y_{10}\neq y_{20}$. For simplicity and without loss of generality,
we assume that the initial beam widths satisfy
$W_{10}^{(x)}=W_{20}^{(x)}\equiv W_{0}^{(x)}$ and
$W_{10}^{(y)}=W_{20}^{(y)}\equiv W_{0}^{(y)}$. Therefore, the initial
condition for the two-beam collision problem is given by:

$\displaystyle\psi_{j}(x,y,0)=A_{j}(0)\exp\left[-\frac{(x-x_{j0})^{2}}{2W^{(x)2}_{0}}-\frac{(y-y_{j0})^{2}}{2W^{(y)2}_{0}}+i\alpha_{j0}\right],$
(55)

for $j=1,2$. Since the initial condition is separable for both beams, we can
use Eq. (38) for calculating the collision-induced amplitude shift. For this
initial condition, $c_{p1}^{(y)}=c_{p2}^{(x)}=\pi^{1/2}$, and the functions
$G_{j}^{(y)}$ are given by Eq. (51) with
$W_{10}^{(y)}=W_{20}^{(y)}=W_{0}^{(y)}$. Substituting these expressions into
Eq. (38) and integrating with respect to $y$, we obtain:

$\displaystyle\Delta
A_{1}^{(c)}=-\frac{(2\pi)^{1/2}\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}$
$\displaystyle\times\frac{W_{0}^{(x)}W_{0}^{(y)2}}{(W_{0}^{(y)4}+4z_{c}^{2})^{1/2}}\exp\left[-\frac{W_{0}^{(y)2}(y_{20}-y_{10})^{2}}{2(W_{0}^{(y)4}+4z_{c}^{2})}\right].$
(56)

Thus, the effect of partial beam overlap on the amplitude shift is contained
in the transverse part of the expression for $\Delta A_{1}^{(c)}$:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{transverse
factor}=\frac{W_{0}^{(y)2}}{2^{1/2}(W_{0}^{(y)4}+4z_{c}^{2})^{1/2}}\exp\left[-\frac{W_{0}^{(y)2}(y_{20}-y_{10})^{2}}{2(W_{0}^{(y)4}+4z_{c}^{2})}\right].$
(57)

We see that $\Delta A_{1}^{(c)}$ is a Gaussian function of the separation
between the beam centers at $z_{c}$, $y_{20}-y_{10}$. The width of the
Gaussian function is equal to $2^{1/2}W^{(y)}(z_{c})$, where
$W^{(y)}(z_{c})=(W_{0}^{(y)2}+4z_{c}^{2}/W_{0}^{(y)2})^{1/2}$ is the width of
both beams in the transverse direction at $z_{c}$.

It is unclear if the approximations used by the perturbation theory hold when
the separation between the beam centers at $z_{c}$ is relatively large. For
this reason, it is important to check the predictions of Eqs. (56) and (57) by
numerical solution of Eq. (1). We take on this important numerical
investigation by carrying out simulations with Eq. (1) and by measuring the
dependence of $\Delta A_{1}^{(c)}$ on $y_{20}-y_{10}$. For brevity, we
describe the results of these simulations briefly without showing the
corresponding figures. The physical parameter values are $\epsilon_{3}=0.01$
and $d_{11}=20$. The initial values of the beam parameters in Eq. (55) are
$A_{j}(0)=1$, $\alpha_{j0}=0$, $x_{10}=0$, $x_{20}=-20$, $y_{10}=0$,
$W_{0}^{(x)}=4$, $W_{0}^{(y)}=5$, and the value of $y_{20}$ is varied in the
interval $-10\leq y_{20}\leq 10$. The final propagation distance is $z_{f}=2$,
and the beam centers are well separated at $z=0$ and at $z=z_{f}$. In addition
to $\Delta A_{1}^{(c)}$, we measure the relative error (in percentage)
$|\Delta A_{1}^{(c)(num)}-\Delta A_{1}^{(c)(th)}|\times 100/|\Delta
A_{1}^{(c)(th)}|$, where $\Delta A_{1}^{(c)(th)}$ is given by Eq. (56). The
agreement between the numerical simulations result and the perturbation theory
prediction is very good. More specifically, the relative error is smaller than
$7.9\%$ in the entire interval $-10\leq y_{20}\leq 10$. Thus, our numerical
simulations confirm the prediction of Eqs. (56) and (57) for Gaussian
dependence of $\Delta A_{1}^{(c)}$ and its transverse factor on
$y_{20}-y_{10}$. Furthermore, these simulations demonstrate that the
perturbation method of subsections II.2 and II.3 is applicable for fast
collisions between partially overlapping beams, even when the beam centers are
relatively far from each other at the collision distance $z_{c}$.

###  III.4 Dependence of the amplitude shift on the orientation angle between
the beams

Another phenomenon that exists only in spatial dimension higher than 1 is
associated with direction dependent collision-induced effects, i.e., with
effects that exist due to some anisotropy in the system. In particular, we are
interested in studying the effects of anisotropy in the initial condition. In
a simple setup, this anisotropy can be characterized by a single angle, e.g.,
the angle $\theta_{0}$ between a “preferred” direction in the initial
condition and the $x$ axis. To illustrate this situation, consider the case
where the initial shape of beam 1 $|\psi_{1}(x,y,0)|$ is wider along one
direction that we denote by $x^{\prime}$, and narrower along the perpendicular
direction that we denote by $y^{\prime}$. We can then define the angle
$\theta_{0}$, as the angle that the $x^{\prime}$ axis forms with the $x$ axis
of our coordinate system. Thus, $\theta_{0}$ is the angle between the relative
velocity vector (between the beam centers) and the $x^{\prime}$ axis of beam
1. In addition, if beam 2 is circularly symmetric, or is elongated along the
$x$ or the $y$ axes, then $\theta_{0}$ can also be regarded as the orientation
angle between the two beams. An important question about the collision
dynamics in this anisotropic setup concerns the dependence of the collision-
induced amplitude shift $\Delta A_{1}^{(c)}$ on the orientation angle
$\theta_{0}$. In the current subsection, we address this important question by
both analytic calculations and numerical simulations with Eq. (1).

We consider the following anisotropic collision setup, which consists of two
initially well separated Gaussian beams. In this setup, beam 1 is elongated
along its $x^{\prime}$ axis, which forms an angle $\theta_{0}$ with the $x$
axis, while beam 2 is circular in the $xy$ plane. Figure 5 shows the contour
plot of $|\psi_{j}(x,y,0)|$ for this initial condition in the case where
$\theta_{0}=\pi/4$. The initial condition can be written as:

$\displaystyle\\!\\!\\!\\!\\!\\!\psi_{1}^{\prime}(x^{\prime},y^{\prime},0)=A_{1}(0)\exp\left[-\frac{x^{\prime
2}}{2W^{(x)2}_{10}}-\frac{y^{\prime 2}}{2W^{(y)2}_{10}}+i\alpha_{10}\right],$
(58)

and

$\displaystyle\\!\\!\\!\\!\\!\\!\psi_{2}(x,y,0)=A_{2}(0)\exp\left[-\frac{(x-x_{20})^{2}}{2W^{2}_{20}}-\frac{y^{2}}{2W^{2}_{20}}+i\alpha_{20}\right],$
(59)

where $\psi_{1}^{\prime}(x^{\prime},y^{\prime},z)$ denotes the electric field
of beam 1 in the $(x^{\prime},y^{\prime},z)$ coordinate system,
$W_{10}^{(x)}>W_{10}^{(y)}$, and

$\displaystyle x^{\prime}=x\cos\theta_{0}+y\sin\theta_{0},$ $\displaystyle
y^{\prime}=y\cos\theta_{0}-x\sin\theta_{0}.$ (60)

Substituting relation (60) into Eq. (58), we obtain:

$\displaystyle\psi_{1}(x,y,0)=A_{1}(0)\exp\left[-B_{1}x^{2}-B_{2}y^{2}-B_{3}xy+i\alpha_{10}\right],$
(61)

where

$\displaystyle
B_{1}=\frac{\cos^{2}\theta_{0}}{2W_{10}^{(x)2}}+\frac{\sin^{2}\theta_{0}}{2W_{10}^{(y)2}},$
$\displaystyle
B_{2}=\frac{\sin^{2}\theta_{0}}{2W_{10}^{(x)2}}+\frac{\cos^{2}\theta_{0}}{2W_{10}^{(y)2}},$

and

$\displaystyle
B_{3}=\left(\frac{1}{W_{10}^{(x)2}}-\frac{1}{W_{10}^{(y)2}}\right)\sin\theta_{0}\cos\theta_{0}.$

Notice that the initial condition for beam 1 is not separable in the $(x,y,z)$
coordinate system. Therefore, the investigation described in the current
subsection also provides an example for collision-induced dynamics in a
collision with a nonseparable initial condition.

Figure 5: (Color online) A contour plot of the initial beam shapes
$|\psi_{j}(x,y,0)|$ for the anisotropic initial condition of Eqs. (59) and
(61). In this example, the initial beam widths are $W_{10}^{(x)}=8$,
$W_{10}^{(y)}=2$, and $W_{20}=2$, and the orientation angle is
$\theta_{0}=\pi/4$.

The initial condition for beam 1 is nonseparable, and therefore we need to
calculate $\Delta A_{1}^{(c)}$ by using the general expression, which is given
by Eq. (19). It is straightforward to show that for the current setup,
$C_{p1}=\pi W_{10}^{(x)}W_{10}^{(y)}$. Furthermore, since the initial
condition for beam 2 is separable, we can use Eq. (34), where
$c_{p2}^{(x)}=\pi^{1/2}$. Substitution of these relations into Eq. (19) yields

$\displaystyle\\!\\!\\!\\!\Delta
A_{1}^{(c)}=-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{\pi^{1/2}|d_{11}|}\frac{W_{20}}{W_{10}^{(x)}W_{10}^{(y)}}$
$\displaystyle\times\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\,G_{2}^{(y)2}(y,z_{c})\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\;\tilde{\Psi}_{10}^{2}(x,y,z_{c}).$
(62)

Since diffraction is isotropic, the unperturbed linear propagation equation
for beam 1 in the $(x^{\prime},y^{\prime},z)$ coordinate system has the same
form as in the $(x,y,z)$ coordinate system. Thus, the unperturbed propagation
equation for $\tilde{\psi}^{\prime}_{10}$ is:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!i\partial_{z}\tilde{\psi}^{\prime}_{10}+\partial_{x^{\prime}}^{2}\tilde{\psi}^{\prime}_{10}+\partial_{y^{\prime}}^{2}\tilde{\psi}^{\prime}_{10}=0.$
(63)

Therefore, we can calculate $\tilde{\Psi}_{10}(x,y,z)$ by solving Eq. (63)
with the initial condition (58) in the $(x^{\prime},y^{\prime},z)$ coordinate
system, and by using Eq. (60) to express the solution in the $(x,y,z)$
coordinate system. The solution of Eq. (63) with the Gaussian initial
condition (58) is described in Appendix D. Using Eqs. (131)-(135) in this
Appendix, we obtain:

$\displaystyle\tilde{\Psi}^{\prime}_{10}(x^{\prime},y^{\prime},z)=\frac{W_{10}^{(x)}W_{10}^{(y)}}{(W_{10}^{(x)4}+4z^{2})^{1/4}(W_{10}^{(y)4}+4z^{2})^{1/4}}$
$\displaystyle\times\exp\left[-\frac{W_{10}^{(x)2}x^{\prime
2}}{2(W_{10}^{(x)4}+4z^{2})}-\frac{W_{10}^{(y)2}y^{\prime
2}}{2(W_{10}^{(y)4}+4z^{2})}\right].$ (64)

Using the transformation relations (60) in Eq. (64), and using Eq. (51) with
$W_{20}^{(y)}=W_{20}$ for $G_{2}^{(y)}(y,z_{c})$, we obtain:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!G_{2}^{(y)2}(y,z_{c})\tilde{\Psi}_{10}^{2}(x,y,z_{c})=\frac{W_{10}^{(x)2}W_{10}^{(y)2}W_{20}^{2}}{(W_{10}^{(x)4}+4z_{c}^{2})^{1/2}(W_{10}^{(y)4}+4z_{c}^{2})^{1/2}(W_{20}^{4}+4z_{c}^{2})^{1/2}}$
$\displaystyle\times\exp\left[-b_{1}^{2}x^{2}-2b_{2}xy-b_{3}^{2}y^{2}\right],$
(65)

where

$\displaystyle
b_{1}=\left(\frac{W_{10}^{(x)2}\cos^{2}\theta_{0}}{W_{10}^{(x)4}+4z_{c}^{2}}+\frac{W_{10}^{(y)2}\sin^{2}\theta_{0}}{W_{10}^{(y)4}+4z_{c}^{2}}\right)^{1/2},$
(66) $\displaystyle
b_{2}=\left(\frac{W_{10}^{(x)2}}{W_{10}^{(x)4}+4z_{c}^{2}}-\frac{W_{10}^{(y)2}}{W_{10}^{(y)4}+4z_{c}^{2}}\right)\sin\theta_{0}\cos\theta_{0},$
(67)

and

$\displaystyle
b_{3}=\left(\frac{W_{20}^{2}}{W_{20}^{4}+4z_{c}^{2}}+\frac{W_{10}^{(x)2}\sin^{2}\theta_{0}}{W_{10}^{(x)4}+4z_{c}^{2}}+\frac{W_{10}^{(y)2}\cos^{2}\theta_{0}}{W_{10}^{(y)4}+4z_{c}^{2}}\right)^{1/2}.$
(68)

Substituting Eq. (65) into Eq. (62) and carrying out the double integration,
we obtain the following expression for the collision-induced amplitude shift:

$\displaystyle\Delta
A_{1}^{(c)}=-\frac{2\pi^{1/2}\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}$
$\displaystyle\times\frac{W_{10}^{(x)}W_{10}^{(y)}W_{20}^{3}}{(W_{10}^{(x)4}+4z_{c}^{2})^{1/2}(W_{10}^{(y)4}+4z_{c}^{2})^{1/2}(W_{20}^{4}+4z_{c}^{2})^{1/2}(b_{1}^{2}b_{3}^{2}-b_{2}^{2})^{1/2}},$
(69)

where

$\displaystyle
b_{1}^{2}b_{3}^{2}-b_{2}^{2}=\frac{W_{20}^{2}}{W_{20}^{4}+4z_{c}^{2}}\left(\frac{W_{10}^{(x)2}\cos^{2}\theta_{0}}{W_{10}^{(x)4}+4z_{c}^{2}}+\frac{W_{10}^{(y)2}\sin^{2}\theta_{0}}{W_{10}^{(y)4}+4z_{c}^{2}}\right)$
$\displaystyle+\frac{W_{10}^{(x)2}W_{10}^{(y)2}}{(W_{10}^{(x)4}+4z_{c}^{2})(W_{10}^{(y)4}+4z_{c}^{2})}.$
(70)

We see that even the relatively simple anisotropic setup of the two-beam
collision considered in the current subsection leads to a nontrivial
dependence of $\Delta A_{1}^{(c)}$ on the orientation angle $\theta_{0}$. This
nontrivial dependence on $\theta_{0}$ can also be associated with the
nonseparable nature of the initial condition for beam 1.

We check the predictions of Eq. (69) for the dependence of the collision-
induced amplitude shift on $\theta_{0}$ by carrying out numerical simulations
with Eq. (1) with the initial condition of Eqs. (59) and (61). The simulations
are carried out for $\theta_{0}$ values in the interval
$0\leq\theta_{0}\leq\pi/2$. The physical parameter values are
$\epsilon_{3}=0.01$ and $d_{11}=20$. The initial values of the beam parameters
in Eqs. (59) and (61) are $A_{j}(0)=1$, $\alpha_{j0}=0$, $x_{20}=-20$,
$W_{10}^{(x)}=8$, $W_{10}^{(y)}=2$, and $W_{20}=2$. The final propagation
distance is $z_{f}=2$, and therefore, the beam centers are well separated at
$z_{f}$. Figure 6 shows the initial beam shapes $|\psi_{j}(x,y,0)|$, and the
beam shapes $|\psi_{j}(x,y,z)|$ obtained in the numerical simulation with
$\theta_{0}=\pi/4$ at the intermediate distance $z_{i}=1.2$, and at $z_{f}=2$.
We observe that both beams experience broadening due to diffraction but do not
develop extended tails. The dependence of $\Delta A_{1}^{(c)}$ on the
orientation angle $\theta_{0}$ is shown in Fig. 7. We observe very good
agreement between the perturbation theory prediction and the numerical
simulations result. In particular, the numerical simulations confirm the
expectation that for the chosen parameter values, the value of $|\Delta
A_{1}^{(c)}|$ would be larger for smaller orientation angles, since in this
case, beam 2 traverses through the wider part of beam 1. Furthermore, the
relative error (in percentage) $|\Delta A_{1}^{(c)(num)}-\Delta
A_{1}^{(c)(th)}|\times 100/|\Delta A_{1}^{(c)(th)}|$ is smaller than $6.3\%$
in the entire interval $0\leq\theta_{0}\leq\pi/2$. Therefore, the numerical
simulations with Eq. (1) confirm the perturbation theory prediction for a
nontrivial dependence of $\Delta A_{1}^{(c)}$ on $\theta_{0}$ due to the
anisotropic (and nonseparable) nature of the initial condition.

Figure 6: (Color online) Contour plots of the beam shapes $|\psi_{j}(x,y,z)|$
at $z=0$ (a), $z=z_{i}=1.2$ (b), and $z=z_{f}=2$ (c) in a fast collision
between two Gaussian beams with the anisotropic initial condition of Eqs. (59)
and (61). The plots represent the beam shapes obtained by numerical solution
of Eq. (1) with parameter values $\epsilon_{3}=0.01$ and $d_{11}=20$. The
orientation angle is $\theta_{0}=\pi/4$. Figure 7: (Color online) The
collision-induced amplitude shift of beam 1 $\Delta A_{1}^{(c)}$ vs the
orientation angle $\theta_{0}$ in a fast collision between two Gaussian beams
with the anisotropic initial condition of Eqs. (59) and (61). The red circles
represent the result obtained by numerical simulations with Eq. (1). The solid
blue curve represents the perturbation theory prediction of Eq. (69).

###  III.5 Collision-induced change in the beam shape

It was shown in subsection II.2 that the two-beam collision in the presence of
weak cubic loss leads to a change of the beam shape in the transverse
direction. In contrast, it was shown in Refs. PNH2017B ; NHP2020  that within
the leading order of the perturbative calculation for the one-dimensional
case, the beam shape is preserved during a fast two-beam collision in the
presence of cubic loss. Thus, the collision-induced change in the beam shape
in the transverse direction is clearly a collisional effect that exists only
in spatial dimension higher than 1. In the current subsection, we investigate
this effect in detail for a concrete two-beam setup by both analytic
calculations and numerical simulations.

To enable a more accurate comparison between the perturbation theory
predictions and the numerical simulations, we assume that the effects of the
optical medium’s cubic loss on single-beam propagation are negligible compared
with cubic loss effects on inter-beam interaction. This situation can be
realized, for example, in certain semiconductors, in which two-photon
absorption (2PA) associated with the simultaneous absorption of two photons
with the same wavelength (degenerate 2PA) is much weaker than 2PA associated
with the simultaneous absorption of two photons with different wavelengths
(nondegenerate 2PA) Hagan2002 ; Hagan2011 ; Rauscher97 . Under this
assumption, the dynamics of the two-beam collision is described by the
following perturbed linear propagation model, in which the perturbation terms
are only due to the effects of weak cubic loss on two-beam interaction:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!i\partial_{z}\psi_{1}+\partial_{x}^{2}\psi_{1}+\partial_{y}^{2}\psi_{1}=-2i\epsilon_{3}|\psi_{2}|^{2}\psi_{1},$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!i\partial_{z}\psi_{2}+id_{11}\partial_{x}\psi_{2}+\partial_{x}^{2}\psi_{2}+\partial_{y}^{2}\psi_{2}=-2i\epsilon_{3}|\psi_{1}|^{2}\psi_{2}.\\!\\!\\!\\!\\!\\!\\!\\!$
(71)

Similar to Eq. (1), we assume that the cubic loss coefficient satisfies
$0<\epsilon_{3}\ll 1$. We consider the change in the beam shape in a collision
between two Gaussian beams as a concrete example. Therefore, the initial
condition for the collision problem is given by Eq. (49). This choice of the
initial condition allows us to obtain an explicit analytic expression for the
collision-induced change in the shape of beam 1 $\phi_{1}$ in the post-
collision interval.

Since the initial condition (49) is separable for both beams, we can calculate
$\phi_{1}$ in the post-collision interval by employing Eq. (48). In addition,
since the effects of cubic loss on single-beam propagation are negligible, we
can replace $A_{j}(z_{c}^{-})$ by $A_{j}(0)$ everywhere in the calculation.
Therefore, the coefficient $\tilde{a}_{1}$, which is defined in Eq. (43),
takes the form a1zc :

$\displaystyle\tilde{a}_{1}=2\pi^{1/2}\epsilon_{3}A_{1}(0)A_{2}^{2}(0)W_{20}^{(x)}/|d_{11}|,$
(72)

where $c_{p2}^{(x)}=\pi^{1/2}$ is used. The function $g_{1}^{(x)}(x,z)$ in Eq.
(48) is obtained by the solution of the unperturbed linear propagation
equation with the initial condition (49). This function is given by Eqs. (127)
and (129) in Appendix D. Additionally, in Appendix C we show that the inverse
Fourier transform of
$\hat{g}_{12}^{(y)}(k_{2},z_{c})\exp[-ik_{2}^{2}(z-z_{c})]$ is given by:

$\displaystyle{\cal
F}^{-1}\left(\hat{g}_{12}^{(y)}(k_{2},z_{c})\exp[-ik_{2}^{2}(z-z_{c})]\right)=$
$\displaystyle\frac{W_{10}^{(y)}W_{20}^{(y)2}\exp\left[-q_{1}(z_{c})y^{2}/R_{1}^{4}(z,z_{c})+i\chi_{1}^{(y)}(y,z)\right]}{(W_{10}^{(y)4}+4z_{c}^{2})^{1/4}(W_{20}^{(y)4}+4z_{c}^{2})^{1/2}R_{1}(z,z_{c})},$
(73)

where $q_{1}(z_{c})$, $R_{1}(z,z_{c})$, and $\chi_{1}^{(y)}(y,z)$ are given by
Eqs. (117), (122), and (123), respectively. Substituting Eqs. (72), (73),
(127), and (129) into Eq. (48), we obtain the following expression for
$\phi_{1}(x,y,z)$ in the post-collision interval:

$\displaystyle\phi_{1}(x,y,z)=\frac{\tilde{a}_{1}W_{10}^{(x)}W_{10}^{(y)}W_{20}^{(y)2}}{(W_{10}^{(x)4}+4z^{2})^{1/4}(W_{10}^{(y)4}+4z_{c}^{2})^{1/4}(W_{20}^{(y)4}+4z_{c}^{2})^{1/2}R_{1}(z,z_{c})}$
$\displaystyle\times\exp\left[-\frac{W_{10}^{(x)2}x^{2}}{2(W_{10}^{(x)4}+4z^{2})}-\frac{q_{1}(z_{c})y^{2}}{R_{1}^{4}(z,z_{c})}+i\chi_{1}^{(tot)}(x,y,z)\right],$
(74)

where the total phase factor $\chi_{1}^{(tot)}$ is given by:

$\chi_{1}^{(tot)}(x,y,z)=\chi_{10}^{(x)}(x,z)+\chi_{1}^{(y)}(y,z)+\alpha_{10}+\pi.$
(75)

Comparing Eqs. (75) and (135), we observe that the beam’s phase factor in the
post-collision interval is different from the phase factor of the unperturbed
beam. We therefore define the collision-induced change in the beam’s phase
factor by:

$\Delta\chi_{1}^{(tot)}(x,y,z)=\chi_{10}(x,y,z)-\chi_{1}^{(tot)}(x,y,z).$ (76)

Using Eqs. (75) and (135), we obtain

$\Delta\chi_{1}^{(tot)}(y,z)=\chi_{10}^{(y)}(y,z)-\chi_{1}^{(y)}(y,z)-\pi,$
(77)

where $\chi_{10}^{(y)}$ and $\chi_{1}^{(y)}$ are given by Eqs. (130) and
(123). Thus, only the dependence of the phase factor on the transverse
coordinate is changed by the collision, while the dependence on the
longitudinal coordinate remains unchanged. This change in the $y$ dependence
of the phase factor is a result of the change in the beam’s shape inside the
collision interval, which leads to a change in the y-dependence of
$\phi_{1}(x,y,z)$ in the post-collision interval [see detailed discussions
following Eqs. (36) and (48)].

The collision-induced change in the beam shape can be characterized by the
fractional intensity reduction factor $\Delta I_{1}^{(r)}$, which is defined
by:

$\displaystyle\Delta
I_{1}^{(r)}(x,y,z)=\frac{\tilde{I}_{1}(x,y,z)-I_{1}(x,y,z)}{\tilde{I}_{1}(x,y,z)}=1-\frac{I_{1}(x,y,z)}{A_{1}^{2}(0)\tilde{\Psi}_{10}^{2}(x,y,z)},$
(78)

where $I_{1}(x,y,z)=|\psi_{1}(x,y,z)|^{2}$ is the intensity of beam 1 at
$(x,y,z)$ in the presence of cubic loss, and
$\tilde{I}_{1}(x,y,z)=A_{1}^{2}(0)\tilde{\Psi}_{10}^{2}(x,y,z)$ is the
intensity of beam 1 in the absence of cubic loss. Thus, $\Delta I_{1}^{(r)}$
measures the ratio between the intensity decrease of beam 1, which is induced
solely by the effects of cubic loss on the collision, and the intensity of
beam 1 in the unperturbed single-beam propagation problem. We obtain the
perturbation theory prediction for the fractional intensity reduction factor
by expressing $\Delta I_{1}^{(r)}$ in terms of $\phi_{1}$. Using Eq. (13) and
keeping terms up to order $\epsilon_{3}/d_{11}$, we arrive at

$\displaystyle I_{1}\simeq
A_{1}^{2}(0)|\tilde{\psi}_{10}|^{2}+A_{1}(0)\left(\tilde{\psi}_{10}\phi_{1}^{*}+\tilde{\psi}_{10}^{*}\phi_{1}\right).$
(79)

Substitution of Eq. (79) together with the relation
$|\tilde{\psi}_{10}|=\tilde{\Psi}_{10}$ into Eq. (78) yields:

$\displaystyle\Delta
I_{1}^{(r)}=-\frac{\tilde{\psi}_{10}\phi_{1}^{*}+\tilde{\psi}_{10}^{*}\phi_{1}}{A_{1}(0)\tilde{\Psi}_{10}^{2}}.$
(80)

Using Eq. (11) and the relation $\phi_{1}=|\phi_{1}|\exp[i\chi_{1}^{(tot)}]$,
we obtain:

$\tilde{\psi}_{10}\phi_{1}^{*}+\tilde{\psi}_{10}^{*}\phi_{1}=2\tilde{\Psi}_{10}|\phi_{1}|\cos\left[\Delta\chi_{1}^{(tot)}\right].$
(81)

Therefore, we can also express $\Delta I_{1}^{(r)}$ as:

$\displaystyle\Delta
I_{1}^{(r)}=-\frac{2|\phi_{1}|\cos\left[\Delta\chi_{1}^{(tot)}\right]}{A_{1}(0)\tilde{\Psi}_{10}}.$
(82)

For a separable initial condition, $\Delta\chi_{1}^{(tot)}$ is given by Eq.
(77). In addition, the $x$ dependences of $\tilde{\Psi}_{10}$ and of the
leading order of $|\phi_{1}|$ are identical. As a result, in this case, the
dependence on $x$ cancels out on the right hand side of Eq. (82), and $\Delta
I_{1}^{(r)}$ becomes independent of $x$. Therefore, for a separable initial
condition, the expression for $\Delta I_{1}^{(r)}$ takes the form:

$\displaystyle\Delta
I_{1}^{(r)}(y,z)=-\frac{2|\phi_{1}(x,y,z)|\cos\left[\chi_{10}^{(y)}(y,z)-\chi_{1}^{(y)}(y,z)-\pi\right]}{A_{1}(0)\tilde{\Psi}_{10}(x,y,z)}.$
(83)

The fractional intensity reduction factor for the collision setup considered
in the current subsection is given by Eq. (83), where $\phi_{1}$,
$\tilde{\Psi}_{10}$, $\chi_{10}^{(y)}$, and $\chi_{1}^{(y)}$ are given by Eqs.
(74), (132), (130), and (123), respectively. We note that the effects of the
collision-induced change in the beam’s phase factor are included in Eq. (83)
via the dependence of the expression on the right hand side of this equation
on $\cos\left[\chi_{10}^{(y)}(y,z)-\chi_{1}^{(y)}(y,z)-\pi\right]$. We will
see in the next paragraphs that these effects can lead to negative values of
$\Delta I_{1}^{(r)}$ in certain intervals of $y$, that is, to a localized
increase in the intensity of beam 1 relative to the intensity in the
unperturbed single-beam propagation problem.

Figure 8: (Color online) Contour plots of the beam shapes $|\psi_{j}(x,y,z)|$
at $z=0$ (a), $z=z_{i}=0.96$ (b), and $z=z_{f}=1.6$ (c) in a fast collision
between two Gaussian beams with parameter values $\epsilon_{3}=0.01$ and
$d_{11}=25$. The plots represent the beam shapes obtained by numerical
solution of Eq. (71) with the initial condition (49).

We check the perturbation theory predictions for the collision-induced change
in the beam shape by extensive numerical simulations with Eq. (71). The
simulations are carried out with $\epsilon_{3}=0.01$ and with $d_{11}$ values
that are varied in the intervals $4\leq|d_{11}|\leq 60$. The parameter values
of the initial condition (49) are $A_{j}(0)=1$, $\alpha_{j0}=0$, $x_{20}=\pm
20$, $W_{10}^{(x)}=3$, $W_{10}^{(y)}=2$, $W_{20}^{(x)}=2$, and
$W_{20}^{(y)}=1$. The final propagation distance is
$z_{f}=2z_{c}=-2x_{20}/d_{11}$, and therefore, the beam centers are well
separated at $z_{f}$. Figure 8 shows the initial beam shapes
$|\psi_{j}(x,y,0)|$, and the beam shapes $|\psi_{j}(x,y,z)|$ obtained in the
simulation with $d_{11}=25$ at $z_{i}=0.96>z_{c}$, and at $z_{f}=1.6$. We
observe that both beams experience significant broadening due to diffraction.
Figure 9 shows the collision-induced change in the shape of beam 1 obtained in
the simulation with $d_{11}=25$ at $z=z_{f}$ $|\phi_{1}^{(num)}(x,y,z_{f})|$.
The perturbation theory prediction $|\phi_{1}^{(th)}(x,y,z_{f})|$, which is
obtained by Eq. (74) is also shown. The agreement between the simulation
result and the perturbation theory prediction is very good both near the
beam’s maximum and at the tails. We quantify the difference between
$|\phi_{1}^{(num)}(x,y,z)|$ and $|\phi_{1}^{(th)}(x,y,z)|$ by defining the
relative error (in percentage) $E_{r}^{(|\phi_{1}|)}(z)$ in the following
manner:

$\displaystyle E_{r}^{(|\phi_{1}|)}(z)=100\times\left[\int dx\int
dy\,|\phi_{1}^{(th)}(x,y,z)|^{2}\right]^{-1/2}$
$\displaystyle\times\left\\{\int dx\int
dy\left[\;\left|\phi_{1}^{(th)}(x,y,z)\right|-\left|\phi_{1}^{(num)}(x,y,z)\right|\;\right]^{2}\right\\}^{1/2},$
(84)

where the integration is carried out over the entire domain in the $xy$ plane,
which is used in the numerical simulation. We find that the value of
$E_{r}^{(|\phi_{1}|)}(z_{f})$ for $d_{11}=25$ is $7.0\%$, in accordance with
the good agreement between simulation and theory that is observed in Fig. 9.

Figure 9: (Color online) The collision-induced change in the shape of beam 1
$|\phi_{1}(x,y,z_{f})|$ at $z_{f}=1.6$ in a fast two-beam collision with
parameter values $\epsilon_{3}=0.01$ and $d_{11}=25$. Top: the perturbation
theory prediction of Eq. (74). Bottom: the result obtained by numerical
solution of Eq. (71).

We now turn to analyze the behavior of the fractional intensity reduction
factor as a function of $y$. Figure 10 shows the $y$ dependence of $\Delta
I_{1}^{(r)}$ at $z=z_{f}$ obtained in the numerical simulation with
$d_{11}=25$ Delta_I_1 . The perturbation theory prediction of Eqs. (83) and
(74) is also shown. The agreement between the simulations result and the
theoretical prediction is very good. Based on this result and on similar
results that are obtained with other choices of the physical parameter values
we conclude that the perturbation theory correctly captures the spatial
distribution of the intensity reduction in fast two-beam collisions in the
presence of weak cubic loss. We also point out that according to Eqs. (83) and
(74), $\Delta I_{1}^{(r)}(y,z_{f})$ attains negative values at intermediate
values of $y$, e.g., in the intervals $3.85\leq|y|\leq 6.25$ (for
$d_{11}=25$). This prediction is confirmed by the numerical simulation. In
particular, the numerically obtained $\Delta I_{1}^{(r)}(y,z_{f})$ attains
negative values in the intervals $3.8\leq|y|\leq 6.1$, in very good agreement
with the result of Eqs. (83) and (74). The negative values of $\Delta
I_{1}^{(r)}(y,z_{f})$ correspond to a localized increase in the intensity of
beam 1 relative to the intensity in the unperturbed single-beam propagation
problem. According to the perturbation theory, these values are a result of
the collision-induced change in the phase factor of beam 1, which affects the
value of $\Delta I_{1}^{(r)}(y,z_{f})$ via its dependence on
$\cos\,[\,\Delta\chi_{1}^{(tot)}(y,z_{f})\,]$ [see Eqs. (83) and (77)]. To
check if this is indeed the case, we compare the numerical simulation result
for the $y$ dependence of $\cos\,[\,\Delta\chi_{1}^{(tot)}(y,z_{f})\,]$ with
the perturbation theory prediction, which is obtained by using Eqs. (77),
(130), and (123). This comparison is shown in Fig. 11. We observe good
agreement between the results of the perturbative calculation and the
numerical simulation for this quantity. In particular, the perturbation theory
result for $\cos\,[\,\Delta\chi_{1}^{(tot)}(y,z_{f})\,]$ attains positive
values in the intervals $3.85\leq|y|\leq 6.25$, while the simulation result
attains positive values in the intervals $3.8\leq|y|\leq 6.1$. These intervals
are the same as the ones, in which the values of the theoretically and
numerically obtained $\Delta I_{1}^{(r)}(y,z_{f})$ are negative. We therefore
conclude that the relative localized intensity increase for beam 1 at
intermediate $y$ values is indeed a result of the collision-induced change in
the phase factor of this beam.

Figure 10: (Color online) The $y$ dependence of the fractional intensity
reduction factor for beam 1 at $z=z_{f}$ $\Delta I_{1}^{(r)}(y,z_{f})$ in a
two-beam collision with parameter values $\epsilon_{3}=0.01$ and $d_{11}=25$.
The solid blue curve corresponds to the perturbation theory prediction of Eqs.
(83) and (74). The other two curves represent the results obtained by
numerical solution of Eq. (71). The dashed red curve is obtained by averaging
$\Delta I_{1}^{(r)}(x,y,z_{f})$ over the $x$-interval $[-2,2]$. The dashed-
dotted green curve is obtained by using the numerically computed value of
$\Delta I_{1}^{(r)}(0,y,z_{f})$. Figure 11: (Color online) The $y$ dependence
of the cosine of the collision-induced change in the phase factor of beam 1 at
$z=z_{f}$ $\cos\,[\,\Delta\chi_{1}^{(tot)}(y,z_{f})\,]$ in a two-beam
collision with parameter values $\epsilon_{3}=0.01$ and $d_{11}=25$. The solid
blue curve corresponds to the perturbation theory prediction, which is
obtained by using Eqs. (77), (130), and (123). The other two curves represent
the results obtained by numerical solution of Eq. (71). The dashed red curve
is obtained by averaging $\cos\,[\,\Delta\chi_{1}^{(tot)}(x,y,z_{f})\,]$ over
the $x$-interval $[-2,2]$. The dashed-dotted green curve is obtained by using
the numerically computed value of
$\cos\,[\,\Delta\chi_{1}^{(tot)}(0,y,z_{f})\,]$.

Finally, we study the dependence of the fractional intensity reduction factor
on the value of the beam-steering coefficient by measuring $\Delta
I_{1}^{(r)}(0,z_{f})$ as a function of $d_{11}$. Figure 12 shows the
dependence of $\Delta I_{1}^{(r)}(0,z_{f})$ on $d_{11}$ obtained in the
numerical simulations together with the theoretical prediction of Eqs. (83)
and (74). The agreement between the simulations result and the perturbation
theory prediction is excellent over the entire interval of $d_{11}$ values.
More specifically, the relative error in the approximation of $\Delta
I_{1}^{(r)}(0,z_{f})$ (in percentage), which is defined by $|\Delta
I_{1}^{(r)(num)}(0,z_{f})-\Delta I_{1}^{(r)(th)}(0,z_{f})|\times 100/|\Delta
I_{1}^{(r)(th)}(0,z_{f})|$, is smaller than $0.6\%$ for $10\leq|d_{11}|\leq
60$, and smaller than $1.2\%$ for $4\leq|d_{11}|<10$. We also checked the
dependence of $\Delta A_{1}^{(c)}$ on $d_{11}$, and obtained very good
agreement between the perturbation theory prediction and the numerical
simulations result (similar to what is shown in Figs. 2 and 4). Based on these
results and on the results shown in Figs. 9 \- 11 we conclude that the
perturbation theory of subsections II.2 and II.3 enables accurate calculation
of both the change in the beam shapes and the dynamics of the beam amplitudes
in fast two-beam collisions in the presence of weak cubic loss.

Figure 12: (Color online) The $d_{11}$ dependence of the fractional intensity
reduction factor for beam 1 at $y=0$ and $z=z_{f}$, $\Delta
I_{1}^{(r)}(0,z_{f})$, in a fast two-beam collision with $\epsilon_{3}=0.01$.
The solid blue curve corresponds to the perturbation theory prediction of Eqs.
(83) and (74). The other two curves are obtained from the numerical solution
of Eq. (71). The red squares represent the result obtained by averaging
$\Delta I_{1}^{(r)}(x,0,z_{f})$ over the $x$-interval $[-2,2]$. The green
circles represent the result obtained by using the numerically computed value
of $\Delta I_{1}^{(r)}(0,0,z_{f})$.

##  IV Fast collisions between pulsed-beams in spatial dimension 3

###  IV.1 Introduction

We consider the dynamics of collisions between two pulsed-beams in a three-
dimensional linear optical medium with weak cubic loss. Similar to sections II
and III, we assume that the pulsed-beams propagate along the $z$ axis, and
that the propagation is accurately described by the paraxial approximation.
For each value of $z$, the distribution of the optical field depends on the
spatial coordinates $x$ and $y$, and on time $t$. Therefore, using the
terminology that was introduced in section II.1, the spatial dimension is 3,
and the propagation is described by a $(3+1)$-dimensional propagation model.

We take into account the effects of first and second-order dispersion,
isotropic diffraction, and weak cubic loss, and neglect beam-steering. Since
beam-steering is neglected, the relative velocity vector between the beam
centers lies along the $t$ axis in the $txy$ space. Therefore, the dynamics of
the collision is described by the following weakly perturbed linear
propagation model:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!i\partial_{z}\psi_{1}+\partial^{2}_{t}\psi_{1}+d_{2}\partial^{2}_{x}\psi_{1}+d_{2}\partial^{2}_{y}\psi_{1}=-i\epsilon_{3}|\psi_{1}|^{2}\psi_{1}-2i\epsilon_{3}|\psi_{2}|^{2}\psi_{1},$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!i\partial_{z}\psi_{2}+id_{11}\partial_{t}\psi_{2}+\partial^{2}_{t}\psi_{2}+d_{2}\partial^{2}_{x}\psi_{2}+d_{2}\partial^{2}_{y}\psi_{2}=-i\epsilon_{3}|\psi_{2}|^{2}\psi_{2}-2i\epsilon_{3}|\psi_{1}|^{2}\psi_{2},$
(85)

where $\psi_{j}$ are proportional to the electric fields of the beams, $x$,
$y$, and $z$ are the spatial coordinates, and $t$ is time Dimensions2 ;
sign_beta_2 . In Eq. (1), $d_{11}$ is the first-order dispersion coefficient,
$d_{2}$ is the diffraction coefficient, and $\epsilon_{3}$ is the cubic loss
coefficient. The term $id_{11}\partial_{t}\psi_{2}$ describes the effects of
first-order dispersion, the terms $\partial^{2}_{t}\psi_{j}$ describe the
effects of second-order dispersion, and $d_{2}\partial_{x}^{2}\psi_{j}$ and
$d_{2}\partial_{y}^{2}\psi_{j}$ describe the effects of isotropic diffraction.
Additionally, the terms $-i\epsilon_{3}|\psi_{j}|^{2}\psi_{j}$ and
$-2i\epsilon_{3}|\psi_{k}|^{2}\psi_{j}$ describe intra-beam and inter-beam
effects due to weak cubic loss.

We consider collisions between pulsed-beams with general initial shapes and
with tails that decrease sufficiently fast, such that the values of the
integrals
$\int_{-\infty}^{\infty}dt\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dy\,|\psi_{j}(t,x,y,0)|^{2}$
are finite. We assume that the initial pulsed-beams can be characterized by
the following parameters. (1) The initial amplitudes $A_{j}(0)$. (2) The
initial widths of the pulsed-beams along the $t$, $x$, and $y$ axes,
$W_{j0}^{(t)}$, $W_{j0}^{(x)}$, and $W_{j0}^{(y)}$. (3) The initial positions
of the beam centers $(t_{j0},x_{j0},y_{j0})$. (4) The initial phases
$\alpha_{j0}$. Thus, the initial electric fields can be written as:

$\displaystyle\psi_{j}(t,x,y,0)=A_{j}(0)h_{j}(t,x,y)\exp(i\alpha_{j0}),$ (86)

where $h_{j}(t,x,y)$ is real-valued. We are equally interested in the
important case, where the initial electric fields of the two pulsed-beams are
separable. i.e., where each of the functions $\psi_{j}(t,x,y,0)$ is a product
of three functions of $t$, $x$, and $y$ laser_modes . In this case, the
initial electric fields can be expressed as:

$\displaystyle\psi_{j}(t,x,y,0)=A_{j}(0)h_{j}^{(t)}[(t-t_{j0})/W_{j0}^{(t)}]h_{j}^{(x)}[(x-x_{j0})/W_{j0}^{(x)}]$
$\displaystyle\times h_{j}^{(y)}[(y-y_{j0})/W_{j0}^{(y)}]\exp(i\alpha_{j0}).$
(87)

We are interested in the collision-induced dynamics of complete fast
collisions. We therefore obtain conditions on the physical parameter values
for these collisions. As stated earlier, the relative velocity vector of the
pulsed-beams centers lies along the $t$ axis in the $txy$ space. As a result,
the conditions for a complete fast collision involve the initial widths of the
pulsed-beams along the $t$ axis $W_{j0}^{(t)}$, as well as the initial and
final values of the $t$ coordinate of the pulsed-beam centers $t_{j0}$ and
$t_{j}(z_{f})$, respectively. The conditions for a complete collision are
obtained by requiring that the pulsed-beams are well-separated in time before
and after the collision. This requirement yields the following inequalities:
$|t_{20}-t_{10}|\gg W_{10}^{(t)}+W_{20}^{(t)}$ and
$|t_{2}(z_{f})-t_{1}(z_{f})|\gg W_{1}^{(t)}(z_{f})+W_{2}^{(t)}(z_{f})$, where
$W_{j}^{(t)}(z_{f})$ are the pulsed-beam widths along the $t$ axis at
$z=z_{f}$. The collision length $\Delta z_{c}$ is given by $\Delta
z_{c}=2(W_{10}^{(t)}+W_{20}^{(t)})/|d_{11}|$. The assumption of a fast
collision means that $\Delta z_{c}$ is much smaller than the length scale
$z_{D}^{(min)}$, which is the smallest dispersion length or diffraction length
in the problem. By definition,
$z_{D}^{(min)}=\min\left\\{z_{d1}^{(t)},z_{d2}^{(t)},z_{D1}^{(x)},z_{D2}^{(x)},z_{D1}^{(y)},z_{D2}^{(y)}\right\\}$,
where $z_{dj}^{(t)}$ are the dispersion lengths of the pulsed-beams, and
$z_{Dj}^{(x)}$ and $z_{Dj}^{(y)}$ are the diffraction lengths along the $x$
and $y$ axes. Requiring that $\Delta z_{c}\ll z_{D}^{(min)}$, we obtain
$2(W_{10}^{(t)}+W_{20}^{(t)})\ll|d_{11}|z_{D}^{(min)}$, as the condition for a
fast collision.

###  IV.2 The perturbation theory predictions for the collision-induced
changes in the pulsed-beam shape and amplitude

The perturbative calculation of the collision-induced changes in the shapes
and amplitudes of the pulsed-beams in spatial dimension 3 is very similar to
the calculation that was presented in section II for the two-dimensional case.
For this reason and for brevity, we do not present the entire derivation of
the equations for the collision-induced dynamics. Instead, we present only the
two main results of the perturbative calculation, namely, the expressions for
the collision-induced changes in the shape and amplitude of pulsed-beam 1,
$\Delta\Phi_{1}(t,x,y,z_{c})$ and $\Delta A_{1}^{(c)}$. The notations for the
physical quantities are the same as the ones used in section II, apart from
the fact that $\psi_{j}$, $\psi_{j0}$, $\tilde{\psi}_{j0}$,
$\tilde{\Psi}_{j0}$, and $\chi_{j0}$ are now functions of $t$, $x$, $y$, and
$z$.

We start by describing the results of the perturbative calculation for the
general initial condition (86). In this case, the collision-induced change in
the shape of pulsed-beam 1 inside the collision interval is given by:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\Delta\Phi_{1}(t,x,y,z_{c})\\!=\\!-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}\tilde{\Psi}_{10}(t,x,y,z_{c})\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!d\tilde{t}\,\bar{\Psi}_{20}^{2}(\tilde{t},x,y,z_{c}),$
(88)

where $\tilde{t}=t-t_{20}-d_{11}z$. Additionally, the collision-induced change
in the amplitude of pulsed-beam 1 is

$\displaystyle\\!\\!\\!\\!\Delta
A_{1}^{(c)}=-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{C_{p1}|d_{11}|}$
$\displaystyle\times\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dt\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\;\tilde{\Psi}_{10}^{2}(t,x,y,z_{c})\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!d\tilde{t}\;\bar{\Psi}_{20}^{2}(\tilde{t},x,y,z_{c}),$
(89)

where $C_{p1}$ is given by:

$\displaystyle
C_{p1}=\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dt\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\;\tilde{\Psi}_{10}^{2}(t,x,y,0).$
(90)

Further insight into the collision dynamics is obtained when the initial
condition is of the form (87), which is completely separable for both beams.
This case is also of special interest for practical reasons laser_modes . In
this case, the collisional change in the shape of pulsed-beam 1 in the
collision interval is given by:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\Delta\Phi_{1}(t,x,y,z_{c})\\!=\\!-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}c_{p2}^{(t)}W_{20}^{(t)}$
$\displaystyle\times
G_{1}^{(t)}(t,z_{c})G_{1}^{(x)}(x,z_{c})G_{1}^{(y)}(y,z_{c})G_{2}^{(x)2}(x,z_{c})G_{2}^{(y)2}(y,z_{c}).$
(91)

The real-valued functions $G_{j}^{(t)}(t,z)$, $G_{j}^{(x)}(x,z)$, and
$G_{j}^{(y)}(y,z)$ in Eq. (91) are defined by:

$\displaystyle\tilde{\psi}_{j0}(t,x,y,z)=G_{j}^{(t)}(t,z)G_{j}^{(x)}(x,z)G_{j}^{(y)}(y,z)$
$\displaystyle\times\exp\left\\{i\left[\chi_{j0}^{(t)}(t,z)+\chi_{j0}^{(x)}(x,z)+\chi_{j0}^{(y)}(y,z)+\alpha_{0}\right]\right\\},$
(92)

where $\tilde{\psi}_{j0}$ are the solutions to the unperturbed linear
propagation equations with the separable initial condition (87), and
$\chi_{j0}^{(t)}(t,z)$, $\chi_{j0}^{(x)}(x,z)$, and $\chi_{j0}^{(y)}(y,z)$ are
the (real-valued) phase factors. In addition, the coefficient $c_{p2}^{(t)}$
is given by the equation

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!dt\,G_{j}^{(t)2}(t,z)=\int_{-\infty}^{\infty}\\!\\!\\!\\!dt\,G_{j}^{(t)2}(t,0)=W_{j0}^{(t)}\int_{-\infty}^{\infty}\\!\\!\\!\\!ds\,h_{j}^{(t)2}(s)=W_{j0}^{(t)}c_{pj}^{(t)}.$
(93)

Similar to the situation in the two-dimensional case, it follows from Eq. (91)
that the $t$ dependence of the pulsed-beam is not changed by the collision at
all (within the leading order of the perturbative calculation).

The collision-induced change in the amplitude of pulsed-beam 1 in the case of
a completely separable initial condition is given by:

$\displaystyle\\!\\!\\!\\!\Delta
A_{1}^{(c)}=-\frac{2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}\frac{c_{p2}^{(t)}W_{20}^{(t)}}{c_{p1}^{(x)}W_{10}^{(x)}c_{p1}^{(y)}W_{10}^{(y)}}$
$\displaystyle\times\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dx\,G_{1}^{(x)2}(x,z_{c})G_{2}^{(x)2}(x,z_{c})\\!\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!dy\,G_{1}^{(y)2}(y,z_{c})G_{2}^{(y)2}(y,z_{c}).$
(94)

We observe that Eq. (94) has the same form as Eq. (39), where the longitudinal
factor is now $c_{p2}^{(t)}W_{20}^{(t)}$ and the overall factor is
$2\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})/|d_{11}|$. This finding is
in accordance with the expectation that the form (39) is valid for a general
spatial dimension when the initial condition is separable in the longitudinal
direction for both beams. In addition, similar to the situation in spatial
dimensions 1 and 2, the longitudinal factor $c_{p2}^{(t)}W_{20}^{(t)}$ is
universal in the sense that it does not depend on the exact details of the
initial pulsed-beam shapes and on the collision distance $z_{c}$.

###  IV.3 Numerical simulations for pulsed-beam collisions

The predictions of the perturbative calculation in subsection IV.2 are based
on several simplifying assumptions. In particular, it is assumed that the
pulsed-beams are sharply peaked in the $txy$ space throughout the propagation,
and that as a result, the integration with respect to $z$ can be extended to
$\pm\infty$. In addition, the explicit conditions for the validity of the
approximations employed by the perturbative calculation are not known. For
these reasons, it is important to check the predictions of the perturbative
calculation by numerical simulations with the perturbed linear propagation
model (85).

As explained in the beginning of subsection IV.2, the perturbative calculation
of the collision-induced dynamics in spatial dimension 3 is very similar to
the calculation for spatial dimension 2. For this reason and for brevity, we
do not present the results of the numerical simulations for all the
collisional setups that were considered in section III in the two-dimensional
case. Instead, we present only the simulations results for the collisional
setup, which is used for checking the theoretical predictions for universality
of the longitudinal part in the expression for $\Delta A_{1}^{(c)}$ (the setup
considered in subsection III.2 for spatial dimension 2).

Similar to subsection III.2, we choose two initial conditions with widely
different pulsed-beam profiles in the longitudinal (temporal) direction. More
specifically, the $t$ dependence of the pulsed-beams in the first initial
condition is Gaussian, i.e., it has rapidly decaying tails. In contrast, the
$t$ dependence of the pulsed-beams in the second initial condition is given by
a generalized Cauchy-Lorentz distribution, i.e., it has slowly decaying tails,
whose decay is characterized by a power-law. The initial profiles of the
pulsed-beams in the transverse direction (that is, the initial dependence on
$x$ and $y$) is taken as Gaussian, since this choice enables the explicit
calculation of the integrals with respect to $x$ and $y$ on the right hand
side of Eq. (94). The numerical simulations with the two types of initial
pulsed-beams, which posses widely different temporal (longitudinal) profiles,
provide a careful test to the perturbation theory prediction for universal
behavior of the longitudinal part in the expression for the collision-induced
amplitude shift.

The initial condition for a collision between two Gaussian pulsed-beams is
given by:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\psi_{1}(t,x,y,0)=A_{1}(0)\exp\left[-\frac{t^{2}}{2W^{(t)2}_{10}}-\frac{x^{2}}{2W^{(x)2}_{10}}-\frac{y^{2}}{2W^{(y)2}_{10}}+i\alpha_{10}\right],$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\psi_{2}(t,x,y,0)=A_{2}(0)\exp\left[-\frac{(t-t_{20})^{2}}{2W^{(t)2}_{20}}-\frac{x^{2}}{2W^{(x)2}_{20}}-\frac{y^{2}}{2W^{(y)2}_{20}}+i\alpha_{20}\right].$
(95)

Additionally, the initial condition for a collision between two Cauchy-
Lorentz-Gaussian pulsed-beams is given by:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\psi_{1}(t,x,y,0)=A_{1}(0)\left[1+\frac{2t^{4}}{W^{(t)4}_{10}}\right]^{-1}\exp\left[-\frac{x^{2}}{2W^{(x)2}_{10}}-\frac{y^{2}}{2W^{(y)2}_{10}}+i\alpha_{10}\right],$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\psi_{2}(t,x,y,0)=A_{2}(0)\left[1+\frac{2(t-t_{20})^{4}}{W^{(t)4}_{20}}\right]^{-1}\exp\left[-\frac{x^{2}}{2W^{(x)2}_{20}}-\frac{y^{2}}{2W^{(y)2}_{20}}+i\alpha_{20}\right].$
(96)

The calculation of the collision-induced amplitude shift by using Eq. (94) is
similar to the one that was presented in subsection III.2 for the two-
dimensional problem. This calculation yields the following expression for
$\Delta A_{1}^{(c)}$:

$\displaystyle\Delta
A_{1}^{(c)}=-\frac{2b\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}\frac{W_{10}^{(x)}W_{10}^{(y)}W_{20}^{(t)}W_{20}^{(x)2}W_{20}^{(y)2}}{(W_{10}^{(x)2}+W_{20}^{(x)2})^{1/2}(W_{10}^{(y)2}+W_{20}^{(y)2})^{1/2}}$
$\displaystyle\times\frac{1}{(W_{10}^{(x)2}W_{20}^{(x)2}+4d_{2}^{2}z_{c}^{2})^{1/2}(W_{10}^{(y)2}W_{20}^{(y)2}+4d_{2}^{2}z_{c}^{2})^{1/2}},$
(97)

where $b=\pi^{1/2}$ for Gaussian pulsed-beams, and $b=3\pi/2^{11/4}$ for
Cauchy-Lorentz-Gaussian pulsed-beams. The longitudinal part of the expression
for $\Delta A_{1}^{(c)}$, $c_{p2}^{(t)}W_{20}^{(t)}=bW_{20}^{(t)}$, is
universal. On the other hand, the transverse part of the expression, which is
given by

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{transverse
factor}=\frac{W_{10}^{(x)}W_{10}^{(y)}W_{20}^{(x)2}W_{20}^{(y)2}}{(W_{10}^{(x)2}+W_{20}^{(x)2})^{1/2}(W_{10}^{(y)2}+W_{20}^{(y)2})^{1/2}}$
$\displaystyle\times\frac{1}{(W_{10}^{(x)2}W_{20}^{(x)2}+4d_{2}^{2}z_{c}^{2})^{1/2}(W_{10}^{(y)2}W_{20}^{(y)2}+4d_{2}^{2}z_{c}^{2})^{1/2}},$
(98)

is clearly not universal. We define the quantity $\Delta A_{1}^{(c)(s)}$,
which is used for measuring the deviation of the $d_{11}$ dependence of
$\Delta A_{1}^{(c)}$ from the $1/|d_{11}|$ scaling, by a simple generalization
of the definition in the two-dimensional case. For this purpose, we note that
the collision distance $z_{c}$ satisfies $z_{c}=(t_{10}-t_{20})/d_{11}$.
$\Delta A_{1}^{(c)(s)}$ is defined as the approximate expression for the
collision-induced amplitude shift that is obtained by neglecting the terms
$4d_{2}^{2}(t_{10}-t_{20})^{2}/d_{11}^{2}$ in the denominator of Eq. (97).
Therefore, $\Delta A_{1}^{(c)(s)}$ is given by:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\Delta
A_{1}^{(c)(s)}=-\frac{2b\epsilon_{3}A_{1}(z_{c}^{-})A_{2}^{2}(z_{c}^{-})}{|d_{11}|}\frac{W_{20}^{(t)}W_{20}^{(x)}W_{20}^{(y)}}{(W_{10}^{(x)2}+W_{20}^{(x)2})^{1/2}(W_{10}^{(y)2}+W_{20}^{(y)2})^{1/2}}.$
(99)

It is clear from the definition of $\Delta A_{1}^{(c)(s)}$ that the difference
$|\Delta A_{1}^{(c)}-\Delta A_{1}^{(c)(s)}|$ measures the deviation of the
$d_{11}$ dependence of $\Delta A_{1}^{(c)}$ from the $1/|d_{11}|$ scaling
observed in spatial dimension 1.

We check the perturbation theory predictions for $\Delta A_{1}^{(c)}$ and for
the universality of the longitudinal part in the expression for $\Delta
A_{1}^{(c)}$ by numerical simulations with Eq. (85) with the initial
conditions (95) and (96). These initial conditions possess very different
temporal (longitudinal) pulsed-beam profiles. In particular, the tails of the
Gaussian pulsed-beams decay rapidly with increasing values of $|t|$ and
$|t-t_{20}|$. For such pulsed-beams, the sharp-peak approximation that is used
in the derivation of Eqs. (89) and (94) is expected to hold. In contrast, the
tails of the Cauchy-Lorentz-Gaussian pulsed-beams are slowly decaying with
increasing values of $|t|$ and $|t-t_{20}|$, and the decay is characterized by
a power-law. It is unclear if the sharp-peak approximation is valid for
pulsed-beams with such slowly decaying tails. Therefore, the numerical
simulations with the initial conditions (95) and (96) provide a careful check
for the validity of the perturbation theory approximations for pulsed-beams
with widely different temporal (longitudinal) distributions. In this way, the
simulations help to establish the universality of the longitudinal part in the
expression for $\Delta A_{1}^{(c)}$.

Equation (85) is numerically solved by the split-step method with periodic
boundary conditions Agrawal2001 ; Yang2010 . Since we are interested in fast
collisions, we perform the simulations for $d_{11}$ values in the intervals
$4\leq|d_{11}|\leq 60$. The other physical parameters values in Eq. (85) are
chosen as $\epsilon_{3}=0.01$ and $d_{2}=1.5$, as an example. The parameter
values of the initial pulsed-beams are $A_{j}(0)=1$, $\alpha_{j0}=0$,
$t_{20}=\pm 15$, $W_{10}^{(t)}=2$, $W_{10}^{(x)}=3$, $W_{10}^{(y)}=4$,
$W_{20}^{(t)}=3$, $W_{20}^{(x)}=4$, and $W_{20}^{(y)}=5$. The final
propagation distance is $z_{f}=2z_{c}=-2t_{20}/d_{11}$. We emphasize that
results similar to the ones presented in the current subsection are obtained
in simulations with other values of the physical parameters. For each initial
condition, we compare the numerically obtained dependence of $\Delta
A_{1}^{(c)}$ on $d_{11}$ with the theoretical predictions of Eqs. (97) and
(99). We also describe the behavior of the relative errors in the
approximation of $\Delta A_{1}^{(c)}$ (in percentage), $E_{r}^{(1)}$ and
$E_{r}^{(2)}$, which were defined in subsection III.2.

We first discuss the simulations results for Gaussian pulsed-beams. Figure 13
shows the values of $|\psi_{j}(t,x,y,z)|$ obtained in the simulation with
$d_{11}=20$ at three specific planes (cross-sections) at the distances $z=0$,
$z=z_{i}=0.9$, and $z=z_{f}=1.5$ zi_values . At each distance, we choose the
three cross-sections (planes) such that the main bodies of the pulsed-beams
are shown clearly cross_sections . In particular, one cross-section, which is
denoted by $CS^{(0)}_{2}$, is located at the plane $x=0$. Additionally, the
other two cross-sections are located at the planes $t=t_{j}(z)$, where
$j=1,2$, $t_{j}(z)$ is the $t$ coordinate of the $j$th pulsed-beam’s center,
and $z$ can take the values $0$, $z_{i}$, or $z_{f}$. The latter cross-
sections are denoted by $CS^{(t_{j}(z))}_{1}$, where $j=1,2$. It is seen that
the pulsed-beams experience broadening due to the effects of second-order
dispersion and diffraction. Additionally, the maximum values of
$|\psi_{j}(t,x,y,z)|$ decrease with increasing distance, mainly due to the
broadening. The dependence of $\Delta A_{1}^{(c)}$ on $d_{11}$ obtained in the
simulations is shown in Fig. 14 together with the two theoretical predictions
of Eqs. (97) and (99). The agreement between the simulations result and the
prediction of Eq. (97) is very good. More specifically, the relative error
$E_{r}^{(1)}$ is smaller than $1.9\%$ for $10\leq|d_{11}|\leq 60$ and smaller
than $4.5\%$ for $4\leq|d_{11}|<10$. In addition, we observe good agreement
between the simulations result and the more crude approximation $\Delta
A_{1}^{(c)(s)}$ for large $|d_{11}|$ values, but there is a noticeable
difference between the results for intermediate $|d_{11}|$ values. In
particular, the relative error $E_{r}^{(2)}$ is smaller than $10.4\%$ for
$10\leq|d_{11}|\leq 60$ and smaller than $39.3\%$ for $4\leq|d_{11}|<10$.
Thus, as expected from Eqs. (97) and (99), the deviation of the $d_{11}$
dependence of $\Delta A_{1}^{(c)}$ from the $1/|d_{11}|$ scaling increases
with decreasing value of $|d_{11}|$.

Figure 13: (Color online) Contour plots of the pulsed-beam shapes
$|\psi_{j}(t,x,y,z)|$ on three planes in the $txy$ space at $z=0$ (a),
$z=z_{i}=0.9$ (b), and $z=z_{f}=1.5$ (c) in a fast collision between two
Gaussian pulsed-beams. The physical parameter values are $\epsilon_{3}=0.01$,
$d_{2}=1.5$, and $d_{11}=20$. The plots represent the beam shapes obtained by
numerical solution of Eq. (85) with the initial condition (95). The plane
$x=0$ is denoted by $CS^{(0)}_{2}$, and the planes $t=t_{j}(z)$ with $j=1,2$
are denoted by $CS^{(t_{j}(z))}_{1}$. Figure 14: (Color online) The collision-
induced amplitude shift of pulsed-beam 1 $\Delta A_{1}^{(c)}$ vs the first-
order dispersion coefficient $d_{11}$ in a fast collision between two Gaussian
pulsed-beams with parameter values $\epsilon_{3}=0.01$ and $d_{2}=1.5$. The
red circles represent the result obtained by numerical simulations with Eq.
(85) with the initial condition (95). The solid blue and dashed green curves
represent the theoretical predictions of Eqs. (97) and (99), respectively.

We now describe the results of the simulations for collisions between Cauchy-
Lorentz-Gaussian pulsed-beams, which serve as an example for pulsed-beams with
tails that exhibit slow (power-law) temporal decay. In this case, it is not
clear if the sharp-peak approximation used in the perturbative calculation is
valid. For brevity, we discuss the simulations results without showing the
corresponding figures. The numerical simulation with $d_{11}=20$ shows that
the pulsed-beams experience considerable broadening and develop extended tails
due to second-order dispersion and diffraction. As a result, the maximum
values of $|\psi_{j}(t,x,y,z)|$ decrease with increasing $z$. Furthermore,
despite the broadening of the pulsed-beams, the agreement between the
numerical simulations result for $\Delta A_{1}^{(c)}$ and the theoretical
prediction of Eq. (97) is very good. In particular, the relative error
$E_{r}^{(1)}$ is less than $1.3\%$ for $10\leq|d_{11}|\leq 60$ and less than
$10.2\%$ for $4\leq|d_{11}|<10$. The values of $E_{r}^{(1)}$ are comparable to
the values obtained for collisions between Gaussian pulsed-beams. Based on
these findings and on similar results obtained with other parameter values, we
conclude that the longitudinal part of the expression for $\Delta A_{1}^{(c)}$
is indeed universal, since it is not sensitive to the details of the initial
pulsed-beams shapes. We also note that the relative error $E_{r}^{(2)}$ in the
approximation of the amplitude shift by $\Delta A_{1}^{(c)(s)}$ is smaller
than $9.8\%$ for $10\leq|d_{11}|\leq 60$ and smaller than $42.9\%$ for
$4\leq|d_{11}|<10$. These values are significantly larger than the
corresponding values of $E_{r}^{(1)}$. Thus, the agreement between the
numerical result and the perturbation theory prediction for $\Delta
A_{1}^{(c)}$ is significantly better than the agreement between the numerical
result and the cruder approximation $\Delta A_{1}^{(c)(s)}$. Additionally, as
in the case of Gaussian pulsed-beams, the deviation of the $d_{11}$ dependence
of $\Delta A_{1}^{(c)}$ from the $1/|d_{11}|$ scaling increases with
decreasing value of $|d_{11}|$.

##  V Conclusions

We studied the dynamics of fast collisions between two optical beams in linear
optical media with weak cubic loss in spatial dimension higher than 1. For
this purpose, we introduced a perturbation method, which generalizes the
perturbation method developed in Refs. PNH2017B ; NHP2020  for the one-
dimensional case in three major ways. First, it extends the perturbative
calculation from spatial dimension 1 to spatial dimension 2, and enables the
extension of the calculation to a general spatial dimension in a
straightforward manner. Second, it provides a perturbative calculation of the
collision-induced dynamics of the beam shape both in the collision interval
and outside of the collision interval. In contrast, the perturbative
calculation in Refs. PNH2017B ; NHP2020  was limited to the collision interval
only. Third, it enables the discovery and analysis of several collision-
induced effects, which exist only in spatial dimension higher than 1.

We used the generalized two-dimensional version of the perturbation method to
obtain formulas for the collision-induced changes in the beam shapes and
amplitudes. This was done both for a general initial condition and for the
important case of an initial condition, which is separable for both beams. We
found that for a general initial condition, the collision leads to a change in
the beam shape in the direction transverse to the relative velocity vector
between the beam centers. Additionally, we found that for a separable initial
condition, the beam shape in the longitudinal direction is not changed by the
collision within the leading order of the perturbative calculation. Moreover,
we showed that for a separable initial condition, the longitudinal part in the
expression for the collision-induced amplitude shift is universal, while the
transverse part is proportional to the integral of the product of the beam
intensities with respect to the transverse coordinate. We also demonstrated
that the same behavior of the longitudinal and transverse parts of the
expression for the amplitude shift exists in collisions between pulsed optical
beams in spatial dimension 3.

We checked these predictions of the generalized perturbation method along with
other predictions concerning the effects on the collision of partial beam
overlap and anisotropy in the initial condition by extensive numerical
simulations with the perturbed linear propagation model in spatial dimensions
2 and 3. The simulations in spatial dimension 2 were carried out for four
different two-beam collision setups. These setups demonstrate four major
collisional effects and properties that either exist only in spatial dimension
higher than 1, or are qualitatively different from their one-dimensional
counterparts. (1) The universality of the longitudinal part in the expression
for the collision-induced amplitude shift. (2) The effect of partial beam
overlap. (3) The effect of anisotropy in the initial condition. (4) The
collision-induced change in the beam shape in the transverse direction. The
prediction for universal behavior of the longitudinal part in the expression
for the amplitude shift was also checked in spatial dimension 3 by numerical
simulations of collisions between pulsed optical beams.

In all the simulation setups we obtained very good agreement between the
perturbation theory and the numerical simulations. In particular, in setup
(1), the simulations showed that the longitudinal part in the expression for
the collision-induced amplitude shift is universal in the sense that it is not
sensitive to the details of the initial beam shapes, and that this is true in
both spatial dimensions 2 and 3. Additionally, the simulations in setup (2)
verified the validity of the transverse part in the expression for the
amplitude shift and demonstrated that the generalized perturbation method can
be employed for fast collisions between partially overlapping beams. The
simulations in the anisotropic setup [setup (3)] verified the complex
dependence of the expression for the amplitude shift on the orientation angle
$\theta_{0}$, which was predicted by the perturbation theory. We attributed
this complex dependence to the nonseparable nature of the initial condition in
the anisotropic setup. Moreover, the numerical simulations in setup (4)
confirmed the perturbation theory predictions for the collision-induced change
in the beam shape in the transverse direction. Based on the results of the
latter simulations we concluded that the generalized perturbation method that
we developed in the current paper enables accurate calculation of both the
change in the beam shape and the dynamics of the beam amplitude in fast two-
beam collisions in the presence of weak cubic loss.

In summary, our study extended the results of the previous works in Refs.
PNH2017B ; QMN2018 ; NHP2020  in two major ways. First, it generalized the
perturbation method of Refs. PNH2017B ; QMN2018 ; NHP2020  from spatial
dimension 1 to a general spatial dimension. Second, it demonstrated a variety
of collision-induced physical effects, which exist only for spatial dimension
higher than 1. We point out that in another study, we developed a similar
perturbation method for analyzing fast two-pulse collisions in systems
described by linear diffusion-advection equations with weak quadratic loss in
spatial dimension higher than 1 PHN2020 . Using the latter perturbation method
and numerical simulations, we found that the collision-induced changes in
pulse shapes and amplitudes in the perturbed linear diffusion-advection
systems exhibit similar behavior to the one reported in the current paper.
Thus, the perturbation methods developed in the current paper and in Ref.
PHN2020  are very valuable tools for analyzing fast-collision dynamics in
linear physical systems with weak nonlinear dissipation. Indeed, as described
in the current paper and in Ref. PHN2020 , these methods enabled deep insight
into many collision-induced effects in the perturbed linear systems in spatial
dimensions 2 and 3. We also comment that detailed analytic results on
collisions between pulse (or beam) solutions of linear or nonlinear evolution
models in the presence of nonlinear dissipation in spatial dimension higher
than 1 are quite scarce. Therefore, the current work and the work in Ref.
PHN2020  also significantly extended the understanding of the more general
high-dimensional problem of two-pulse (or two-beam) collisions in the presence
of nonlinear dissipation.

## Acknowledgments

Q.M.N. and T.T.H. are supported by the Vietnam National Foundation for Science
and Technology Development (NAFOSTED) under Grant No. 107.99-2019.340.

##  Appendix A The relation between $\Delta A_{1}^{(c)}$ and
$\Delta\Phi_{1}(x,y,z_{c})$

In this Appendix, we present the derivation of the relation (17) between the
collision-induced amplitude shift $\Delta A_{1}^{(c)}$ and the collision-
induced change in the beam shape $\Delta\Phi_{1}(x,y,z_{c})$. This relation is
used in subsection II.2 to obtain Eq. (19) from Eq. (16). The derivation is
carried out for a collision in spatial dimension 2, but it can be generalized
in a straightforward manner to spatial dimension $n$.

We first recall that the amplitude dynamics of a single beam propagating in
the presence of diffraction and linear or nonlinear loss can be approximately
determined by an energy balance equation of the form
$\partial_{z}\int_{-\infty}^{\infty}\\!\\!\\!dx\int_{-\infty}^{\infty}\\!\\!\\!dy\,|\psi_{1}(x,y,z)|^{2}=...$,
where the right hand side of the equation is determined by the type of the
loss perturbation. [See, for example, Eqs. (109) and (112) in Appendix B]. A
fast collision in the presence of nonlinear loss at the distance $z=z_{c}$
leads to a jump in the value of
$\int_{-\infty}^{\infty}\\!\\!\\!dx\int_{-\infty}^{\infty}\\!\\!\\!dy\,|\psi_{1}(x,y,z)|^{2}$
at $z=z_{c}$. Therefore, for a fast collision, the term
$\partial_{z}\int_{-\infty}^{\infty}\\!\\!\\!dx\int_{-\infty}^{\infty}\\!\\!\\!dy\,|\psi_{1}(x,y,z)|^{2}$
in the equation that determines the amplitude dynamics can be replaced by:

$\displaystyle\Delta_{P}=\int_{-\infty}^{\infty}\\!\\!\\!dx\int_{-\infty}^{\infty}\\!\\!\\!dy\,|\psi_{1}(x,y,z_{c}^{+})|^{2}$