You can use Rfast::Dist() instead of sf::st_distance() to fast calculate the Euclidean distance matrix between points. Alternatively, you can use parallelDist::parDist(), which is multithreaded, but returns a triangle instead of a matrix.

Code
library("sf")
library("Rfast")
set.seed(1)
bbox = st_bbox(c(xmin = 0, xmax = 1, ymax = 0, ymin = 1))
bbox = st_as_sfc(bbox)
n = 10000
pts = st_sample(bbox, n)
system.time({
  m1 = st_distance(pts, which = "Euclidean")
})
##    user  system elapsed 
##    1.18    0.19    1.36
system.time({
  df = st_coordinates(pts)
  m2 = Rfast::Dist(df, method = "euclidean")
})
##    user  system elapsed 
##    0.55    0.18    0.72
all(m1 == m2)
## [1] TRUE
m1[1:5, 1:5]
##            1         2         3         4          5
## 1 0.00000000 0.6211186 0.7377293 0.6443854 0.06633203
## 2 0.62111859 0.0000000 0.2091530 0.7790792 0.65260806
## 3 0.73772930 0.2091530 0.0000000 0.7084691 0.78236784
## 4 0.64438539 0.7790792 0.7084691 0.0000000 0.70947711
## 5 0.06633203 0.6526081 0.7823678 0.7094771 0.00000000