This paper presents a search for pair production of higgsinos, the supersymmetric partners of the Higgs bosons, in scenarios with gauge-mediated supersymmetry breaking. Each higgsino is assumed to decay into a Higgs boson and a nearly massless gravitino. The search targets events where each Higgs boson decays into <a:math xmlns:a="http://www.w3.org/1998/Math/MathML" display="inline"><a:mi>b</a:mi><a:mover accent="true"><a:mi>b</a:mi><a:mo stretchy="false">¯</a:mo></a:mover></a:math>, leading to a reconstructed final state with at least three energetic <e:math xmlns:e="http://www.w3.org/1998/Math/MathML" display="inline"><e:mi>b</e:mi></e:math>-jets and missing transverse momentum. Two complementary analysis channels are used, with each channel specifically targeting either low or high values of the higgsino mass. The low-mass (high-mass) channel exploits <g:math xmlns:g="http://www.w3.org/1998/Math/MathML" display="inline"><g:mrow><g:mn>126</g:mn><g:mtext> </g:mtext><g:mo stretchy="false">(</g:mo><g:mn>139</g:mn><g:mo stretchy="false">)</g:mo><g:mtext> </g:mtext><g:mtext> </g:mtext><g:msup><g:mrow><g:mi>fb</g:mi></g:mrow><g:mrow><g:mo>−</g:mo><g:mn>1</g:mn></g:mrow></g:msup></g:mrow></g:math> of <k:math xmlns:k="http://www.w3.org/1998/Math/MathML" display="inline"><k:msqrt><k:mi>s</k:mi></k:msqrt><k:mo>=</k:mo><k:mn>13</k:mn><k:mtext> </k:mtext><k:mtext> </k:mtext><k:mi>TeV</k:mi></k:math> data collected by the ATLAS detector during Run 2 of the Large Hadron Collider. No significant excess above the Standard Model prediction is found. At 95% confidence level, masses between 130 GeV and 940 GeV are excluded for higgsinos decaying exclusively into Higgs bosons and gravitinos. Exclusion limits as a function of the higgsino decay branching ratio to a Higgs boson are also reported. © 2024 CERN, for the ATLAS Collaboration 2024 CERN