The two ways Feynman diagrams appear in mathematics

If I understand the question correctly, the search is for a calculation of the asymptotic expansion of Gaussian integrals using concepts and techniques from category theory. Here is one such calculation:

Feynman diagrams via graphical calculus (2001)

There is a very close connection between the graphical formalism for ribbon categories and Feynman diagrams. Although this correspondence is frequently implied, we know of no systematic exposition in existing literature; the aim of this paper is to provide such an account. In particular, in deriving Feynman diagrams expansion of Gaussian integrals as an application of the graphical formalism for symmetric monoidal categories, we discuss in detail how different kinds of interactions give rise to different families of graphs and show how symmetric and cyclic interactions lead to “ordinary” and “ribbon” graphs respectively.

Before interpreting them in more advanced language like "string diagrams" or "monoidal closed categories" it might be good to stress that Feynman diagrams are very elementary combinatorial objects which encode contractions of tensors. By tensor I mean an array of numbers like $A=(A_{a,b,c})$ with indices $a,b,c$ running over some finite sets which need to be specified. If you have such objects say $B_{abcd}$ and $C_{ab}$ you can construct new ones like $$ T_{a,b,c,d}:=\sum_{e,f,\ldots,l} C_{ae}C_{bg} B_{eghf} C_{fi}C_{hk}B_{iklj}C_{jc}C_{ld} $$ Obviously one can produce tons of similar examples of increasing complexity and it is desirable to have a way of encoding precisely such constructions. A natural way of doing that is basically to use pictures or graphs which is what Feynman diagrams are.

Linear algebra "done wrong", i.e., matrix algebra is the particular case where tensors have one (vectors) or two indices (matrices) only. Although, the $n$-dimensional determinant introduces a higher (Levi-Civita) tensor $\epsilon_{i_1\ldots i_n}$ given by the sign of the permutation $i_1\ldots i_n$ (and zero if indices are repeated).

A significant part of functional analysis is about studying what happens when discrete summation indices $a,b,\ldots$ become continuous variables that are integrated over. Then, matrices $C_{a,b}$ become kernels $C(x,y)$ which one can make sense of, e.g, as distributions, by invoking the Schwartz Kernel Theorem.

Feynman diagrams, as they are used in quantum field theory, typically correspond to this infinite-dimensional generalization. For instance, the diagram for the expression $T_{abcd}$ above becomes the main second order contribution to the four-point function of the $\phi^4$ model in $d$ dimensions if one decides that the tensor $C_{ab}$ now becomes the kernel $C(x,y)$ of the operator $-\Delta+m^2 I$ in $\mathbb{R}^d$ and one lets the tensor $B_{abcd}$ become the kernel $B(x,y,z,u)$ of the distribution in $S'(\mathbb{R}^{4d})$ given by the action $$ f\longmapsto\ \int_{\mathbb{R}^d} f(x,x,x,x)\ d^dx $$ on test functions $f\in S(\mathbb{R}^{4d})$.

As for why this should have to do with the Laplace/stationary phase method, the reason is because Gaussian integration is an "algebraic" operation. Namely, it can be expressed as a differential operator (albeit of infinite order). For example if $\mu$ is the centered Gaussian measure on $\mathbb{R}^n$ with covariance $C_{a,b}$, then for any polynomial $P\in \mathbb{R}[x_1,\ldots,x_n]$, $$ \int_{\mathbb{R}^n} P(x)\ d\mu(x)=\left.\exp\left(\frac{1}{2}\sum_{a,b=1}^n C_{a,b} \frac{\partial^2}{\partial x_a\partial x_b}\right)\ P(x)\ \right|_{x=0}\ . $$

Note that Haar integration on $SU(n)$ can also be expressed as an infinite order differential operator (see my two answers to How to constructively/combinatorially prove Schur-Weyl duality? ). So Feynman diagrams also appear in invariant theory/representation theory (see my answer to Who invented diagrammatic algebra? for some examples and pictures by Kempe in the case of the invariants of the binary quintic that are given explicitly in nondiagrammatic fashion in my answer to Explicit formulas for invariants of binary quintic forms ).