# Numerical Simulation of Noise in Pulsed Brillouin Scattering

Oscar A. Nieves  School of Mathematical and Physical Sciences, University of
Technology Sydney, 15 Broadway, Ultimo NSW 2007, Australia.  Corresponding
author<EMAIL_ADDRESS>Matthew D. Arnold  School
of Mathematical and Physical Sciences, University of Technology Sydney, 15
Broadway, Ultimo NSW 2007, Australia.  Michael. J. Steel  Macquarie University
Research Centre for Quantum Engineering (MQCQE), MQ Photonics Research Centre,
Department of Physics and Astronomy, Macquarie University, NSW 2109,
Australia.  Mikołaj K. Schmidt  Macquarie University Research Centre for
Quantum Engineering (MQCQE), MQ Photonics Research Centre, Department of
Physics and Astronomy, Macquarie University, NSW 2109, Australia.  Christopher
G. Poulton  School of Mathematical and Physical Sciences, University of
Technology Sydney, 15 Broadway, Ultimo NSW 2007, Australia.

###### Abstract

We present a numerical method for modelling noise in Stimulated Brillouin
Scattering (SBS). The model applies to dynamic cases such as optical pulses,
and accounts for both thermal noise and phase noise from the input lasers.
Using this model, we compute the statistical properties of the optical and
acoustic power in the pulsed spontaneous and stimulated Brillouin cases, and
investigate the effects of gain and pulse width on noise levels. We find that
thermal noise plays an important role in the statistical properties of the
fields, and that laser phase noise impacts the SBS interaction when the laser
coherence time is close to the time-scale of the optical pulses. This
algorithm is applicable to arbitrary waveguide geometries and material
properties, and thus presents a versatile way of performing noise-based SBS
numerical simulations, which are important in signal processing, sensing,
microwave photonics and opto-acoustic memory storage.

††journal: josab

##  1 Introduction

Stimulated Brillouin Scattering (SBS) is an opto-acoustic process that results
from the interaction between two counter-propagating optical fields, the pump
and the Stokes, as well as an acoustic wave inside a dielectric medium [1, 2,
3, 4, 5]. This interaction has been used for applications including narrow-
band radio-frequency (RF) and optical signal filtering [6, 7], phase
conjugation and precision spectroscopy [1], novel laser sources [8, 9], and in
recent experiments in opto-acoustic memory storage [10]. One of the key
challenges of simulating the SBS interaction is modelling of thermal noise,
which is present in all real systems and which can significantly affect
performance [11, 12, 13]. Simulating noise in the SBS equations is complicated
because of the nonlinear coupling between the envelope fields: beyond the
undepleted pump regime the noise is multiplicative and can only be understood
in the context of statistical moments using multiple independent realizations
[14]. Thermal noise in SBS has been simulated numerically in earlier studies
[11, 12], with these investigations concentrating on the noise properties of
the Stokes signal that arises spontaneously in response to a strong,
continuous-wave (CW) pump. More recent simulations [15] have incorporated both
thermal and laser noise in the SBS interaction, but have focused on single-
mode structures such as micro-ring resonators in steady-state laser
conditions. A numerical method for solving the transient SBS equations with
laser and thermal noise is needed for accurately predicting and characterizing
the noise in modern integrated SBS waveguide experiments [16, 2, 10].

In this paper, we present a numerical method by which the transient SBS
equations with thermal noise can be solved for pulses of arbitrary shape and
size, in arbitrary waveguide geometries. The method allows for the inclusion
of input laser noise in the form of stochastic boundary conditions. We apply
this method to the case of a short chalcogenide waveguide and use the model to
compute the statistics of the output envelope fields. We examine the dynamics
of the noise when the Stokes arises spontaneously from the thermal field, and
for the case when it is seeded with an input pulse at the far end of the
waveguide. We demonstrate the transition from the low-gain, short pulse case,
in which noise is amplified by the pump, to the high gain, long pulse regime
in which coherent amplification occurs. In this latter situation, we show that
while the output pulses remain smooth, significant fluctuations in the peak
powers arising from the thermal field can persist. We also show that, within
the framework of this model, phase noise from the pump only has a significant
impact on Stokes noise when the laser coherence time matches the time scales
of the pulses involved in the interaction. Finally, we investigate the
convergence of this numerical method, and find that it yields linear
convergence in both the average power and variance of the power for three
fields in the SBS interaction, which is in agreement with the Euler-Mayurama
scheme for solving stochastic ordinary differential equations.

##  2 Method

###  2.1 The SBS equations

We consider backward SBS interactions in a waveguide of finite length $L$
along the $z$-axis, in which a pump pulse with angular frequency $\omega_{1}$
is injected into the waveguide at $z=0$ and propagates in the positive
$z$-direction, while a signal pulse is injected at $z=L$ and propagates in the
negative $z$-direction, as shown in Fig. 1. The spectrum of the signal pulse
is centered around the Brillouin Stokes frequency
$\omega_{2}=\omega_{1}-\Omega$, which is down-shifted from the pump by the
Brillouin shift $\Omega$, and its spectral extent lies entirely within the
Brillouin linewidth $\Delta\nu_{B}$. When these two pulses interact, energy is
transferred from the pump to the signal via the acoustic field, resulting in
coherent amplification of the signal around the Brillouin frequency. At the
same time, as the pump moves through the waveguide, it interacts with the
thermal phonon field and generates an incoherent contribution to the Stokes
field which also propagates in the negative $z$-direction. This noisy Stokes
field combines with the coherent signal to form a noisy amplified output field
centered around the Stokes frequency. The interaction can be modelled using
three envelope fields for the pump ($a_{1}(z,t)$), Stokes ($a_{2}(z,t)$) and
acoustic field ($b(z,t)$), according to the equations [14]

$\displaystyle\frac{\partial a_{1}}{\partial z}+\frac{1}{v}\frac{\partial
a_{1}}{\partial t}+\frac{1}{2}\alpha a_{1}$
$\displaystyle=i\omega_{1}Q_{1}a_{2}b^{*},$ (1a) $\displaystyle\frac{\partial
a_{2}}{\partial z}-\frac{1}{v}\frac{\partial a_{2}}{\partial
t}-\frac{1}{2}\alpha a_{2}$ $\displaystyle=i\omega_{2}Q_{2}a_{1}b,$ (1b)
$\displaystyle\frac{\partial b}{\partial z}+\frac{1}{v_{a}}\frac{\partial
b}{\partial t}+\frac{1}{2}\alpha_{\text{ac}}b$ $\displaystyle=i\Omega
Q_{a}a_{1}^{*}a_{2}+\sqrt{\sigma}R(z,t).$ (1c)

Figure 1: Illustration of the SBS interaction, showing the pump, Stokes and
acoustic powers on a photonic chip waveguide. In (a), the pump and Stokes
pulses are injected into opposite ends of the waveguide, and the acoustic
field is made up of random thermal fluctuations. In (b), the optical fields
have interacted inside the waveguide, the Stokes depletes the pump to gain
some energy, and the rest of the energy goes to the acoustic field.

Here $\alpha$ and $\alpha_{\text{ac}}$ are the optical and acoustic loss
coefficients respectively (in units of m-1), along with optical group velocity
$v>0$ and acoustic group velocity $v_{a}>0$. The envelope fields $a_{1,2}$ and
$b$ have units of W1/2. The coefficients $Q_{1,2,a}$ represent the coupling
strength of the SBS interaction, which depend on the optical and acoustic
modes of the waveguide [17]; from local conservation of energy, we have
$Q_{2}=Q_{1}^{*}$ and $Q_{a}=Q_{1}$ [18]. Here we focus on the single
acoustic-mode case, which we can choose by tuning the laser frequencies and
relying on the large free spectral range of the acoustic modes. This model can
further be extended by including additional acoustic fields with their own
opto-acoustic coupling constants and potentially different noise properties
[18].

The boundary conditions for the pump and signal fields are applied by
specifying the input values $a_{1}(0,t)$ and $a_{2}(L,t)$ respectively. These
boundary conditions depend on the laser properties, such as the linewidth, and
may contain noise. Thermal noise in the waveguide is introduced through the
complex-valued space-time white noise function $R(z,t)$, which has mean
$\left\langle R(z,t)\right\rangle=0$ and auto-correlation function
$\left\langle
R(z,t)R^{*}(z^{\prime},t^{\prime})\right\rangle=\delta(z-z^{\prime})\delta(t-t^{\prime})$.
The noise strength is derived by analytically solving (1c) in the absence of
any optical fields [14], and is $\sigma=k_{B}T\alpha_{\text{ac}}$ where
$k_{B}$ is the Boltzmann constant and $T$ is the temperature of the waveguide.

We begin with the observation that the propagation distance of the acoustic
wave over the time-scale of the interaction is very small [11]. We therefore
apply the limit $\partial_{z}b\rightarrow 0$ in (1c), which simplifies to

$\frac{1}{v_{a}}\frac{\partial b}{\partial
t}+\frac{1}{2}\alpha_{\text{ac}}b=i\Omega
Q_{a}a_{1}^{*}a_{2}+\sqrt{\sigma}R(z,t).$ (2)

This has the formal solution

