J. Eur. Opt. Society-Rapid Publ. 21, 23( 2025) 233
68 Massaro LM, et al., Heterogeneous random laser with switching activity visualized by replica symmetry breaking maps, ACS Photonics 8, 376( 2021). https:// doi. org / 10.1021 / acsphotonics. 0c01803.
69 Ducci S, Ramazza PL, González-Viñas W, Arecchi FT, Order parameter fragmentation after a symmetry-breaking transition, Phys Rev Lett. 83, 5210( 1999). https:// doi. org / 10.1103 / PhysRevLett. 83.5210.
70 Hendry I, et al., Spontaneous symmetry breaking and trapping of temporal Kerr cavity solitons by pulsed or amplitude-modulated driving fields, Phys. Rev. A. 97, 053834( 2018). https:// doi. org / 10.1103 / physreva. 97.053834.
71 Butcher JC, A history of Runge-Kutta methods, Appl. Numer. Math. 20, 247( 1996). https:// doi. org / 10.1016 / 0168-9274( 95) 00108-5.
72 Bezanson J, Edelman A, Karpinski S, Shah VB, Julia: A fresh approach to numerical computing, SIAM Rev. 59, 65( 2017). https:// doi. org / 10.1137 / 141000671.
Appendix: Numerical implementation
A. 1 Normalization
In the numerical implementation of the dynamics of equations( 3) and( 4), it is important to adapt a proper normalization. In order to minimize the computational cost and to maximize the accuracy, good prescription is that the key parameters and variables assume values close to the unity. In our simulations, we normalized the system according to the following substitutions:
a 0 1 ¼ a 1 = a 1 ¼ 1; a 0 2 ¼ a 2 = a 1; d = a 1; qffiffiffiffiffi qffiffiffi t 0 ¼ a 1 Lt = t r; s 0 2a
¼ 1 jb 2 s; A 0; B 0 c
¼ 1 j a 1
A; B; b 0 FFðSHÞ ¼ b FFðSHÞ a 1
; b 0 1; FFðSHÞ ¼ q ffiffiffiffiffiffiffiffiffiffiffiffi
2 a 1jb 2; FFj b 1; FFðSHÞ; b 0 2; FFðSHÞ ¼ b 2; FFðSHÞ = jb 2; FF j c 0 1 ¼ c 1 = c 1 ¼ 1 qffiffiffiffiffiffi c 0 2; 12; 21 ¼ c 2; 12; 21 = c 1; j 0 1
¼ c 1 a 1 j; q S 0 ¼ ffiffiffiffiffiffiffi S;
c 1 a 3 1 L2 ðA1Þ
where we have defined the normalized detuning D. The resulting coupled equations are formally identically to the system of equations( 3) and( 4), except for the prefactor t R / L on the LHS, which vanishes. However, subject to this substitution, the system is dimensionless and typical solutions are of the order 1. It follows that the machine virtual memory used is minimized, and both the CPU efficiency and accuracy are optimized.
A. 2 Numerical algorithm
In our simulations we make used of a standard split-step Fourier algorithm [ 49, 54 ]. In order to illustrate the procedure, we rewrite the system according to:
@ U
¼ LU þ NU þ S; @ t ðA2Þ
where the linear L, nonlinear N and pump S operators are defined, in normalized units, by comparison with equations( 3) and( 4):
U ¼ A0
; ðA3Þ
B 0
0 �a 0 1 � i � i b0 2; FF @ 2
2 LU ¼
@ s; 0 2
B @
�a 0 2 � i2 � @ b0 1 @ s � i b0 2; SH 0 2
@ 2
@ s 0 2
1 C A U;
0 h i
1 ic 0 1
A0 j 2 þ 2ic 0
12B 0 j 2 A 0 þ ij 0 B 0 A 0
B C
NU ¼ @ h i
A; ic 0 2 jB0 j 2 þ 2ic 0 21 jA0 j 2 B 0 þ ij 0 A 02
S ¼ S0
0 ðA4Þ ðA5Þ
: ðA6Þ
The two-steps algorithm works in both the time and Fourier domain in a single discretization step Dt, which connects two consecutive instants t k and t k + 1 = t k + Dt. The linear block @ U L =@ t ¼ LU þ S is solved analytically in the Fourier domain. The field U is first transformed to ~ U by a standard fast Fourier transform( FFT) numerical operator, i. e. ~ Uðt k; xÞ FFTfUðt k; sÞg. Then, the linear contribution to the evolution of the jth component of the field at the kth step is computed by:
U~ L; j ðt kþ1; xÞ ¼ U~ L j ðt k; xÞe ~ S jt þ ~
j L e ~ jt � 1: ðA7Þ
L~ j
This expression is valid for the two-field components j = 1,2, and we can note that the second term on the RHS is null for j = 2. Next, we transform back ~ U L; j ðt kþ1; xÞ to the time domain by means of an inverse fast Fourier transform operation( IFFT), i. e. U L; j ðt kþ1; sÞ ¼IFFTð ~ U L; j ðt kþ1; xÞÞ. At this point, we calculate numerically the nonlinear contribution by means of a 4th order Runge-Kutta solver( RGK4) [ 71 ]:
U N; j ðt kþ1; sÞ ¼RGK4ðNU L; j ðt kþ1; sÞÞ: ðA8Þ
The two-steps( A7) and( A8) are subject to the condition that the linear L and nonlinear N operators act separately and independently on the two waves field U. This condition is mathematically fulfilled if and only if their commutator vanishes, i. e. ½ L; N ŠU ¼ 0 [ 49 ]. This is true in the limit Dt? 0, and thus the accuracy of procedure crucially relies on the choice of small enough discretization steps Dt.
A. 3 Dispersion effects
The contribution of chromatic dispersion effects is evaluated in the Fourier domain and embedded in the linear operator LðxÞ ~. The system of equations( 3) and( 4) neglects leading order terms with respect to the GVD. Within this approximation, the operator LðxÞ ~ reads as:
�a 0 1 LðxÞ ~ � i þ! ib0 2; FF x2
�a 0 2 � 2i � b0 1 x þ: ðA9Þ ib0 2; SH x2
This approximation is optimal in a narrow band spectral domain, while for broadband dynamics can be detrimental. In our case, we have a doubly resonant system which owns spectral components covering the whole octave spanning FF? SH domain, with important dispersive effects. For this reason, the use of a full dispersion profile is of crucial importance for the implementation of reliable and realistic models. In multi-envelope systems [ 39 ], we can expand