Curl of a second-order tensor

Assuming that by $e_{ijk}$ you mean the totally anti-symmetric tensor $\epsilon_{ijk}$, the expression you cite only is valid in three dimensions (since only in three dimensions $\epsilon_{ijk}$ has three indices). With the above assumption, the equation you provide can be implemented as follows

twotensorCurl3D[S_List] := Module[{},
 Id = IdentityMatrix[3];
 eps = LeviCivitaTensor[3];
 var = Table[x[i], {i, 1, 3}];
 Sum[eps[[i, j, k]] D[S[[m, j]], var[[i]]] TensorProduct[Id[[All, k]],Id[[All, m]]], {i, 1, 3}, {j, 1, 3}, {k, 1, 3}, {m, 1, 3}]
]

The bits of code mean the following. twotensorCurl3D[S_List] defines a function name and makes sure that only a list can be passed to the function. := as opposed to simply = makes sure that the right hand side is only evaluated once the function is actually called with a specific input. Module[{},...] is simply a wrapper that allows several computational steps to be done within the function before the results are put out. You could define local variables in the {} like {a,b,c}. We write the three dimensional identity matrix into the variable Id. Then the unit vectors $e_i$ can be accessed through the columns of this matrix: Id[[All,i]]. The position variables, in respect to which the nabla is taking derivatives, are defined as x[i] with i=1,2,3. The remaining code is self explanatory. The result of the last line in the Module is given back by the function after evaluation because it is not suppressed by use of ; at the end of the line.

Recently, I communicated with wolfram support about this issue and they send me back an illustrative notebook which compares the definition that you got from wikipedia with that of Mathematica. I will put the notebook here so you can download it and hope it will be useful for future readers of this post. Here is also some parts of the notebook

(*The constant vector c in wikipedia definition*)

c = {c1, c2, c3}

(*A function which constructs tensor fields of rank n*)

A[n_] := Array[Subscript[a, ##][x, y, z] &, ConstantArray[3, n]]

(*The curl definition for a 3 dimensional space with the use Mathematica's Curl command*)

Curl3D[array_, vars_] := 
 With[{n = ArrayDepth[array]}, 
  Transpose[Map[Curl[#, vars] &, array, {n - 1}], Reverse[Range[n]]]]

(*Checking the recursive identities that wikipedia used for defining curl*)

Expand[Curl3D[A[1], {x, y, z}].c] === Div[Cross[A[1], c], {x, y, z}]

Expand[Curl3D[A[2], {x, y, z}].c] === Expand[Curl3D[c.A[2], {x, y, z}]]

Expand[Curl3D[A[3], {x, y, z}].c] === Expand[Curl3D[c.A[3], {x, y, z}]]

(*Curl of a second order tensor in Cartesian coordinates*)

Expand[Curl3D[A[2], {x, y, z}]] // TableForm

The curl of a tensor field

Alas, the example for the curl of a tensor field in Mathematica 12 is not what is defined in the literature for continuum mechanics. The scalar result they have is the negative of half the trace of the curl.

Mathematica has sufficient functions to correctly compute the curl of a vector or tensor if the definitions given in the attached file are followed. For a second-order tensor, a single line command:

Transpose[Div[Dot[T[x,y,z], LeviCivitaTensor[3]], {x, y, z}]]

is all you need after defining the tensor T,

For example:

T[x_,y_,z_] := {{x y, x y^2, x y^3}, {x^2  y, x^2 y^2, x^2 y^3}, {x^3  y, x^3 y^2, x^3  y^3}}

using the example in Mathematica 12.