$b(z,t)=iv_{a}\Omega
Q_{a}\int_{-\infty}^{t}e^{-\frac{\Gamma}{2}(t-s)}a_{1}^{*}(z,s)a_{2}(z,s)\text{d}s+D(z,t),$
(3)

where $\Gamma=v_{a}\alpha_{\text{ac}}$ is decay rate of the acoustic field,
namely $\Gamma=1/\tau_{a}$, and is related to the Brillouin linewidth via
$\Gamma=2\pi\Delta\nu_{B}$. The thermal noise enters through the function

$D(z,t)=v_{a}\sqrt{\sigma}\int_{-\infty}^{t}\mathrm{e}^{-\frac{\Gamma}{2}(t-s)}R(z,s)\text{d}s.$
(4)

This function $D$ is a stochastic integral with zero mean $\left\langle
D(z,t)\right\rangle=0$ since the function $R(z,s)$ is itself a zero-mean
stochastic process. The auto-correlation function of $D$ at two times and two
points in space is found by following the derivation in [14], which uses the
stochastic Fubini theorem [19] to obtain the expression:

$\left\langle
D(z,t)D^{*}(z^{\prime},t^{\prime})\right\rangle=\frac{v_{a}\sigma}{\alpha_{\text{ac}}}\delta(z-z^{\prime})\exp\left\\{-\frac{\Gamma}{2}\left|t-t^{\prime}\right|\right\\}.$
(5)

Upon substitution of (3) into (1a) and (1b), and assuming that the fields
$a_{1,2}$ are everywhere zero for $t<0$, we obtain the pair of equations

$\frac{\partial a_{1}}{\partial z}+\frac{1}{v}\frac{\partial a_{1}}{\partial
t}+\frac{1}{2}\alpha a_{1}=i\omega_{1}Q_{1}a_{2}(z,t)D^{*}(z,t)\\\
-\frac{1}{4}g_{1}\Gamma
a_{2}(z,t)\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}a_{1}(z,s)a_{2}^{*}(z,s)\text{d}s,$
(6) $\frac{\partial a_{2}}{\partial z}-\frac{1}{v}\frac{\partial
a_{2}}{\partial t}-\frac{1}{2}\alpha a_{2}=i\omega_{2}Q_{2}a_{1}(z,t)D(z,t)\\\
-\frac{1}{4}g_{2}\Gamma
a_{1}(z,t)\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}a_{1}^{*}(z,s)a_{2}(z,s)\text{d}s,$
(7)

where $g_{1}=g_{0}\omega_{1}/\omega_{2}$, $g_{2}=g_{0}$, and the SBS gain
parameter $g_{0}=4v_{a}\omega_{2}\Omega|Q_{2}|^{2}/\Gamma$ (with units of
m-1W-1) [14].

The approach of the numerical method is to solve (6) and (7) in a stepwise
manner to find the optical fields; the optical fields at each time step are
then substituted into (3) to obtain the acoustic envelope field, and the
process is repeated. At each time step the solution requires calculation of
the thermal noise function $D(z,t)$ which behaves as a random walk in time
while remaining white in space. The optical equations are solved with the
input boundary conditions $a(0,t)=a_{p}(t)$ and $a(L,t)=a_{s}(t)$; in general,
these boundary conditions may be stochastic to account for noise in the input
lasers. In the following we first describe the approach taken to compute the
thermal noise function, then discuss the inclusion of noise into the boundary
conditions, before describing the iterative algorithm itself.

It should be noted that it is also possible to solve (1a), (1b) and (1c)
directly without integrating the acoustic envelope field in time first (as in
(3)), and this procedure would yield the same results. However, since the
thermal background field is assumed to be in an equilibrium state by $t=0$,
this alternative method would require simulating the acoustic envelope field
for a very long time $t<0$. This is computationally less efficient and poses
no advantages over the present method.

###  2.2 Computing the thermal noise function

The function $D(z,t)$ contains all the thermal noise information about the
system. To model $D(z,t)$ numerically, we note that its evolution in time
corresponds to an Ornstein-Uhlenbeck process [20]. Equation (4) can be written
in Itô differential form [21] as

$\text{d}D(z_{j},t)=-\frac{1}{2}\Gamma
D(z_{j},t)\text{d}t+v_{a}\sqrt{\sigma}R(z_{j},t)\text{d}t,$ (8)

where the $z$ axis is discretized on the equally spaced grid $z_{j}$ with
spacing $\Delta z$. We know that
$R(z_{j},t_{n})\text{d}t=\frac{1}{\sqrt{\Delta z}}\text{d}W_{j}(t_{n})$ where
$\text{d}W_{j}(t_{n})$ is the standard complex-valued Wiener increment in
time, and the scaling factor arises from the Dirac-delta nature of the
continuous-space auto-correlation function of $D(z,t)$. The complex increment
$\text{d}W_{j}(t)$ is a linear combination of two independent real Wiener
processes

$\text{d}W_{j}(t)=\frac{1}{\sqrt{2}}\left[\text{d}W_{j}^{(1)}(t)+i\text{d}W_{j}^{(2)}(t)\right],$
(9)

where
$\left\langle\text{d}W_{j}^{(p)}(t)\text{d}W_{j}^{(q)}(t)\right\rangle=\delta_{pq}\text{d}t$,
where $\delta_{pq}$ is the Kronecker delta. Integrating (8) from 0 to $t$
yields the analytic solution

$D(z_{j},t)=e^{-\frac{1}{2}\Gamma
t}D_{0}(z_{j})+v_{a}\sqrt{\frac{\sigma}{\Delta
z}}\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}\text{d}W_{j}(s),$ (10)

where $D_{0}(z_{j})$ is the cumulative random walk from $t=-\infty$ up to
$t=0$. This quantity is calculated using

$D_{0}(z_{j})=\frac{1}{\sqrt{2}}\left[\mathcal{N}_{z_{j}}^{(1)}\left(0,\frac{v_{a}\sigma}{\alpha_{\text{ac}}\Delta
z}\right)+i\mathcal{N}_{z_{j}}^{(2)}\left(0,\frac{v_{a}\sigma}{\alpha_{\text{ac}}\Delta
z}\right)\right],$ (11)

where $\mathcal{N}_{z_{j}}^{(1,2)}\left(0,v_{a}\sigma/\alpha_{\text{ac}}\Delta
z\right)$ are normal random variables with zero mean and variance
$v_{a}\sigma/(\alpha_{\text{ac}}\Delta z)$, independently sampled at each
$z_{j}$. Numerically, we can compute the integral in (10) following the
procedure in Appendix A. Thus, we simulate (10) as a random walk using
discrete increments $\Delta t$

$D(z_{j},t_{n+1})=e^{-\frac{1}{2}\Gamma\Delta t}D(z_{j},t_{n})\ \\\
+\gamma(\Delta
t)\left[\mathcal{N}_{z_{j},t_{n}}^{(1)}\left(0,1\right)+i\mathcal{N}_{z_{j},t_{n}}^{(2)}\left(0,1\right)\right],$
(12)

where

$\gamma(\Delta t)=v_{a}\sqrt{\frac{\sigma\left(1-e^{-\Gamma\Delta
t}\right)}{2\Delta z\Gamma}},$ (13)

and setting the initial value as $D(z_{j},t_{0})=D_{0}(z_{j})$. The random
numbers $\mathcal{N}_{z_{j},t_{n}}^{(1,2)}\left(0,1\right)$ are independently
sampled at each point $(z_{j},t_{n})$. Figure 2 shows multiple realizations of
$D(z_{j},t)$ at an arbitrary point $z_{j}$ and its ensemble average.

Figure 2: Multiple independent realizations (dashed grey) of the modulus
squared of the thermal function $|D(z_{j},t)|^{2}$ at an arbitrary position
$z_{j}$, the numerical ensemble average over 20 realizations (red) and the
analytic ensemble average (blue). We use a temperature of $T=300$ K,
$v_{a}=2500$ m/s, $\tau_{a}=5.3$ ns, $\Delta z=0.79$ mm and $\Delta t=6.43$
ps.

###  2.3 Noisy boundary conditions

Input laser noise can be an important feature in SBS experiments. In the
context of the SBS envelope equations, it enters in the form of random phase
fluctuations at the inputs of the waveguide, namely $z=0$ for the pump field
and $z=L$ for the Stokes field. We simulate this laser phase noise in the
input fields by modeling the boundary conditions as

$\displaystyle a_{1}(0,t_{n})$
$\displaystyle=a_{p}(t_{n})=\sqrt{P_{1}^{\text{in}}(t_{n})}e^{i\phi_{1}(t_{n})},$
(14) $\displaystyle a_{2}(L,t_{n})$
$\displaystyle=a_{s}(t_{n})=\sqrt{P_{2}^{\text{in}}(t_{n})}e^{i\phi_{2}(t_{n})},$
(15)

where $P_{1}^{\text{in}}(t)$ and $P_{2}^{\text{in}}(t)$ are deterministic
envelope shape functions for the pump and Stokes fields respectively
representing input power from the lasers. The variables $\phi_{1}(t)$ and
$\phi_{2}(t)$ are stochastic phase functions modeled as zero-mean independent
Brownian motions. The variation in the phase $\phi(t)$ is related to the
laser’s intrinsic linewidth $\Delta\nu_{L}$, or conversely the coherence time
$\tau_{\text{coh}}=1/(\pi\Delta\nu_{L})$, via the expression
$\left\langle\left[\phi(t+\tau)-\phi(t)\right]^{2}\right\rangle=2\pi\Delta\nu_{L}|\tau|$,
where $\tau=t^{\prime}-t$ for the two times $t^{\prime}$ and $t$ [22, 23, 24].
Following a similar numerical procedure to [25], we compute $\phi_{j}(t)$ as

