  An amplitude analysis of $B^0 \to J/\psi \phi K_S^0$ decays is performed
using proton-proton collision data corresponding to an integrated luminosity of
9 ${\rm fb}^{-1}$, collected with the LHCb detector at centre-of-mass energies
of 7, 8 and 13 TeV. Evidence with a significance of $4.0$ standard deviations
of a structure in the $J/\psi K_S^0$ system, named $T_{\psi
s1}^{\theta}(4000)^0$, is seen, with its mass and width measured to be
$3991~_{-10}^{+12}~^{+\phantom{0}9}_{-17}$ MeV and
$105~_{-25}^{+29}~^{+17}_{-23}$ MeV, respectively, where the first uncertainty
is statistical and the second systematic. The $T_{\psi s1}^{\theta}(4000)^0$
state is likely to be the isospin partner of the $T_{\psi s1}^{\theta}(4000)^+$
state, previously observed in the $J/\psi K^+$ system of the $B^+ \to J/\psi
\phi K^+$ decay. When isospin symmetry for the charged and neutral $T_{\psi
s1}^{\theta}(4000)$ states, is assumed, the signal significance increases to
$5.4$ standard deviations.
