The polarization of $\Lambda$ and $\bar{\Lambda}$ hyperons along the beam direction has been measured relative to the second and third harmonic event planes in isobar Ru+Ru and Zr+Zr collisions at $\sqrt{s_{NN}}$ = 200 GeV. This is the first experimental evidence of the hyperon polarization by the triangular flow originating from the initial density fluctuations. The amplitudes of the sine modulation for the second and third harmonic results are comparable in magnitude, increase from central to peripheral collisions, and show a mild $p_T$ dependence. The azimuthal angle dependence of the polarization follows the vorticity pattern expected due to elliptic and triangular anisotropic flow, and qualitatively disagree with most hydrodynamic model calculations based on thermal vorticity and shear induced contributions. The model results based on one of existing implementations of the shear contribution lead to a correct azimuthal angle dependence, but predict centrality and $p_T$ dependence that still disagree with experimental measurements. Thus, our results provide stringent constraints on the thermal vorticity and shear-induced contributions to hyperon polarization. Comparison to previous measurements at RHIC and the LHC for the second-order harmonic results shows little dependence on the collision system size and collision energy.
$sgn(\alpha_H)\langle\cos(\theta_p^{\ast})\rangle$ of $\Lambda$ and $\bar{\Lambda}$ as a function of hyperon azimuthal angle relative to the second-order event plane in isobar collisions at 200 GeV.
$sgn(\alpha_H)\langle\cos(\theta_p^{\ast})\rangle$ of $\Lambda$ and $\bar{\Lambda}$ as a function of hyperon azimuthal angle relative to the third-order event plane in isobar collisions at 200 GeV.
$P_z$ sine coefficients of $\Lambda+\bar{\Lambda}$ as a function of centrality in isobar collisions at 200 GeV.