$\phi_{j}(t)=\sqrt{2\pi\Delta\nu_{L}}\int_{0}^{t}\text{d}W_{j}(s),$ (16)

where $\text{d}W(s)$ is a real-valued Wiener process increment in time. To
generate the random walk numerically, we cast this integral as an Itô
differential equation
$\text{d}\phi_{j}(t)=\sqrt{2\pi\Delta\nu_{L}}\text{d}W_{j}(t)$, which is
discretized using an Euler-Mayurama [26] scheme as

$\phi_{j}(t_{n+1})=\phi_{j}(t_{n})+\sqrt{2\pi\Delta\nu_{L}}\sqrt{\Delta t}\
\mathcal{N}_{t_{n}}(0,1),$ (17)

where $\mathcal{N}_{t_{n}}(0,1)$ is a standard normally distributed random
number sampled at each $t_{n}$. A simulation of a single realization of the
noisy boundary conditions is shown in Figure 3.

Figure 3: Single realization of the noisy boundary conditions. The plots show
(a) pump power, (b) pump phase, (c) Stokes power and (d) Stokes phase. Both
pulses are Gaussian with FWHM of 2 ns. The laser linewidth used here is
$\Delta\nu_{L}=100$ MHz.

###  2.4 The numerical algorithm

We now present the main numerical algorithm of this paper. The algorithm
consists of two consecutive steps: first, we solve Eqs. (1a) and (1b) in the
absence of optical loss or nonlinear interactions. In other words, we solve
the following pair of advection equations

$\displaystyle\frac{\partial a_{1}}{\partial z}+\frac{1}{v}\frac{\partial
a_{1}}{\partial t}$ $\displaystyle=0,$ (18) $\displaystyle\frac{\partial
a_{2}}{\partial z}-\frac{1}{v}\frac{\partial a_{2}}{\partial t}$
$\displaystyle=0.$ (19)

With the boundary conditions $a_{1}(0,t)=a_{p}(t)$ and $a_{2}(L,t)=a_{s}(t)$,
these have the elementary solutions

$\displaystyle a_{1}(z,t)$ $\displaystyle=a_{p}\left(t-\frac{z}{v}\right),$
(20) $\displaystyle a_{2}(z,t)$
$\displaystyle=a_{s}\left(t-\frac{L-z}{v}\right).$ (21)

Setting the numerical grid parameter $\Delta z=v\Delta t$ further simplifies
(20) and 21 to

$\displaystyle a_{1}(z_{j},t_{n})\ $ $\displaystyle\leftarrow\
a_{1}(z_{j-1},t_{n-1}),$ (22) $\displaystyle a_{2}(z_{j},t_{n})\ $
$\displaystyle\leftarrow\ a_{2}(z_{j+1},t_{n-1}),$ (23)

such that the optical fields are shifted in space by exactly $\Delta z$ during
each time iteration. The envelope field $b(z,t)$ is assumed to remain
stationary in space during each time step, as is typical in the context of SBS
experiments involving pulses [14]. After the fields are shifted across the
waveguide, we solve the time evolution equations at each point $z_{j}$
independently; i.e. we solve

$\frac{1}{v}\frac{\partial a_{1}(z_{j},t)}{\partial t}=-\frac{1}{2}\alpha
a_{1}(z_{j},t)-\frac{1}{4}g_{1}\Gamma a_{2}(z_{j},t)I_{1,2}^{*}(z_{j},t)\\\
+i\omega_{1}Q_{1}a_{2}(z_{j},t)D^{*}(z_{j},t),$ (24)
$\frac{1}{v}\frac{\partial a_{2}(z_{j},t)}{\partial t}=-\frac{1}{2}\alpha
a_{2}(z_{j},t)+\frac{1}{4}g_{2}\Gamma a_{1}(z_{j},t)I_{1,2}(z_{j},t)\\\
-i\omega_{2}Q_{2}a_{1}(z_{j},t)D(z_{j},t),$ (25)

where the interaction integral $I_{1,2}(z_{j},t)$ is computed as

$I_{1,2}(z_{j},t_{n})=\frac{\Delta t}{2}e^{-\frac{\Gamma}{2}n\Delta
t}\bigg{[}I_{1,2}(z_{j},t_{n-1})\ \\\
+a_{1}^{*}(z_{j},t_{n-1})a_{2}(z_{j},t_{n-1})e^{\frac{\Gamma}{2}(n-1)\Delta
t}\ \\\ +a_{1}^{*}(z_{j},t_{n})a_{2}(z_{j},t_{n})e^{\frac{\Gamma}{2}n\Delta
t}\bigg{]}\leavevmode\nobreak\ .$ (26)

To integrate the envelope fields $a_{1}$ and $a_{2}$ in time, we use an Euler-
Mayurama scheme [27], which yields the following finite-difference equations

$a_{1}(z_{j},t_{n+1})=\left[1-\frac{v\alpha\Delta
t}{2}\right]a_{1}(z_{j},t_{n})\\\ -v\Delta t\left[\frac{g_{1}\Gamma
I_{1,2}^{*}(z_{j},t_{n})}{4}-i\omega_{1}Q_{1}D^{*}(z_{j},t_{n})\right]a_{2}(z_{j},t_{n}),$
(27) $a_{2}(z_{j},t_{n+1})=\left[1-\frac{v\alpha\Delta
t}{2}\right]a_{2}(z_{j},t_{n})\\\ +v\Delta t\left[\frac{g_{2}\Gamma
I_{1,2}(z_{j},t_{n})}{4}-i\omega_{2}Q_{2}D(z_{j},t_{n})\right]a_{1}(z_{j},t_{n}).$
(28)

The acoustic field is computed at each $z_{j}$ and $t_{n+1}$ after computing
$a_{1,2}$, via the equation

$b(z_{j},t_{n+1})=iv_{a}\Omega Q_{a}I_{12}(z_{j},t_{n+1})+\sqrt{\Delta
z}D(z_{j},t_{n+1}).$ (29)

The $\sqrt{\Delta z}$ factor in front of $D(z_{j},t_{n+1})$ ensures that the
variance of $b$ is independent of the numerical grid resolution.

Once all the fields are computed at $t_{n+1}$, we repeat the drift steps in
(22) and (23) and the entire process is iterated until the optical fields have
propagated across the waveguide. The steps of this numerical method are given
in Algorithm 1.

Algorithm 1 Numerical algorithm

1:Compute $D(z_{j},t_{n})$ for all $t_{n}$

2:Compute $\phi_{1,2}(t_{n})$ for all $t_{n}$

3:Set $a_{1,2}=0$ inside $z\in[0,L]$

4:for $n=1$ to $N_{t}-1$ do$\triangleright$ $N_{t}=$ size of time grid

5: Insert noisy boundary conditions in $a_{1,2}$ at $t_{n}$

6: Shift optical fields $a_{1,2}$ in space by $\Delta z$

7: Compute interaction integral $I_{1,2}(z_{j},t_{n})$

8: Compute $a_{1,2}(z_{j},t_{n+1})$ from $a_{1,2}(z_{j},t_{n})$

9: Compute $b(z_{j},t_{n+1})$ end for

###  2.5 Statistical properties of the fields

The iterative scheme in Algorithm 1 computes a single realization of the SBS
interaction given a specific set of input parameters. We must repeat this
process $M$ times with the same input parameters to build an ensemble of $M$
independent simulations, from which statistical properties may be calculated.
For instance, the true average of the power for all three fields ($P_{1,2}$
for the optical fields and $P_{a}$ for the acoustic field) may be calculated
as

$\left\langle
P_{1,2}(z_{j},t_{n})\right\rangle=\left\langle\left|a_{1,2}(z_{j},t_{n})\right|^{2}\right\rangle\approx\frac{1}{M}\sum_{m=1}^{M}\left|a_{1,2}^{(m)}(z_{j},t_{n})\right|^{2},$
(30) $\left\langle
P_{a}(z_{j},t_{n})\right\rangle=\left\langle\left|b(z_{j},t_{n})\right|^{2}\right\rangle\\\
\approx\frac{1}{M}\sum_{m=1}^{M}\left|b^{(m)}(z_{j},t_{n})\right|^{2},$ (31)

where $m$ refers to a specific realization of each process. Similarly, we
compute the standard deviation in the power at each point $(z_{j},t_{n})$ as

$\displaystyle\text{std}\left[P_{1,2}(z_{j},t_{n})\right]$
$\displaystyle=\sqrt{\left\langle\left[P_{1,2}(z_{j},t_{n})\right]^{2}\right\rangle-\left\langle
P_{1,2}(z_{j},t_{n})\right\rangle^{2}},$ (32)
$\displaystyle\text{std}\left[P_{a}(z_{j},t_{n})\right]$
$\displaystyle=\sqrt{\left\langle\left[P_{a}(z_{j},t_{n})\right]^{2}\right\rangle-\left\langle
P_{a}(z_{j},t_{n})\right\rangle^{2}}.$ (33)

