5 where λ$_{x}$ ( z ) and λ$_{y}$ ( z ) the normalized diffraction coeffi- cients corresponding to the two transverse directions: Next, we perform the variational analysis using the three dimensional Gaussian Ansatz in the following form: Following similar procedure as we used in the two- dimensional case [Sec. II C], we arrive at the set of vari- ational equations λ$_{x}$ ( z ) = cos( B sin z ) , λ$_{y}$ ( z ) = cos( B cos z ) , B = − x$_{0}$$_{m}$ 4 π $^{2}$n$_{0}$ dx 2 s λ$_{d}$z$_{s}$ . (28) u = A ( z ) exp { i φ ( z ) − 1 2 [ x 2 W $^{2}$( z ) + y 2 V $^{2}$( z ) + t 2 T $^{2}$( z ) ] + + i 2 [ b ( z ) x 2 + c ( z ) y 2 + β ( z ) t $^{2}$] } , (29) W ′ = λ$_{1}$bW, V ′ = λ$_{2}$cV, b ′ = λ$_{1}$ W 4 − λ$_{1}$b 2 − E 2 3 / $^{2}$W $^{3}$V T , (30) c ′ = λ$_{2}$ V 4 − λ$_{2}$c 2 − E 2 3 / $^{2}$W V $^{3}$T , T ′′ = 1 T 3 − E 2 3 / $^{2}$W V T $_{2}$, FIG. 3. An example of an unstable solution. Parameters are the same as in Fig. 2 except E = 0 . 6. along the y axis where E ≡ A $^{2}$W V T . We perform numerical analysis to identify periodic so- lutions of the variational equations (30). Such solutions correspond to three-dimensional light bullets, where the light-bullet extension along all dimensions ( x, y, t ) is re- stored after each bending period. We present an example of such solution in Fig. 4(b), which illustrates a stable and mobile three-dimensional light bullet. Indeed, the light bullet spatial and temporal localization widths are re- stored to the same values after each bending period.The solid line at the bottom shows the transverse position x$_{2}$ of the pulse center as a function of the propagation distance z for a finite initial beam tilt. The light bullet wobbles in the transverse direction during the propaga- tion, and it is clear that on average the pulse moves across the array. where y = y$_{d}$/y$_{s}$ is the normalized coordinate, and we assume the same scaling of the transverse spatial dimen- sions, y$_{s}$ = x$_{s}$ . Then, the coupling phase coefficients are defined similar to Eq. (3), where x$_{0}$ ( z ) and y$_{0}$ ( z ) define the transverse bending pro- file. We further consider the helical waveguide bending pro- file [see Fig. 4(a)], where x$_{0}$$_{m}$ = y$_{0}$$_{m}$ is the bending amplitude (normalized to x$_{s}$ ) and L is the normalized bending period. i ∂a$_{n,m}$ ∂z + C ( e − ip$_{x}$ ( z $^{)}$a$_{n}$$_{+1}$$_{,m}$ + e ip$_{x}$ ( z $^{)}$a$_{n}$$_{−}$$_{1}$$_{,m}$ ) + + C e − ip$_{y}$ ( z $^{)}$a$_{n,m}$$_{+1}$ + e ip$_{y}$ ( z $^{)}$a$_{n,m}$$_{−}$$_{1}$ + (24) + D 2 ∂ $^{2}$a$_{n,m}$ ∂t 2 + γ | a$_{n,m}$ | $^{2}$a$_{n,m}$ = 0 , p$_{x}$ ( z ) = 2 πn$_{0}$dx 2 s λ$_{d}$z$_{s}$ x ′ $_{0}$, p$_{y}$ ( z ) = 2 πn$_{0}$dx 2 s λ$_{d}$z$_{s}$ y ′ $_{0}$, (25) x$_{0}$ ( z ) = x$_{0}$$_{m}$ cos(2 πz/L ) , y$_{0}$ ( z ) = y$_{0}$$_{m}$ sin(2 πz/L ) , (26) i ∂u ∂z + λ$_{x}$ ( z ) 2 ∂ $^{2}$u ∂x 2 + λ$_{y}$ ( z ) 2 ∂ $^{2}$u ∂y 2 + σ$_{D}$ 2 ∂ $^{2}$u ∂t 2 + σ$_{γ}$ | u | $^{2}$u = 0 , (27) We summarize the existence and stability properties of three-dimensional light bullets in Figs. 4(c,d). We show the dependence between the energy E and the diffrac- tion map strength S on the phase increment Δ φ for one bending period. In comparison with the two-dimensional case [Fig. 1(d)], this dependence is much more complex, forming a loop-like structure. Interestingly, stable so- lutions, depicted with a solid line, are present for both positive and negative curve slopes of the E (Δ φ ) depen- dence, which indicates that the Vakhitov-Kolokolov cri- terion [36] is not applicable here. Following the approach described above in Sec. II B, we consider wavepackets that are slowly varying in space and describe their profiles with a continuous function, a$_{n,m}$ ( t, z ) = ˜ u ( x = nd, y = md, t, z ). Then, we apply the scaling coefficients according to Eq. (9), and derive the approximate continuous equation [cf. Eq. (10)]: