High-order differentiation of forms with external operators

There is an issue concerning the high-order differentiation of forms that contain `ExternalOperator`-s.

We start with a standard form and its first derivative
```python
N = ExternalOperator(u, function_space=V_ufl) # N(u; v')
E = N * dx # N(u; v') * dx
F = derivative(E, u, u_test)
F_expanded = expand_derivatives(F) # Action( v1 * dx; dN(u; v1, v'))
```
Then we want to differentiate `F` or `F_expanded`, but unfortunately, there is an error in `formoperators` when we differentiate `Action` (see the error below)
```python
# Differentiation expansion -> Second derivative -> Second expansion
J_expanded = derivative(F_expanded, u, u_trial)
J_expanded_expanded = expand_derivatives(J_expanded) # Error
```

Without the `F` expansion, this works, but this is worse since one gets 0, which is not correct:
```python
# Double differentiation
J = derivative(F, u, u_trial)
J_expanded = expand_derivatives(J) # == 0 <--- Works, but incorrectly
```

It seems that `Action` is not properly differentiated at some point. The AD algorithm cannot handle the situation where the right part has to be differentiated as well as the left one.

An obvious workaround is to manually define a new form `F_impl`, where we explicitly use the right part of `F_expanded`, which is the derivative of `N`. Basically, we replace the action `F_expanded` with a normal form:
```python
# Replace Action with a form explicitly
dN = F_expanded.right() # dN(u; v1, v')
F_expl = dN * u_test * dx
J_expl = derivative(F_expl, u, u_trial)
J_expl_expanded = expand_derivatives(J_expl) # Action( v1 * v2 * dx; ddN(u; v2, v1, v')) <--- Works correctly
```

How should one handle the issue? Introduce a replacement mechanism for `Action` before differentiation, or is it possible to augment the current algorithm recursively?

Full MWE:
```python
from ufl.core.external_operator import ExternalOperator
from ufl import Coefficient, Argument, derivative, dx
from ufl.algorithms import expand_derivatives
from mpi4py import MPI
from dolfinx.mesh import create_rectangle, CellType
from dolfinx.fem import functionspace

domain = create_rectangle(
MPI.COMM_WORLD,
[[0.0, 0.0], [1.0, 1.0]],
[10, 10],
cell_type=CellType.triangle,
)
V = functionspace(domain, ("Lagrange", 1))

V_ufl = V.ufl_function_space()
u = Coefficient(V_ufl)
u_test = Argument(V_ufl, 1)
u_trial = Argument(V_ufl, 2)
N = ExternalOperator(u, function_space=V_ufl) # N(u; v')
E = N * dx # N(u; v') * dx
F = derivative(E, u, u_test)
F_expanded = expand_derivatives(F) # Action( v1 * dx; dN(u; v1, v'))

# Replace Action with a form explicitly
dN = F_expanded.right() # dN(u; v1, v')
F_expl = dN * u_test * dx
J_expl = derivative(F_expl, u, u_trial)
J_expl_expanded = expand_derivatives(J_expl) # Action( v1 * v2 * dx; ddN(u; v2, v1, v')) <--- Works correctly

# Double differentiation
J = derivative(F, u, u_trial)
J_expanded = expand_derivatives(J) # == 0 <--- Works, but incorrectly

# Differentiation expansion -> Second derivative -> Second expansion
J_expanded = derivative(F_expanded, u, u_trial)
J_expanded_expanded = expand_derivatives(J_expanded) # Error
```

Error:
```python
---> [36](vscode-notebook-cell:?execution_count=18&line=36) J_expanded = derivative(F_expanded, u, u_trial)
37 J_expanded_expanded = expand_derivatives(J_expanded) # Error

File /dolfinx-env/lib/python3.12/site-packages/ufl/formoperators.py:411, in derivative(form, coefficient, argument, coefficient_derivatives)
409 dright = derivative(right, coefficient, argument, coefficient_derivatives)
410 # Leibniz formula
--> [411](https://vscode-remote+attached-002dcontainer-002b7b22636f6e7461696e65724e616d65223a2266656e696373782d302e31302e30227d.vscode-resource.vscode-cdn.net/dolfinx-env/lib/python3.12/site-packages/ufl/formoperators.py:411) return action(
412 adjoint(dleft, derivatives_expanded=True), right, derivatives_expanded=True
413 ) + action(left, dright, derivatives_expanded=True)
414 else:
415 raise NotImplementedError(
416 "Action derivative not supported when the left argument is not a 1-form."
417 )

File /dolfinx-env/lib/python3.12/site-packages/ufl/formoperators.py:205, in action(form, coefficient, derivatives_expanded)
203 if isinstance(form, Form):
204 return compute_form_action(form, coefficient)
--> [205](https://vscode-remote+attached-002dcontainer-002b7b22636f6e7461696e65724e616d65223a2266656e696373782d302e31302e30227d.vscode-resource.vscode-cdn.net/dolfinx-env/lib/python3.12/site-packages/ufl/formoperators.py:205) return Action(form, coefficient)

File /dolfinx-env/lib/python3.12/site-packages/ufl/action.py:91, in Action.__new__(cls, *args, **kw)
...
--> [217](https://vscode-remote+attached-002dcontainer-002b7b22636f6e7461696e65724e616d65223a2266656e696373782d302e31302e30227d.vscode-resource.vscode-cdn.net/dolfinx-env/lib/python3.12/site-packages/ufl/action.py:217) V_left = left.arguments()[-1].ufl_function_space().dual()
218 else:
219 raise TypeError("Action left argument must be either Coefficient or BaseForm")

IndexError: tuple index out of range
```