The standard deviation is useful when comparing with experiments, since it
gives a quantitative measure of the size of the power fluctuations in the
measured optical fields.

##  3 Results and Discussion

We demonstrate the numerical method by simulating the SBS interaction of the
three fields with both thermal noise ($T=300$ K, $\Delta\nu_{B}=30$ MHz) and
laser noise ($\Delta\nu_{L}=100$ kHz), using a chalcogenide waveguide of
length 50 cm, with the properties in Table 1. Although our formalism includes
optical loss through the factor $\alpha$, we have chosen $\alpha=0$ in the
simulations to focus on the effect of net SBS gain and pulse properties on the
noise. Here we study the noisy SBS interaction in two different cases:
spontaneous scattering and stimulated scattering, and investigate the effects
of pump width and SBS gain on the noise properties of the Stokes field.

Parameter | Value
---|---
Waveguide length $L$ | 50 cm
Waveguide temperature $T$ | 300 K
Refractive index $n$ | 2.44
Acoustic velocity $v_{a}$ | 2500 m/s
Brillouin linewidth $\Delta\nu_{B}$ | 30 MHz
Brillouin shift $\Omega/2\pi$ | 7.7 GHz
Brillouin gain parameter $g_{0}$ | 423 m-1W-1
Optical wavelength $\lambda$ | 1550 nm
Laser linewidth $\Delta\nu_{L}$ | 100 kHz
Peak pump power | 1 W
Peak Stokes power | 0$-$1 mW
Simulation time $t_{f}$ | up to 80 ns
Pump pulse FWHM | 0.5$-$5 ns
Stokes pulse FWHM | 1 ns
Grid size (space) $N_{z}$ | 1001
Grid size (time) $N_{t}$ | 2601
Step-size $\Delta t$ | 4.07 ps
Table 1: Simulation parameters using a chalcogenide waveguide of the type
shown in [28].

###  3.1 The spontaneous Brillouin scattering case

We first consider the situation in which there is no input Stokes field from
an external laser source, and the Stokes arises purely from the interaction
between the pump and the thermal field — this situation is customarily
referred to as spontaneous or spontaneously-seeded Brillouin scattering. We
specify a Gaussian pump pulse of varying widths and constant peak power, with
input phase noise ($\Delta\nu_{L}=100$ kHz). Setting the waveguide temperature
at 300 K and the pump FWHM of 2 ns, in Fig. 4(a)$-$(c) we see that the thermal
acoustic field interacts with the pump to generate an output Stokes signal. At
the same time, the Stokes field depletes some of the pump and amplifies the
acoustic field, which leads to more Stokes energy being generated. The noisy
character of the Stokes field in Fig. 4(b) is due to the incoherent thermal
acoustic background, which generates multiple random Stokes frequencies. In
this short-pump regime, the SBS amplification is small, and the generated
Stokes field remains incoherent.

Figure 4: Waterfall plots for a single numerical realization of (a) pump
power, (b) Stokes power and (c) acoustic power in the spontaneous scattering
case, using a Gaussian pump of FWHM 2 ns and peak power of 1 W. Plots
(d)$-$(f) show the standard deviation of the field powers at each point
$(z,t)$, calculated from 100 independent realizations of the SBS interaction.

As we increase the width of the pump to 5 ns, the net SBS gain in the
waveguide also increases. In this long-pump regime, the (spontaneously-
generated) Stokes field is amplified coherently, as shown in Fig. 5(b).
However, it should be noted that, although the Stokes output becomes smooth,
there is significant variation in the peak Stokes power from one independent
realization to the next, as illustrated in Fig. 6(a) and (b). The standard
deviation of the Stokes power over multiple independent realizations increases
with longer pump pulses, as shown in Fig. 5(e).

Figure 5: Waterfall plots for a single numerical realization of (a) pump
power, (b) Stokes power and (c) acoustic power in the spontaneous scattering
case, using a Gaussian pump of FWHM 5 ns and peak power of 1 W. Plots
(d)$-$(f) show the standard deviation of the field powers at each point
$(z,t)$, calculated from 100 independent realizations of the SBS interaction.
Figure 6: Multiple independent realizations of the spontaneously generated
Stokes power across the waveguide for (a) 2 ns wide pump and (b) 5 ns wide
pump. These snapshots are taken at the time when the peak of the pump pulse
reaches $z=50$ cm.

As the pump becomes very long we approach the CW regime, in which the pump
power ramps up quickly at $z=0$ and is kept at a constant value. If the
waveguide is sufficiently long, the spontaneously generated Stokes field is
amplified coherently until pump depletion begins to take effect, initially at
$z=0$ and then throughout the length of the waveguide, until both Stokes and
pump fields relax into the steady-state configuration in which the pump
decreases exponentially, as shown in Fig. 7(a)$-$(b). When such a steady state
is reached, the depletion induced by the spontaneously-seeded Stokes may
inhibit Brillouin scattering from an input Stokes pulse injected at $z=L$.

Figure 7: Waterfall plots for a single numerical realization of (a) pump power
and (b) Stokes power in the spontaneous scattering case, using a CW pump with
1 W peak power. Figure 8: Computations of the spontaneously generated Stokes
field at $z=0$ over 500 independent realizations, using a 5 ns Gaussian pump
pulse, with varying input peak pump power. Plot (a) shows the output ensemble
averaged Stokes power at $z=0$, normalized by the maximum power at each input
pump, (b) shows the FWHM of the Stokes in time domain. Plot (c) shows the
normalized power spectral density (PSD) of the Stokes field, and (d) is the
FWHM of the Stokes in frequency domain.

Returning to the pulsed case, we investigate the effect of increasing the peak
pump power, and therefore the overall SBS gain, on the amplification of the
spontaneous Stokes field. Figure 8 shows how the Stokes spectral linewidth
increases for input pump powers between 0.1$-$2 W for a Gaussian pump pulse
with fixed FWHM of 5 ns. The increase in linewidth occurs due to the
transition from linear to nonlinear SBS amplification: in the linear
amplification regime, the spontaneously generated Stokes field retains a
constant temporal width while its peak power increases with input pump power.
In the nonlinear amplification regime, the Stokes field undergoes temporal
compression as a result of the central peak of the pulse being amplified
faster than the tails. Beyond 2 W of peak pump power, the spectral linewidth
of the Stokes field narrows as pump depletion becomes significant, because the
Stokes field is prevented from uniformly experiencing exponential gain
throughout the waveguide, an effect which is also observed in the CW pump case
[29].

###  3.2 The effect of laser phase noise

Our previous simulations included laser phase noise corresponding to a laser
linewidth of 100 kHz in the pump. This is equivalent to a coherence time of
$\tau_{\text{coh}}=3.2$ $\mu$s, which is at least 100 times larger than the
characteristic time of the SBS interaction in Fig. 4$-$7. For this reason it
is understandable that no contribution from the laser phase noise to the
optical or acoustic fields was observed. The contribution of laser phase noise
can however be observed if the linewidth of the pump is suffiently broad. We
therefore consider the CW-pump regime with zero Stokes input power, with a
laser linewidth of 100 MHz, which corresponds to a coherence time of 3.2 ns
(Fig. 9). We see a significant contribution from the laser phase noise in the
form of amplitude fluctuations, which are completely absent in the 100 kHz
linewidth case (Fig. 7). From this we infer that when the laser coherence time
$\tau_{\text{coh}}$ is comparable to the pulse widths $\tau_{p,s}$, the
fluctuations in the phase are fast enough to be transferred to the envelope of
the pulse. However, when $\tau_{\text{coh}}\gg\tau_{p,s}$, the noisy character
of the envelope fields will vanish. This has important implications for the
case of pulsed SBS: phase noise can only play a significant role in the
interaction if $\tau_{\text{coh}}\leq\tau_{p,s}$. For lasers with a relatively
small linewidth, such as in the kHz range, phase noise will only become a
significant effect when operating in the long-pulse or CW regime.

Figure 9: Waterfall plots for a single numerical realization of (a) pump power
and (b) Stokes power in the spontaneous scattering case, using a CW pump with
1 W peak power and a laser linewidth of 100 MHz.

###  3.3 The stimulated Brillouin scattering case

We now examine the case of seeded Brillouin scattering, in which a Stokes
signal is injected at $z=L$. We first consider a 1 mW peak power Stokes pulse
of FWHM 1 ns in the same chalcogenide waveguide as before. The pump is a
Gaussian pulse of constant peak power of 1 W, with a width of 2 ns. As can be
seen in Fig. 10, the Stokes pulse remains smooth throughout the interaction,
and although the standard deviation over 100 independent realizations is
approximately 1.4% of the peak value, there are no visible fluctuations in the
power across space or time in Fig. 10(b). A closer look at multiple individual
realizations in Fig. 11(a) reveals that there is a measurable level of
variation in the Stokes power, although each individual realization of the
Stokes field is smooth. By increasing the pump width to 5 ns as shown in Fig.
11(b), we also increase the standard deviation in the Stokes, however each
independent realization appears smoother compared to Fig. 11(a). This further
demonstrates how in the longer pump, high SBS gain regime, the amplification
of the Stokes is sufficient to cancel random phase differences in the Stokes
field, as we observed in the spontaneous scattering case in Fig. 5.

