How to fix Mathematica crashing in the following case?

Question

I have a problem with Mathematica when trying to glue two tables together.

The tables have the form

tab2={{d1,f1,g1,...},{d2,f2,g2,...},...}
tab1 = {{a1,b1,c1},{a1,b1,c2},...,{a1,b2,c1},{a1,b2,c2},...{a2,b1,c1},{c2,b1,c1},...}

The resulting table would be

tab3 = {{a1,b1,c1,d1,f1,g1,...},{a1,b1,c1,d2,f2,g2,...},...}

My problem is that for the case of specific tab1, tab2 Mathematica crashes when trying to compile tab3, Most likely this is due to the total length of tab3 (it just quits kernel without any messages). In particular, it works if reducing the dimension of tab1 by, for instance, 50 times. My question is how to fix this problem. My machine: Mathematica 12.1@Windows10, Acer ConceptD7 Ezel, 16 Gb RAM, i7-10875H.

Please find the code defining an example tab2 (in the code, it is called TableDistributionNtoYZtempReduced) below. Note that it is complicated, as in reality tab2 is computed much less trivially. The dimensionality is controlled by the parameter isimulval (set to 10000):

MassParticleTest = 1.5;
mProdTomXratio2body = 0.3;
DecayType = "2-body";
mProd3Body1 = 0;
mProd3Body2 = 0;
mProd3Body3 = 0;
mproduct1 = 
  If[DecayType == "2-body", mProdTomXratio2body*MassParticleTest, 
   If[DecayType == "3-body", mProd3Body1]];
mproduct2 = 
  If[DecayType == "2-body", mProdTomXratio2body*MassParticleTest, 
   If[DecayType == "3-body", mProd3Body2]];
(*Number of rows in the table*)
isimulval = 10000;
(*Position of the point in coordinate space in terms of spherical \
coordinates*)
phVal[px_, py_] = 
  If[py > 0, 
   ArcCos[px/Sqrt[px^2 + py^2]], -ArcCos[px/Sqrt[px^2 + py^2]]];
thVal[px_, py_, pz_] = ArcCos[pz/Sqrt[px^2 + py^2 + pz^2]];
TableDistributionNtoYZtempReducedTemp = 
  Hold@Compile[{{tab, _Real, 2}, imax}, 
       Table[{thVal[tab[[i]][[1]], tab[[i]][[2]], tab[[i]][[3]]], 
         phVal[tab[[i]][[1]], tab[[i]][[2]]], tab[[i]][[4]], 
         thVal[tab[[i]][[5]], tab[[i]][[6]], tab[[i]][[7]]], 
         phVal[tab[[i]][[5]], tab[[i]][[6]]], tab[[i]][[8]]}, {i, 1, 
         imax, 1}]] /. DownValues@phVal /. DownValues@thVal // 
   ReleaseHold;
imax = 10000;
Block2BodyDistrAtRest := Block[{},
  nVectorParticle1 = Table[RandomPoint[Sphere[]], {i, 1, imax, 1}];
  nVectorParticles := 
   Hold@Compile[{{nVectorParticle1, _Real, 2}, mN, mX, mY, imax}, 
      Table[{Sqrt[((mN^2 + mX^2 - mY^2)/(2*mN))^2 - mX^2]
          nVectorParticle1[[i]][[1]], 
        Sqrt[((mN^2 + mX^2 - mY^2)/(2*mN))^2 - mX^2]
          nVectorParticle1[[i]][[2]], 
        Sqrt[((mN^2 + mX^2 - mY^2)/(2*mN))^2 - mX^2]
          nVectorParticle1[[i]][[3]], ((mN^2 + mX^2 - mY^2)/(
         2*mN)), -Sqrt[((mN^2 - mX^2 + mY^2)/(2*mN))^2 - mY^2]
           nVectorParticle1[[i]][[1]], -Sqrt[((mN^2 - mX^2 + mY^2)/(
            2*mN))^2 - mY^2]
           nVectorParticle1[[i]][[2]], -Sqrt[((mN^2 - mX^2 + mY^2)/(
            2*mN))^2 - mY^2] nVectorParticle1[[i]][[3]], ((
         mN^2 - mX^2 + mY^2)/(2*mN))}, {i, 1, imax, 1}]] // 
    ReleaseHold;
  nVectorParticlesmXmY = 
   nVectorParticles[nVectorParticle1, MassParticleTest, mproduct1, 
    mproduct2, imax]]
TableDistributionNtoYZtempReduced = 
  If[DecayType == "2-body", 
   TableDistributionNtoYZtempReducedTemp[Block2BodyDistrAtRest, 
    isimulval], 
   If[DecayType == "3-body", 
    TableDistributionNtoYZtempReducedTemp[Block3BodyDistrAtRest[imax],
      isimulval], 0]];

The same is true about tab1 (called TableHNLgridWithEN). Its dimensionality is fixed and equal to 229190:

AzimuthalAcceptanceMotherParticle[thN_, xLongN_] = 0.16;
thGivenExperimentDetMin = 0.44;
thGivenExperimentDetMax = 0.88;
xLongToDecVolGivenExperiment = 60;
xLongNMaxGivenExperiment = 85;
NumberphN = 10;
RandomphNvalues[thN_, xLongN_] := RandomReal[{-Pi, Pi}, NumberphN]
Needs["CCompilerDriver`"]
CCompilers[]
DefaultCCompiler[]
TableHNLgridTemp1 = 
   Flatten[ParallelTable[{thN, xLongN, 
      AzimuthalAcceptanceMotherParticle[thN, xLongN], 
      RandomphNvalues[thN, xLongN]}, {thN, 
      1.001 thGivenExperimentDetMin + 10^-5., 
      0.99 thGivenExperimentDetMax, (0.99 thGivenExperimentDetMax - \
(1.001 thGivenExperimentDetMin + 10^-5.))/12.}, {xLongN, 
      1.005 xLongToDecVolGivenExperiment, 
      xLongNMaxGivenExperiment, (xLongNMaxGivenExperiment - 
         1.005 xLongToDecVolGivenExperiment)/40.}], {1, 
     2}]; // AbsoluteTiming
TableHNLgridTemp12 = TableHNLgridTemp1[[All, {1, 2, 3}]];
TableHNLgridTemp13 = TableHNLgridTemp1[[All, 4]];
ENvalues = Table[10^EN, {EN, Log10[1.1*MassParticleTest], 3.2, 0.07}];
TableHNLgridWithEN = 
   Flatten[Table[{ENvalues[[k]], TableHNLgridTemp12[[i]][[1]], 
      TableHNLgridTemp12[[i]][[2]], TableHNLgridTemp12[[i]][[3]], 
      TableHNLgridTemp13[[i]][[j]]}, {k, 1, Length[ENvalues], 1}, {i, 
      1, Length[TableHNLgridTemp12], 1}, {j, 1, NumberphN, 1}], {1, 2,
      3}]; // AbsoluteTiming

The code which computes tab3 (which has the dimensionality isimulval*229190 in the example reproducing the problem) and fails:

cf = Compile[{{x, _Real, 1}, {y, _Real, 2}}, 
   Table[Join[x, Compile`GetElement[y, i]], {i, 1, Length[y]}], 
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
   Parallelization -> True, RuntimeOptions -> "Speed"];
TableHNLgridWithENfinal = 
   Flatten[cf[TableHNLgridWithEN, 
     TableDistributionNtoYZtempReduced], {1, 2}]; // AbsoluteTiming