Is it true that $\prod_{i=1}^n\left(3 + \frac{1}{3i}\right) < \left(3 + \frac{1}{n}\right)^n$?

If you know the gamma function$$P_n=\prod\limits_{i=1}^n\left(3 + \dfrac{1}{3i}\right)=\frac{3^n\, \Gamma \left(n+\frac{10}{9}\right)}{\Gamma \left(\frac{10}{9}\right) \Gamma (n+1)}$$ Taking its logarithm, using Stirling approximation and continuing with Taylor expansion, we have $$\log(P_n)=n \log(3)+\frac 19\log(n)-\log \left(\Gamma \left(\frac{10}{9}\right)\right)+\frac 5{81n}+O\left(\frac 1{n^2}\right)$$

Doing the same for

$$Q_n=\left(3 + \dfrac{1}{n}\right)^n$$

$$\log(Q_n)=n \log (3)+\frac{1}{3}-\frac{1}{18 n}+O\left(\frac{1}{n^2}\right)$$

$$\log(Q_n)-\log(P_n)=-\frac{\log (n)}{9}+\frac{1}{3}+\log \left(\Gamma \left(\frac{10}{9}\right)\right)-\frac{19}{162 n}+O\left(\frac{1}{n^2}\right)$$ $$\frac {Q_n}{P_n}=e^{\log(Q_n)-\log(P_n)}=\frac{\sqrt[3]{e} (162 n-19) \Gamma \left(\frac{10}{9}\right)}{162 n^{10/9}}$$ and, according to this approximation $\frac {Q_n}{P_n} < 1$ as soon as $n>9$.

It is quite inaccurate and we need more terms for a better bound; Adding two terms to the expansion shows that it could be true for $n > 11$.

Now, it is time for numerical checks (the numbers are rounded to the closest integer) $$\left( \begin{array}{ccc} n & P_n & Q_n \\ 10 & 81025 & 81963 \\ 11 & 245529 & 246008 \\ 12 & 743407 & 738322 \\ 13 & 2249282 & 2215727 \\ 14 & 6801399 & 6649146 \\ 15 & 20555340 & 19952558 \end{array} \right)$$

Here is an approach which uses a couple of properties of the prime numbers. First, as discussed in the question comments, e.g., with Professor Vector's comment, the right hand side of your inequality is

$$f(n) = \left(3 + \frac{1}{n}\right)^n = 3^n\left(1 + \frac{\left(\frac{1}{3}\right)}{n}\right)^n \tag{1}\label{eq1A}$$

As shown at the end of the Formal definition section of Wikipedia's "Exponential function" article,

$$e^x = \lim_{n \to \infty}\left(1 + \frac{x}{n}\right)^n \tag{2}\label{eq2A}$$

In your case, $x = \frac{1}{3}$, so for large $n$,

$$f(n) \sim 3^ne^{\frac{1}{3}} \tag{3}\label{eq3A}$$

In particular, the limit existing means $\left(1 + \frac{1}{3n}\right)^n$ is bounded, so there's a constant $C$ where, for all positive integers $n$,

$$f(n) \lt C(3^n) \tag{4}\label{eq4A}$$

The left side of your inequality is

$$g(n) = \prod_{i=1}^n\left(3 + \frac{1}{3i}\right) = 3^n\left(\prod_{i=1}^n\left(1 + \frac{1}{9i}\right)\right) \tag{5}\label{eq5A}$$

Note

$$(1 + x_1)(1 + x_2) = 1 + x_1 + x_2 + x_1 x_2 \tag{6}\label{eq6A}$$

$$(1 + x_1)(1 + x_2)(1 + x_3) = 1 + x_1 + x_2 + x_3 + x_1 x_2 + x_1 x_3 + x_2 x_3 + x_1 x_2 x_3 \tag{7}\label{eq7A}$$

In general, the product of $n$ factors $1 + x_i$ is the sum of all of the $x_i$ taken $0$ at a time (i.e., just $1$), $1$ at a time, $2$ at a time, etc., up to $n$ at a time. Thus, if all $x_i \ge 0$, they'll all be non-negative terms, so for all integers $n \ge 1$,

$$\prod_{i = 1}^{n}(1 + x_i)^n \gt \sum_{i = 1}^{n}x_i \tag{8}\label{eq8A}$$

Using this with the right side of \eqref{eq5A} gives

$$\prod_{i=1}^n\left(1 + \frac{1}{9i}\right) \gt \sum_{i=1}^n\left(\frac{1}{9i}\right) \tag{9}\label{eq9A}$$

The Prime number theorem shows for large integers $m$,

$$p_m \sim m\ln(m) \tag{10}\label{eq10A}$$

This means there's an integer $m_0$ where, for all $m \ge m_0$,

$$p_m \gt 9m \tag{11}\label{eq11A}$$

Thus, for any $n \ge m_0$,

$$\sum_{i=m_0}^n\left(\frac{1}{9i}\right) \gt \sum_{i=m_0}^n\left(\frac{1}{p_i}\right) \tag{12}\label{eq12A}$$

However, such as discussed in Does the sum of reciprocals of primes converge?, and shown in Wikipedia's Divergence of the sum of the reciprocals of the primes article, the right side of \eqref{eq12A} diverges. Thus, the left side must also diverge, so it'll eventually become larger than the $C$ in \eqref{eq4A}. This means there's an $n_0$ where, for all $n \ge n_0$,

$$g(n) \gt f(n) \tag{13}\label{eq13A}$$

i.e., your inequality will not be true. Note this also shows that if you change the $\frac{1}{3i}$ in $g(n)$ defined in \eqref{eq5A} to instead be $\frac{1}{ki}$ for some positive constant $k$, then no matter how large $k$ is, your inequality will still always eventually no longer be true.

Update: The right side of \eqref{eq9A} is $\frac{1}{9}$ of the Harmonic series which is known to diverge. Using this is simpler and more direct than utilizing the properties of the prime numbers I used above.