Figure 10: Waterfall plots for a single numerical realization of (a) pump
power, (b) Stokes power and (c) acoustic power in the stimulated scattering
case, using a Gaussian pump pulse of width 2 ns and peak power 1 W. The input
Stokes pulse has width 1 ns and peak power 1 mW. Plots (d)$-$(f) show the
standard deviation in the fields at each point $(z,t)$ for 100 independent
realizations of the SBS interaction. Figure 11: Multiple independent
realizations of the Stokes power across the waveguide for (a) 2 ns wide pump
and (b) 5 ns wide pump. These snapshots are taken at the time in which the
peaks of the pump and Stokes meet in the middle of the waveguide.

###  3.4 Convergence of the method

We now study the convergence of the numerical method by looking at the
statistical properties of the power in each field at fixed points on $(z,t)$.
We use a default minimum step-size in time $\Delta t_{\text{min}}=40.7$ fs
against which we compare the results for larger step-sizes $\Delta t$. We
compute the relative error in the power and variance of the power, taken over
1,000 independent realizations. These results correspond to what is known as
weak convergence in stochastic differential equations [26], where the mean
value of a random quantity, in our case the power, converges at a specific
rate with respect to the step-size used.

The results for the convergence computations are shown in Fig. 12. As expected
from the Euler-Mayurama scheme [26], the convergence rate is at most linear
for the mean power of all three fields. A similar rate of convergence is
recorded for the variance in each power, showing a one-to-one error reduction
with step-size. Although some higher order methods exist which implement
higher order Taylor expansions and Runge-Kutta schemes [26, 30, 31, 32], these
methods only work with ordinary stochastic differential equations; numerical
methods for partial stochastic differential equations are an active area of
research in applied mathematics [33].

Figure 12: Convergence plots showing the relative error in the ensemble
averaged powers (a)$-$(c) and in the variance of the powers (d)$-$(f), as a
function of the step-size $\Delta t$ used in the numerical grid. The reference
step-size used is $\Delta t_{\text{min}}=40.7$ fs. The calculations are based
on a sample size of 1,000 independent simulations of the fields. The test
problem consists of two optical Gaussian pulses for the pump and Stokes of
width 1 ns, with peak powers pump 100 mW (pump) and 10 $\mu$W (Stokes). The
statistical properties of $P_{1}$, $P_{2}$ and $P_{a}$ are calculated from
$P_{1}(L,t_{\text{max}})$, $P_{2}(0,t_{\text{max}})$ and
$P_{a}(L/2,t_{\text{max}})$ respectively, where $t_{\text{max}}$ is the time
at which the peaks of the optical pulses reach the opposite ends of the
waveguide. The computations include thermal noise in the waveguide at
temperature 300 K, and input laser phase noise with linewidth 100 kHz. The
waveguide properties are given in Table 1.

##  4 Conclusion

We have presented a numerical method by which the fully-dynamic coupled SBS
equations in both CW and pulsed scenarios with thermal and laser noise can be
solved. The method offers linear convergence in both the average power and
variance of the power of the optical and acoustic fields, with variances that
do not depend on step-size. From our simulations, we find that the noise
properties of the fields rely on the length of the optical pulses involved as
well as on the net SBS gain in the waveguide. For short-pump, low gain
regimes, the spontaneous Stokes field is incoherently amplified and exhibits
large spatial and temporal fluctuations, whereas for the long-pump, high gain
regime the field is amplified coherently, resulting in a smooth field but with
large variations in peak power between independent realizations. Similar
observations are made for the stimulated scattering case using a Stokes
signal. We also find that laser phase noise does not play a significant role
in the SBS interaction unless the laser coherence time is comparable to the
characteristic time-scales of the SBS interaction.

##  Appendix A Appendix A

The integral term in (10) can be evaluated using the properties of Itô
integrals. Firstly, since the integrand is a deterministic function of time,
and $\text{d}W_{j}(s)$ is a normally distributed stochastic process, the
integral is also a normally distributed stochastic process. Secondly,
$\text{d}W_{j}(s)$ is a complex-valued process, so the integral can be split
into two statistically independent real-valued integrals

$\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}\text{d}W_{j}(s)=\frac{1}{\sqrt{2}}\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}\text{d}W_{j}^{(1)}(s)\\\
+i\frac{1}{\sqrt{2}}\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}\text{d}W_{j}^{(2)}(s),$
(34)

each of these real integrals will have the same statistical properties, namely

$\left\langle\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}\text{d}W_{j}^{(q)}(s)\right\rangle=0.$
(35)

The variance is derived using the Itô isometry property for a stochastic
process $X(t)$ [34]

$\left\langle\left(\int_{0}^{t}X(s)\text{d}W(s)\right)^{2}\right\rangle=\left\langle\int_{0}^{t}X^{2}(s)\text{d}s\right\rangle.$
(36)

Using this property, we write

$\left\langle\left(\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}\text{d}W_{j}^{(q)}(s)\right)^{2}\right\rangle=\frac{1}{\Gamma}\left(1-e^{-\Gamma
t}\right),$ (37)

which leads to the result for the variance

$\text{Var}\left[\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}\text{d}W_{j}^{(q)}(s)\right]=\frac{1}{\Gamma}\left(1-e^{-\Gamma
t}\right).$ (38)

This means the integral can be computed as a normal random variable as

$\int_{0}^{t}e^{-\frac{\Gamma}{2}(t-s)}\text{d}W_{j}(s)\sim\sqrt{\frac{1-e^{-\Gamma
t}}{2\Gamma}}\left[\mathcal{N}_{z_{j},t}^{(1)}(0,1)+i{N}_{z_{j},t}^{(1)}(0,1)\right],$
(39)

which leads to (12).

## Disclosures

The authors declare no conflicts of interest.

## Acknowledgements

The authors acknowledge funding from the Australian Research Council (ARC)
(Discovery Project DP200101893), the Macquarie University Research Fellowship
Scheme (MQRF0001036) and the UTS Australian Government Research Training
Program Scholarship (00099F). Part of the numerical calculations were
performed on the UTS Interactive High Performance Computing (iHPC) facility.

## Data Availability Statement

Data underlying the results presented in this paper are not publicly available
at this time, but the data and accompanying code used to generate it can be
obtained from the authors on reasonable request.

