L cholesky factor gets NA for fit

I have a STAN model that runs without divergences but that generates NA for the upper triangular matrix entries in the L // Cholesky factor of the correlation matrix (which are supposed to just be 0). This prevents bridgesampling from working for example.

any idea what it could be or what i am doing wrong?

data {

    // Metadata
    int<lower=1>  N;                   // Number of total observations
    int<lower=1>  J[N];                // Subject-indicator per observation
    int<lower=1>  K[N];                // Bandit-indicator per observation
    
    // Data
    int        Y[N];                   // Response (risky = 1, safe = 0)
    real       R[N];                   // Outcome (win = 1, loss = 0)

}
transformed data {

    int  NJ = max(J);                  // Number of total subjects
    int  NK = max(K);                  // Number of total bandits

}
parameters {

    // Subject parameters
    vector[4]     theta_mu;            // Population-level effects
    matrix[4,NJ]  theta_pr;            // Standardized subject-level effects
    
    // Subject variances
    cholesky_factor_corr[4] L;         // Cholesky factor of correlation matrix
    vector<lower=0>[4] sigma;          // Subject-level standard deviations

}
transformed parameters {

    vector[NJ]  b1;                    // Inverse temperature
    vector[NJ]  a1;                    // Learning rate (positive)
    vector[NJ]  a2;                    // Learning rate (negative)
    vector[NJ]  q0;                    // Initial Q-value

    // Construction block
    {
    
    // Rotate random effects
    matrix[NJ,4] theta = transpose(diag_pre_multiply(sigma, L) * theta_pr);
    
    // Construct random effects
    b1 = (theta_mu[1] + theta[,1]) * 8;
    a1 = Phi_approx(theta_mu[2] + theta[,2]);
    a2 = Phi_approx(theta_mu[3] + theta[,3]);
    q0 = Phi_approx(theta_mu[4] + theta[,4]);
    
    }

}
model {
   real delta;
   real eta;

    // Initialize Q-values
    real  Q[NJ,NK] = to_array_2d(rep_matrix(q0, NK));

    // Construct linear predictor
    vector[N] mu;
    for (n in 1:N) {
    
        // Compute (weighted) difference in expected values
        mu[n] = b1[J[n]] * (Q[J[n],K[n]] - 0.5);
        
        // Compute prediction error
         delta = R[n] - Q[J[n],K[n]];
        
        // Assign trial-level learning rate
         eta = R[n] * a1[J[n]] + (1-R[n]) * a2[J[n]];
        
        // Update state-action values
        Q[J[n],K[n]] += eta * delta;
        
    }
    
    // Likelihood
    target += bernoulli_logit_lpmf(Y | mu); 
    
    // Priors
    target += normal_lpdf(theta_mu | 0, 2);
    target += std_normal_lpdf(to_vector(theta_pr));
    target += student_t_lpdf(sigma | 3, 0, 1);
    target += lkj_corr_cholesky_lpdf(L | 1);

}
generated quantities {

    matrix[4,4] Omega = multiply_lower_tri_self_transpose(L);
    vector[NJ]  ratio = (a1 - a2) ./ (a1 + a2);
    
}