I would like to simulate a response using posterior_predict() with a model where sample_prior='only' was specified. As written in the docs, this gives me a matrix with n_samples rows and N columns.

Is it possible to somehow also know which exact samples from the prior were used in order to simulate each row of size N?

How does posterior_predict pick a concrete sample from the prior distributions, randomly or just picking the already sampled values from the prior one by one?

sample_prior='only' does not store the priors, is this intended?

The sample_prior='only' actually does not involve direct sampling from the prior - it just runs the model through Stan while adding only the contributions from the prior to the target density and not contributions from the likelihood. So the prior is sampled via Stan’s MCMC and there is no immediate way to trace the samples back to univariate draws from the relevant prior distributions. So I don’t think you can get the prior distributions in the way you want.

However, it is usually not difficult to write an R script to simulate from the prior directly - and I’ve often find writing such a simulator very helpful in understanding what the model is actually doing.