## References

  * [1] B. J. Eggleton, C. G. Poulton, P. T. Rakich, M. J. Steel, and G. Bahl, “Brillouin integrated photonics,” Nature Photonics 13, 664–677 (2019).
  * [2] R. Pant, D. Marpaung, I. V. Kabakova, B. Morrison, C. G. Poulton, and B. J. Eggleton, “On-chip stimulated Brillouin scattering for microwave signal processing and generation,” Laser Photonics Rev. 8, 653–666 (2014).
  * [3] L. Brillouin, “Diffusion de la lumière et des rayons X par un corps transparent homogène,” AnPh 9, 88–122 (1922).
  * [4] R. W. Boyd, _Nonlinear optics_ (Elsevier, 2003).
  * [5] A. Kobyakov, M. Sauer, and D. Chowdhury, “Stimulated Brillouin scattering in optical fibers,” Adv. Opt. Photonics 2, 1–59 (2010).
  * [6] A. Choudhary, Y. Liu, D. Marpaung, and B. J. Eggleton, “On-chip Brillouin filtering of RF and optical signals,” IEEE J. Sel. Top. Quantum Electron. 24, 1–11 (2018).
  * [7] H. Jiang, L. Yan, W. Pan, B. Luo, and X. Zou, “Ultra-high speed RF filtering switch based on stimulated Brillouin scattering,” Opt. Lett. 43, 279–282 (2018).
  * [8] B. Steinhausser, A. Brignon, E. Lallier, J.-P. Huignard, and P. Georges, “High energy, single-mode, narrow-linewidth fiber laser source using stimulated Brillouin scattering beam cleanup,” Opt. Express 15, 6464–6469 (2007).
  * [9] W. Loh, A. A. Green, F. N. Baynes, D. C. Cole, F. J. Quinlan, H. Lee, K. J. Vahala, S. B. Papp, and S. A. Diddams, “Dual-microcavity narrow-linewidth Brillouin laser,” Optica 2, 225–232 (2015).
  * [10] M. Merklein, B. Stiller, K. Vu, S. J. Madden, and B. J. Eggleton, “A chip-integrated coherent photonic-phononic memory,” Nat. Commun. 8, 1–7 (2017).
  * [11] R. W. Boyd, K. Rzazewski, and P. Narum, “Noise initiation of stimulated Brillouin scattering,” Phys. Rev. A 42, 5514 (1990).
  * [12] A. L. Gaeta and R. W. Boyd, “Stochastic dynamics of stimulated Brillouin scattering in an optical fiber,” Phys. Rev. A 44, 3205 (1991).
  * [13] M. Ferreira, J. Rocha, and J. Pinto, “Analysis of the gain and noise characteristics of fiber Brillouin amplifiers,” Opt. Quantum Electron. 26, 35–44 (1994).
  * [14] O. A. Nieves, M. D. Arnold, M. Steel, M. K. Schmidt, and C. G. Poulton, “Noise and pulse dynamics in backward stimulated Brillouin scattering,” Optics Express 29, 3132–3146 (2021).
  * [15] R. O. Behunin, N. T. Otterstrom, P. T. Rakich, S. Gundavarapu, and D. J. Blumenthal, “Fundamental noise dynamics in cascaded-order Brillouin lasers,” Phys. Rev. A 98, 023832 (2018).
  * [16] W. Zhang and R. A. Minasian, “Widely tunable single-passband microwave photonic filter based on stimulated Brillouin scattering,” IEEE Photonics Technology Letters 23, 1775 (2011).
  * [17] B. C. Sturmberg, K. B. Dossou, M. J. Smith, B. Morrison, C. G. Poulton, and M. J. Steel, “Finite element analysis of stimulated Brillouin scattering in integrated photonic waveguides,” J. Light. Technol. 37, 3791–3804 (2019).
  * [18] C. Wolff, M. J. Steel, B. J. Eggleton, and C. G. Poulton, “Stimulated Brillouin scattering in integrated photonic waveguides: Forces, scattering mechanisms, and coupled-mode analysis,” Phys. Rev. A 92, 013836 (2015).
  * [19] M. Veraar, “The stochastic Fubini theorem revisited,” Stochastics 84, 543–551 (2012).
  * [20] G. E. Uhlenbeck and L. S. Ornstein, “On the theory of the Brownian motion,” Physical review 36, 823 (1930).
  * [21] N. G. Van Kampen, “Stochastic differential equations,” Physics reports 24, 171–228 (1976).
  * [22] B. Moslehi, “Analysis of optical phase noise in fiber-optic systems employing a laser source with arbitrary coherence time,” Journal of lightwave technology 4, 1334–1351 (1986).
  * [23] A. Debut, S. Randoux, and J. Zemmouri, “Linewidth narrowing in Brillouin lasers: Theoretical analysis,” Physical Review A 62, 023803 (2000).
  * [24] C. Wei, M. Zhou, Z. Hui-Juan, and L. Hong, “Stimulated Brillouin scattering-induced phase noise in an interferometric fiber sensing system,” Chinese Physics B 21, 034212 (2012).
  * [25] Y. Atzmon and M. Nazarathy, “Laser phase noise in coherent and differential optical transmission revisited in the polar domain,” Journal of lightwave technology 27, 19–29 (2009).
  * [26] P. Kloeden and E. Platen, _Numerical Solution of Stochastic Differential Equations_ (Springer-Verlag Berlin Heidelberg, 1992), 1st ed.
  * [27] X. Wang and S. Gan, “The tamed Milstein method for commutative stochastic differential equations with non-globally Lipschitz continuous coefficients,” Journal of Difference Equations and Applications 19, 466–490 (2013).
  * [28] Y. Xie, A. Choudhary, Y. Liu, D. Marpaung, K. Vu, P. Ma, D.-Y. Choi, S. Madden, and B. J. Eggleton, “System-level performance of chip-based Brillouin microwave photonic bandpass filters,” J. Light. Technol. 37, 5246–5258 (2019).
  * [29] A. L. Gaeta and R. W. Boyd, “Stochastic dynamics of stimulated brillouin scattering in an optical fiber,” Physical Review A 44, 3205 (1991).
  * [30] R. L. Honeycutt, “Stochastic Runge-Kutta algorithms. i. white noise,” Physical Review A 45, 600 (1992).
  * [31] R. L. Honeycutt, “Stochastic Runge-Kutta algorithms. ii. colored noise,” Physical Review A 45, 604 (1992).
  * [32] A. Tocino and R. Ardanuy, “Runge–kutta methods for numerical solution of stochastic differential equations,” Journal of Computational and Applied Mathematics 138, 219–241 (2002).
  * [33] Z. Zhang and G. Karniadakis, _Numerical methods for stochastic partial differential equations with white noise_ (Springer-Verlag, 2017), 1st ed.
  * [34] B. Øksendal, “Stochastic differential equations,” in _Stochastic differential equations,_ (Springer, 2003), pp. 29–33.


Liam Parker<EMAIL_ADDRESS>
Princeton University, Princeton, NJ, USA and Emre Onal<EMAIL_ADDRESS>
Princeton University, Princeton, NJ, USA and Anton Stengel
<EMAIL_ADDRESS>
Princeton University, Princeton, NJ, USA and Jacob Intrater
<EMAIL_ADDRESS>
Princeton University, Princeton, NJ, USA

# Neural Collapse in the Intermediate Hidden Layers of Classification Neural
Networks

###### Abstract

Neural Collapse ($\mathcal{NC}$) gives a precise description of the
representations of classes in the final hidden layer of classification neural
networks. This description provides insights into how these networks learn
features and generalize well when trained past zero training error. However,
to date, $\mathcal{NC}$ has only been studied in the final layer of these
networks. In the present paper, we provide the first comprehensive empirical
analysis of the emergence of $\mathcal{NC}$ in the intermediate hidden layers
of these classifiers. We examine a variety of network architectures,
activations, and datasets, and demonstrate that some degree of $\mathcal{NC}$
emerges in most of the intermediate hidden layers of the network, where the
degree of collapse in any given layer is typically positively correlated with
the depth of that layer in the neural network. Moreover, we remark that: (1)
almost all of the reduction in intra-class variance in the samples occurs in
the shallower layers of the networks, (2) the angular separation between class
means increases consistently with hidden layer depth, and (3) simple datasets
require only the shallower layers of the networks to fully learn them, whereas
more difficult ones require the entire network. Ultimately, these results
provide granular insights into the structural propagation of features through
classification neural networks.

##  1 Introduction

Modern, highly-overparameterized deep neural networks have exceeded human
performance on a variety of computer vision tasks [1, 2, 3]. However, despite
their many successes, it remains unclear how these overparameterized networks
converge to solutions which generalize well. In a bid to demystify neural
networks’ performance, a recent line of inquiry has explored the internally
represented ‘features’ of these networks during the Terminal Phase of Training
(TPT), i.e. when networks are trained past the point of zero error on the
training data [4, 5, 6].

The phenomenon of Neural Collapse ($\mathcal{NC}$), introduced by Papyan, Han,
and Donoho [7, 8], represents one such avenue of inquiry. $\mathcal{NC}$
refers to the emergence of a simple geometry present in neural network
classifiers that appears during TPT. Specifically, $\mathcal{NC}$ describes
the phenomena by which neural network classifiers converge to learning
maximally separated, negligible-variance representations of classes in their
last layer activation maps. However, despite the extensive documentation of
$\mathcal{NC}$ in the last-layer representations of classification neural
networks [9, 10, 11, 12], there has been no exploration of its presence
throughout the intermediate hidden layers of these networks.

In the present study, we provide a detailed account of the emergence of
$\mathcal{NC}$ in the intermediate hidden layers of neural networks across a
range of different settings. Specifically, we investigate the impact of
varying architectures, datasets, and activation functions on the degree of
$\mathcal{NC}$ present in the intermediate layers of classification networks.
Our results show that some level of $\mathcal{NC}$ typically occurs in these
intermediate layers in all explored settings, where the strength of
$\mathcal{NC}$ in a given layer increases as the depth of that layer within
the network increases. By examining the presence of $\mathcal{NC}$ not only in
the final hidden layer but also the intermediate hidden layers of
classification networks, we gain a more nuanced understanding of the
mechanisms that drive the behavior of these networks.

##  2 Methodology

###  2.1 Network Architecture and Training

We examine the degree of $\mathcal{NC}$ present in the hidden layers of three
different neural network classifiers. Two of these models are popular in the
computer vision community, and have been extensively studied and widely
adopted: VGG11 [13] and ResNet18 [14]. We also train a fully-connected (FC)
classification network, MLP6, with network depth of $\ell=6$ and layer width
of $d=4096$ for each of its hidden layers. MLP6 serves as a toy model in which
to more easily explore $\mathcal{NC}$. In addition to varying network
architecture, we also vary the activation functions. Specifically, we explore
the effects on $\mathcal{NC}$ of the ReLU, Tanh, and LeakyReLU activation
functions. We train these classification neural networks on five popular
computer vision classification datasets: MNIST [15], CIFAR10 and CIFAR100
[16], SVHN [17], and FashionMNIST [18]. To rebalance the datasets, MNIST and
SVHN were subsampled to $N=5,000$, $N=4,600$, and $N=600$ images per class,
respectively. We normalize all datasets but do not perform any data
augmentation. We use stochastic gradient descent with 0.9 momentum, Mean
Square Error Loss (MSE) 111We use MSE loss rather than Cross Entropy loss as
it has been shown to exhibit a greater degree of $\mathcal{NC}$ in
classification neural networks while maintaining a greater degree of
mathematical clarity [7, 19]., $\lambda=10^{-5}$ weight decay, and the one-
cycle learning rate scheduler for all training [20].

###  2.2 Intermediate Layer Analysis

