  Multiplex networks consist of a fixed set of nodes connected by several sets
of edges which are generated separately and correspond to different networks
("layers"). Here, a simple variant of the Ising model on multiplex networks
with two layers is considered, with spins located in the nodes and edges
corresponding to ferromagnetic interactions between them. Critical temperatures
for the ferromagnetic transition are evaluated for the layers in the form of
random Erd\"os-R\'enyi graphs or heterogeneous scale-free networks using the
mean-field approximation and the replica method, from the replica symmetric
solution. Both methods require the use of different "partial" magnetizations,
associated with different layers of the multiplex network, and yield
qualitatively similar results. If the layers are strongly heterogeneous the
critical temperature differs noticeably from that for the Ising model on a
network being a superposition of the two layers, evaluated in the mean-field
approximation neglecting the effect of the underlying multiplex structure on
the correlations between the degrees of nodes. The critical temperature
evaluated from the replica symmetric solution depends sensitively on the
correlations between the degrees of nodes in different layers and shows
satisfactory quantitative agreement with that obtained from Monte Carlo
simulations. The critical behavior of the magnetization for the model with
strongly heterogeneous layers can depend on the distributions of the degrees of
nodes and is then determined by the properties of the most heterogeneous layer.
