In the multi-urn Ehrenfest model $N$ balls are distributed among $d + 1 (d \geqq 2)$ urns. If we label the urns $0, 1, \cdots, d,$ then the system is said to be in state $\mathbf{i} = (i_1, i_2, \cdots, i_d)$ when there are $i_j$ balls in urn $j (j = 1, 2, \cdots, d)$ and $N - \mathbf{1}\cdot\mathbf{i}$ balls in urn 0. (The vector $\mathbf{1}$ has all its components equal to 1 and $\mathbf{x} \cdot \mathbf{y}$ is the usual scalar product.) At discrete epochs a ball is chosen at random from one of the $d + 1$ urns; each of the $N$ balls has probability $1/N$ of being selected. The ball chosen is removed from its urn and placed in urn $i (i = 0, 1, \cdots, d)$ with probability $p^i,$ where the $p^i$'s are elements of a given vector, $(p^0, \mathbf{p})$, satisfying $p^i > 0$ and $\sum^d_{i=0} p^i = 1.$ We shall let $\mathbf{X}_N(k)$ denote the state of the system after the $k$th such rearrangement of balls. Our interest in this paper is to obtain limit theorems for the sequence of processes $\{\mathbf{X}_N(k): k = 0, \cdots, N\}$ as $N$ tends to infinity. For the classical Ehrenfest model $(d = 1, p^0 = p^1 = \frac{1}{2})$ Kac [7] showed that the distribution of $(X_N(\lbrack Nt\rbrack) - N/2)/(N/2)^{\frac{1}{2}}$ converges as $N \rightarrow \infty$ to the distribution of the Ornstein-Uhlenbeck process at time $t$ having started at $y_0$ at $t = 0,$ provided $X_N(0) = \lbrack(N/2)^{\frac{1}{2}}y_0 + N/2\rbrack.$ (The symbol $\lbrack x\rbrack$ denotes the integer part of $x$.) Recently, Karlin and McGregor [8] obtained a similar result for the continuous time version of the model with $d = 2$; in this version the random selection of balls is done at the occurrence of events of an independent Poisson process. In addition, they obtained a local limit theorem for the transition function. The proof in [7] depended on the continuity theorem of characteristic functions. On the other hand, the proof in [8] used the properties of the spectral representation of the $n$-step transition probabilities which is available for these processes. These results suggested the direction we shall follow in this paper. A preliminary calculation indicates that the process $\{\mathbf{X}_N(k): k = 0, \cdots, N\}$ is attracted to the pseudo-equilibrium state $N\mathbf{p}$ and that states far from $B\mathbf{p}$ will only occur rarely. Thus it is natural to consider the fluctuations of $\mathbf{X}_N(k)$ about $N\mathbf{p}$ measured in an appropriate scale. For our purposes the appropriate processes to consider are $\{\mathbf{Y}_N(k): k = 0, \cdots, N\}$, where $\mathbf{Y}_N(k) = (\mathbf{X}_N(k) - N\mathbf{P})/N^{\frac{1}{2}}.$ Next we define a sequence of stochastic processes $\{\mathbf{y}_N(t): 0 \leqq t \leqq 1\}$ which are continuous, linear on the intervals $((k - 1)N^{-1}, kN^{-1})$, and satisfy $\mathbf{y}_N(kN^{-1}) = \mathbf{Y}_N(k)$ for $k = 0, 1, \cdots, N.$ In other words we let $\mathbf{y}_N(t) = \mathbf{Y}_N(k) + (Nt - k)(\mathbf{Y}_N(k + 1) - \mathbf{Y}_N(k))$ if $kN^{-1} \leqq t \leqq (k + 1)N^{-1}$. Throughout this paper we shall let $X_N^i(0) = \lbrack N^{\frac{1}{2}}y_0^i + Np^i\rbrack,$ where $y_0 = (y_0^1, \cdots, y_0^d)$ is an arbitrary, but fixed, element of $R^d.$ (It will always be understood that $N$ is sufficiently large so that $0 \leqq X_N^i(0) \leqq N$ for all $i = 1, 2, \cdots, d$, where $X_N^i(\cdot)$ is the $i$th component of the vector $\mathbf{X}_N(\cdot). R^d$ is $d$-dimensional Euclidean space.) Observe that this initial condition implies that $|Y_N^i (0) - y_0^i| \leqq N^{-\frac{1}{2}}.$ With this initial condition and the Markov structure of the model, the processes $\{\mathbf{X}_N(k): k = 0, \cdots, N\}$ for $N = 1, 2, \cdots$ can be defined on a probability triple $(\Omega_N, \mathscr{F}_N, P_N).$ We shall let $C_d\lbrack 0, 1\rbrack$ denote the product space of $d$ copies of $C\lbrack 0, 1\rbrack$, the space of continuous functions on $\lbrack 0, 1\rbrack$ with the topology of uniform convergence, and endow $C_d\lbrack 0, 1\rbrack$ with the product topology. The topological Borel field of $C_d\lbrack 0, 1\rbrack$ will be denoted by $\mathscr{C}_d$. Clearly, the transformation taking the sequence $\{\mathbf{X}_N(k): k = 0, \cdots, N\}$ into $\{\mathbf{y}_N(t): 0 \leqq t \leqq 1\}$ is measurable and induces a probability measure on $\mathscr{C}_d$. We shall denote this induced measure by $\mu_N(\cdot; \mathbf{y}_0)$. The general notion of weak convergence of a sequence of probability measures is defined as follows. Let $S$ be a metric space and $\mathscr{S}$ be the Borel field generated by the open sets of $S$. If $\nu_N$ and $\nu$ are probability measures on $\mathscr{S}$ and if the $\lim_{N\rightarrow\infty} \int _sf d\nu_N = \int _sf d\nu$ for every bounded, continuous function $f$ on $S$, then we say that $\nu_N$ converges weakly to $\nu$ and write $\nu_N\Longrightarrow \nu$. The principal result of this paper is that $\mu_N(\cdot; \mathbf{y}_0) \Longrightarrow \mu (\cdot; \mathbf{y}_0)$ as $N \rightarrow \infty$, where $\mu(\cdot; y_0)$ is the probability measure on $\mathscr{C}_d$ of a $d$-dimensional diffusion process, $\mathbf{y}(\cdot)$, starting at the point $\mathbf{y}_0$. The limit process $\mathbf{y}(\cdot)$ is a $d$-dimensional analog of the Ornstein-Uhlenbeck process whose distribution at time $t$ is a multi-variate normal with mean vector $e^{-t}\mathbf{y}_0$ and covariance matrix $\Sigma$, where the elements of $\Sigma$ are \begin{equation*}\begin{split}\sigma_{ij} = (1 - e^{-2t})p^i(1 - p^i),\quad i = j, \\ = -(1 - e^{-2t})p^ip^j,\quad i \neq j.\end{split}\end{equation*} For applications it is useful to note that $\mu_N(\cdot; \mathbf{y}_0)\Longrightarrow\mu(\cdot;\mathbf{y}_0)$ is equivalent to the statement that $\lim_{N\rightarrow\infty} \mu_N(\{f(\mathbf{y}(\cdot)) \leqq \alpha\}; \mathbf{y}_0) = \mu(\{f(\mathbf{y}(\cdot)) \leqq \alpha\}; \mathbf{y}_0)$ for all functionals $f$ on $C_d\lbrack 0, 1\rbrack$ which are continuous almost everywhere with respect to $\mu(\cdot; \mathbf{y}_0)$. To establish weak convergence two steps are usually required. First, the convergence of the finite-dimensional distributions (fdd) of the approximating processes, $\{\mathbf{y}_N(t): t \geqq 0, N = 1, 2, \cdots\}$ in our case, to the corresponding fdd of the limiting process must be obtained. Second, the probability that the approximating processes can have large fluctuations between points at which these processes are determined by their fdd must be shown to be small. The notion of weak convergence is intimately related to the so-called invariance principles. An invariance principle was first introduced for the case of sums of independent, identically distributed random variables by Erdos and Kac [4] and generalized by Donsker [3]. The method of proof used by Erdos-Kac and Donsker was later modified by Billingsley [1] for dependent random variables. In carrying out the second step outlined above we shall follow Billingsley's argument in Theorems 2.3 and 3.1 of [1]. A program similar in nature to ours was carried out by Lamperti [10] for a particular class of Markov processes. This paper is organized into the following sections: Section 2 is devoted to our analog of the central limit theorem (clt), namely, that the distribution of $\mathbf{y}_N(t)$ converges to the distribution of $\mathbf{y}(t)$ at a single fixed value of $t$. In Section 3 the limit process $\mathbf{y}(t)$ is identified and the properties of the process needed here are discussed. Section 4 completes the proof of the convergence of the fdd of $\{\mathbf{y}_N(t)\}$ to those of $\{\mathbf{y}(t)\}$. The proofs in both Sections 2 and 4 are carried out using the Levy continuity theorem for characteristic functions. Section 5 provides the proof required to show weak convergence. The main tool here, in addition to Billingsley's theorems mentioned above, is the result of Stone (1961) on the weak convergence of random walks. Finally, in Section 6 applications are mentioned along with a suggestion as to how the multi-urn Ehrenfest model might be used to study certain problems in statistical mechanics, networks of queues, and epidemic theory.