To assess the extent of $\mathcal{NC}$ in the hidden layers of our
classifiers, we perform a step of “$\mathcal{NC}$ analysis” at various points
during training. This analysis involves freezing the network and passing all
the training samples through it. We then collect the network’s post-activation
representation of each sample after every hidden layer of interest.
Specifically, we collect post-activations after each FC layer in MLP6, each
convolutional layer in VGG11, and each convolutional layer in ResNet18. We
then flatten these post-activation representations into vectors
$\bm{h}^{j}_{i,c}$, where $j$ is the hidden layer, $i$ is the training sample,
and $c$ is the class. We then compute four quantities in each of these hidden
layer post-activation vectors to track $\mathcal{NC}$: intra-class variance
collapse $(\mathcal{NC}1)$, intra-class norm equality $(\mathcal{NC}2)$,
inter-class maximal angular separation $(\mathcal{NC}2)$, and simplification
to nearest-class center classifier $(\mathcal{NC}4)$ following the general
outline provided in [7]. The specifics of these quantities are provided in
Appendix A.

##  3 Results

We present the results for each of the $\mathcal{NC}$ conditions for MLP6,
VGG11, and ResNet18 with ReLU activation functions trained to TPT on
FashionMNIST, MNIST, SVHN, CIFAR10, and CIFAR100 in Figure 1 to Figure 4.

Figure 1: Intra-class Variance Collapse ($\mathcal{NC}1$) in the ReLU
classifiers’ intermediate hidden layers. The results are generated at various
points in training, where the blue dotted line indicates $\mathcal{NC}1$ at
initialization and the red solid line indicates $\mathcal{NC}1$ after TPT.

$\mathcal{NC}1$: The within-class variability decreases linearly with layer
depth in the shallower layers of the classification networks, indicating that
the earlier fully-connected/convolutional layers are equally effective in
clustering samples of the same class. However, in the deeper layers of the
networks, the within-class variability plateaus, suggesting that the network
has already maximally reduced the variance between same-class samples even
when they have only partially propagated through the network. This behavior is
observed in most tested network architectures and datasets, except for
CIFAR100. Ultimately, the earlier layers of the classifiers primarily group
same-class samples together and contribute to the generation of
$\mathcal{NC}1$.

Figure 2: Equal Norms ($\mathcal{NC}2$) in the ReLU classifiers’ intermediate
hidden layers. The results are generated at various points in training, where
the blue dotted line indicates $\mathcal{NC}2$ at initialization and the red
solid line indicates $\mathcal{NC}2$ after TPT. Figure 3: Maximal Angles
($\mathcal{NC}2$) in the ReLU classifiers’ intermediate hidden layers. The
results are generated at various points in training, where the blue dotted
line indicates $\mathcal{NC}2$ at initialization and the red solid line
indicates $\mathcal{NC}2$ after TPT.

$\mathcal{NC}2$: The two phenomena related to the emergence of the simplex ETF
structure in the class means, i.e. the convergence of class means to equal
norms and to maximal equal angles, also exhibit a somewhat linear relationship
between the degree of collapse in any given hidden layer and that layer’s
depth in the network. However, unlike $\mathcal{NC}1$, the collapse continues
to strengthen even in the deeper layers of the network, rather than plateauing
after the first few layers; this is more prevalent for the angular separation
between different class means than for the similarity between their norms.
This phenomenon also persists across most architectures and datasets, and
suggests that the network continues to separate classes as it feeds samples
forward through its full depth. This makes sense, as features extracted in the
shallower layers of the network can be used to learn more and more effective
ways of separating different-class samples in deeper layers, leading to an
increase in the recorded strength of $\mathcal{NC}2$ over layer depth.

Figure 4: Simplification to NCC ($\mathcal{NC}4$) in the ReLU classifiers’
intermediate hidden layers. The results are generated at various points in
training, where the blue dotted line indicates $\mathcal{NC}4$ at
initialization and the red solid line indicates $\mathcal{NC}4$ after TPT.

$\mathcal{NC}4$: The degree of $\mathcal{NC}4$ in any given layer during
training seems to be influenced by the presence of $\mathcal{NC}1$ and
$\mathcal{NC}2$ in that layer. For the nearest class mean to accurately
predict network output, samples need to be close to their class mean
($\mathcal{NC}1$) and class means should be well-separated ($\mathcal{NC}2$).
In most experiments, $\mathcal{NC}4$ decreases linearly in the shallower
layers and plateaus in the deeper layers. However, the plateauing occurs later
than $\mathcal{NC}1$ due to the additional angular separation between class
means in these deeper layers, which contributes valuable information for
classification. The degree of $\mathcal{NC}4$ in the $j$-th hidden layer
indicates how much of the network’s classification ability is captured by its
$\leq j$ hidden layers. If the mismatch between the nearest neighbor
classification in the $j$-th layer and the total network classification is
zero, then all of the network’s classification ability is already present by
the $j$-th layer. For simpler datasets like MNIST and SVHN, the NCC mismatch
between shallower layers and the network output reaches zero, while for more
complex datasets like CIFAR100, the NCC mismatch only reaches zero in the
final layers. This observation aligns with the notion that complex datasets
require the deeper networks for complete classification, while simpler
datasets can achieve it with the shallower layers; however, it will be
important for future studies to observe how this generalizes to
test/validation datasets.

However, despite these general trends, we note that the models trained on
CIFAR100 exhibit a number of unique behaviors not consistent with our broader
observations on $\mathcal{NC}$, most striking of which is no significant
decrease in within-class variability over training. Instead, the data seems
largely noisy for this $\mathcal{NC}$ condition. This merits future
investigation across other challenging datasets such as TinyImageNet, as well
as more exploratory analysis.

In addition to the experiments performed on the classification networks with
ReLU activation functions above, we also perform the same set of experiments
on Tanh and LeakyReLU classifiers in Appendix B and Appendix C respectively.
These experiments largely demonstrate the same characteristics as the ReLU
experiments.

##  4 Conclusions

Our work demonstrates that the Neural Collapse ($\mathcal{NC}$) phenomenon
occurs throughout most of the hidden layers of classification neural networks
trained through the terminal phase of training (TPT). We demonstrate these
findings across a variety of settings, including varying network
architectures, classification datasets, and network activations functions, and
find that the degree of $\mathcal{NC}$ present in any hidden layer is
typically correlated with the depth of that layer in the network. In
particular, we make the following specific observations: (1) almost all of the
reduction in intra-class variance in the samples occurs in the shallower
layers of the classification networks, (2) the angular separation between
class means is increased consistently as samples propagate through the entire
network, and (3) simpler datasets require only the shallower layers of the
networks to fully learn them, whereas more difficult ones require the entire
network. Ultimately, these results provide a granular view of the structural
propagation of features through classification networks. In future work, it
will be important to analyze how these results generalize to held-out
validation data. For example, do our observations of $\mathcal{NC}4$ serve as
a proxy for the $\leq j$ network’s ability to classify data extend to
validation data? Moreover, is there a broader relationship between
$\mathcal{NC}$ and network generalization/over-training?

##  5 Contributions

Parker: Initiated problem concept; led experiments and analysis design;
collaborated on computational work; collaborated on results analysis; wrote
the paper. Onal: Collaborated on experiments and analysis design; collaborated
on computational work; collaborated on results analysis. Stengel: Collaborated
on background research and initial experiment design; collaborated on
computational work. Intrater: Collaborated on problem concept and background
research.

##  6 Acknowledgements

We would like to thank Professor Boris Hanin for his generous guidance and
support in the completion of this paper as well as for his insightful and
exciting class, Deep Learning Theory, which ultimately led to the creation of
this project.

## References

  * [1] Ian Goodfellow, Yoshua Bengio and Aaron Courville  “Deep learning”  MIT press, 2016
  * [2] Alex Krizhevsky, Ilya Sutskever and Geoffrey E Hinton  “Imagenet classification with deep convolutional neural networks”  In _Communications of the ACM_ 60.6 AcM New York, NY, USA, 2017, pp. 84–90
  * [3] Yann LeCun, Yoshua Bengio and Geoffrey Hinton  “Deep learning”  In _nature_ 521.7553 Nature Publishing Group, 2015, pp. 436–444
  * [4] Tolga Ergen and Mert Pilanci  “Revealing the structure of deep neural networks via convex duality”  In _International Conference on Machine Learning_ , 2021, pp. 3004–3014  PMLR
  * [5] Cong Fang, Hangfeng He, Qi Long and Weijie J Su  “Exploring deep neural networks via layer-peeled model: Minority collapse in imbalanced training”  In _Proceedings of the National Academy of Sciences_ 118.43 National Acad Sciences, 2021, pp. e2103091118
  * [6] Zhihui Zhu et al.  “A geometric analysis of neural collapse with unconstrained features”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 29820–29834
  * [7] Vardan Papyan, XY Han and David L Donoho  “Prevalence of neural collapse during the terminal phase of deep learning training”  In _Proceedings of the National Academy of Sciences_ 117.40 National Acad Sciences, 2020, pp. 24652–24663
  * [8] XY Han, Vardan Papyan and David L Donoho  “Neural collapse under mse loss: Proximity to and dynamics on the central path”  In _arXiv preprint arXiv:2106.02073_ , 2021
  * [9] Dustin G. Mixon, Hans Parshall and Jianzong Pi  “Neural collapse with unconstrained features”  In _arXiv_ , 2020  DOI: 10.48550/ARXIV.2011.11619
  * [10] Jianfeng Lu and Stefan Steinerberger  “Neural Collapse with Cross-Entropy Loss”  In _arXiv_ , 2020  DOI: 10.48550/ARXIV.2012.08465
  * [11] Tom Tirer and Joan Bruna  “Extended Unconstrained Features Model for Exploring Deep Neural Collapse”  In _arXiv_ , 2022  DOI: 10.48550/ARXIV.2202.08087
  * [12] Jinxin Zhou et al.  “Are All Losses Created Equal: A Neural Collapse Perspective”  In _arXiv_ , 2022  DOI: 10.48550/ARXIV.2210.02192
  * [13] Karen Simonyan and Andrew Zisserman  “Very Deep Convolutional Networks for Large-Scale Image Recognition”  arXiv, 2014  DOI: 10.48550/ARXIV.1409.1556
  * [14] Karen Simonyan and Andrew Zisserman  “Very deep convolutional networks for large-scale image recognition”  In _arXiv preprint arXiv:1409.1556_ , 2014
  * [15] Yann LeCun  “The MNIST database of handwritten digits”  In _http://yann. lecun. com/exdb/mnist/_ , 1998
  * [16] Alex Krizhevsky and Geoffrey Hinton  “Learning multiple layers of features from tiny images”  Toronto, ON, Canada, 2009
  * [17] Yuval Netzer et al.  “Reading digits in natural images with unsupervised feature learning”, 2011
  * [18] Han Xiao, Kashif Rasul and Roland Vollgraf  “Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms”  In _arXiv preprint arXiv:1708.07747_ , 2017
  * [19] Vardan Papyan  “Measurements of three-level hierarchical structure in the outliers in the spectrum of deepnet hessians”  In _arXiv preprint arXiv:1901.08244_ , 2019  URL: https://arxiv.org/pdf/1901.08244.pdf
  * [20] Leslie N. Smith and Nicholay Topin  “Super-Convergence: Very Fast Training of Neural Networks Using Large Learning Rates”, 2018  arXiv:1708.07120 [cs.LG]

