I just extended @Michael-E2 answer by including FunctionCompile:

lx = Range[40000, 59999];
bx = RandomSample[Range[20000], 10000];
ares = MapAt[#^2 &, lx, {#} & /@ bx]; // AbsoluteTiming
bres = SubsetMap[#^2 &, lx, {#} & /@ bx]; // AbsoluteTiming
(cres = lx; Do[cres[[k]] = cres[[k]]^2, {k, bx}]); // AbsoluteTiming
dres = Compile[{{lx, _Integer, 1}, {bx, _Integer, 1}}, 
     Block[{dres = lx}, Do[dres[[k]] = dres[[k]]^2, {k, bx}];
      dres]][lx, bx]; // AbsoluteTiming
(eres = lx; eres[[bx]] = eres[[bx]]^2); // AbsoluteTiming
ffun = FunctionCompile[
   Function[
    {Typed[lx, TypeSpecifier["PackedArray"]["MachineInteger", 1]],
     Typed[bx, TypeSpecifier["PackedArray"]["MachineInteger", 1]]},
    Block[{dres = lx}, Do[dres[[k]] = dres[[k]]^2, {k, bx}]; dres]]];
fres = ffun[lx, bx]; // AbsoluteTiming

ares == bres == cres == dres == eres == fres

(*
  {1.36477,Null}
  {0.0527128,Null}
  {0.0252494,Null}
  {0.0004774,Null}
  {0.0001841,Null}
  {0.0000768,Null}
  True
*)

I just extended @Michael-E2 answer by including FunctionCompile:

lx = Range[40000, 59999];
bx = RandomSample[Range[20000], 10000];
ares = MapAt[#^2 &, lx, {#} & /@ bx]; // AbsoluteTiming
bres = SubsetMap[#^2 &, lx, {#} & /@ bx]; // AbsoluteTiming
(cres = lx; Do[cres[[k]] = cres[[k]]^2, {k, bx}]); // AbsoluteTiming
dres = Compile[{{lx, _Integer, 1}, {bx, _Integer, 1}}, 
     Block[{dres = lx}, Do[dres[[k]] = dres[[k]]^2, {k, bx}];
      dres]][lx, bx]; // AbsoluteTiming
(eres = lx; eres[[bx]] = eres[[bx]]^2); // AbsoluteTiming
ffun = FunctionCompile[
   Function[
    {Typed[lx, TypeSpecifier["PackedArray"]["MachineInteger", 1]], 
     Typed[bx, TypeSpecifier["PackedArray"]["MachineInteger", 1]]}, 
    Block[{fres = lx}, Do[fres[[k]] = fres[[k]]^2, {k, bx}]; fres]]];
fres = ffun[lx, bx]; // AbsoluteTiming

ares == bres == cres == dres == eres == fres

(*
  {1.29098,Null}
  {0.0542107,Null}
  {0.0229712,Null}
  {0.0005062,Null}
  {0.0001768,Null}
  {0.0000713,Null}
  True
*)