  We consider the structure-dependent amplitude of the decay $B_s \to l^+ l^-
\gamma$ $(l=e,\mu)$ in a model based on the effective Hamiltonian for $b
\bar{s} \to l^+ l^-$ containing the Wilson coefficients $C_7,C_9$ and $C_{10}$.
The form factors characterising the matrix elements $< \gamma | \bar{s}
\gamma_\mu (1 \mp \gamma_5) b | \bar{B}_s>$ and $< \gamma | \bar{s}
\sigma_{\mu\nu} (1 \mp \gamma_5) b | \bar{B}_s>$ are taken to have the
universal form $f_V \approx f_A \approx f_T \approx f_{B_s} M_{B_s} R_s / (3
E_\gamma)$ suggested by recent work in QCD, where $R_s$ is a parameter related
to the light cone wave function of the $B_s$ meson. Simple expressions are
obtained for the charge asymmetry $A(x_\gamma)$ and the photon energy spectrum
$d \Gamma/ d x_\gamma (x_\gamma = 2 E_\gamma/M_{B_s})$. The decay rates are
calculated in terms of the decay rate of $B_s \to \gamma \gamma$. The branching
ratios are estimated to be $Br(B_s \to e^+ e^- \gamma) = 2.0 \times 10^{-8}$
and $Br(B_s \to \mu^+ \mu^- \gamma) = 1.2 \times 10^{-8}$, somewhat higher than
earlier estimates.
