code: def _expected_hessian_sphere_immersion(point):
      theta = point[..., 0]
      phi = point[..., 1]
      hessian_immersion_x = gs.array(
          [
              [
                  -gs.sin(theta) * gs.cos(phi),
                  -gs.cos(theta) * gs.sin(phi),
              ],
              [
                  -gs.cos(theta) * gs.sin(phi),
                  -gs.sin(theta) * gs.cos(phi),
              ],
          ]
      )
      hessian_immersion_y = gs.array(
          [
              [
                  -gs.sin(theta) * gs.sin(phi),
                  gs.cos(theta) * gs.cos(phi),
              ],
              [
                  gs.cos(theta) * gs.cos(phi),
                  -gs.sin(theta) * gs.sin(phi),
              ],
          ]
      )
      hessian_immersion_z = gs.array(
          [[-gs.cos(theta), 0.0], [0.0, 0.0]]
      )
      hessian_immersion = gs.stack(
          [
              hessian_immersion_x,
              hessian_immersion_y,
              hessian_immersion_z,
          ],
          axis=0,
      )
      return hessian_immersion
  
imports: import geomstats.backend
  from geomstats.test.data import TestData
  from .riemannian_metric import RiemannianMetricComparisonTestData, RiemannianMetricTestData
label_desc: Use image node_id 3 for calling a global function with example usage: _expected_hessian_sphere_immersion(point) and returns: hessian_immersion