Dear (more advanced) users of Mathematica,

I'm still a beginner and currently trying to solve a system of delayed partial differential equations. Two of them describe the evolution of of two optical fields ($s_\text{P}$ and $s_\text{St}$) along a waveguide (position z):

(1) $\frac{\partial s_\text{P}\left(z,t\right)}{\partial z}=i \kappa_\text{S} R\left(z,t\right) s_\text{St}(z,t)$, and (2) $\frac{\partial s_\text{St}\left(z,t\right)}{\partial z}=i \kappa_\text{S} R^*\left(z,t\right) s_\text{P}\left(z,t\right)$.

These equation are coupled by the amplitude $R\left(z,t\right)$ of vibrations driven by the optical fields, which itself depends on time:

(3) $\frac{\partial R\left(z,t\right)}{\partial t} + \frac{\Gamma}{2} R\left(z,t\right) = i \kappa_\text{R} s_\text{P}\left(z,t\right) s_\text{St}^* \left(z,t\right)$.

I can obtain the steady-state solution ($\partial R/\partial t =0$) very easily by putting (3) into (1),(2) and applying NDSolve and the boundary conditions ($s_0\left(0,t\right) = \sqrt{0.9}$ and $s_\text{-1}(0,t)=\sqrt{0.1}$) to these two equations. But what I'd be interested in is the build-up of R over time. If I simply integrate (3) with $R(z,0) = 0$, I end up with an "delayed PDE" error.

Does anyone have an idea how I could proceed?

Thanks for your help in advance!

EDIT: This is my kind of naive code to solve this problem:

    (*DEFINE ALL THE CONSTANTS*)

Clear[ΩR, αS, κR, κS, γR, sar, c0];
ω0 = 2*Pi*2.99792458*10^8/(1.55*10^(-6));   (*optical frequency in Hz*)
c0 = 2.99792458*10^8;                       (*speed of light in vacuum*)
ΩR = 2*Pi*5.611*10^6;                       (*mechanical frequency in Hz*)
γR = 2*Pi*400;                              (*damping rate/linewidth in Hz*)
κS = -219.666*0.1;                          (*coupling rate*)
κR = κS*((c0*ΩR)/ω0)                        (*coupling rate*)  

(*SET UP PDEs TO SOLVE*)

Clear[sP, sSt, eqs, vars, bcs];
eqs = {{D[sP[z, t], z] == I κS*(I κR Exp[-γR*t/2]*
    Integrate[sP[z, d]*Conjugate[sSt[z, d]]*Exp[γR*d/2], {d, 0, t}])*sSt[z, t]},
       {D[sSt[z, t], z] == I κS*Conjugate[I κR Exp[-γR*t/2]*
    Integrate[sP[z, d]*Conjugate[sSt[z, d]]*Exp[γR*d/2], {d, 0, t}]]*sP[z, t]}}
vars = {sP[z, t], sSt[z, t]};      (*variables to solve for*)
bcs = {sP[0, t] == Sqrt[0.9], sSt[0, t] == Sqrt[0.1]};   (*boundary conditions at z=0*);
NDSolve[{eqs, bcs}, vars, {t, 0, 0.02}, {z, 0, 0.2}]

EDIT 2: Defining the integration of R as a function outside the NDsolve-block also seems to work.

(*DEFINE ALL THE CONSTANTS*)

Clear[ΩR, αS, κR, κS, γR, sar, c0];
ω0 = 2*Pi*2.99792458*10^8/(1.55*10^(-6));   (*optical frequency in Hz*)
c0 = 2.99792458*10^8;                       (*speed of light in vacuum*)
ΩR = 2*Pi*5.611*10^6;                       (*mechanical frequency in Hz*)
γR = 2*Pi*400;                              (*damping rate/linewidth in Hz*)
κS = -219.666*0.1;                          (*coupling rate*)
κR = κS*((c0*ΩR)/ω0)                        (*coupling rate*)  


R[t_, P_, S_] := I κR Exp[-γR*t/2]* Integrate[P*Conjugate[S]*Exp[γR*d/2], {d, 0, t}]

(*SET UP PDEs TO SOLVE*)

Clear[sP, sSt, eqs, vars, bcs];
eqs = {{D[sP[z, t], z] == 
    I κS*(R[t, sP[z, t], sSt[z, t]])*sSt[z, t]}, {D[sSt[z, t], 
     z] == I κS*Conjugate[R[t, sP[z, t], sSt[z, t]]]*sP[z, t]}}
vars = {sP[z, t],sSt[z, t]};(*variables to solve for*)
bcs = {sP[0, t] == Sqrt[0.9], sSt[0, t] == Sqrt[0.1]};(*boundary conditions for input fields at z=0*);
sol = NDSolve[{eqs, bcs}, vars, {t, 0, 0.003}, {z, 0, 0.1}]

(* GENERATE PLOTS *)
 
DensityPlot[sP[z, t]^2 /. sol, {t, 0, 0.003}, {z, 0, 0.1}, 
 ColorFunction -> "Rainbow", 
 FrameLabel -> {"Time (s)", "Propagation distance (m)"}, 
 PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 80]

DensityPlot[sSt[z, t]^2 /. sol, {t, 0, 0.003}, {z, 0, 0.1}, 
 ColorFunction -> "Rainbow", 
 FrameLabel -> {"Time (s)", "Propagation distance (m)"}, 
 PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 80]

DensityPlot[
 Abs[R[t, sP[z, t], sSt[z, t]]] /. sol, {t, 0, 0.003}, {z, 0, 0.1}, 
 ColorFunction -> "Rainbow", 
 FrameLabel -> {"Time (s)", "Propagation distance (m)"}, 
 PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 40]