##  Appendix A Intermediate Layer Analysis

In this section, we provide greater details on how the various $\mathcal{NC}$
conditions are computed from the flattened $\bm{h}^{j}$ post-activation
vectors after they have been randomly sampled (if they have been randomly
sampled at all). Specifically, from these post-activation vectors, we first
compute four preliminary quantities for each layer $j$, which we will use
later to calculate the $\mathcal{NC}$ metrics:

  1. 1.

$j$-th layer global mean:
$\bm{\mu}_{G}^{j}=\operatorname{Ave}^{j}_{i,c}\bm{h}_{i,c}^{j}$,

  2. 2.

$j$-th layer class means:
$\bm{\mu}_{c}^{j}=\operatorname{Ave}^{j}_{i}\bm{h}_{i,c}^{j}$,

  3. 3.

$j$-th layer within-class covariance:
$\Sigma_{W}^{j}=\operatorname{Ave}^{j}_{i}(\bm{h}_{i,c}^{j}-\bm{\mu}_{c}^{j})(\bm{h}_{i,c}^{j}-\bm{\mu}_{c}^{j})^{T}\in\mathbb{R}^{p_{j}\times
p_{j}}$,

  4. 4.

$j$-th layer between-class covariance:
$\Sigma_{B}^{j}=\operatorname{Ave}_{c}^{j}(\bm{\mu}_{c}^{j}-\bm{\mu}_{G}^{j})(\bm{\mu}_{c}^{j}-\bm{\mu}_{G}^{j})^{T}\in\mathbb{R}^{p_{j}\times
p_{j}}$,

From these quantities, we compute the following for each of the hidden layers
across the entire training set, following the general outline provided by [7].
Importantly, we do not track $\mathcal{NC}3$, as it does not yield any
meaningfully intuitive insights in the context of intermediate layer analysis:

$(\mathcal{NC}1)$ Intra-Class Variance Collapse:

$\operatorname{Tr}(\bm{\Sigma}_{B}^{j+}\bm{\Sigma}_{W}^{j}/C).$

Here $+$ denotes the pseudoinverse calculated with singular value
decomposition (SVD). Ignoring the first term, the variation becomes negligible
as the weights approach their class means; thus, the vanishing of this
quantity implies the vanishing of intra-class variance in the intermediate
activation maps. This quantity is a more precise version of the one introduced
in [7].

$(\mathcal{NC}2)$ Convergence to Simplex ETF: The first quantity, ”Equal
Norms”, is

$\operatorname{std}^{j}_{c}(\|\bm{\mu}_{c}^{j}-\bm{\mu}_{G}^{j}\|_{2})/\operatorname{avg}^{j}_{c}(\|\bm{\mu}_{c}^{j}-\bm{\mu}_{G}^{j}\|_{2}).$

As this quantity converges to zero, the coefficients of variation of class
means similarly vanish, implying that the class-means become equinormal. The
second quantity, ”Maximal Angles”, is

$\operatorname{avg}_{c\neq
c^{\prime}}\left(\left|\frac{\langle\bm{\mu}_{c}^{j}-\mu_{G}^{j}\bm{\mu}_{c^{\prime}}^{j}-\bm{\mu}_{G}^{j}\rangle}{\|\bm{\mu}_{c}^{j}-\bm{\mu}_{G}^{j}\|_{2}\|\bm{\mu}_{c^{\prime}}^{j}-\bm{\mu}_{G}^{j}\|_{2}}+\frac{1}{C-1}\right|\right).$

The first component in the average denotes $\cos_{\mu}(c^{j},c^{\prime j})$,
i.e. the angle between the average activation map of class $c$ and that of
class $c^{\prime}$ in layer $j$. Thus this represents the average difference
between these angles and $\frac{-1}{C-1}$, and therefore the convergence to
zero corresponds to maximum separation for globally centered vectors.
Together, these terms represent convergence of the class means to a Simplex
ETF.

$(\mathcal{NC}4)$ Simplification to Nearest Class-Center:

$1-\operatorname{avg}_{i,c}\mathbb{1}\\{f_{\theta}(x_{i,c})=\arg\min_{c^{\prime}}\|\bm{h}^{j}_{i,c}-\bm{\mu}_{c^{\prime}}^{j}\|_{2}\\}$

Here $1\\{.\\}$ is the indicator function. This quantity represents the
proportion of times that the entire network $f_{\theta}$ disagrees with the
result that would have been obtained by simply taking the nearest class mean
for an activation map in the $j$th layer. As this term goes to zero, it
implies that the classifier’s behavior simplifies to a nearest class-mean
decision rule. Importantly, as this quantity vanishes in intermediate layers,
it implies that the features necessary for the classifier to correctly
classify

##  Appendix B Tanh Activations

We present the results for each of the $\mathcal{NC}$ conditions for the
various architectures with Tanh activation functions trained to TPT on the
various datasets in Figure 9 to Figure 12.

Figure 5: Intra-class Variance Collapse ($\mathcal{NC}1$) in the Tanh
classifiers’ intermediate hidden layers. The results are generated at various
points in training, where the blue dotted line indicates $\mathcal{NC}1$ at
initialization and the red solid line indicates $\mathcal{NC}1$ after TPT.
Figure 6: Equal Norms ($\mathcal{NC}2$) in the Tanh classifiers’ intermediate
hidden layers. The results are generated at various points in training, where
the blue dotted line indicates $\mathcal{NC}2$ at initialization and the red
solid line indicates $\mathcal{NC}2$ after TPT. Figure 7: Maximal Angles
($\mathcal{NC}2$) in the Tanh classifiers’ intermediate hidden layers. The
results are generated at various points in training, where the blue dotted
line indicates $\mathcal{NC}2$ at initialization and the red solid line
indicates $\mathcal{NC}2$ after TPT. Figure 8: Simplification to NCC
($\mathcal{NC}4$) in the Tanh classifiers’ intermediate hidden layers. The
results are generated at various points in training, where the blue dotted
line indicates $\mathcal{NC}4$ at initialization and the red solid line
indicates $\mathcal{NC}4$ after TPT.

##  Appendix C LeakyReLU Activations

We present the results for each of the $\mathcal{NC}$ conditions for the
various architectures with LeakyReLU activation functions trained to TPT on
the various datasets in Figure 9 to Figure 12.

Figure 9: Intra-class Variance Collapse ($\mathcal{NC}1$) in the LeakyReLU
classifiers’ intermediate hidden layers. The results are generated at various
points in training, where the blue dotted line indicates $\mathcal{NC}1$ at
initialization and the red solid line indicates $\mathcal{NC}1$ after TPT.
Figure 10: Equal Norms ($\mathcal{NC}2$) in the LeakyReLU classifiers’
intermediate hidden layers. The results are generated at various points in
training, where the blue dotted line indicates $\mathcal{NC}2$ at
initialization and the red solid line indicates $\mathcal{NC}2$ after TPT.
Figure 11: Maximal Angles ($\mathcal{NC}2$) in the LeakyReLU classifiers’
intermediate hidden layers. The results are generated at various points in
training, where the blue dotted line indicates $\mathcal{NC}2$ at
initialization and the red solid line indicates $\mathcal{NC}2$ after TPT.
Figure 12: Simplification to NCC ($\mathcal{NC}4$) in the LeakyReLU
classifiers’ intermediate hidden layers. The results are generated at various
points in training, where the blue dotted line indicates $\mathcal{NC}4$ at
initialization and the red solid line indicates $\mathcal{NC}4$ after TPT.