This is related to the Coupon collector's problem as noted in the comments.

Building off of [this][1] post, the probability of observing $k$ unique letters in $m$ random uniform samples from an alphabet of size $n$ is:

$$\bigg\{\!{m\!\atop{k}}\bigg\}\binom{n}{k}\frac{k!}{n^m}=\bigg\{\!{m\!\atop{k}}\bigg\}\frac{n!}{n^m(n-k)!}$$

Where $\big\{\!{m\!\atop{k}}\big\}$ is the Stirling number of the second kind.

For large $m$, $\ln\!\big(\big\{\!{m\!\atop{k}}\big\}\big)$ can be [approximated][2].

Here is an R function that returns the probability of every $k$:

    library(copula)
    
    coupons1 <- function(n, m) {
      l <- min(m, n)
      k <- 1:l
      if (m < 200) {
        logS <- log(Stirling2.all(m)[k])
      } else {
        # estimate the log Stirling numbers
        v <- m/k
        G <- 1/v
        vexpv <- v/exp(v)
        for (i in 1:5) G <- G - (G - (vexpG <- vexpv*exp(G)))/(1 - vexpG) # Newton's method
        logS <- (log(v - 1) - log(v*(1 - G)))/2 + (m - k)*(log(v - 1) - log(v - G)) + m*log(k) - k*log(m) + k*(1 - G) + lgamma(m + 1) - lgamma(k + 1) - lgamma(m - k + 1)
        if (l == m) logS[m] <- 0
      }
      exp(logS + lgamma(n + 1) - lgamma(n - k + 1) - m*log(n))
    }

The probability of $k=1,2...43000$ for $n=m=43000$:

    system.time(k1 <- coupons1(43e3, 43e3))
    #>    user  system elapsed 
    #>    0.02    0.00    0.01

    plot(26900:27500, k1[26900:27500], xlab = "k", ylab = "p(k)", col = "blue")

[![enter image description here][3]][3]

Comparing that result to a brute-force approach:

    Rcpp::cppFunction("
      NumericVector coupons2(const int& n, const int& m) {
        int maxk;
        int n1 = n - 1;
        if (n > m) {
          maxk = m;
        } else {
          maxk = n;
        }
        
        NumericVector k (maxk);
        k(0) = 1;
        
        for (int i = 1; i < maxk; i++) {
          for (int j = i; j > 0; j--) {
            k(j) = (k(j - 1)*(n - j) + k(j)*(j + 1))/n;
          }
          k(0) = k(0)/n;
        }
        
        for (int i = maxk; i < m; i++) {
          for (int j = n1; j > 0; j--) {
            k(j) = (k(j - 1)*(n - j) + k(j)*(j + 1))/n;
          }
          k(0) = k(0)/n;
        }
        
        return k;
      }
    ")
    
    system.time(k2 <- coupons2(43e3, 43e3))
    #>    user  system elapsed 
    #>   12.29    0.00   12.31

The relative error from using the Stirling number approximation is small for a large $m$.

    max(abs(k1[26900:27500] - k2[26900:27500])/k2[26900:27500])
    #> [1] 8.288009e-07

    points(26900:27500, k2[26900:27500], col = "orange", pch = 20)
    legend("topright", legend = c("k1", "k2"), col = c("blue", "orange"), pch = c(1, 20))

[![enter image description here][4]][4]

A note on the `occupancy` R package:
------------------------------------

The `occupancy` package features the distribution for this problem, with an option for exact or approximate calculations. `docc` returns the PMF of `k`. However, it is much slower than `coupons1` and `coupons2` above, with the computation time growing approximately geometrically for $m=n$. Additionally, the approximation provided is poor compared to the one provided by `coupons1`:

    n <- 43e3L
    system.time(k3 <- occupancy::docc(1:n, n, n, approx = TRUE))
    #>    user  system elapsed 
    #>    1.61    2.00    3.61
    max(abs(k3[26900:27500] - k2[26900:27500])/k2[26900:27500])
    #> [1] 0.01156602
    
Time `occupancy::docc` for $m=n=4300$ with `approx = FALSE`:
    
    n <- 4300L
    system.time(occupancy::docc(1:n, n, n))
    #>    user  system elapsed 
    #>   35.58    0.14   35.75

  [1]: https://math.stackexchange.com/a/693254/697491
  [2]: https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind#Asymptotic_approximation
  [3]: https://i.sstatic.net/gZYAZ.png
  [4]: https://i.sstatic.net/1wGdh.png