Generation of sponteneous Brillouin scattering density fluctuation

I am trying to simulate the density fluctuation w.r.t time using FD method for stimulated brillouin scattering process originating from noise given by equation:

for a specific length at time 't' .

The noise term is a langevin noise source with zero mean Gaussian random variable. where .

However upon simulating this using finite difference (FD) for given 'z' the ρ is going to infinity which I am sure should not be the situation. The link for the paper I am reffering to: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.42.5514

Kindly suggest what is the problem. How to correct it?

I code is as follows:

clc
clear
close all

%% Parameters

kB=1.38e-23;    % Boltzmann Const.
rho0=2210;      % mean density in kg/m^3
Aeff=55.41e-12; % effective area of the medium in m^2
velo_aco=5960;  % acoustic velosity in m/s
Gamma=1e8;      % inverse of sound decay rate (sec^-1)
T=330;          % temperature in Kelvin (K)

%% defining the noise term

Q=(2*kB*T*rho0*Gamma)/(velo_aco^2*Aeff); % variance
dz=1;                                    % length in meter
Qt=(2*kB*T*rho0*Gamma)/(velo_aco^2*Aeff*dz); % % variance/meter
f=sqrt(Qt)*randn(1,1000)+2210;               % value of noise term    
% f=1e9*randn(1,1000);

j=sqrt(-1);
q=(2*pi/1550e-9)+(2*pi/1550.08e-9);
L_aco=10e9;
t=linspace(0,1,1000);
dt=t(2)-t(1);

%rho=0.5*2210+7e-4*(exp(j*(q*dz-L_aco.*t))+exp(-j*(q*dz-L_aco.*t)));
rho(1)=2210+7e-4;%*(exp(j*(q*dz-L_aco*0))+exp(-j*(q*dz-L_aco*0)));
rho1=zeros(1,1000);

%% claculation of density fluctuation with time

for i=1:length(t)-1
   rho(i+1)= (f(i)-0.5*Gamma*rho(i))*dt+rho(i);
   rho(i);
   t1=(f(i)-0.5*Gamma*rho(i))*dt
   rho1(i)=rho(i+1)-rho(i);
end

%% plots

figure(1)
plot(t,f);

figure(2)
plot(t,rho)

figure(3)
histogram(f,'Normal','pdf');
M = linspace(min(f),max(f),1000);                                     % range of samples to compute the theoretical pdf
fx_theory = pdf('normal',M,2210,sqrt(Qt));                            % theoretical normal probability density
hold on; plot(M,fx_theory,'r'); 

[David Goodmanson]

Hi Aritra,

things are a bit different with noise, but still I am assuming that the units of f are the same as the units of gamma*rho. But it doesn't seem to be working out that way.

[Aritra Paul]

Hi David,

Actually the value f is generated using expression given in the mentioned article (link). However I have noticed your point.

Thanks,

Aritra

[David Goodmanson]

Hi Aritra,

I don't have access to the paper. But even if I did, if everything was copied in verbatim and all the constants have units as you have listed them, the units still do not agree. For myself, until that was resolved I would always be asking if the units and the bad results are merely coincidental.

[Aritra Paul]

Thanks David,

I am looking into that and let you know as fast as possible.

[David Goodmanson]

HI Aritra,

no hurry, I am interested to know what you eventually find out